Quantitative laser speckle flowmetry of the in vivo microcirculation using sidestream dark field microscopy

We present integrated Laser Speckle Contrast Imaging (LSCI) and Sidestream Dark Field (SDF) flowmetry to provide real-time, noninvasive and quantitative measurements of speckle decorrelation times related to microcirculatory flow. Using a multi exposure acquisition scheme, precise speckle decorrelation times were obtained. Applying SDFLSCI in vitro and in vivo allows direct comparison between speckle contrast decorrelation and flow velocities, while imaging the phantom and microcirculation architecture. This resulted in a novel analysis approach that distinguishes decorrelation due to flow from other additive decorrelation sources. ©2013 Optical Society of America OCIS codes: (170.2655) Functional monitoring and imaging; (170.3880) Medical and biological imaging; (170.6480) Spectroscopy, speckle; (170.0180) Microscopy. References and links 1. B. Fagrell and M. Intaglietta, “Microcirculation: its significance in clinical and molecular medicine,” J. Intern. Med. 241(5), 349–362 (1997). 2. K. R. Mathura, G. J. Bouma, and C. Ince, “Abnormal microcirculation in brain tumours during surgery,” Lancet 358(9294), 1698–1699 (2001). 3. Y. Sakr, M. J. Dubois, D. De Backer, J. Creteur, and J. L. Vincent, “Persistent microcirculatory alterations are associated with organ failure and death in patients with septic shock,” Crit. Care Med. 32(9), 1825–1831 (2004). 4. P. T. Goedhart, M. Khalilzada, R. Bezemer, J. Merza, and C. Ince, “Sidestream Dark Field (SDF) imaging: a novel stroboscopic LED ring-based imaging modality for clinical assessment of the microcirculation,” Opt. Express 15(23), 15101–15114 (2007). 5. J. G. G. Dobbe, G. J. Streekstra, B. Atasever, R. van Zijderveld, and C. Ince, “Measurement of functional microcirculatory geometry and velocity distributions using automated image analysis,” Med. Biol. Eng. Comput. 46(7), 659–670 (2008). 6. A. Fercher and J. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). 7. Z. Wang, S. Hughes, S. Dayasundara, and R. S. Menon, “Theoretical and experimental optimization of laser speckle contrast imaging for high specificity to brain microcirculation,” J. Cereb. Blood Flow Metab. 27(2), 258–269 (2007). 8. H. Cheng, Q. Luo, Q. Liu, Q. Lu, H. Gong, and S. Zeng, “Laser speckle imaging of blood flow in microcirculation,” Phys. Med. Biol. 49(7), 1347–1357 (2004). 9. M. Draijer, E. Hondebrink, T. Leeuwen, and W. Steenbergen, “Review of laser speckle contrast techniques for visualizing tissue perfusion,” Lasers Med. Sci. 24(4), 639–651 (2009). 10. D. D. Duncan and S. J. Kirkpatrick, “Can laser speckle flowmetry be made a quantitative tool?” J. Opt. Soc. Am. A 25(8), 2088–2094 (2008). 11. A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multiexposure speckle imaging,” Opt. Express 16(3), 1975–1989 (2008). 12. R. Bezemer, E. Klijn, M. Khalilzada, A. Lima, M. Heger, J. van Bommel, and C. Ince, “Validation of nearinfrared laser speckle imaging for assessing microvascular (re) perfusion,” Microvasc. Res. 79(2), 139–143 (2010). 13. T. Durduran, R. Choe, W. Baker, and A. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). 14. D. Bicout and G. Maret, “Multiple light scattering in Taylor-Couette flow,” Physica A 210(1-2), 87–112 (1994). #194675 $15.00 USD Received 25 Jul 2013; revised 12 Sep 2013; accepted 27 Sep 2013; published 7 Oct 2013 (C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002347 | BIOMEDICAL OPTICS EXPRESS 2347 15. R. Bonner and R. Nossal, “Model for laser Doppler measurements of blood flow in tissue,” Appl. Opt. 20(12), 2097–2107 (1981). 16. A. B. Parthasarathy, S. M. Kazmi, and A. K. Dunn, “Quantitative imaging of ischemic stroke through thinned skull in mice with Multi Exposure Speckle Imaging,” Biomed. Opt. Express 1(1), 246–259 (2010). 17. J. C. Ramirez-San-Juan, R. Ramos-García, I. Guizar-Iturbide, G. Martínez-Niconoff, and B. Choi, “Impact of velocity distribution assumption on simplified laser speckle imaging equation,” Opt. Express 16(5), 3197–3203 (2008). 18. D. A. Boas and A. K. Dunn, “Laser speckle contrast imaging in biomedical optics,” J. Biomed. Opt. 15, 011109 (2010). 19. P. Zakharov, A. Völker, A. Buck, B. Weber, and F. Scheffold, “Quantitative modeling of laser speckle imaging,” Opt. Lett. 31(23), 3465–3467 (2006). 20. J. W. Goodman, Speckle phenomena in optics: theory and applications (Roberts and Company Publishers, Greenwood Village, CO, 2007). 21. J. D. Briers, G. Richards, and X. W. He, “Capillary blood flow monitoring using laser speckle contrast analysis (LASCA),” J. Biomed. Opt. 4(1), 164–175 (1999). 22. J. W. Goodman, “Some fundamental properties of speckle,” J. Opt. Soc. Am. A 66(11), 1145–1150 (1976). 23. A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic imaging of cerebral blood flow using laser speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). 24. R. Bandyopadhyay, A. Gittings, S. Suh, P. Dixon, and D. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. sci. Instrum. 76, 093110 (2005). 25. T. Yoshimura, “Statistical properties of dynamic speckles,” J. Opt. Soc. Am. A 3(7), 1032–1054 (1986). 26. H. Cheng, Q. Luo, S. Zeng, S. Chen, J. Cen, and H. Gong, “Modified laser speckle imaging method with improved spatial resolution,” J. Biomed. Opt. 8(3), 559–564 (2003). 27. D. D. Duncan, S. J. Kirkpatrick, and R. K. Wang, “Statistics of local speckle contrast,” J. Opt. Soc. Am. A 25(1), 9–15 (2008). 28. J. Qiu, P. Li, W. Luo, J. Wang, H. Zhang, and Q. Luo, “Spatiotemporal laser speckle contrast analysis for blood flow imaging with maximized speckle contrast,” J. Biomed. Opt.15, 016003 (2010). 29. D. M. de Bruin, R. H. Bremmer, V. M. Kodach, R. de Kinkelder, J. van Marle, T. G. van Leeuwen, and D. J. Faber, “Optical phantoms of varying geometry based on thin building blocks with controlled optical properties,” J. Biomed. Opt. 15(2), 025001 (2010). 30. T. L. Alexander, J. E. Harvey, and A. R. Weeks, “Average speckle size as a function of intensity threshold level: comparisonof experimental measurements with theory,” Appl. Opt. 33(35), 8240–8250 (1994). 31. S. J. Kirkpatrick, D. D. Duncan, and E. M. Wells-Gray, “Detrimental effects of speckle-pixel size matching in laser speckle contrast imaging,” Opt. Lett. 33(24), 2886–2888 (2008). 32. O. Thompson, M. Andrews, and E. Hirst, “Correction for spatial averaging in laser speckle contrast analysis,” Biomed. Opt. Express 2(4), 1021–1029 (2011). 33. J. M. Brown and A. J. Giaccia, “The unique physiology of solid tumors: opportunities (and problems) for cancer therapy,” Cancer Res. 58(7), 1408–1416 (1998). 34. C. W. Song, “Effect of local hyperthermia on blood flow and microenvironment: a review,” Cancer Res. 44(10 Suppl), 4721s–4730s (1984). 35. V. M. Kodach, D. J. Faber, J. van Marle, T. G. van Leeuwen, and J. Kalkman, “Determination of the scattering anisotropy with optical coherence tomography,” Opt. Express 19(7), 6131–6140 (2011). 36. H. Ullah, A. Mariampillai, M. Ikram, and I. Vitkin, “Can temporal analysis of optical coherence tomography statistics report on dextrorotatory-glucose levels in blood?” Laser Phys. 21(11), 1962–1971 (2011).


