Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Particle streak velocimetry-optical coherence tomography: a novel method for multidimensional imaging of microscale fluid flows

Open Access Open Access

Abstract

We present a new OCT method for flow speed quantification and directional velocimetry: particle streak velocimetry-OCT (PSV-OCT). PSV-OCT generates two-dimensional, 2.5-vector component (vx,|vy|,vz) maps of microscale flow velocity fields. Knowledge of 2.5-vector components also enables the estimation of total flow speed. The enabling insight behind PSV-OCT is that tracer particles in sparsely-seeded fluid flow trace out streaks in (x,z,t)-space. The streak orientations in x-t and z-t yield vx and vz, respectively. The in-plane (x-z plane) residence time yields the out-of-plane speed |vy|. Vector component values are generated by fitting streaks to a model of image formation that incorporates equations of motion in 3D space. We demonstrate cross-sectional estimation of (vx,|vy|,vz) in two important animal models in ciliary biology: Xenopus embryos (tadpoles) and mouse trachea.

© 2016 Optical Society of America

Corrections

Kevin C. Zhou, Brendan K. Huang, Ute A. Gamm, Vineet Bhandari, Mustafa K. Khokha, and Michael A. Choma, "Particle streak velocimetry-optical coherence tomography: a novel method for multidimensional imaging of microscale fluid flows: erratum," Biomed. Opt. Express 7, 2360-2361 (2016)
https://opg.optica.org/boe/abstract.cfm?uri=boe-7-6-2360

1. Introduction

Optical methods are indispensable in the study of cilia-driven fluid flow [1]. Quantifying cilia-driven fluid flow has relevance in respiratory diseases such as asthma, chronic obstructive pulmonary disease, cystic fibrosis, and primary ciliary dyskinesia. Cilia-driven fluid flow also plays a role in cerebrospinal fluid flow in the central nervous system and ova transport in the oviducts. Recently, there has been growing interest in exploiting the cross-sectional nature of OCT to quantify ciliary physiology [2–10]. Cross-sectional imaging of cilia-driven fluid flow is important for several reasons. First, since cilia-driven fluid flow is a surface-driven flow, it is important to localize flow velocity measurements with respect to the surface that is generating those velocities. Second, since cilia-driven fluid flow is not amenable to simplifying assumptions (e.g. parabolic velocity profile for Poiseuille flow), there are not straightforward models to which data can be fit. For example, while three vector component information about parabolic flow can be inferred from, for example, one vector component measurements of vascular blood flow by modeling the vessels as a cylindrical tube [11–15], such model fitting is, in general, not possible for cilia-driven fluid flow.

Previous velocimetry techniques developed in OCT include Doppler-OCT [16–18], multiangle Doppler-OCT [19–21], particle tracking velocimetry (PTV) [2], digital particle image velocimetry (DPIV) [9, 22–24], speckle tracking [25], dynamic light scattering (DLS) [26], and directional DLS [27, 28]. Doppler-OCT is able to estimate only the axial component of the velocity by analyzing the rate of phase evolution over time. PTV-OCT tracks the center-of-mass location of sparsely distributed tracer particles in 2D space in between frames and hence gets the two in-plane velocity components. DPIV-OCT also estimates in-plane velocity components using frame-to-frame displacement, but operates in an intermediate scattering regime by analyzing spatial correlation at displaced spatial locations. Speckle tracking is a similar idea, but operates in a densely scattering regime in which particles are unresolvable. DLS-OCT looks at temporal decorrelation of densely seeded, sub-resolution particles, with the idea that faster moving particles will cause the signal to decorrelate faster. Assuming that the effects of diffusion are accounted for or that diffusion is small compared to the translational motion, the total speed can be estimated. Directional DLS-OCT is an extension that imparts a directional scan bias that allows recovery of an additional lateral velocity component. In both methods, Doppler can be used to get the axial component. Another recently reported method, called OCT micro-particle image velocimetry [29], analyzes the intensity trace as a function of time as a single particle enters and exits the beam, which, like DLS, enables the estimation of total speed. Table 1 summarizes these methods.

Tables Icon

Table 1. A comparison of velocimetry techniques in OCT

We observed that these techniques can be broadly categorized as either spatial- or temporal-correlation techniques. In particular, techniques such as Doppler and DLS can estimate velocity based on how the signal changes in one particular location over time, while techniques such as PTV, DPIV, and speckle tracking examine correlations at neighboring spatial locations. Herein, we present a novel OCT-based velocimetry method called particle streak velocimetry-OCT (PSV-OCT), which combines elements of both spatial (PTV) and temporal (DLS or OCT micro PIV) techniques. PSV-OCT generates 2D, 2.5-vector component flow velocity fields v(x,z) = (vx(x,z),|vy(x,z)|,vz(x,z)) from time series of 2D B-scans. x and z are the in-plane (lateral and axial, respectively) image axes and y is out-of-plane axis. Like PTV-OCT, PSV-OCT generates flow velocity estimates from images of a sparsely seeded flow field. While PTV tracks in-plane center-of-mass particle position [2, 30], PSV-OCT fits a spatiotemporal volumetric streak generated by particle motion in (x,z,t)-space to a model of the streak that is parameterized by (vx,|vy|,vz). Because PSV-OCT does not require identification of individual particles, it avoids the particle-pairing issues and simplifies segmentation. As such, our method can operate in a potentially denser particle regime than PTV-OCT. Our technique also enables the extraction of the out-of-plane speed |vy|. Also note that (vx,|vy|,vz) is sufficient for estimating total flow speed |v|. PTV-OCT falls short in this respect, because it does not account for differing intensities of the particles nor the lateral beam waist and hence does not obtain quantitative out-of-plane speed estimates [2]. Another feature of our method, which will be clarified below, is that our method does not require calibration of the OCT beam waist—effectively, this means that we do not need to measure the beam width, nor do we have to account for the depth-dependent beam width in Fourier-domain OCT. Overall, modeling of the spatiotemporal OCT intensity signal in (x,z,t)-space enables detailed multidimensional flow velocity imaging of cilia-driven fluid flow.

