Non-Invasive Blood Flow Speed Measurement Using Optics

Non-invasive measurement of the arterial blood speed gives important health information such as cardio output and blood supplies to vital organs. The magnitude and change in arterial blood speed are key indicators of the health conditions and development and progression of diseases. We demonstrated a simple technique to directly measure the blood flow speed in main arteries based on the diffused light model. The concept is demonstrated with a phantom that uses intralipid hydrogel to model the biological tissue and an embedded glass tube with flowing human blood to model the blood vessel. The correlation function of the measured photocurrent was used to find the electrical field correlation function via the Siegert relation. We have shown that the characteristic decorrelation rate (i.e., the inverse of the decoherent time) is linearly proportional to the blood speed and independent of the tube diameter. This striking property can be explained by an approximate analytic solution for the diffused light equation in the regime where the convective flow is the dominating factor for decorrelation. As a result, we have demonstrated a non-invasive method of measuring arterial blood speed without any prior knowledge or assumption about the geometric or mechanic properties of the blood vessels.


Introduction
Oxygenated red blood cells in the blood flow deliver essential nutrients and oxygen to organs and limbs to maintain their homeostatic conditions and proper functions. Oximeter detects the oxygen saturation level of red blood cells by using wavelength dependent absorption for oxygenated and deoxygenated blood [1]. Such measurements can be done with wearable devices such as smart watches. On the other hand, there has been no convenient and non-invasive method to measure the amount of blood supply or blood speed in major arteries.
Reduced blood supply is usually a sign of diseases and medical conditions. Dehydration, thromboses formation, and other physiological and pathological conditions can cause short-term or long-term changes in the blood speed at a given location of the artery. The information helps early diagnosis and monitoring of diseases with impaired blood flow [2][3][4]. Therefore, it is desirable to directly measure the traveling velocity of red blood cells at well-defined locations such as the carotidal artery to the head, femoral artery to the lower limb, brachial artery to the upper arm, spinal arteries to the spinal cord, renal arteries to the kidneys, hepatic artery to the liver, and pulmonary artery to the lungs, etc. The arterial blood speed at the well-defined positions provide direct and unambiguous information about blood supply to the sites of health concerns.
Doppler frequency shift of ultrasonic waves offers a non-invasive measurement of the blood flow speed at the well-defined position. However, acoustic measurements require close physical contacts of the ultrasound transducers with the tissue to allow efficient coupling of the acoustic wave, and gel is needed to ensure the quality and reliability of the contact. Furthermore, the measured signal is highly dependent on the angle between the probe and the blood flow direction. The most convenient and ergonomic direction is to have the transducer perpendicular to the blood vessel. However, this arrangement will produce no Doppler shift signal, thus making the measurement less convenient in a homecare or nursing home setting for self-administered operation [5]. In this regard, an optical method is potentially more desirable because a light wave can enter the tissue via free space coupling and once entering the tissue, as we will discuss later, it travels diffusively. In addition, the optical hardware including semiconductor light sources and detectors are more compact than acoustic transmitters and receivers.
People have used dynamic scattering of laser light by red blood cells to measure the speed for blood vessels near the tissue surface [6,7]. Similar to the Doppler effect of sound propagation, the frequency of the scattered light by moving particles is also shifted. Although the operation principle seems to be straightforward, the physical model of the "optical Doppler device" is based on a key assumption that the light is only scattered once by a moving object (e.g., a traveling red blood cell) and the rest of the scattering is produced by quasi static objects such as skin and fat. This assumption may be reasonable for probing near surface tissue blood flow, but not for major arteries where around 40% of the volume of the blood vessel with a diameter of a few millimeters is occupied by traveling red blood cells. Under the approximation of single light scattering by moving particles, the detected scattered light will display a frequency shift proportional to the velocity of the particle and the amount of frequency shift can be detected using the homodyne coherent detection technique [8]. However, the above technique cannot be applied to multiple scatterings by a large number of moving particles. In such situations, the frequency spectrum of the light wave is broadened, making the interpretation of signal difficult. Because of this limit, the current technique of optical Doppler shift has been limited to measuring the blood flow in microflow channels near the skin although the blood supply by major arteries produces more valuable information for health and disease conditions. In addition, the laser Doppler flowmetry can only measure relative speed [9].
To measure the blood flow speed in major arteries, we extend the technique of optical scattering method by taking into account the effect of multiple scattering by red blood cells. Instead of measuring the Doppler shift in the optical frequency, we measure the decorrelation time of the reflected light. The decorrelation time and spectral broadening can be easily related by the Wiener-Khinchin Theorem, we find that the key information can be better extracted and understood in the time domain. The key results of our work are the following: (a) we have shown experimentally and theoretically that one can measure the arterial blood flow velocity using a simple optical setup comprised of a single semiconductor laser and a single detector, both coupled to a multimode fiber bundle, and (b) we discovered that for main arteries where the diffused light model applies, the decorrelation time is inversely proportional to the RBC speed and we can measure its value without any prior knowledge about the anatomy and tissue mechanic properties of the artery.
Once entering the tissue, light has a very short scattering mean free path and the light becomes diffusive [10]. Modeling of light propagation in this regime led to the development of the oximeter back in 1970s [10] and many researchers have utilized the diffused light to construct images through highly diffusive biological media, such as reconstructed breast cancer images [11]. The standard configuration for diffused light correlation spectroscopy includes a point source and a point detector, and the measured correlation function depends on the relative position between the source and the detector [12][13][14]. The method of diffused light correlation spectroscopy has been used in blood perfusion measurement [15], including brain circulations [16][17][18]. The setup uses a single-mode fiber coupled detector to allow single-mode transmission of two orthogonal light polarizations [19]. The single spatial mode yields the best signal-to-noise ratio since it detects one speckle. However, this setup produces very low light intensity and requires a single photon detector (e.g., a single-photon avalanche detector (SPAD) or a photomultiplier tube (PMT)), making the setup quite sophisticated and subject to interference and stray light.
To make the system robust and easy to operate, in this paper we use multimode fiber and a regular photodetector to detect temporal fluctuations of the light intensity. In the following, we first present our experimental setup and measurement results, and then describe the physical model to show how the measured results can be used to find the red blood cell velocity. We will show that in a simple and robust setup that is relatively insensitive to the optical alignment between the fiber bundle and the tissue, where we can reliably measure the absolute value of the traveling speed of red blood cells in a blood vessel that has a millimeter diameter embedded in a layer of tissue. The results offer a promising solution for non-invasive measurement of the blood supply into organs and limbs.