Introduction
The microcirculation comprises the smallest functional vessels of the vascular network, where blood and tissue interact to create an environment necessary for cell and tissue survival [1].It provides the link between cellular processes and system physiology.Microcirculation imaging has been successfully applied in internal, surgical and intensive care medicine where its significance has become particularly apparent in the diagnosis and treatment of septic shock [2,3].Good image contrast and determination of functional parameters are key requirements for microcirculation imagers.Recently, sidestream dark field (SDF) video imaging has evolved from these clinical requirements [4] resulting in a handheld device which captures light reflected from the subsurface tissue bed, focused onto a camera to visualize movement of red blood cells (RBCs) in the capillaries.High quality images can be acquired at locations where the microcirculation is superficial, such as mucosal tissue (e.g.sublingual), conjunctiva and organ surfaces.
Functional microcirculation parameters such as vessel geometry and -density can be straightforwardly obtained from SDF images.The RBCs can be tracked to quantify the blood flow velocities within the vessels semi-automatically [5].However, the maximum measurable flow velocities are limited by the frame rate of the camera (e.g.~1 mm/s at 30 frames per second) [4,5].An alternative technique to quantify microcirculation dynamics is based on laser speckle photography [6], but only after digital cameras became readily available Laser Speckle Contrast Imaging (LSCI) increased in popularity [7][8][9].Speckles are formed when coherent light is randomly reflected (scattered) and constructively or destructively interferes at the detector to form bright respectively dark spots ('speckles').When the scatterers are in motion the speckle pattern fluctuates, resulting in blurring of speckles at the detector within a finite integration time.LSCI links the reduced speckle contrast to the dynamics of the scatterers (flowing RBC's).The speckle contrast depends on the combination of integration time and the decorrelation time of the speckle pattern.However, the quantitative relation between flow velocity and speckle contrast/decorrelation time remains topic of debate [10].On the other hand, the advantage of LSCI is the sensitivity to considerably higher flow velocities than SDF imaging (typically 10 mm/s for LSCI [11]).
Combining conventional SDF imaging with LSCI enhances the dynamic range of microcirculation perfusion measurements, and improves the quantification algorithm of perfusion.LSCI and SDF imaging have been shown to be sensitive to the same microcirculation perfusion dynamics, using the handheld SDF imaging system and a separate commercial LSCI system with different optics, magnification and resolution [12].Here, we present an integrated system where both SDF and speckle contrast images can be obtained using one hand held SDF-LSCI system.This enabled us to acquire true RBC velocities (conventional SDF mode) and speckle contrast values (SDF-LSCI mode) of the same microcirculation region in vivo.We applied a multi exposure acquisition scheme which has been shown to result in reliable estimates for speckle decorrelation times [11].This dualmode approach therefore provided fundamental insight in the quantitative relation between RBC flow velocities and time-integrated speckle contrast decorrelation.This resulted in the development of a novel speckle contrast analysis scheme that distinguishes decorrelation due to blood flow from other additive decorrelation sources, and is a first step towards absolute quantitative laser speckle flowmetry in vivo.

Theory
After interaction of coherent light with tissue, each pixel on the camera receives light that has been scattered from various positions within the sample.The associated path length distribution leads to an interference amplitude that depends on the arrangement of the scattering medium with respect to the pixel, so that an image with randomly varying intensity (speckle pattern) is produced.Moving scatterers will naturally lead to fluctuations in the detected speckle intensity.Quantification of this fluctuation by calculating either temporal or spatial statistics of the speckle pattern provides information on the motion of the scatterers.

Temporal assessment of speckle dynamics
The rate of change of the speckle pattern is quantified through the characteristic decorrelation time(s) of the scattered electric field τ C , which parameterizes the electric field autocorrelation function g 1 (τ) = E(t + τ)E(t)/E(t)E(t).This cannot be measured directly in most experiments so that in practice, the temporal autocorrelation of the detected intensity g 2 (τ) = I(t + τ)I(t)/I(t)I(t) is recorded.Both are related through the Siegert relation: where β M is a measurement-geometry specific constant discussed in Section 2.2.The characteristic time(s) of g 1 (τ), can thus be obtained from measurements of g 2 (τ) after proper assessment of β M .For independent particles with isotropic dynamics, g 1 (τ) is given in [13]: where q = 2k sin(θ/2), with k the wave-vector k = 2π/λ; λ the wavelength; and θ the angle between incident and reflected light.For low-NA detection optics in backscattering, θ ≈180° and q ≈2k.The <Δr 2 (τ)> is the mean-square displacement of the scatterer in time interval τ.
Multiple sources of decorrelation can be present at the same time [10].For example, particles in flow can also undergo Brownian motion due to finite temperatures.If Brownian motion is not affected by the velocity gradient, the total mean square displacement can be separated in a diffuse and convective part.The decorrelation of the scattered fields due to these two independent motions can then be treated separately [14].For Brownian motion, <Δr 2 (τ)> = 6D b τ [13] where D b is the Brownian motion diffusion coefficient.For flows with a Gaussian velocity distribution, Bonner and Nossal [15] showed that <Δr 2 (τ)> = <V 2 >τ 2 , where <V 2 > is the second moment of the velocity distribution.Neglecting other sources of decorrelation (e.g.rotations of the particles), we have: where Both Lorentzian and Gaussian form is encountered in the LSCI literature, often discussing which one is more appropriate [7,9,10,16,17].However, the presence of multiple sources of decorrelation validates the use of a combinational model such as Eq. ( 3).Additionally, a sampled biological tissue volume likely consist of stationary structures and moving blood.The influence of the static component on g 1 (τ) can be taken into account using a modified Siegert relation [11,18,19] introducing a dependence on / ( ) , with I f the detected intensity of the fluctuating scattered light, I s the detected intensity of the light scattered by static components.

Spatial assessment of speckle dynamics
Direct measurement of temporal statistics requires long acquisition times, or a small field of view so that research has focused on spatial assessment of speckle dynamics.This approach relies on quantification of loss in speckle contrast due to temporal averaging of speckle fluctuations during the integration time of a camera.Speckle contrast K is defined as the ratio of the standard deviation (σ i ) to the mean (<I>) of the intensity: .
Contrast K is related to the time integrated g 1 (τ) [20]: ( ) where T is the exposure time of the camera and β M is a factor propagating from Eq. ( 1).
Although τ c can be obtained from a measurement using a single T, it has been shown [11] that more robust estimates of τ C can be obtained using a multi exposure scheme, e.g.measuring K(T) vs. T followed by non-linear curve fitting.The probability density function of the measured intensity is a gamma-PDF [20] with shape parameter M = 1/β M , i.e. the maximum attainable speckle contrast K max = 1/√M.A fully developed speckle pattern gives K max = 1 (β M = 1).In practice K max < 1 due to depolarization; incoherence of the light source; vibrations in the experimental setup; etc [7,21,22].Another cause for reduced K max , is failure to meet the Nyquist sampling criterion i.e. pixel size ≤ ½ speckle size (optical resolution in the camera plane).
An often used approximation is 'Brownian motion only' [6,8,11,23], i.e. the first term in Eq. (3).Using the modified Siegert relation as before, Eq. ( 5) evaluates to: where x = T/τ c and C noise an added noise term for measurement noise.
In case of directional dynamics, the second term of Eq. ( 3) is retained, and speckle contrast is given by [16,18]: .
The model choice has consequences for the retrieved τ c and is a complicating factor in speckle contrast analysis [10].A combinational model, as leading to Eq. ( 3) can be evaluated using e.g.Mathematica, but is not reiterated here because of the arguments in Section 2.4.

Relation between τ C and V
Equation (3) suggests an inverse linear relation between τ C and V for both Brownian motion and directional flow and has been observed in many independent experiments [8,11,24].A common relation is V = λ/(2πτ C ), which reduces to approximately V = 0.1/τ C [μm/ms] for visible laser light [21].This relation differs by 1-2 orders of magnitude [9] from the analysis of Bonner [15] that takes into account the form factor of red blood cells, leading to V = 3.5/τ 1/2 [μm/ms].Also suggested is V = ω/τ C ; where ω is a measure of the point spread function of the system [10].In the camera plane, ω corresponds to the minimal speckle size, and τ C to the time it takes to "replace" a speckle pattern with an uncorrelated pattern.Yoshimura [25] gave an extensive overview of statistical properties of dynamic speckle patterns in different imaging geometries, demonstrating (again) an inverse relation between τ C and V, but proportionality was highly dependent on the exact experimental layout.Imaging a speckle pattern on a camera, yields V = r s /Mτ C , where r s is the speckle size in the object plane and M is the magnification of the imaging system.The relation between τ C and V may thus depend on a combination of scatterer properties and imaging geometry, but remains unsettled.Another factor hampering quantification is multiple scattering by moving objects which reduces τ C because of an apparent increase of <Δr 2 (τ)> in a given time interval τ.

Practical considerations
Our preliminary in vivo measurements showed that 'static tissue' also exhibit a non-negligible speckle decorrelation.We hypothesize that this is due to 1) dynamic processes encountered along the photon path before arriving at the image plane (tissue matrix, other vessels); 2) movement of muscles in the sublingual region (subject) and hand (operator) while holding the Microscan.This 'offset' decorrelation is thus also present in data obtained from vessels in addition to the sought decorrelation due to flowing RBCs.Since, analogous to the derivation leading to Eq. ( 3), the total autocorrelation function is given by the product of the two (or more) autocorrelation functions we write: where flow refers to decorrelation due to blood flow and offset due to photon path dynamics and muscle movements.Likely, g 1,offset is best modeled as random ("Brownian").For simplicity we assume that a Gaussian form describes both processes, such that flow dynamics are characterized by a convoluted time constant: ( ) , where τ c,total can be derived from a multi exposure fit of K measured at the vessel region, while τ c,offset can be measured at a tissue region adjacent to the vessel.LSCI suffers from a trade-off between K-accuracy and spatial resolution, as K must be determined over a sufficient number of pixels.A practical approach [26][27][28] is to determine K in a small spatial region from subsequently acquired images, assuming ergodicity -ensemble statistics equal temporal statistics.Even if ergodicity is violated (e.g. the presence of a static fraction (1-ρ)) correct K-values are obtained as long as ρ is constant between frames, and the time between acquisitions is large compared to τ C .It has been shown that a multi exposure acquisition scheme and fit of Eq. ( 6) or (7) adequately corrects for differences in spatial and temporal sampling [11], and in this report we show similar assurance for spatiotemporal sampling.