2. Mathematical description of particle streak velocimetry-OCT (PSV-OCT)

Particles in sparsely-seeded flow trace out streaks in (x,y,t)-space (Fig. 1). More formally, a 2D intensity OCT image of N particles taken in the xz-plane at a position yn along the y-axis (out-of-plane direction) can be modeled as:

image(r)n=1N[psfy(y)δy(y)]y=yn[psfr(r)δr(rrn)].
r = (x,z) is particle position in the xz-plane, psfr(r) is the intensity point-spread function (PSF) in the xz-plane (not necessarily isotropic in OCT), psfy(y) is the PSF along the y-axis, (xn,zn,yn) = (rn,yn) is the location of the nth particle, and δy(y) and δr(r) are the 1D the 2D Dirac delta functions in their subscripted domains. If a series of images are acquired over time, and if each particle moves at a constant velocity over the acquisition period, then Eq. (1) is modified to account for this in-plane trajectory. In particular, rn and yn are no longer a constant, but rather a linear function in time:
rn(t)=r0,n+vr,n×(tt0)
yn(t)=|vy,n|(tt0).
Here, vr,n is the is the vectorial velocity of the nth particle in the xz plane. Each particle traces out a streak in (x,z,t)-space as it moves in time. The width of the streak is given by the PSF width. The characteristic length (e.g. 1/e2 length, or the σ parameter of a Gaussian) of the streak is inversely proportional to the out-of-plane speed. Hence, the spatial-temporal orientation and length of a streak gives us three pieces of information: the in-plane velocity vector components (vx,vz) and the out-of-plane speed |vy|. In-plane particle motion traces out the shape of the in-plane PSF (psfx,z(x,z)). The characteristic length of the streak is inversely proportional to the out-of-plane speed and directly proportional to the point-spread function width along the y-axis. The slope of the line segment generated by translation of the nth particle contains information about the in-plane speed. Given these pieces of information, we can write an expression that models each streak as an oriented 3D Gaussian function. We can fit this model to the streaks in the image (Fig. 1) to extract vx, vz, and |vy| for each streak. Thus, assuming an anisotropic Gaussian PSF, the streak is an oriented oblong Gaussian:
I(x,z,t)=Ioexp((tto)22σt2)×exp((x[xo+vx(tto)])22σxy2(z[zo+vz(tto)])22σz2)+Iback
Here, σxy is the lateral point spread function width, σz is the axial point spread function width, Iback is the background intensity and I0 + Iback is the peak intensity of the streak, which occurs at the coordinate (x0,z0,t0). σt is the streak characteristic length and |vy| = σxy / σt. There is sign ambiguity in vy because estimation of vy depends not on streak orientation in (x,z,t)-space but rather streak length: particles moving at + vy generate the same streak length as particles moving at -vy. For fitting, vx was expanded to vx = tan(a)cos(b) and vz was expanded to vz = tan(a)sin(b). For a (x’,z’,t’) coordinate system centered on an individual streak, a is the angle that the streak makes with the t’ axis and b is the angle that the projection of the streak onto the x’z’-plane makes with the x’-axis. tan(a) is the total in-plane (x’z’-plane) speed (Fig. 2). The cos(b) and sin(b) operations project that total speed onto the x’ and z’ axes, respectively. We used the lsqcurvefit() nonlinear least squares regression function in MATLAB (Natick, MA, USA) to fit each streak to the model Eq. (3). Because the fitting procedure estimates σxy, the method does not require prior knowledge of the point spread function width. In other words, our method is calibration-free.

 figure: Fig. 1

Fig. 1 3D (x,z,t) particle streaks generated by the ciliated skin of a Xenopus embryo. a) A single frame from a time series of 2D images (Visualization 1). The bright dots are particles to be tracked. b-d) The entire 3D stack (Visualization 2).

Download Full Size | PDF

 figure: Fig. 2

Fig. 2 Two-dimensional (z,t) streak model. The streak is the product of the in-plane trajectory and the in-plane residence time, σt = σxy / |vy|. If σt is small (top row), the in-plane residence is small and thus the out-of-plane speed is fast. If σt is large (bottom row), the speed is slow. Sign ambiguity in vy arises because σt is insensitive to the sign of vy.

Download Full Size | PDF

3. Materials and methods

3.1 Materials