Experimental Setup
The experimental setup and measured results are first presented here, followed by the theoretical analysis in the later sections.

Tissue Phantom Fabrication
We used intralipid as the scattering agent to mimic the tissue by making its optical properties of tissue to make the scattering coefficient µ s and anisotropy factor g similar to the the values of real tissue. Because of its similar optical properties to the bilipid membrane of cells, intralipid is commonly used to simulate tissue scattering [20][21][22]. In addition, the absorption coefficient of intralipid is low and its refractive index is close to that of soft tissue [23]. The tissue phantom is created using intralipid (Sigma-Aldrich 20%) in gelatin gel. One percent concentration of intralipid was chosen to achieve a scattering coefficient of around 10 cm −1 at 784 nm [20]. The phantom was initially prepared at 90 • C, and then poured into a mold where a glass tube was pre-inserted at a depth of 5 mm from the phantom surface (measured from the center of the tube). The liquid phantom was immediately placed into a freezer at −18 • C for 30 min for rapid solidification to prevent sedimentation of intralipid to achieve a uniform scattering property. Then the sample was placed in a refrigerator at 6 • C for 30 min to further solidify. The subsequent experiment was usually completed within 30 min after the phantom fabrication to prevent evaporation induced dehydration of the phantom, which could reduce the surface height. Glass tubes of three different inner diameters were used to mimic the arteries. Tygon tubing was connected to both ends of the glass tube, with one end connected to a 10 mL syringe mounted on a syringe pump and the other end connected to a reservoir (See Figure 1). Figure 1. Experimental Setup. The laser light illuminates the tissue phantom surface at a 45-degree angle and the scattered diffused light is collected through the reflection probe fiber bundle which is coupled to a photoreceiver. The glass tube that emulates a blood vessel is embedded in the tissue phantom 5 mm below the surface.