Sidestream dark field and speckle microscope system
An SDF-imaging system (Microscan, Microvision Medical, The Netherlands) was modified to allow flexible coupling between light sources for tissue illumination.The SDF-principle is based on a central imaging pathway surrounded by illumination light.Specular reflection is blocked while multiple scattered light illuminates the subsurface vessels and is focused on the camera.Absorption of green light by hemoglobin in red blood cells (RBCs) results in image contrast between tissue and flowing RBCs, an example can be viewed in Media 1.In our system, four large multimode plastic optical fibers (POF ESKA, fiber core 980 µm, NA 0.5) surrounding the central lens tube were used for light delivery.For conventional SDF microcirculation imaging, green light provided by a band pass filtered (550 ± 20 nm) broadband light source (Xenon light, Karl Storz Endoskope, Germany) was coupled into the fibers.For SDF-LSCI imaging, red laser light (632.8nm He/Ne, Spectra Physics, US) was coupled into the fibers using a flip mirror which simultaneously blocked the green light, see Fig. 1.To prevent overexposure a variable neutral density filter was placed in front of the laser.The tip of the Microscan including the four optical fibers is covered with a disposable sterile cap, and can be placed safely on organ and tissue surfaces.The camera plane can be axially translated with respect to the lens system within the tube (not shown, 5x magnification) which results in translation of the focal plane from tissue surface to a maximal depth of 400 µm.Light was detected with a monochrome camera (IEEE 1394, Guppy F-080B, Allied Vision Technologies, Germany), software controlled by self-written scripts (LabVIEW, National Instruments, US).By flipping the mirror between image sequences, conventional 'SDF-absorption mode' and 'LSCI mode' videos could be recorded sequentially.The principal imaging mechanisms are illustrated in the enlargement in Fig. 1: absorption of green light results in contrast between RBCs and tissue, while reflection (and minimal absorption) of red coherent light results in constructive and destructive interference patterns of which the temporal properties reflect the dynamic contrast between RBCs and tissue.