All OCT measurements were made using a spectral domain OCT system (Thorlabs Telesto) with 1325 nm center wavelength. Using subresolution particles, we estimated the lateral full width at half max (FWHM) to be 8-9 μm near the focus, and the axial FWHM to be 12 μm in air (9 μm in water). The B-scan images were composed of 2,048 A-scans, scanned over a lateral range of 3-3.2 mm (~1.5 μm per A-scan) at a rate of 28,000 A-scans per second (13.7 frames per second). A-scans were 512 pixels over a depth range (in air) of 2.494 mm. Typically, 100 frames were acquired for the videos. This duration enabled imaging a sufficient number of particles to generate a reasonably dense flow velocity field. Likewise, lower particle concentrations require longer OCT video acquisition times. We used 5 μm polystyrene microspheres (Bangs labs) as flow tracers in all experiments. 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.

3.2 PSV-OCT algorithm and workflow

Our PSV-OCT procedure is as follows. First, stationary objects such as flow phantom walls and animal bodies are removed based in a minimum projection across time and mathematical morphological operations to remove false regions. Next, a 3D (x,z,t) image is subject to a liberal intensity threshold to identify 3D streaks, so as to capture as many of the particle streaks as possible. To filter out falsely identified streaks due to noise, we used one of the following methods: 1) the temporal dimension of the streak must be at least n pixels, where n is user-defined (e.g., 10), or 2) we calculated the orientation of each streak and kept the middle m% (e.g., 95%). The procedure for estimating 3D orientation is described in the next paragraph. For method 1, n and the intensity threshold should be jointly selected, depending on the background noise. Other strategies for filtering out false streaks, such as mathematical morphological operations that mitigate noise, may be used as well. Further filtering can be done after the fitting based on fitting diagnostics, which we describe later in this section.

We use principal component analysis (PCA) to estimate the orientation of individual binary streaks segmented by the intensity threshold. This estimate will be crude because of the arbitrariness of the intensity threshold. This involved obtaining the 3D coordinates (x,z,t) of the region of connected pixels and finding the eigenvectors of the 3-by-3 covariance matrix. The eigenvector corresponding to the largest eigenvalue is the one that points in the direction of the streak. We do not use the other principal components. This orientation from PCA serves as the initial guess for our fitting procedure. Using the orientation parameters, we define a local 3D neighborhood around the streak that defines the domain over which nonlinear fitting is to occur (Fig. 3). This region should be on the liberal side so as not to cut off part of the streak with the tradeoff being that, the larger the region, the more time the fitting will take. Lastly, using information from the PCA to inform initial estimates, each streak is fit to the oblong Gaussian in Eq. (3) in order to estimate (vx,|vy|,vz). Using PCA results as initial estimates facilitates faster convergence of the nonlinear fitting process. All fitting was performed on the linear intensity data; that is, we do not use the logarithm of the intensity values as is typically used when displaying OCT images, because doing so invalidates the Gaussian model.

 figure: Fig. 3

Fig. 3 A sample streak fit. The three panels show the three possible 2D projections from the 3D streak (x-z, x-t, and z-t). The smaller, background circles are the image pixels, whose color represents the pixel intensity. As such, the yellow-white circles are where the streak is located. The large region around the streak was the region over which the fitting was done; the orientation of this region is derived from PCA. The fit itself is the continuous region in the background. The turquoise circles represent particle tracking based on center of mass per frame.

Download Full Size | PDF

Figure 3 shows a sample streak. A region was isolated around the streak based on the PCA-derived orientation, and the fitted oblong Gaussian curve is shown to match up quite closely to the intensity values. For comparison, we also show the particle tracking result.

4. Results

4.1 Validation in a flow phantom

As an initial demonstration, we used PSV-OCT to quantify flow in a cylindrical flow phantom. The image plane orientation with respect to the channel ranged from ψ = 0° (perpendicular to the channel axis) to 90° (parallel to the channel axis; Fig. 4(a)). We tilted the axis of the cylindrical channel (φ) so that there would be a non-zero vertical velocity component. If our technique works, then at ψ = 0° the out-of-plane speed should be the greatest while the in-plane (xz) velocity should be approximately 0, while at ψ = 90° the in-plane velocity should be maximum while the out-of-plane speed should be approximately 0. The axial velocity (z) should be approximately constant across all angles because the image plane is being rotated about the z-axis. Lastly, the total speed should remain constant as a function of ψ. Figure 4(a)-4(c) confirms these expected patterns, validating our method. Peak velocities measured using PSV-OCT for ψ = 90° (~1.7 mm/s) also are consistent with peak velocities predicted from the tube geometry (~1.75 mm/s) and those obtained using our PTV-OCT approach [2] (~1.5 mm/s). We typically use spatial binning in PTV-OCT, which will underestimate the peak value of a parabolic profile. Also, shape of the PSV-OCT flow profiles were slightly deviated off-center, which we attribute to strong refraction at air/glass and glass/water interfaces. Lateral and axial point-spread function widths (σxy and σz, respectively) were consistent with width estimates made by imaging 1 μm particles. Lastly, using PSV-OCT estimates of σxy as a function of axial position z (Fig. 4(d)) were used to estimate an imaging center wavelength of l = 1312 nm, which is quite close to the Thorlabs reported value of 1325 nm. We used [31]

w(z)=wo1+[λ(zzf)πwo2]2
for estimating l using nonlinear fitting. Here, w is the 1/e2 beam radius (w = 2σxy), wo is the waist, and zf is the axial location of the beam waist.

 figure: Fig. 4