Blood Preparation
12 mL of human whole blood collected on the same day of the experiment was purchased from the San Diego Blood Bank. EDTA was added to the blood sample to prevent coagulation.

Optical and Electronics Setup
We used an off-the-shelf multicore reflection fiber bundle probe (Model:RP21 from Thorlabs) that consists of one center core surrounded by 6 peripheral cores to couple the input and reflected light. The center core was used to deliver the laser light to the blood vessel and the 6 peripheral cores that are merged into a single output were coupled to the photodetector (see Figure 1). The fiber bundle was mounted about 45 degrees relative to the surface of the tissue phantom to prevent specular refection from the tissue surface. The fiber bundle was coarsely aligned to the position of the glass tube in the tissue phantom. The laser used in the experiment was an off-the-shelf fiber coupled laser diode at 785 nm, biased slightly above its threshold current to produce an output of 5 mW. The output of the detection fiber was coupled to a silicon APD with an integrated transimpedance amplifier (Thorlabs APD410A) to achieve a transimpedance gain of 500 kV/A at 10 MHz bandwidth. The APD multiplication gain was set to be the lowest (M ∼ 10). The output signal from the APD photoreceiver entered a data acquisition board (Advantech PCI-E 1840) at a sampling rate of 20 Msps, yielding a Nyquist bandwidth of 10 MHz. The total data acquisition time for the measurement was 25 s. We treated each 5 ms duration as one recording section, so 25 s of measurement produced 5000 sections for analysis and noise cancellation. In all the data presented next, we used the average data over 25 s for the mean and over 1 s to determine the variations shown in error bars. We tested blood flows with three tube inner diameters: 0.93 mm, 1.14 mm, and 1.68 mm. The tube diameters were measured by a digital microscope to assure high accuracy. For each tube diameter, we flow the blood at different speeds.

Experimental Results
We have used the setup in Figure 1 to measure the photocurrent correlation defined as is the photocurrent at time"t" and is the ensemble average. τ is the time delay between the instantaneous photocurrents, and is the variable for the correlation function.
One important relation is to convert the photocurrent correlation into the normalized electrical field correlation function g 1 represented by Equation (1): Using the Siegert Relation [24], we can obtain a relation between the magnitude of g 1 (τ) and the photocurrent correlation as shown in Equation (2): where i(t) is the measured photocurrent, N is the number of spatial modes coupled into the detection fiber bundle and collected by the detector, δ(t) is a delta function, and e is the electron charge. The last term in Equation (2) represents the shot noise. For a 25 s measurement that can be divided into 5000 5 ms sections, the ensemble average is the average of these 5000 sections. With the 5000 ensemble averaging, we achieved a decent signal-to-noise ratio and a reasonable data acquisition time. It has been shown that there is significant SNR improvements in detection of phantom optical property and tissue blood oxygenation by using multimode fibers [25]. The magnitude of electrical field correlation g 1 under different flow speeds of blood is shown in semi-log plots in Figures 2a, 3a and 4a, where each figure shows measurements performed with a different tube diameter. The x-axis of each plot is the logarithmic of τ in the correlation function. From these plots we made two interesting observations: (a) At relatively low blood flow speed, the curve shows a characteristic of logistic function; (b) at high blood speed, the curve behaves like a superposition of two logistic functions. These characteristics become more apparent if we take the derivative of g 1 with ln(τ), showing that ln(T F ) occurs at the inflection point where the magnitude of slope reaches the maximum.
We will discuss the physics of such characteristics in a later section. By examining the interesting features of the measured electrical field correlation g 1 , we can extract the ln(T F ). In the low blood speed regime where the|g 1 | plot shows two superimposed logistic functions, we chose the inflection point in the regime of smaller ln(τ) values. We will elucidate the reasons for such a choice in the next section. Essentially, each of the two superimposed logistic functions represents a corresponding regime of light scattering mechanisms. Figures 2b, 3b and 4b show the plot of 1/T F versus the blood flow velocity with different tube diameters. Amazingly, we obtain a simple linear relation between 1/T F and the blood speed. More interestingly, we have found that the 1/T F versus speed curve is independent of the tube diameter. As shown in Figure 5, the three curves of 1/T F versus speed measured from different tube diameters completely overlap and can be represented by one simple relation independent of the tube diameter. This result suggests that our measurement setup can potentially measure the blood speed in different artery sizes without having to know the exact dimension of the blood vessel diameter. This discovery is very important in practical applications because it shows that from the g 1 (τ) curve, which can be obtained from the correlation of the photocurrent, we can obtain the speed of the blood directly for different arteries at a given position without any prior knowledge of the anatomy of the blood vessel. The discovery of this important relation requires a sound physical foundation, to rule out the possibility of being simply coincidental. The physical model and mathematical analysis will be discussed next.