In vitro flow phantom imaging
We designed a sublingual tissue simulating optical phantom based on a protocol by De Bruin et al. [29].Silicone elastomer (Sylgard® 184 Silicone Elastomer DOW/Corning, US) was mixed with titanium dioxide (TiO 2 ) particles (anatase form, Sigma Aldrich, US) to mimic tissue scattering ([TiO 2 ] = 1 mg/ml, corresponding to a reduced scattering coefficient µ' s = 1.5 mm −1 at 632.8 nm) [29].A small glass tube (diameter 300 µm) was embedded in the silicone and connected to a reservoir of Intralipid® (2.5%), a scattering solution of fat particles, subject to a constant pressure to create a stable flow through the channel.The resulting flow velocity was estimated by measuring the throughput volume over an extended period of time.Figure 2 shows an optical coherence tomography image of a cross-section of the flow phantom, together with a schematic drawing.

In vivo sublingual tissue imaging
In vivo sublingual microcirculation SDF and speckle images were obtained with the integrated SDF-LSCI device from a healthy volunteer.The handheld device was gently pressed against sublingual tissue ensuring undisruptive blood flow.The SDF and SDF-LSCI modes were alternated to record data from the same microcirculation region.

Data analysis
The SDF-LSCI device recorded 8-bit monochrome image frames, 1024 x 768 pixels at 30 frames-per-second (except for exposure times > 34 ms).Image and data analysis was performed using self-written software (Labview, National Instruments, USA and Matlab, MathWorks, USA).A prerecorded bias frame (average dark frame) was subtracted from each image frame before analysis, which made the noise-term in Eq. ( 6) and ( 7) negligible.The speckle contrast values recorded at multiple exposure times in the range 0.5 -80 ms were calculated according to Eq. ( 4) using spatial and spatiotemporal sampling (see also Section 2.4), and were fitted to Eq. ( 6) or ( 7) using a nonlinear fitting procedure based on the Levenberg-Marquardt algorithm.For our in vitro flow phantom measurements, β M was apriori estimated from the static part of the phantom, at short integration time T. The same approach is not possible for sublingual tissue measurements, since completely static regions are absent in vivo due to background tissue and/or probe movement.Therefore, for each individual vessel or tissue region β M was estimated from an initial fit with fit parameters β M , ρ and τ c .We note that β M is expected to be different in in vitro and in vivo situations, due to differences in optical properties affecting multiple scattering, depolarization, etc. which are likely to be different in both experiments.In subsequent evaluations of a given phantom or vessel/tissue region, β M was kept fixed and ρ was estimated at T >> τ c , since the first two terms in Eqs. ( 6), ( 7) will be negligible (K 2 = β M (1-ρ) 2 ).Finally, either Eq. ( 6) or ( 7) was fitted to the speckle contrast data with τ c the only fitting parameter.Equation ( 6) was used when no flow, but only Brownian motion was present in phantom experiments; Eq. ( 7) in all other cases.In vivo, determination of τ c corresponding to flow regions and to non-flow regions allowed extraction of the speckle decorrelation time associated with flow only according to Eq. ( 9).Image frames acquired using conventional SDF-mode could be used to measure the blood flow velocities, using commercially available software (Automated Vascular Analysis, AVA3.0,Microvision Medical, The Netherlands).

Validation of SDF-LSCI device
Following Section 2.2 the intensity PDF of a static speckle pattern determines the maximum attainable speckle contrast.Figure 3(a) shows the intensity PDF of the speckle pattern resulting from the static part of the flow phantom imaged with SDF-LSCI (green triangles), and the black line represents a gamma-PDF with M = 2.2 [20,30]; the maximum contrast is K max = 0.67 and β M = 0.45.Spatial integration and under-sampling of speckles can also reduce speckle contrast, therefore we estimate our speckle size using the second order statics of the intensity distribution [31,32], such as the power spectral density (PSD) and the normalized autocovariance of the intensity in the spatial and the frequency domain respectively (Fig. 3(b)  and 3c).The measured PSD (Fig. 3(b)) straightforwardly demonstrates that the speckles are sampled above the Nyquist frequency because there is no signal content at the maximum spatial frequency.More quantitatively, Alexander et al. calculate the average speckle diameter as half the full width of the autocovariance function [30], while Goodman relates the square-root of the 'equivalent area' of the 2D autocovariance function to the width of a speckle [20], and Kirkpatrick et al examine the FWHM of the autocovariance function to estimate the minimum speckle size [31].Figure 3(c) shows the normalized autocovariance of the speckle intensity with an equivalent area of ~6 pixels and a full width of ~5 pixels; resulting in a speckle size of 2.4 respectively 2.5 pixels.We estimated the FWHM of the PSD by fitting a Gaussian and calculated the corresponding FWHM of the (also Gaussian) autocovariance function according to cov.

8ln(2) /
auto PSD FWHM FWHM = = 2.2 pixels.It is thus assured that the speckles are adequately sampled above the Nyquist frequency, with an estimated minimum speckle size ~2.2pixels (10.2 µm) and average speckle size ~2.5 pixels (11.6 µm).Lastly, the temporal behavior of the light source and light delivery by multimode fibers is investigated.The Spearman correlation coefficient of subsequent image frames (T = 1 s), recorded from a static phantom region, showed a strong correlation (r s > 0.5) up to 5 seconds.In view of the expected speckle decorrelation times in vivo this is sufficiently long.

Validation of speckle contrast measurement
Accurate estimation of K benefits from a large number of sampled speckles and thus a large spatial region size for calculation of <I> and σ i .The estimated K-value is reduced with respect to its true value for a small region-size to speckle-size ratio and the uncertainty in K is larger [27,28,32].Large spatial regions result in a reduced spatial resolution.Spatiotemporal sampling can alleviate this resolution loss by (also) averaging over temporal pixels.In microcirculation imaging both spatial and temporal resolution are important, therefore we define a spatiotemporal local region [27,28].
To validate K-estimation we recorded 20 in vitro and in vivo images (T = 1 ms) and a large area representing the center of the flow channel respectively vessel was selected for which constant flow (thus constant K) was assumed.The flow velocity in the channel was 7 mm/s while that in the vessel was <2 mm/s resulting in lower absolute speckle contrast values for the channel compared to the vessel.To determine the smallest local region size that still gives an accurate estimation of K we varied the number of pixels in the spatial (N s × N s ) and temporal (N t ) dimension.For the varying local region sizes (N s × N s × N t pixels) K-values were calculated for the selected channel and vessel area (see Fig. 4(a)).Next, the mean <K> and coefficient of variance (CV = σ K /<K>) were determined, where σ K is the standard deviation in K.For the largest local region size a total of 45 K-values could be obtained, thus, for fair comparison 45 K-values were selected for each local region size to calculate <K> and CV.The global K-value was calculated from all pixels comprising the area representing the channel respectively vessel in the raw speckle image, and is regarded as the 'true' K-value.In Fig. 4(b) the results for <K> in vivo and in vitro are plotted as a dependency on N s for 4 different N t 's, with global K represented by the dashed line.In Fig. 4(c) and 4(d), the obtained CV is plotted for in vitro respectively in vivo K-values.As expected, the estimated <K> converges to the global value with increasing region size, and the CV decreases.Figure 4 also shows that for a scenario with a large fraction of static scatterers, like the flow phantom, the estimated <K> is mostly dependent on the spatial dimension N s , although the CV of <K> depends on both N s and N t .Based on Fig. 4 we defined our local region for subsequent analysis to be 7 x 7 x 20 (N s × N s × N t ) pixels, which means a negligible underestimation of K and a CV below 8%, while keeping the spatial resolution at 7 µm.

SDF-LSCI in vitro: flow phantom
The flow phantom enabled recording of raw speckle images for a large range of Intralipid® flows (in the range V = 0 -20 mm/s) and exposure times (0.5 -80 ms).The raw speckle images were converted to speckle contrast images using a local region of 7 × 7 × 20 (N s × N s × N t ) pixels.In Fig. 5(a) the multi exposure K-values are plotted for five different flow velocities and Brownian motion, as well as for a static (non-channel) region of the phantom.The flow data points are fitted to Eq. ( 7) (Gaussian autocorrelation function), while the Brownian motion data is fitted to Eq. ( 6) (Lorentzian autocorrelation function), <R adj 2 > averaged over all 9 fits = 0.86 ± 0.08.β M was estimated at 0.33, and ρ at 0.6 ± 0.02.To verify that spatiotemporal sampling resulted in a comparable fitting outcome compared to spatial sampling, we also estimated K-values using a spatial local region consisting of an equal amount of pixels (31 × 31 instead of 7 × 7 × 20) and repeated the same fitting procedure.Again, the goodness of fit was high (<R adj 2 > averaged over all 9 fits = 0.94 ± 0.03).The obtained fitting parameter τ c from both sampling procedures is plotted in Fig. 5(b) as a dependency of flow velocity.(as in the flow phantom) diverges the estimate of K when sampled spatially compared to spatiotemporally [18].As can be seen in Fig. 5(b), the absolute τ c values from fitting the spatially and spatiotemporally sampled multi exposure curves are similar (difference in slope linear fit <10%).Thus, the adopted model based on the modified Siegert relation properly corrects for the presence of a static component in the speckle images in the SDF-LSCI geometry, as was shown for other geometries before [11].7) for the flow data, while Eq. ( 6) is fitted to the Brownian motion data.(b)Fit parameter τ c (plotted as 1/τ c ) from multi exposure speckle contrast fits for several flow velocities between 0 -20 µm/ms (mm/s), applying a spatiotemporal (purple squares) and a spatial (orange triangles) local region, consisting of 7 × 7 × 20 and 31 × 31 pixels respectively, and its linear fit (dashed line).