Fig. 4 PSV-OCT in a capillary flow phantom. (a) The intersection of the imaging plane (x0,z0) with a cylinder of radius r’ leads to an ellipse with major and minor axes a and b, determined by the angle of imaging with respect to the tube. The tube was rotated about the global x-axis by an angle φ = 9°. The imaging plane was then rotated about the global z-axis by an angle ψ, from 0° to 90°. Intensity images at three ψ are shown. Flow was estimated using PSV-OCT at each of these ψ. (b) Predicted and estimated in-plane speed (vx2 + vz2)1/2, out-of-plane speed |vy|, and total speed v = (vx2 + vy2 + vz2)1/2 at three angles with respect to flow. When ψ = 0°, the imaging plane is orthogonal to the direction of flow except for a small component due to vertical tilt φ. At ψ = 90°, flow is predicted to be completely in plane. (c) Peak flow speed of the parabolic flow profile, estimated by non-linear least squared fitting to Poiseuille flow in a tube, as a function of ψ. The error here, a 95% confidence interval (CI), is estimated by calculating the Jacobian during non-linear fitting. (d) Lateral point spread function (PSF) width σxy as a function of axial location (ψ = 0).

Download Full Size | PDF

4.2 PSV-OCT in cilia-driven fluid flow in Xenopus tropicalis embryo

After demonstration in a calibrated flow phantom, we demonstrated PSV-OCT in two important animal models of ciliary disease: Xenopus and mouse. Figure 5 visualizes 2-dimensional, 2.5-vector component cilia-driven fluid flow fields generated by a ciliated Xenopus tropicalis embryo. Xenopus (frog) embryos have a ciliated skin that generates a vigorous, head-to-tail flow, making them a versatile and readily manipulated model in ciliary biology [32]. For imaging, embryos were pharmacologically immobilized using benzocaine, so that the only motion in the images is fluid flow due to the cilia. In Fig. 5, in-plane speed is encoded by both quiver size and quiver color. Out-of-plane speed is indicated by the size and color of a circle at the base of each quiver. Both in-plane and out-of-plane speeds decrease with increasing distance away from the ciliated embryo surface. Having access to out-of-plane speed confirms that, for most of the flow field near the surface, the flow is dominated by predominantly head-to-tail (i.e. in-plane flow). Interestingly, however, the out-of-plane flow speeds begin to approach in-plane flow speeds at the boundary between the ciliated embryo and the base of the well in which the embryo sits. That is, the presence of a large out-of-plane flow component (and recirculatory flow patterns) is a geometrical consequence of flow near the edges of a finite-sized ciliated surface.

 figure: Fig. 5

Fig. 5 Quantification of cilia-driven flow in X. tropicalis tadpole epithelium. (a) Maximum intensity projection image showing streaks traced out by particles, with primary direction of flow from head (h) to tail (t). Overlay of vector flow field estimated with PSV-OCT showing in-plane velocity (vx,vz) and direction indicated by length and direction of arrow, and out-of-plane speed |vy| indicated by size and color of circle placed at base of each arrow. Axes: h- head; t- tail; l-left, r-right. (b,c) inset of location of region of higher out-of-plane flow near the (a) head of the embryo and (b) tail of the embryo.

Download Full Size | PDF

4.3 PSV-OCT in cilia-driven fluid flow in mouse trachea

Lastly, we demonstrated PSV-OCT in postnatal day (PND) 21 ex vivo mouse trachea. We used PND 21 mice since, at that age, 1325 nm OCT can visualize cilia-driven fluid flow directly through the tracheal wall [8, 9]. Tracheae were kept at 37°C. PSV-OCT was particularly useful in quantifying the degree of out-of-plane flow when complex recirculatory flow profiles are present. We minimize recirculation in the context of experimental diagnostic tests [8] 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. For example, we found here that significant out-of-plane motion occurs near the surface of the trachea, while flow is more directed down the axis of the trachea in the center of the lumen (Fig. 6(b), 6(c)). Such findings would not be possible with standard particle tracking velocimetry.

 figure: Fig. 6

Fig. 6 Quantification of cilia-driven flow in mouse trachea. (a) Maximum intensity projection of flow in tracheal lumen showing tail-head flow directed along the bottom and top surface, and recirculatory flow in the center of the trachea due to closed outflow. Overlay of vector flow field estimated with PSV-OCT showing in-plane velocity (vx,vz) and direction indicated by length and direction of arrow, and out-of-plane speed |vy| indicated by size and color of circle placed at base of each arrow. Axes: h- head; t- tail; d- dorsal; v- ventral. (b,c) Quantification of relative values of out-of-plane flow |vy| versus in-plane flow (vx2 + vz2)1/2. In two regions of interest (ROI), one near the surface of the tracheal wall (ROI 1, red) and one in the center of the lumen (ROI 2, green). (c) Histogram of ratio of values of |vy|/(vx2 + vz2)1/2 . Note the logarithmic scale. Values less than one indicate flows that are primarily in-plane, while values greater than 1 indicate flows that are primarily out-of-plane.

Download Full Size | PDF

4.4 Nonlinear least squares fit diagnostics