Diffused Light Equation for Moving Scatters
In order to understand the underlining physics of the relation between the inverse of characteristic decorrelation time, 1/T F , and the flow speed of the blood, we describe the physical model and mathematical formulation for our experiment in this section. We will outline the key steps and, through approximations, produce an analytical relation between the blood speed and the decorrelation time. The detailed numerical results for more general analyses without approximations will be presented in a separate paper.
People usually take two approaches to model light propagation in a strongly scattered biological medium. In one approach, we start with the wave equation and introduce the scattering and absorption characteristics along the optical path. Despite its rigor in the mathematical formulation, to make the result useful, approximations have to be made to make the problem solvable. Twersky's theory and Dyson's equation fall into this category [10]. In another approach, we can formulate the problem in a transport equation that deals with photon energy transport. These two methods eventually give rise to the same result. In this paper we take the second approach because of its relative simplicity.
The governing equation is a radiative transfer equation (RTE) in the following form: in which I is the light specific intensity with a unit Wm −2 sr −1 , ρ is the particle concentration, σ s is the scattering cross section, σ t is the total scattering cross section which is the summation of scattering cross section and absorption cross section. p(ŝ,ŝ ) is the normalized differential scattering cross section which is sometimes called the phase function and it is unitless, and S( r,ŝ) is the source intensity which has a unit of Wm −3 sr −1 . Figure 6 illustrates the orientations and coordinates for an infinitesimal scattering volume. Equation (3) contains 5 coordinates: x, y, z, θ, φ where x, y, z define the position vector r of the light intensity, and θ, φ represent the beam's propagation direction. The unit vectorsŝ andŝ represent the direction of the incident light and the propagation direction after scattering. Solving Equation (3) would produce the full solution for the light scattering problem for any scatter concentration. However, Equation (3) is very complicated to solve. Fortunately, in most biological samples, the scatter density is so high that photons quickly lose the memory of their path histories after multiple scatterings and it can be justified to assume the light intensity depends on its position (x, y, z) with a slight flux flow in the direction of propagation (θ, φ). This would lead to a simplification of the problem, and the average intensity could be described by the diffusion equation which is only dependent on the position (x, y, z). Mathematically, it means that we can expand the light intensity in spherical harmonics by keeping the zero-order and first-order term.
where the first term is the average intensity and the second term is the small photon flux in the direction of propagation. Applying the diffusion approximation to the RTE, we can obtain a steady state diffusion Equation (5) to model light propagation in a diffusive medium when the scatters are not moving. As light needs to be scattered multiple times to become diffusive, an important requirement for the diffusion approximation is that the cross section for light scattering is much stronger than the cross section for light absorption. A problem would arise when dealing with light intensity close to the boundary since the light is highly directional at the boundary, violating the diffusive condition. Different boundary conditions have been explored to resolve this problem, including adoption of an extended boundary condition which uses Taylor expansion to convert the Robin boundary condition into a Dirichlet type boundary condition to simplify the solution for Equation (5) [26].
where D = 1/3µ s is the photon diffusivity with a unit of meter and µ s is the reduced scattering coefficient. µ a is the absorption coefficient with a unit m −1 . S( r) is the source (unit of Wm −3 ) that depends only on the location but not on the propagation direction, different from S( r,ŝ) in Equation (3). When the scatters exhibit motions, the scattered intensity would include time as a parameter. For statistical optics, it is natural to use a field correlation function to capture this dynamic process. Here, we are still talking about the steady state response of the system so that the time dependence can be represented as a parameter in the differential equation.
We start with the case of a single scattering event by a moving object. The normalized single scattering function can be written as: where g s 1 is the single scattering normalized correlation function, E is the scalar electrical field of light, q is the photon momentum transfer by each scattering, ∆r 2 i (τ) is the RMS displacement of the scatters in a duration of τ. The above equation describes the electrical field correlation function due to single scattering. This single scattering function can be incorporated into the radiative transfer equation, yielding the so-called correlation transfer equation which is the dynamic counter part of the static radiative transfer equation [27]: Similarly, the diffusion approximation can also be applied to obtain the correlation diffusion equation which is again the counterpart of the static diffusion equation described previously. Equation (8) is the correlation diffusion equation. The diffusion equation can then be solved with appropriate boundary conditions to be compared with experimental results.
where n is the refactive index of scattering media, k 0 is the wave vector in vacuum, µ s is the reduced scattering coefficient which is the scattering coefficient modified due to scattering anisotropy. For isotropic scattering, the scattering coefficient would be the same reduced scattering coefficient.