SDF-LSCI in vivo: sublingual microcirculation
The in vivo data set consists of videos of several regions of the sublingual microcirculation acquired in SDF and SDF-LSCI mode sequentially, recorded with multiple exposure times (in the range 0.5 -80 ms) for the SDF-LSCI mode.To illustrate the recording and analysis procedure and resulting processed images, a video is shown in Media 2 (1.6 Mb), Media 3, (7.9 Mb), and Fig. 6. Figure 6 shows four panels each representing the same area of sublingual microcirculation.Figures 6(a) and 6(b) (top left and top right panel in both media files) represent unprocessed conventional SDF and raw speckle images respectively.The corresponding processed K-image at exposure time 10 ms is shown in Fig. 6(c).In Media 3 (large file) raw speckle images and processed SDF-LSCI images are shown for T = 1, 2, 10, 20 and 50 ms consecutively, each as 2 second videos, and Media 2 only for T = 10 ms.Note that to enhance the contrast in Fig. 6(c) the scale bar for speckle contrast K ranges from 0.1 -0.4,while in both media files (bottom left panel) the K-scale bar ranges from 0 -1 to accommodate for K-values from short and long exposure times.To increase the spatial and temporal visualization, the raw SDF-LSCI videos in both media files (top right panel) were processed for each pixel (with local region N s = 7, N t = 20), thus only every 7th K-value in space and every 20th K-value in time are independent.For each exposure time 150 frames were recorded, consequently each pixel represented an average of 7 K-values to perform a multi exposure fit (Eq.( 7)), resulting in a 1/τ c map of the sublingual microcirculation region as shown in Fig. 6(d) and bottom right panel of both media files.
To average the physiological variation in flow, SDF videos of 5 seconds were captured on average 5 times (at least twice) for each vessel before, during and after multi exposure laser speckle imaging.SDF and SDF-LSCI images were recorded for 5 different sublingual regions in total.The minimum blood flow recoded was 0.2 mm/s and the maximum was higher than the reliable measurement range (>1.2 mm/s) using the RBC tracking algorithm in AVA3.0, resulting in reliably measurements for 15 out of 19 analyzed vessels.A multi exposure speckle acquisition scheme (0.5 -80 ms) was performed twice for each vessel, and for each exposure time 150 frames were recorded (≥ 5 sec) again averaging for physiological heterogeneity in flow.In these processed K-images for each vessel at least 25 (in space) x 15 (in time) independent K-values were selected to calculate mean <K>.The standard deviation in K-values was between 1 and 10% with an average of 4%.In Fig. 7(a) three examples of multi exposure K-curves and their corresponding fit result (Eq.( 7)) are shown for three blood vessels, ranging from 0.2 -1.1 mm/s.We note that the maximum K-value of the multi exposure speckle contrast curves in vivo (Fig. 7(a)) is higher than the static data points in Fig. 5a, likely due to differences in β M as discussed in Section 3.4.The goodness of fit was high, <R adj 2 > averaged over all 19 fits = 0.98 ± 0.01.The average value of β M was 0.44 ± 0.06, and ρ 0.85 ± 0.02.
A bar plot of τ c,flow according to Eq. ( 9) versus in vivo blood flow velocity is shown in Fig. 7(b).Due to physiological variance in blood flow velocity we grouped the vessels, where each group consisted of minimally 3 vessels.Subscribing to our hypothesis, the linear correlation between blood flow velocity and 1/τ c was higher for the corrected τ c,flow (r 2 = 0.95) than the uncorrected τ c,total (r 2 = 0.40) for the grouped vessels.The slope and intercept [95% upper CI -lower CI] of the linear fit were for 1/ τ c,flow versus flow: 0.32 [0.19 -0.44] [µm −1 ] (slope) and 0.37 [0.29 -0.44] [ms −1 ] (intercept).In addition, τ c,total (uncorrected vessel τ c ) correlated highly to τ c,offset (adjacent tissue region) (r 2 = 0.99).