Because we are fitting a model to the image data, we have access to fit diagnostics that can tell us how well the data match the model. This simplest approach is to analyze the residuals in a bulk manner; that is, we can look at the r value, the correlation coefficient between the fit and the pixel values, or the root mean-square error (RMSE). Low values for these metrics indicates either that the fitting algorithm did not converge to a good local minimum, and/or that the streak identified may not have been a streak. In either case, we would want to omit these points or perhaps reanalyze them. In our presented data, we simply omitted based on user-defined thresholds. Figure 7 shows the the correlation coefficients and the square-root of the MSEs plotted jointly for the mouse trachea data. The fact that most of the fits are of high correlation and relatively low MSEs is indicative that our fitting was accurate. Furthermore, as reported above, the lateral beam waist values also led to a good estimate of the imaging wavelength, a good reality check that the fits were reasonable (Fig. 4(d)). We drew thresholds on the MSE and correlation coefficient to exclude poor fits.

 figure: Fig. 7

Fig. 7 The correlation coefficients and root mean-square errors (with arbitrary intensity units) of all streak fits (prior to any goodness-of-fit filtering) for the tadpole data and the mouse trachea data. Ideally, the fits will be of low mean-square error and correlation coefficients close to 1. It is thus reassuring that the points in the figure are clustered close to the lower right corner.

Download Full Size | PDF

5. Discussion and conclusion

In sum, PSV-OCT is a novel method for multi-vector component velocimetry of microscale biological fluid flows. Our method is distinct compared to prior streak velocimetry-type methods [33, 34] since our estimation process fits the streak data to a function with continuously varying intensity. That is, prior methods neglect image intensity and PSF information but rather identify the beginning and end of a streak and use those endpoints, which are arbitrarily defined, as an interval displacement over some time interval. PSV-OCT is thus able to obtain out-of-plane speed quantitatively because it does not acknowledge hard cutoffs of the streaks, as in principle the PSF is infinite, and it accounts for the depth-dependent beam waist in OCT.

It should be noted that Mujal et al. [29] observed that an out-of-plane speed estimate can be derived from particle duration time within the PSF. The particle duration time, in turn, can be estimated from a time-varying intensity trace. To summarize, their method fits to the parabolic log-intensity (equivalent to our Gaussian model of the linear intensity) profile as a particle traverses the imaging PSF. In this way, their duration estimate is independent of absolute levels of backscatter intensity. However, what differentiates our technique from their method is that their technique was demonstrated in M-mode and does not exploit in-plane trajectory estimation (e.g. PTV, PSV). PSV-OCT, on the other hand, looks at the intensity profile over its spatiotemporal trajectory because the method simultaneously tracks in-plane particle position as well as the rate at which the particle is leaving the PSF. In other words, as mentioned in the introduction, our technique is both temporal- and spatial-correlative while previous techniques use one or the other.

An important assumption of our work is that the trajectories of the particles are at least locally linear. Indeed, we are exploiting the finiteness of the OCT beam waist to sample local windows of a potentially curved trajectory, while techniques that use projection images may exhibit curved trajectories. Such curves are apparent in the projection images in Figs. 5 and 6, from which the blending of multiple particles as a result of projecting gives deceptively curved appearances. However, the supplemental videos, which show the full x-z-t data set, show that the streaks are actually quite linear. Hence, PSV-OCT is advantageous over a technique previously termed particle streak velocimetry [14], which uses projection images from a camera with a long integration time. In particular, there is a loss of time-resolution and a limit on how dense the particle distribution can be to avoid overlapping trajectories. Also, unlike PSV-OCT, they did not obtain out-of-plane speed because they did not account for the PSF. Another note is that they fit parabolas to the streaks, noting that they were not linear. In principle, we could do the same thing with PSV-OCT by modifying Eqs. (2a) and (2b) to be parabolic rather than linear and use a model selection criterion to account for potentially curved trajectories (one can think of this as a second-order Taylor expansion). However, we found this to be unnecessary.

Finally, we acknowledge that there are upper and lower limits on the speed of particles measurable by PSV-OCT. An implicit limit is that the particle translations are sub-resolution; therefore, if the particles move too quickly, the streaks would not be contiguous and would hence be missed by segmentation criteria. However, if the discontiguous streaks could be identified, the fitting would still be able to proceed. The lower bound is limited by Brownian motion, which violates our assumption of dominant translational motion. The lower bound also is limited by imaging noise (e.g. shot noise), since noise will result in a particle image profile that deviates from a smooth Gaussian curve. In general, the longer the particle remains in-plane (the slower the out-of-plane speed), the better both in-plane and out-of-plane velocity estimates are. Another concern is that streaks may be clipped at the beginning or end of the videos. In these cases, as long as the streak is long enough to pass the initial filtering criteria, the streak will still be fit. Even if it is clipped, the remaining parts should still follow the Gaussian intensity profile. The fit may be noisier than that from a complete streak, however.

In sum, PSV-OCT is a novel contribution to a growing body of work that approaches transverse velocimetry and total speed estimation by incorporating equations of motion into models of OCT image formation (see, for example [9, 26, 27, 35, 36],). We expect this technique to be applicable in a wide range of biological settings in which discrete scattering objects (e.g. red blood cells) can be identified.

References and links