Analytic Relation between Decorrelation Time and RBC Flow Speed
We assume a plane wave illumination at the surface. Despite the simplification, the approximation can yield a satisfactory result in good agreement with the experiment. To simplify the practical application, we adopted an experimental design where a single detector covers an area defined by the numerical aperture of 6 hexagonally arranged multimode fiber surrounding the light source (see Figure 1). Assuming the blood vessel is cylindrically symmetric, we used a 2D model, as shown in Figure 7.
We can reduce Equation (8) into the following diffusion equation, in which k = jW(τ) Since the scattering parameters are different inside and outside the blood vessel, two sets of parameters are used. The parameters inside the blood vessel are represented using subscript 1, i.e., k 1 = jW 1 (τ). The parameters outside the blood vessel are represented using subscript 0, i.e., k 0 = jW 0 (τ).
To solve Equation (9), we need to first find the mean squared displacement ∆r 2 (τ) . The mean square displacement for RBCs includes Brownian, shear induced diffusion, and convective motion. We can represent the position variation of RBCs over a given time interval as [28]: where the first term in Equation (13) is caused by Brownian and shear induced diffusions and the second term by convection. For RBCs, the diffusion by Brownian motions is much smaller than the shear induced diffusion [29]. The flow of RBCs inside arteries can be modeled by a laminar flow and the diffusion coefficient can be represented as [30] D α = α s ∂v RBC ∂r α s is a parameter describing the interaction strength among blood cells due to shear and its value has been measured experimentally [31][32][33][34][35][36][37]. The radial dependent shear rate is usually replaced by its average value as multiple scattering will lead to an ensemble averaging across all radial locations. For tissue blood perfusion, the slow flow speed and small blood vessel diameter make the diffusive motion the dominant effect compared to the convective motion [38,39]. The situation is reversed, however, for main arteries and our phantom model where the vessel inner diameter is of millimeter size and the blood flow speed is several centimeters per second [40]. In such cases, the convective motion becomes the dominant effect. Hence, in our analysis, we have ignored the diffusion contribution and included only the convective flow contribution in Equation (13).
We can then write the general solution of the correlation function in Equation (9) as follows, The first term is the inhomogeneous solution and the second term is the homogeneous solution for Equation (9) in the absence of the source. The diffusion equation under a given boundary condition can be solved in polar coordinates and the method of separation of variables. For an approximate analytic solution, we kept only the zeroth-order term since all higher-order terms are insignificant compared with the zeroth order term. A continuity boundary condition is applied at the cylinder interface between the blood vessel and the surrounding tissue. The air-tissue surface is ignored for simplicity in the solution. After solving the differential equation with some approximations based on the numerical value of actual tissue and blood cell scattering properties, the normalized correlation function g 1 can be written in the form of Equation (16). The detailed derivation can be found in the supplementary materials.
where V is the flow velocity of blood, n is the refractive index of blood, k 0 is the wavevector of 784 nm light in vacuum, µ in s is the reduced scattering coefficient of blood (m −1 ), and µ in a is the absorption coefficient of blood (m −1 ). If we define T F as the characteristic decorrelation time, and 1/T F as the characteristic decorrelation rate, we obtain the key relation that the characteristic decorrelation rate is only a function of flow speed and blood optical parameters, as shown in Equation (17). In Equation (16), g 1 (r, a, θ, ∞) is the asymptotic value of g 1 when τ approaches infinity, and its value is a function of detector position and blood vessel diameter. In a separate study that solves Equation (9) numerically, we have validated the approximations that led to the analytic solution for g 1 in Equation (16) and, most importantly, the linear relation between 1/T F and the blood flow velocity.
From the approximated relation in Equation (17), we find that the characteristic decorrelation rate is proportional to the flow speed. The proportional constant has a square root dependence on the ratio between the scattering coefficient and the absorption coefficient. This is intuitive since scattering events cause dephasing, and light absorption would terminate the scattering process.  . Geometry used for theoretical calculation. The circle represents the cross section of blood vessel with parameters: D in , W 1 , k 1 , µ in a , µ in s ; The area outside the circle represents the static scattering media, i.e., tissue with parameters D out , W 0 , k 0 , µ out a , µ out s . We adopt cylindrical coordinates and define the center of the vessel as the origin.