Discussion
We validated sidestream dark field laser speckle contrast imaging (SDF-LSCI), for functional microcirculation imaging in vivo based on contrast due to RBC-concentration (absorption of green light, conventional SDF-mode) and on contrast due to RBC-flow (reflection of coherent light, SDF-LSCI mode).SDF-LSCI is sensitive to a large range of flows, from sub-mm/s to effectively unlimited velocities, making it a suitable tool to image low flows in angiogenic microvessels (wound healing, tumor neovascularization), or for high blood flow imaging during hyperthermia [33,34].A linear relationship between 1/τ c and flows up to 12 mm/s (equivalent to τ c > 0.03 ms) was found in vitro, suggesting that SDF-LSCI can be used to quantify changes in the microcirculation during the development of sepsis [3].9).The vessels were grouped in five different velocity ranges, with minimally 3 vessels per group.
1. "a persistent erroneous formula expressing contrast as a function of integrated instantaneous covariance of intensity" -As pointed out by Bandyopadhyay [24], for robust absolute flow measurements the correct derivation is essential and thus applied in this study (Eq.( 5)).
2. "the inappropriate use of the Lorentzian field correlation relationship" -Differences in τ c fitted with Eq. ( 7) or (8) (Gaussian or Lorentzian) were small in vitro (data not shown), which was previously also found in vivo [16].In [16] a multi exposure scheme was used as well, which seems to reduce the impact of Lorentzian or Gaussian autocorrelation model on the estimate of τ c .However, since the microcirculatory flows are generally low where the error induced by incorrect autocorrelation model is more significant (small T/τ c ratio) [17], we applied the Gaussian model appropriate for directional flow for the subsequent in vivo analysis.
3. "a tendency to operate in the long-exposure asymptotic regime and the subsequent lack of sensitivity" -Fig.5(a) shows that high flow velocities in the phantom experiment allowed only fitting the tail of the multi exposure curve.The reduced linearity of 1/ τ c vs. V for higher velocities (Fig. 5(b)) could thus be caused by a large T/τ c ratio in the applied exposure range.Increasing the intensity to allow for shorter T is expected to result in an increased linear regime.For our in vivo data at T = 0.5 ms, T/τ c,total is well below 1 and the speckle contrast is close to K max (high sensitivity [10]).
4. "a common assumption that the requirements of quasi-elastic light scatter (QLS) and LSCA measurements are the same" -In the optimal QLS measurement, the speckle size matches the pixel size [7,23].However, in LSCI the speckle contrast (thus <I> and σ i ) is highest when speckles are truthfully imaged.Figure 3 shows that for the SDF-LSCI instrument we assure that the speckles are sampled above the Nyquist frequency, with an average speckle size of ~2.5 pixels.pixels/speckle [32], but that any reduction in K is properly corrected for by a 'linear system factor', i.e. β M .
5. "the oft-cited, nonphysical association between the decorrelation time and its associated flow velocity" -The key feature of the modified SDF-LSCI tool is the accessibility to speckle decorrelation times and actual blood flow velocities simultaneously.Having preliminary quantitative in vivo 1/τ c and absolute flow data in hand, investigating this relationship was plausible.The aforementioned association refers to τ c = λ/2πV ≈0.1/V [ms].Interestingly, the (uncorrected) τ c,total found in vivo was between 0.8 -3 ms which according to this relation corresponds to a flow of 0.03 -0.13 mm/s.Without knowledge of actual RBC velocities this seems a plausible range for microvessels, however the results of Section 4.4 show that using the uncorrected τ c,total there is low correlation with flow velocity (r 2 = 0.40).
After correction for offset decorrelation an improved relation between flow and τ c could be established (r 2 = 0.95), Fig. , significantly deviating from this model.The decrease in proportionality constant may be due a larger influence of photons that are multiply scattered from moving scatters, caused by a larger channel (Ø = 300 µm vs. Ø = 30 ± 16 µm in vivo) and smaller scatterers (<700 nm fat particles vs. >10 × larger RBCs [35]).The third suggested model [15], taking into account the size of scatterers yields τ 1/2 ≈3.5/V [ms] for RBCs (close to our values); our proportionality constant in vitro corresponds to a particle of Ø = 0.5 µm, in excellent agreement with our recent determination of the effective scattering size for Intralipid® particles [35].The different physical principles underlying the proposed models confound quantitative flow measurements, further complicated by possible multiple scattering and heterogenic flow profiles.However, we believe our approach addresses all of the five issues discussed above.Most importantly, the dual-mode operating SDF-speckle imager is a valuable method to increase the knowledge on the relationship between τ c and V in vivo.