1. B. K. Huang and M. A. Choma, “Microscale imaging of cilia-driven fluid flow,” Cell. Mol. Life Sci. 72(6), 1095–1113 (2015). [CrossRef]   [PubMed]  

2. S. Jonas, D. Bhattacharya, M. K. Khokha, and M. A. Choma, “Microfluidic characterization of cilia-driven fluid flow using optical coherence tomography-based particle tracking velocimetry,” Biomed. Opt. Express 2(7), 2022–2034 (2011). [CrossRef]   [PubMed]  

3. A. L. Oldenburg, R. K. Chhetri, D. B. Hill, and B. Button, “Monitoring airway mucus flow and ciliary activity with optical coherence tomography,” Biomed. Opt. Express 3(9), 1978–1992 (2012). [CrossRef]   [PubMed]  

4. L. Liu, K. K. Chu, G. H. Houser, B. J. Diephuis, Y. Li, E. J. Wilsterman, S. Shastry, G. Dierksen, S. E. Birket, M. Mazur, S. Byan-Parker, W. E. Grizzle, E. J. Sorscher, S. M. Rowe, and G. J. Tearney, “Method for quantitative study of airway functional microanatomy using micro-optical coherence tomography,” PLoS One 8(1), e54473 (2013). [CrossRef]   [PubMed]  

5. L. Liu, S. Shastry, S. Byan-Parker, G. Houser, K. K Chu, S. E. Birket, C. M. Fernandez, J. A. Gardecki, W. E. Grizzle, E. J. Wilsterman, E. J. Sorscher, S. M. Rowe, and G. J. Tearney, “An autoregulatory mechanism governing mucociliary transport is sensitive to mucus load,” Am. J. Respir. Cell Mol. Biol. 51(4), 485–493 (2014). [CrossRef]   [PubMed]  

6. S. E. Birket, K. K. Chu, L. Liu, G. H. Houser, B. J. Diephuis, E. J. Wilsterman, G. Dierksen, M. Mazur, S. Shastry, Y. Li, J. D. Watson, A. T. Smith, B. S. Schuster, J. Hanes, W. E. Grizzle, E. J. Sorscher, G. J. Tearney, and S. M. Rowe, “A functional anatomic defect of the cystic fibrosis airway,” Am. J. Respir. Crit. Care Med. 190(4), 421–432 (2014). [CrossRef]   [PubMed]  

7. B. K. Huang, U. A. Gamm, S. Jonas, M. K. Khokha, and M. A. Choma, “Quantitative optical coherence tomography imaging of intermediate flow defect phenotypes in ciliary physiology and pathophysiology,” J. Biomed. Opt. 20(3), 030502 (2015). [CrossRef]   [PubMed]  

8. U. A. Gamm, B. K. Huang, M. Syed, X. Zhang, V. Bhandari, and M. A. Choma, “Quantifying hyperoxia-mediated damage to mammalian respiratory cilia-driven fluid flow using particle tracking velocimetry optical coherence tomography,” J. Biomed. Opt. 20(8), 080505 (2015). [CrossRef]   [PubMed]  

9. B. K. Huang, U. A. Gamm, V. Bhandari, M. K. Khokha, and M. A. Choma, “Three-dimensional, three-vector-component velocimetry of cilia-driven fluid flow using correlation-based approaches in optical coherence tomography,” Biomed. Opt. Express 6(9), 3515–3538 (2015). [CrossRef]   [PubMed]  

10. R. Ansari, C. Buj, M. Pieper, P. König, A. Schweikard, and G. Huttmann, “Micro-anatomical and functional assessment of ciliated epithelium in mouse trachea using optical coherence phase microscopy,” Opt. Express 23(18), 23217–23224 (2015). [CrossRef]   [PubMed]  

11. A. Davis, J. Izatt, and F. Rothenberg, “Quantitative Measurement of Blood Flow Dynamics in Embryonic Vasculature Using Spectral Doppler Velocimetry,” Anat. Rec. (Hoboken) 292(3), 311–319 (2009). [CrossRef]   [PubMed]  

12. S. Makita, T. Fabritius, and Y. Yasuno, “Quantitative retinal-blood flow measurement with three-dimensional vessel geometry determination using ultrahigh-resolution Doppler optical coherence angiography,” Opt. Lett. 33(8), 836–838 (2008). [CrossRef]   [PubMed]  

13. Z. Ma, A. Liu, X. Yin, A. Troyer, K. Thornburg, R. K. Wang, and S. Rugonyi, “Measurement of absolute blood flow velocity in outflow tract of HH18 chicken embryo based on 4D reconstruction using spectral domain optical coherence tomography,” Biomed. Opt. Express 1(3), 798–811 (2010). [CrossRef]   [PubMed]  

14. A. S. G. Singh, C. Kolbitsch, T. Schmoll, and R. A. Leitgeb, “Stable absolute flow estimation with Doppler OCT based on virtual circumpapillary scans,” Biomed. Opt. Express 1(4), 1047–1058 (2010). [CrossRef]   [PubMed]  

15. L. Qi, J. Zhu, A. M. Hancock, C. Dai, X. Zhang, R. D. Frostig, and Z. Chen, “Fully distributed absolute blood flow velocity measurement for middle cerebral arteries using Doppler optical coherence tomography,” Biomed. Opt. Express 7(2), 601–615 (2016). [CrossRef]   [PubMed]  