Discussion
The measured characteristic decorrelation rate (1/T F ) is obtained by calculating the first inflection point of the curves in Figures 2-4. The behavior of the curves at a longer time is more complicated, and can be explained by other slower decorrelation processes than convection, such as shear induced diffusion, thus yielding two superimposed logistic functions, as shown in Figures 2-4. To use the model described in the previous section, we need to focus on the short time behavior driven by convection. Figure 8 shows the calculated and measured characteristic decorrelation rate under different flow speeds and tube diameters. An excellent agreement between theory and experiment was achieved, confirming that in the regime where convective flow is the dominating factor for decorrelation, the characteristic decorrelation rate is proportional to the blood speed and independent of the vessel diameter.  (17) and experimentally measured decorrelation rates under different blood flow speeds and tube diameters. The scattered data set with the error bar represents the 95% confidence interval over 25 measurements. Each measurement took 1 s. The following parameters are used in the calculations: n = 1.36, µ in s = 1600 m −1 and µ in a = 1000 m −1 for deoxygenated blood [41][42][43][44].

Conclusions
In this work, we demonstrated a simple technique to directly measure the blood flow speed in main arteries based on the diffused light model. The device uses a single fiber bundle, a diode laser, and a photoreceiver. The concept is demonstrated with a phantom that uses intralipid hydrogel to model the biological tissue and an embedded glass tube with flowing human blood to model the blood vessel. The correlation function of the measured photocurrent was used to find the electrical field correlation function via the Siegert Relation. Interestingly, the measured electric field correlation function g 1 (τ) and ln(T F ) shows a relation similar to the logistic function, allowing us to define the ln(τ), occurring at the first inflection point of g 1 (τ). Surprisingly, the value 1/T F , which we call characteristic decorrelation rate, is found to be linearly proportional to the blood speed and is independent of the diameter of the tube diameter over the size and speed ranges for major arteries. This striking property can be explained by an approximate analytic solution for the diffused light equation in the regime where the convective flow would dominate the decorrelation. This discovery is highly significant because, for the first time, we can use a simple device to directly measure the blood speed in major arteries without any prior knowledge or assumption about the geometry or mechanical properties of the blood vessels. A non-invasive method of measuring arterial blood speed produces important information about health conditions. Although the device and setup can be further optimized (e.g., adding different light wavelengths, different depth of vessel in tissue phantom) and the physical model can be expanded to acquire more information from the measurements, the work has paved its way to a new promising modality for measurements of blood supplies to vital organs.