Comments on hypothesis
Our proposed separation of τ c,offset and τ c,flow for in vivo measurements in Section 2.4 improved the correlation between actual flow velocity and 1/τ c,flow , but needs further verification and explanation of the remaining intercept of 1/τ c at zero flow (Fig. 7(b)).This model was applied to our in vitro data, Fig. 5(b).Here, the intercept of the linear fit of 1/τ c versus V was 1.8 ms −1 , with 95% confidence range (−0.7 − 4.2).This result nearly matches the value of 1/τ c found for Brownian motion of Intralipid® at V = 0 in the flow phantom channel: 1/τ c = 0.4 (0.32 -0.53) ms −1 .Following Eq. ( 3), we interpret the dynamics in flowing Intralipid® as a combination of Brownian motion and directional flow, where the total autocorrelation function is given by the product of the individual functions [10].Applying our analysis of Eq. ( 9) (all processes have Gaussian autocorrelation functions), we interpret decorrelation due to Brownian motion as an offset, e.g.τ c,offset = τ c,Brown .As there were no other sources of decorrelation in our phantom setup, this analysis indeed yields 1/τ c,flow = 0 at 0 mm/s.We again note that the assumption of two Gaussian autocorrelation functions is not strictly correct (since Brownian motion is associated with Lorentzian autocorrelation), but to gain mathematical simplicity the resulting small deviation is acceptable.
This also suggests that the remaining intercept in the in vivo data, Fig. 7(b) (τ c,flow = 3 ms), is partially due to the Brownian motion of RBCs in the vessels that is not corrected for using in vivo τ c,offset (next to the vessel).A τ c,Brown,RBC in the order of 10 ms was measured for rat RBCs in Brownian motion, which have similar size to human RBCs [36].