16. J. A. Izatt, M. D. Kulkarni, S. Yazdanfar, J. K. Barton, and A. J. Welch, “In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomography,” Opt. Lett. 22(18), 1439–1441 (1997). [CrossRef]   [PubMed]  

17. W. Drexler, M. Liu, A. Kumar, T. Kamali, A. Unterhuber, and R. A. Leitgeb, “Optical coherence tomography today: speed, contrast, and multimodality,” J. Biomed. Opt. 19(7), 071412 (2014). [CrossRef]   [PubMed]  

18. Z. Chen, T. E. Milner, D. Dave, and J. S. Nelson, “Optical Doppler tomographic imaging of fluid flow velocity in highly scattering media,” Opt. Lett. 22(1), 64–66 (1997). [CrossRef]   [PubMed]  

19. A. Røyset, T. Støren, F. Stabo-Eeg, and T. Lindmo, “Quantitative measurements of flow velocity and direction using transversal Doppler optical coherence tomography,” Proc. SPIE 6079, 607925 (2006).

20. Y.-C. Ahn, W. Jung, and Z. Chen, “Quantification of a three-dimensional velocity vector using spectral-domain Doppler optical coherence tomography,” Opt. Lett. 32(11), 1587–1589 (2007). [CrossRef]   [PubMed]  

21. R. Haindl, W. Trasischker, A. Wartak, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Total retinal blood flow measurement by three beam Doppler optical coherence tomography,” Biomed. Opt. Express 7(2), 287–301 (2016). [CrossRef]   [PubMed]  

22. C. E. Willert and M. Gharib, “Digital particle image velocimetry,” Exp. Fluids 10(4), 181–193 (1991). [CrossRef]  

23. A. Buchsbaum, M. Egger, I. Burzic, T. Koepplmayr, M. Aigner, J. Miethlinger, and M. Leitner, “Optical coherence tomography based particle image velocimetry (OCT-PIV) of polymer flows,” Opt. Lasers Eng. 69, 40–48 (2015). [CrossRef]  

24. C.-Y. Chen, P. G. Menon, W. Kowalski, and K. Pekkan, “Time-resolved OCT-μPIV: a new microscopic PIV technique for noninvasive depth-resolved pulsatile flow profile acquisition,” Exp. Fluids 54(1), 1–9 (2013). [CrossRef]  

25. C. Sun, B. Standish, and V. X. D. Yang, “Optical coherence elastography: current status and future applications,” J. Biomed. Opt. 16(4), 043001 (2011).

26. J. Lee, W. Wu, J. Y. Jiang, B. Zhu, and D. A. Boas, “Dynamic light scattering optical coherence tomography,” Opt. Express 20(20), 22262–22277 (2012). [CrossRef]   [PubMed]  

27. B. K. Huang and M. A. Choma, “Resolving directional ambiguity in dynamic light scattering-based transverse motion velocimetry in optical coherence tomography,” Opt. Lett. 39(3), 521–524 (2014). [CrossRef]   [PubMed]  

28. K. C. Zhou, B. K. Huang, H. Tagare, and M. A. Choma, “Improved velocimetry in optical coherence tomography using Bayesian analysis,” Biomed. Opt. Express 6(12), 4796–4811 (2015). [CrossRef]   [PubMed]  

29. M. Mujat, R. D. Ferguson, N. Iftimia, D. X. Hammer, I. Nedyalkov, M. Wosnik, and H. Legner, “Optical coherence tomography-based micro-particle image velocimetry,” Opt. Lett. 38(22), 4558–4561 (2013). [CrossRef]   [PubMed]  

30. R. J. Adrian, “Particle-Imaging Techniques for Experimental Fluid-Mechanics,” Annu. Rev. Fluid Mech. 23(1), 261–304 (1991). [CrossRef]  

31. B. E. A. Saleh and M. C. Teich, “Beam Optics,” in Fundamentals of Photonics, 1st ed. (John Wiley & Sons, Inc., New York, 1991).

32. M. E. Werner and B. J. Mitchell, “Understanding ciliated epithelia: the power of Xenopus,” Genesis 50(3), 176–185 (2012). [CrossRef]   [PubMed]  

33. P. E. Dimotakis, F. D. Debussy, and M. M. Koochesfahani, “Particle Streak Velocity-Field Measurements in a Two-Dimensional Mixing Layer,” Phys. Fluids 24(6), 995–999 (1981). [CrossRef]  

34. E. A. V. Jones, M. H. Baron, S. E. Fraser, and M. E. Dickinson, “Measuring hemodynamic changes during mammalian development,” Am. J. Physiol. Heart Circ. Physiol. 287(4), H1561–H1569 (2004). [CrossRef]   [PubMed]  

35. V. J. Srinivasan, H. Radhakrishnan, E. H. Lo, E. T. Mandeville, J. Y. Jiang, S. Barry, and A. E. Cable, “OCT methods for capillary velocimetry,” Biomed. Opt. Express 3(3), 612–629 (2012). [CrossRef]   [PubMed]  

36. N. Weiss, T. G. van Leeuwen, and J. Kalkman, “Localized measurement of longitudinal and transverse flow velocities in colloidal suspensions using optical coherence tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 88(4), 042312 (2013). [CrossRef]   [PubMed]  

Supplementary Material (2)