Conclusion
We studied the combination of two non-invasive in vivo flowmetry techniques, calibrated flowmetry (SDF microscopy) and laser speckle flowmetry (LSCI).This enabled us to compare in vivo microcirculation blood flow to speckle decorrelation characteristics.We found that decorrelation sources other than blood flow confound the relationship between flow and τ c .We proposed an analysis scheme that reduces the contribution of additive decorrelation sources to advance towards absolute quantitative laser speckle flowmetry.The preliminary in vivo quantitative relationship of τ c with blood flow velocity is promising in view of existing models, but warrants further study.We thoroughly validated the integrated instrument to ensure valid speckle sampling and speckle contrast calculations.Using the exposure range of 0.5 -80 ms, the absolute value of fitting parameter τ c was robust to sampling method (spatial or spatiotemporal sampling of K) and 1/τ c was found to scale linearly with flow velocities up to 12 mm/s in vitro.We have shown that SDF-LSCI has the potential to become a quantitative, automatic, non-invasive and affordable microcirculation monitoring tool.

Fig. 1 .
Fig. 1.Modified SDF-imaging system for non-invasive imaging of subsurface microcirculation in dual mode (SDF and SDF-LSCI).Broadband green light (SDF) is highly absorbed by flowing RBCs, resulting in contrast between vessels and tissue.Red coherent light (SDF-LSCI) is scattered by tissue and moving RBCs, resulting in contrasting speckle dynamics between vessels and surrounding tissue.The optical set-up for sequential coupling of SDF light and SDF-LSCI light into the four optical illumination fibers involved 3 lenses (L1 = L2, f = 32 mm, L3, f = 10mm); two filters (F1, band pass interference filter, 550 ± 20 nm and F2, variable neutral density filter, OD 0.04 -2.00); and a flip mirror M to swop between modes.A 5x magnifying lens system in the lens tube focuses the subsurface microcirculation image onto a CCD camera (a video can be viewed in Media 1).Details of imaging modes are illustrated in the inset (green arrows: SDF, absorption; red arrows: SDF-LSCI, speckle).

Fig. 2 .
Fig. 2. Flow phantom design.Optical coherence tomography image of cross-section of flow phantom (left), and schematic drawing of flow phantom (right).

Fig. 3 .
Fig. 3. First and second order statistics of the speckle intensity pattern recorded using the static part of the flow phantom and an exposure time T = 1 ms.(a) Intensity probability density function (PDF) of static speckle pattern (green triangles) and plot of gamma-PDF with M = 2.2 (black line) (b) Power spectral density (PSD) of the static speckle intensity pattern showing sampling above the Nyquist frequency.FWHM PSD ≈2.5 pixel −1 .(c) Normalized autocovariance of static speckle intensity pattern, full width ≈5 pixels, equivalent area ≈6 pixels, FWHM autocov.≈2.2 pixel.

Fig. 4 .
Fig. 4. Validation of local region size for accurate measurement of speckle contrast K. (a) Schematic representation of local region size, where N s is the number of pixels in the spatial dimension and N t in the temporal dimension (total number of pixels: N s × N s × N t ) for which K is calculated.(b) Mean <K>, as calculated from 45 K-values selected in the processed speckle contrast image from the channel (in vitro) or vessel (in vivo) area.The dashed line represents the global K value.(c)Coefficient of variance (CV = σ K /<K>), where σ K is the standard deviation in the 45 K-values.Both <K> and CV are plotted versus N s , for 4 different temporal dimensions N t .All images were recorded at exposure time 1 ms.

Fig. 5 .
Fig. 5. Flow phantom speckle decorrelation results.(a) Multi exposure speckle contrast values (data points) and corresponding fit (lines) for several different flow velocities.The K values (N s = 7, N t = 20) are fitted to Eq. (7) for the flow data, while Eq.(6) is fitted to the Brownian motion data.(b)Fit parameter τ c (plotted as 1/τ c ) from multi exposure speckle contrast fits for several flow velocities between 0 -20 µm/ms (mm/s), applying a spatiotemporal (purple squares) and a spatial (orange triangles) local region, consisting of 7 × 7 × 20 and 31 × 31 pixels respectively, and its linear fit (dashed line).

Fig. 6 .
Fig. 6.In vivo SDF-LSCI recording and analysis procedure.(a) Typical SDF image of sublingual microcirculation.(b) Raw speckle image and (c) processed K-image (N s = 7, N t = 20) of same region as in (a), recorded at T = 10 ms.For each pixel and exposure time, K was estimated to enable a pixel wise multi exposure curve fit (Eq.(7)), resulting in a 1/τ c map of the same region as in (a), shown in (d).In Media 2 and Media 3 the four panels are presented as a video, showing flowing RBC's (top left), raw speckle (top right) and speckle contrast (bottom left) images, and a still frame (bottom right) of the corresponding 1/τ c map. Media 2 (1.6 Mb) represents SDF-LSCI at T = 10 ms, while Media 3 (7.9Mb) represents SDF-LSCI at T = 1, 2, 10, 20 and 50 ms consecutively.

Fig. 7 .
Fig. 7. a.In vivo SDF-LSCI results.(a) Multi exposure K-values (data points) and corresponding fit (lines) for several different flow velocities in vivo (upper curves).For comparison two in vitro multi exposure K-curves are also plotted (lower curves).Equation (7) is used as fitting model.(b) Bar plot of fit parameter τ c,flow (represented as 1/τ c,flow ) of 19 in vivo blood vessels, corrected according to Eq. (9).The vessels were grouped in five different velocity ranges, with minimally 3 vessels per group.