NameDescription
Visualization 1: MOV (7905 KB)      video of polystyrene beads being pushed by the ciliated surface of a tadpole embryo
Visualization 2: MOV (10067 KB)      the same video presented in (x,z,t)-space

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1
Fig. 1 3D (x,z,t) particle streaks generated by the ciliated skin of a Xenopus embryo. a) A single frame from a time series of 2D images (Visualization 1). The bright dots are particles to be tracked. b-d) The entire 3D stack (Visualization 2).
Fig. 2
Fig. 2 Two-dimensional (z,t) streak model. The streak is the product of the in-plane trajectory and the in-plane residence time, σt = σxy / |vy|. If σt is small (top row), the in-plane residence is small and thus the out-of-plane speed is fast. If σt is large (bottom row), the speed is slow. Sign ambiguity in vy arises because σt is insensitive to the sign of vy.
Fig. 3
Fig. 3 A sample streak fit. The three panels show the three possible 2D projections from the 3D streak (x-z, x-t, and z-t). The smaller, background circles are the image pixels, whose color represents the pixel intensity. As such, the yellow-white circles are where the streak is located. The large region around the streak was the region over which the fitting was done; the orientation of this region is derived from PCA. The fit itself is the continuous region in the background. The turquoise circles represent particle tracking based on center of mass per frame.
Fig. 4
Fig. 4 PSV-OCT in a capillary flow phantom. (a) The intersection of the imaging plane (x0,z0) with a cylinder of radius r’ leads to an ellipse with major and minor axes a and b, determined by the angle of imaging with respect to the tube. The tube was rotated about the global x-axis by an angle φ = 9°. The imaging plane was then rotated about the global z-axis by an angle ψ, from 0° to 90°. Intensity images at three ψ are shown. Flow was estimated using PSV-OCT at each of these ψ. (b) Predicted and estimated in-plane speed (vx2 + vz2)1/2, out-of-plane speed |vy|, and total speed v = (vx2 + vy2 + vz2)1/2 at three angles with respect to flow. When ψ = 0°, the imaging plane is orthogonal to the direction of flow except for a small component due to vertical tilt φ. At ψ = 90°, flow is predicted to be completely in plane. (c) Peak flow speed of the parabolic flow profile, estimated by non-linear least squared fitting to Poiseuille flow in a tube, as a function of ψ. The error here, a 95% confidence interval (CI), is estimated by calculating the Jacobian during non-linear fitting. (d) Lateral point spread function (PSF) width σxy as a function of axial location (ψ = 0).
Fig. 5
Fig. 5 Quantification of cilia-driven flow in X. tropicalis tadpole epithelium. (a) Maximum intensity projection image showing streaks traced out by particles, with primary direction of flow from head (h) to tail (t). Overlay of vector flow field estimated with PSV-OCT showing in-plane velocity (vx,vz) and direction indicated by length and direction of arrow, and out-of-plane speed |vy| indicated by size and color of circle placed at base of each arrow. Axes: h- head; t- tail; l-left, r-right. (b,c) inset of location of region of higher out-of-plane flow near the (a) head of the embryo and (b) tail of the embryo.
Fig. 6
Fig. 6 Quantification of cilia-driven flow in mouse trachea. (a) Maximum intensity projection of flow in tracheal lumen showing tail-head flow directed along the bottom and top surface, and recirculatory flow in the center of the trachea due to closed outflow. Overlay of vector flow field estimated with PSV-OCT showing in-plane velocity (vx,vz) and direction indicated by length and direction of arrow, and out-of-plane speed |vy| indicated by size and color of circle placed at base of each arrow. Axes: h- head; t- tail; d- dorsal; v- ventral. (b,c) Quantification of relative values of out-of-plane flow |vy| versus in-plane flow (vx2 + vz2)1/2. In two regions of interest (ROI), one near the surface of the tracheal wall (ROI 1, red) and one in the center of the lumen (ROI 2, green). (c) Histogram of ratio of values of |vy|/(vx2 + vz2)1/2 . Note the logarithmic scale. Values less than one indicate flows that are primarily in-plane, while values greater than 1 indicate flows that are primarily out-of-plane.
Fig. 7
Fig. 7 The correlation coefficients and root mean-square errors (with arbitrary intensity units) of all streak fits (prior to any goodness-of-fit filtering) for the tadpole data and the mouse trachea data. Ideally, the fits will be of low mean-square error and correlation coefficients close to 1. It is thus reassuring that the points in the figure are clustered close to the lower right corner.

Tables (1)

Tables Icon

Table 1 A comparison of velocimetry techniques in OCT

Equations (5)

Equations on this page are rendered with MathJax. Learn more.

image( r ) n=1 N [ ps f y ( y ) δ y ( y ) ] y= y n [ ps f r ( r ) δ r ( r r n ) ] .
r n ( t )= r 0,n + v r,n ×(t t 0 )
y n ( t )=| v y,n |(t t 0 ).
I( x,z,t )= I o exp( ( t t o ) 2 2 σ t 2 )×exp( ( x[ x o + v x (t t o ) ] ) 2 2 σ xy 2 ( z[ z o + v z (t t o ) ] ) 2 2 σ z 2 )+ I back
w( z )= w o 1+ [ λ( z z f ) π w o 2 ] 2
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.