## 1.

## Introduction

Laser Doppler flowmetry (LDF) is a noninvasive technique for estimating the microcirculatory blood flow. The technique utilizes the small frequency shift (i.e., Doppler shift) that arises when light is scattered by moving particles, e.g., red blood cells (RBCs). In conventional LDF (cLDF), a perfusion estimate that increases with both the tissue fraction of red blood cells and RBC flow speed is achieved by calculating the first moment of the detected Doppler power spectrum.^{1, 2} Similarly, the tissue fraction of moving RBCs may be estimated by calculating the zero-order moment of the Doppler power spectrum. Normalizing the perfusion by the RBC tissue fraction gives an estimate of the average blood flow speed.

Being able to monitor the microcirculation in a noninvasive manner makes LDF unique. However, the conventional perfusion estimate is highly affected by the tissue optical properties and the structure of the tissue at the measurement site.^{3} Hence, the cLDF measure is given in a nonabsolute perfusion unit, which complicates inter- and intra-individual comparisons. It is also a measure that responds nonlinearly to the RBC tissue fraction, even at low RBC tissue fractions.^{4} Furthermore, the measurement volume is unknown and affected not only by the optical properties, tissue structure, and measurement setup, but also by the RBC tissue fraction itself.^{5} The cLDF RBC tissue fraction estimate also becomes uncertain due to large fluctuations at low frequencies in the Doppler power spectrum, which also makes the cLDF speed estimate instable. Last, the perfusion estimate cannot be used to separate different flow speeds. Together, these shortcomings may have limited the clinical use of cLDF.

By decomposing the Doppler power spectrum, the LDF measurements can be resolved into flow speeds.^{6, 7, 8, 9} The methods in Refs. 6, 7 were limited to single Doppler shifts, which is an unlikely situation in real measurements. The method presented in this article is an extension of Ref. 8, which handles multiple Doppler shifts in a flow phantom. In Ref. 9, the method in Ref. 7 was extended for multiple Doppler shifts, utilizing the principle of decomposing the Doppler power spectrum into a linear combination of Doppler shift probability distributions. However, with that approach, decomposition was not feasible at realistic noise levels.

The aim of this study was to present and evaluate a model-based quantitative LDF (qLDF) method that eliminates most of the shortcomings mentioned earlier. The method is based on model optimization of an adaptive three-layered Monte Carlo model, which is a simplification and generalization of the six-layered skin model that has previously been presented by our group.^{10} Model parameter estimation is done by comparing measured Doppler power spectra, at two different source–detector separations (0.25 and
$1.2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
), to calculated Doppler power spectra from an iteratively updated model. When the optimal fit is found, the perfusion is calculated from the model. We present the perfusion value in the absolute unit
$\mathrm{g}\phantom{\rule{0.3em}{0ex}}\mathrm{RBC}\u2215100\phantom{\rule{0.3em}{0ex}}\mathrm{g}$
$\text{tissue}\times \mathrm{mm}\u2215\mathrm{s}$
, resolved in three different speed regions, below
$1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
,
$1$
to
$10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
, and above
$10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
. The perfusion is given for a well-defined output volume of
$3.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{3}$
(half-sphere with radius
$1.13\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
). The model is developed for skin relevant optical properties at
$780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
but can be adapted to other wavelengths. The model can also adapt to homogenous tissues by adjusting its geometrical and optical properties.

To test the validity of the method, it is evaluated using simulations where, for example, the tissue optical properties are varied and discrete blood vessels are introduced. The goodness-of-fit for the adaptive tissue model is used to identify measurements that cannot be modeled accurately. An example of the method output during heat stimuli of the skin is also given.

## 2.

## Method Description

## 2.1.

### Data Collection and Calibration

A two-channel Periflux 5000 (Perimed AB, Järfälla, Sweden) modified for $50\text{-}\mathrm{kHz}$ digital recording of the detected signals (time varying part—ac; static light intensity—dc) was used for data collection. The wavelength of the light source (diode laser) was $780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . The probe was custom built by Perimed AB, with one light-emitting fiber placed in the center, one light-collecting fiber placed $0.25\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ (center-to-center distance) from the light-emitting fiber, and another 11 light-collecting fibers placed in a circle $1.2\pm 0.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ from the light-emitting fiber. All fibers were made of silica with a $0.125\text{-}\mathrm{mm}$ diameter and $\mathrm{NA}=0.37$ .

The recorded Doppler power spectra were noise-compensated by subtracting noise spectra compiled from static measurements at corresponding dc levels. The spectra were then
${\mathrm{dc}}^{2}$
-normalized to make them independent of the average light intensity. Additional compensation for the spectral characteristics of the system was also performed.^{10} Unlike conventional LDF, the latter step is important, especially at low frequencies below about
$100\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
.

In conventional LDF, the perfusion signal is calibrated using a measurement in a motility standard of a certain size distribution, concentration, and temperature. In order to be able to compare the measured and simulated spectra, another calibration procedure has been used in this study. The measured spectra were normalized by the total power of a spectrum from a measurement in motility standard (in principle, any scattering liquid would be appropriate), where all photons have been Doppler shifted due to Brownian motion with frequency shifts within the signal bandwidth. Thus, the total power for frequencies above $0\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ in the normalized spectra corresponded to the level of Doppler-shifted photons, and could thus be compared to simulated or forward-calculated spectra in an absolute manner. It should be noted that the calibration measurement had to be averaged over a relatively long time interval, about $30\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ , in order to cancel out high-power, low-frequency fluctuations in the calibration spectra.

## 2.2.

### Skin Model

The basic idea of the proposed method is to iteratively adjust an optical skin model until its output data matches measured LDF data. The model used is a three-layered structure consisting of one epidermal layer and two dermal layers. This is a simplification and generalization of a previously presented six-layered model.^{10} The epidermal layer has a variable thickness and contains no blood, whereas the dermal layers contain variable amounts of blood of variable speed distribution (Sec.
2.3.3). The superior dermal layer is
$0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
thick, and the second dermal layer has an infinite thickness. The optical properties of epidermis and dermis at
$780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
are given in Table 1
. (Ref. 10) Note that the refractive index
$n$
is set to 1.40 to avoid Fresnel reflections at the distinct borders between the layers.

## Table 1

Optical properties of bloodless skin at 780nm .

Tissue type | n (−) | μa (mm−1) | μs′ (mm−1) | μs (mm−1) | g (−) |
---|---|---|---|---|---|

Epidermis (2% eumelanin) | 1.40 | 0.4 | 3.5 | 23 | 0.85 |

Dermis | 1.40 | 0.1 | 2.0 | 13 | 0.85 |

The blood in the model was assumed to have 80% hemoglobin oxygenation, 43% hematocrit, a hemoglobin value of $145\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ Hb/liter blood, and mean cell hemoglobin concentration (MCHC) of $345\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ Hb/liter RBC. Given those values, the absorption coefficient of blood at $780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ is $0.43\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ (Ref. 11). The scattering coefficient ${\mu}_{\mathrm{s}}$ was set to $222\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , and the phase function was set to a Gegenbauer kernel phase function with ${\alpha}_{\mathrm{Gk}}=1.0$ and ${g}_{\mathrm{Gk}}=0.948$ , resulting in an anisotropy factor $g$ of 0.991 and a reduced scattering coefficient ${\mu}_{\mathrm{s}}^{\prime}={\mu}_{\mathrm{s}}(1-g)=2.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ (Ref. 10).

## 2.3.

### Forward Problem

In the forward problem, the Doppler power spectrum for a given RBC tissue fraction, speed distribution, and epidermal thickness is calculated. The forward problem is based on Monte Carlo simulations, but as these are time consuming, a fast post-processing technique is proposed that allows for an approximate, yet accurate, solution to the forward problem. For calculating the Doppler power spectrum using this post-processing technique, knowledge is needed about the intensity of backscattered photons from each layer and the path-length that the photons have propagated in each layer. That knowledge is used for calculating the distribution of the number of Doppler shifts that the photons have accumulated in each layer. For all layers, these distributions are combined with single-shifted optical Doppler spectra, calculated using the assumed speed distribution, to produce a simulated Doppler power spectrum. The first part of the proposed method for post-processing photon transport in a multilayered structure not only is limited to LDF simulations but could also be used for studying tissue absorption effects. The following subsections describe the proposed method in detail. A flowchart giving an overview of the forward problem is shown in Fig. 1 .

## 2.3.1.

#### Adding the effect of absorption

For a single-layer medium where the path-length distribution of the detected photons is known, the effect of absorption on the detected light intensity can be expressed using the Beer-Lambert law:

## 1

$$I={I}_{0}\sum _{j=1}^{D}{p}_{\mathrm{pl}}\left({d}_{j}\right)\mathrm{exp}(-{\mu}_{\mathrm{a}}{d}_{j}),$$For two or more layers, the same principle can be adopted, although it becomes more complex. Now, the path-length distribution ${p}_{\mathrm{pl}}$ is separated into the distributions, ${p}_{\mathrm{pl},m,n}$ , where $m$ is the layer number, and $n$ is the deepest layer that each contributing photon has propagated. For an $N$ -layer model, in total $N(N+1)\u22152$ distributions are thus required. In the first step, all these distributions are normalized to unity, i.e., ${\sum}_{j}{p}_{\mathrm{pl},m,n}\left({d}_{j}\right)=1$ , and then they are modified using the Beer-Lambert law:

## 2

$${p}_{\mathrm{pl},m,n}^{\prime}\left(d\right)={p}_{\mathrm{pl},m,n}\left(d\right)\mathrm{exp}(-{\mu}_{\mathrm{a},m}d).$$We also need another $N$ parameter, ${I}_{0,n}$ , corresponding to the detected intensity for photons that have been propagating in all layers down to layer $n$ in the nonabsorbing medium. Based on this, it is possible to calculate the detected backscattered intensity for photons that have been propagating in all layers down to layer $n$ in the medium:

## 3

$${I}_{n}={I}_{0,n}\prod _{m=1}^{n}\sum _{j=1}^{D}{p}_{\mathrm{pl},m,n}^{\prime}\left({d}_{j}\right).$$The parameters ${p}_{\mathrm{pl},m,n}$ and ${I}_{0,n}$ are given from a Monte Carlo simulation. To account for epidermal thickness variations in the model, ${p}_{\mathrm{pl},m,n}$ and ${I}_{0,n}$ are interpolated from a series of Monte Carlo simulations covering a range of epidermal thicknesses.

## 2.3.2.

#### Doppler power spectrum

For calculating the Doppler power spectrum, the shift distribution, i.e., the distribution of the number of Doppler shifts each photon undergoes, is needed. This can be calculated for each layer using the normalized modified path-length distribution ${p}_{\mathrm{pl},m,n}^{\u2033}\left(d\right)={p}_{\mathrm{pl},m,n}^{\prime}\left(d\right)\u2215{\sum}_{j}{p}_{\mathrm{pl},m,n}^{\prime}\left({d}_{j}\right)$ , the blood tissue fraction in the layer ${\zeta}_{\text{blood},m}$ , and the scattering coefficient of blood, ${\mu}_{\mathrm{s},\text{blood}}$ :

## 5

$${p}_{\text{Shifts},m,n}\left(s\right)=\sum _{j=1}^{D}{p}_{\mathrm{pl},m,n}^{\u2033}\left({d}_{j}\right)\mathrm{Po}(s,\mu ),$$For a known speed distribution in each layer, a single-shifted optical Doppler spectrum ${G}_{1}\left(\beta \right)$ can be generated with knowledge on the scattering phase function of the RBCs, assuming an isotropic relationship between the direction of the photon and the direction of the RBCs. In ${G}_{1}\left(\beta \right)$ , $\beta $ denotes the frequency. This spectrum is calculated as a histogram from $2.5\phantom{\rule{0.3em}{0ex}}\text{million}$ Doppler-shifted photons, where each Doppler shift ${f}_{\mathrm{D}}$ is stochastically generated using

## 6

$${f}_{\mathrm{D}}=-\mathbf{v}\cdot \mathbf{q}=-2\frac{v}{\lambda}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\frac{\theta}{2}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\phi ,$$The optical Doppler spectrum for any $s$ number of shifts can be calculated using the cross-correlation:

${G}_{0}\left(\beta \right)$ , i.e., the non-Doppler-shifted spectrum, equals the Dirac delta function. Since ${G}_{s}\left(\beta \right)$ can be assumed to be symmetric around $\beta =0$ (equal amount of positive and negative Doppler shifts), the cross-correlation can be replaced by a convolution. By Fourier transforming $\left(\mathcal{F}\{\dots \}\right)$ the optical Doppler spectrum, we can thus write:## 8

$$Y\left(\gamma \right)\equiv \mathcal{F}\left\{{G}_{1}\left(\beta \right)\right\}\Rightarrow \mathcal{F}\left\{{G}_{s}\left(\beta \right)\right\}={Y}^{s}\left(\gamma \right).$$## 9

$${Y}_{\mathrm{ms},m,n}\left(\gamma \right)=\sum _{s=0}^{S}{p}_{\text{Shifts},m,n}\left(s\right){Y}_{m,n}^{s}\left(\gamma \right).$$## 10

$${Y}_{\mathrm{ml}}\left(\gamma \right)=\sum _{n=1}^{N}{I}_{n}^{\prime}\prod _{m=1}^{n}{Y}_{\mathrm{ms},m,n}\left(\gamma \right),$$## 11

$${P}_{\mathrm{D}}\left(f\right)={\mathcal{F}}^{-1}\left\{{Y}_{\mathrm{ml}}^{2}\left(\gamma \right)\right\},$$^{13}.

An example of the forward calculation is given in Fig. 2 , where the light thin spectra originate from calculations of the three-layered model, and the dark thick spectra originate from a conventional Monte Carlo simulation of the same setup, i.e., based on the cumulative Doppler shift of each detected photon. The conventional Monte Carlo simulation lasted for several hours, while the forward calculation took on the order of one millisecond per spectrum using an ordinary PC and utilizing its graphics card for hyper-parallel computation.

## 2.3.3.

#### Speed distribution

The speed distribution in each layer is controlled by two parameters that control the slopes, i.e., the derivative, of the distribution for low and high speeds, respectively, in the loglog plane. Furthermore, the slope for speeds below $0.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ is zero, and the slope at $200\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ is $-50$ . The transition between the various slopes is smooth. An example showing the speed distribution with slopes of $-2$ and $-1$ for low and high speeds, respectively, is shown in Fig. 3 . The slopes for the low and high speeds are also drawn in that figure, as well as the zero slope for speeds below $0.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ and the slope of $-50$ at $200\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ . Those four slopes and smooth transitions between them constitute the speed distribution.

## 2.4.

### Inverse problem

To find the optimal skin model for a given measurement, the nonlinear inverse problem of matching calculated spectra to the measured spectra is solved using the Levenberg-Marquardt method. In each iteration, the forward problem is solved for the two source–detector separations 0.25 and $1.2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ using the same skin model parameters. The resulting calculated Doppler power spectra are then compared to the corresponding measured Doppler power spectra. Global convergence of the algorithm is ensured by using multiple start positions. Solving the inverse problem from a starting point relatively close to the solution (for example, the solution in the most recent measurement) takes about $0.5\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ with our implementation of the method using the graphics card on an ordinary desktop computer. A flowchart giving an overview of the inverse problem is shown in Fig. 4 .

## 2.4.1.

#### Model parameters

The applied model has six free parameters. The first parameter controls the epidermal thickness, given in mm. The second and third parameters control the blood concentrations in the two dermal layers. Parameters four and five control the slopes of the speed distributions in the dermis layers at low and high speeds, respectively. Last, the sixth parameter controls the relative difference between the speed distribution slopes of the two dermal layers.

## 2.4.2.

#### Difference between calculated and measured spectra

Both calculated and measured spectra are given in ${2}^{11}$ points evenly distributed between 0 and $25\phantom{\rule{0.3em}{0ex}}\mathrm{kHz}$ . Before calculating the difference between the spectra, they are log-transformed and reduced in length according to

## 12

$${P}_{\mathrm{red}}\left(n\right)={\mathrm{log}}_{10}{\u27e8P\left({f}_{i}\right)\u27e9}_{{j}_{n}\u2a7di<{j}_{n+1}},$$A ${\chi}^{2}$ -error is used for comparing the calculated and measured spectra. This error is calculated as the sum of the squared difference between the calculated and measured reduced spectra plus the sum of the squared difference of the estimated derivative of the reduced spectra. The ${\chi}^{2}$ -error is thus related both to the absolute levels of the reduced spectra and to their shapes, where the latter has been found important for good convergence properties of the least square problem. Mathematically, the ${\chi}^{2}$ -error is expressed as

## 14.

## 2.4.3.

#### Global convergence

The error function is minimized in a least-squares sense using the Levenberg-Marquardt algorithm with six free parameters. However, it is not guaranteed that the minimum found by the Levenberg-Marquardt method is the global minimum of the
${\chi}^{2}$
-error function. Therefore, as the
${\chi}^{2}$
-error function generally contains a number of local minima, other methods have to be used to find the area of attraction to the global minimum, i.e., a point from which the Levenberg-Marquardt method will find the local minimum that equals the global minimum. To find this area of attraction, a combination of a pure random search method and a multiple start method is used.^{14} First, an initial point is evaluated. If the
${\chi}^{2}$
-error of that point is small enough, the Levenberg-Marquardt method starts at that point. This is repeated until the same best local minimum is found three times, when it is assumed that that local minimum is equal to the global minimum. When randomly choosing a new starting point, it is assured that it differs considerably from the previous starting points.

## 2.5.

### Derived quantities

The estimated quantities in qLDF are given from the optimal best-fit model. The presented quantities are given for a half-sphere placed immediately beneath the light-emitting fiber of the probe, with a volume of
$3\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{3}$
(radius
$1.13\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
); see Fig. 5
. This volume is called the output volume, and the quantities from the three layers are weighted depending on the relative part of the half-sphere that each layer occupies. The RBC tissue fraction, average RBC speed, and perfusion (the RBC tissue fraction times their average speed), are given for three speed intervals: below
$1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
,
$1$
to
$10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
, and above
$10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
. The RBC tissue fractions and perfusions in these intervals are simply summarized to give the total RBC tissue fraction and perfusion, respectively, and the mean total speed is calculated by dividing the perfusion with the RBC tissue fraction. Furthermore, the RBC tissue fraction is given in the unit
$\mathrm{g}\phantom{\rule{0.3em}{0ex}}\mathrm{RBC}\u2215100\phantom{\rule{0.3em}{0ex}}\mathrm{g}$
tissue, the speed in mm/s, and the perfusion in
$\mathrm{g}\phantom{\rule{0.3em}{0ex}}\mathrm{RBC}\u2215100\phantom{\rule{0.3em}{0ex}}\mathrm{g}$
$\text{tissue}\times \mathrm{mm}\u2215\mathrm{s}$
. The RBC tissue fraction is given from the blood tissue fraction of the model, where the blood was set to a hematocrit level of 43%, hemoglobin value of
$145\phantom{\rule{0.3em}{0ex}}\mathrm{g}$
Hb/liter blood, and MCHC of
$345\phantom{\rule{0.3em}{0ex}}\mathrm{g}$
Hb/liter RBC. A tissue density of
$1.04\phantom{\rule{0.3em}{0ex}}\mathrm{kg}\u2215\text{liter}$
is assumed.^{15}

## 3.

## Robustness Analysis

The robustness of the presented method was evaluated by applying it on spectra generated directly by conventional Monte Carlo simulations, i.e., when the optical Doppler spectra are given directly by the distribution of cumulated Doppler shifts for the detected photons. Since the blood flow in the simulated models is known, that can be used to evaluate the accuracy of the method for arbitrary complex models. The results from a number of setups are presented in the following. If nothing else is stated, the optical properties of the models equaled the optical properties of the model presented in Sec. 2.1.

## 3.1.

### Three-Layered Models

The qLDF method was evaluated using simulations of three-layered models that could be exactly mimicked by the qLDF model itself. In total, six different setups with various epidermal thicknesses, blood tissue fractions, and speed distributions were evaluated. One extreme case of the three-layered model is when the epidermal layer thickness is zero and the two dermal layers are identical—i.e., a homogeneous model. Such a model constituted one of the six models.

The deviation between the estimated and true perfusion in these six models was between $-9.2\mathrm{\%}$ and 3.4% for the total perfusion. For the separate perfusion in the three speed regions, the deviation was between $-8.1\mathrm{\%}$ and 7.7% relative to the total perfusion. For all six models, the perfusion in the lowest speed region was slightly underestimated, 3.8% on average and maximum 5.2%, whereas the perfusion in the highest speed was overestimated, 4.5% on average and maximum 7.7%.

## 3.2.

### Random Models

For a more objective robustness analysis, the method was also tested on random setups of a complex seven-layered skin model. The top layer of that model mimicked an epidermal layer with variable thickness. Below the epidermis, six dermal layers were located, with variable thicknesses, blood tissue fractions, and speed distributions of the type described in Sec. 2.3.3. The thickness of the epidermal layer and the first five dermal layers was randomly chosen between 0 and $0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , whereas the last dermal layer extended to an infinite depth. For each setup, the blood tissue fraction, mean speed, and perfusion (blood tissue fraction times mean speed), were not allowed to differ more than a factor of 50 between the six dermal layers.

For comparison, both the qLDF and the cLDF method were evaluated on 50 of these random models. For the cLDF method, the perfusion was estimated separately at the two source–detector separations by calculating the first moment of the Doppler power spectra.

The total true perfusion of the simulated model, sorted in ascending order, is plotted against the estimated perfusion in Fig. 6 , and against the conventional perfusion estimates for the two source–detector separations in Fig. 6. The conventional perfusion estimates are normalized to attain a median of the 10 lowest points that equals the median of the true perfusion. For these 10 values, the perfusion is assumed to be in the linear region. The median of the absolute (i.e., ignoring the sign) relative difference between the true and estimated perfusion was 12% using the quantitative method. The corresponding difference for the normalized conventional perfusion estimates was 25% for the short source–detector separation and 26% for the long one. It can be noted that the conventional perfusion estimate, especially for the long source–detector separation, becomes rather nonlinear at high perfusion levels, due to multiple Doppler shifts. Apart from this nonlinearity, the perfusion estimate for the long separation shows a higher robustness (low level of fluctuations) than for the short separation.

The median difference between the perfusion estimates and the true perfusion for the three individual speed regions did not exceed 6% for any of the speed intervals, relative to the total true perfusion. The relative difference for the total estimated perfusion was between $-40$ and 61%, $-14$ and 74% for speeds below $1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ , $-16$ and 21% in the interval $1$ to $10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ , and $-34$ and 41% for speeds above $10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ . The median relative difference between the estimated and true perfusions relative to the total perfusion for the lowest speed region was 2%, for the middle speed region $-3\mathrm{\%}$ , for the highest speed region $-3\mathrm{\%}$ , and for the total perfusion $-5\mathrm{\%}$ .

## 3.3.

### Optical Properties of Static Media

The models presented in this and the coming sections (i.e., Secs. 3.3, 3.4, 3.5, 3.6) are based on the six-layered skin model presented in Ref. 10. Unlike that model, the speed distribution was smooth, as described in Sec. 2.3.3. For the standard six-layered skin model, the deviation between the total estimated perfusion and true perfusion was 0.0%, with $-6.7\mathrm{\%}$ for the lowest speed region, $-1.1\mathrm{\%}$ for the middle speed region, and 7.8% for the highest speed region.

The sensitivity to variations in absorption coefficient of the epidermal layer was evaluated in four different setups. This coefficient is foremost dependent on the tissue fraction of eumelanin. The evaluated absorption coefficients of 0.1, 0.4, 1.6, and
$7.7\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
correspond to a melanin tissue fraction of 0, 2, 10, and 50%, which can be compared to fractions of about 1 to 3% for fair skin and up to about 43% for black skin.^{16, 17}

The absorption coefficient of dermis as well as the scattering coefficient for both epidermis and dermis was also evaluated. These optical properties were varied between $-2$ and 2 standard errors from the mean values presented in Ref. 18, corresponding to $18.7$ to $27.3\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ for ${\mu}_{\mathrm{s}}$ in epidermis, $0.064$ to $0.136\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ for ${\mu}_{\mathrm{a}}$ in dermis, and $10.4$ to $15.6\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ for ${\mu}_{\mathrm{s}}$ in dermis. The effect of changing these optical properties, including the eumelanin tissue fraction, was negligible except for the absorption coefficient of dermis. The findings are summarized in Table 2 , where the relative deviations between the estimated and true perfusion are given, both for the total perfusion and for the perfusion in the three different speed regions relative to the total true perfusion.

## Table 2

Relative deviation (%) between the estimated and true perfusions when changing the optical properties of epidermis and dermis.

Speed region (mm/s) | ||||||||
---|---|---|---|---|---|---|---|---|

Total | v<1 | 1<v<10 | v> 10 | |||||

Min | Max | Min | Max | Min | Max | Min | Max | |

Epidermis ${\mu}_{\mathrm{a}}$ | $-2.3$ | 0.6 | $-9.5$ | $-8.1$ | $-1.6$ | $-0.3$ | 7.9 | 8.9 |

Epidermis ${\mu}_{\mathrm{s}}$ | $-2.4$ | 0.9 | $-8.5$ | $-4.1$ | $-2.3$ | $-0.1$ | 7.3 | 8.3 |

Dermis ${\mu}_{\mathrm{a}}$ | $-7.0$ | 8.4 | $-8.8$ | $-4.2$ | $-4.2$ | 0.9 | 5.0 | 11.7 |

Dermis ${\mu}_{\mathrm{s}}$ | $-2.7$ | 0.9 | $-8.6$ | $-6.7$ | $-2.6$ | 0.0 | 7.4 | 8.5 |

## 3.4.

### Optical Properties of Blood

A change in the hemoglobin oxygen saturation will change the absorption coefficient of blood. Totally reduced hemoglobin increases the absorption coefficient from $0.43$ to $0.58\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ $(+35\mathrm{\%})$ , whereas totally oxygenized hemoglobin decreases it from $0.43$ to $0.39\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ $(-9\mathrm{\%})$ . Furthermore, the absorption coefficient of blood is directly proportional to the mean cell hemoglobin concentration (MCHC) and the hematocrit. Therefore, the effect of a changed ( $\pm 40\mathrm{\%}$ , $0.26$ to $0.6\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ ) absorption coefficient of blood was evaluated. In analogy, the effect of a changed ( $\pm 40\mathrm{\%}$ , $133$ to $311\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ ) scattering coefficient of blood was evaluated.

The RBC shape and formation (aggregation) have an impact on the scattering phase function.^{19} The effect of this has been evaluated by varying the anisotropy factor
$g$
of the scattering phase function of the RBCs from the standard value of 0.991. As the Doppler shift is dependent on
$1-g$
rather than on
$g$
,
$1-g$
was varied
$\pm 40\mathrm{\%}$
, resulting in
$g$
from 0.9874 to 0.9946. For the applied Gegenbauer kernel phase function using
${\alpha}_{\mathrm{Gk}}=1.0,$
this interval corresponds to a
${g}_{\mathrm{Gk}}$
between 0.9365 and 0.9610.

The effect on the perfusion estimate when changing the absorption coefficient of blood is negligible. The situation is the contrary when changing the scattering properties of blood. This is summarized in Table 3 .

## Table 3

Relative deviations (%) between the estimated and true perfusions when changing the optical properties of blood.

Speed region (mm/s) | ||||||||
---|---|---|---|---|---|---|---|---|

Total | v<1 | 1<v<10 | v> 10 | |||||

Min | Max | Min | Max | Min | Max | Min | Max | |

Absorption coefficient | $-2.4$ | $-1.2$ | $-8.8$ | $-7.2$ | $-1.9$ | $-0.9$ | 7.4 | 7.9 |

Scattering coefficient | $-32$ | 34 | $-16$ | 5.8 | $-16$ | 12 | $-0.2$ | 16 |

Anisotropy factor | $-27$ | 22 | $-16$ | $-1.0$ | $-11$ | 9.1 | 0.2 | 14 |

## 3.5.

### Discrete Blood Vessels

For the three-layered model used in the presented method, the blood is homogeneously distributed in the layers. The effect of this simplification was also evaluated by simulating models containing discrete blood vessels. The speed distribution in these blood vessels was parabolic, with zero speed at the vessel wall and twice the average speed in the center.

The effect of locating most of the low-speed flow into discrete capillaries was negligible (within 5% for the total perfusion as well as for the perfusion in the different flow speed regions), even when the capillaries were located directly beneath the source or detector fiber. This was true when the capillaries were placed both in parallel with the skin surface or perpendicular to the skin surface. In these simulations, capillaries with an inner diameter of $8\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ were distributed with a density corresponding to approximately 20 capillary loops per ${\mathrm{mm}}^{2}$ .

For larger blood vessels, on the other hand, the result is very much dependent on the exact location of the vessels. When a large vessel is located just outside the $3\text{-}{\mathrm{mm}}^{3}$ half-sphere volume, it does not affect the true perfusion within the volume, but it still affects the Doppler power spectra and thus the estimated perfusion. In such a situation, the estimated perfusion will be overestimated. However, when a relatively large vessel is located deep in the $3\text{-}{\mathrm{mm}}^{3}$ volume, the perfusion is instead underestimated, but when it is placed shallowly just below the emitting fiber, it is again overestimated, since the Doppler power spectra are largely affected by the vessel. To test this, simulations were performed where a $50\text{-}\mu \mathrm{m}$ vessel with mean speed of $3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ (parabolic flow profile) was successively moved downward in the six-layered skin model, between 0.15 and $1.7\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ below the surface. The vessel was placed in parallel with the skin surface and directly beneath both the light-emitting fiber and the light-collecting fiber at $0.25\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . In Fig. 7 , the true and estimated total perfusion are plotted as a function of vessel depth, whereas the relative difference between the estimated and true total perfusion is plotted in Fig. 7, together with the ${\chi}^{2}$ -error between the simulated Doppler power spectra and the Doppler power spectra calculated from the best-fit model. Two examples of the spectral fit (depths 0.5 and $1.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) are shown in Fig. 8 .

## 3.6.

### Hair and Sweat Glands

The effect of a hair and a sweat gland in the skin was evaluated. The root of the hair shaft was placed straight below the illuminating fiber, in the center of the reticular dermis. It was modeled as a
$0.4\text{-}\mathrm{mm}$
-long cylinder, with a diameter of
$170\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
, tilting away from the closest detector fiber.^{20} The hair shaft was modeled as a
$0.7\text{-}\mathrm{mm}$
-long cylinder, with a diameter of
$40\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
, that extended the hair root and penetrated the epidermis surface about
$0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
away from the illuminating fiber. The optical properties of the hair were approximated to
${\mu}_{\mathrm{a}}=2.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
,
${\mu}_{\mathrm{s}}=13\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
,
$g=0.85,$
and
$n=1.55$
.

The sweat gland was modeled as a $20\text{-}\mu \mathrm{m}$ -diam, $0.1\text{-}\mathrm{mm}$ cylinder placed between the illuminating fiber and the closest light-collecting fiber, perpendicular to the tissue surface. The optical properties of the sweat gland were set to ${\mu}_{\mathrm{a}}=0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , ${\mu}_{\mathrm{s}}=2\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , $g=0.85,$ and $n=1.33$ . The refractive index is equal to the refractive index of water, which is most likely an underestimation.

The results showed a negligible effect for both the hair and the sweat gland, within 5% for both the total perfusion and for the perfusions in the different flow speed regions.

## 4.

## Example Measurement

The presented method has been tested *in vivo* on skin during a heat provocation on the dorsal side of the forearm. The measurements have previously been presented in Ref. 10 and were approved by the local ethical committee (D. No. 15-04). The measurements lasted for
$15\phantom{\rule{0.3em}{0ex}}\mathrm{min}$
. After the first minute, the measurement site was heated to
$38\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$
using the thermostatic probe holder PF 450 (Perimed AB), and after an additional
$5\phantom{\rule{0.3em}{0ex}}\mathrm{min}$
, the skin was heated to
$44\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$
. Details about the measurement system and signal preprocessing can be found in Sec.
2.1 and Refs. 6, 10. In Figs.
9, 10, 11, 12
, an example of the results is shown. This is the same measurement as in the worst-case example in Ref. 10.

In Fig. 9, the change of the RBC tissue fraction and perfusion in the three speed regions below $1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ , $1$ to $10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ , and above $10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ during the heat provocation is shown. The spectral difference ${\chi}^{2}$ is plotted in Fig. 10.

In Fig. 11, the total estimated perfusion—i.e., summarized over all speeds—is shown together with the conventional perfusion estimates for the two source–detector separations. The conventional measures are normalized to the same level as the quantitative measures during the first minute. The shape of the quantitative perfusion measure and the conventional perfusion measures is similar. However, the nonlinear behavior of the conventional perfusion measure can be observed, as the increase is not as high as for the quantitative perfusion. As expected, the nonlinearity is even more severe for the large source–detector separation.

Examples of the spectral fit between the calculated and measured spectra are shown in Fig. 12, both in the beginning (at $30\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ ) and at the end (at $14\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ ) of the measurement.

## 5.

## Discussion

## 5.1.

### Methodological Considerations

The presented qLDF measures are separated into three speed regions; below
$1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
,
$1$
to
$10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
, and above
$10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
. These measures are given for a half-sphere with a volume of
$3\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{3}$
. Additional volumes of 1, 2, 4, and
$5\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{3}$
were also tested. However, since the
$3\text{-}{\mathrm{mm}}^{3}$
volume gave the smallest deviations for the random models presented in Sec.
3.2, and also for the six-layered skin model, that volume was chosen. This volume covers almost equal parts of the two dermis layers of the three-layered model (depending on the epidermal thickness) and is between the measurement volumes normally probed in skin at
$780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
using the 0.25- and
$1.2\text{-}\mathrm{mm}$
source–detector separations.^{5} The volume is arbitrary and could be chosen without any distinct borders to better represent the probed volume, and that volume could be dependent on the best-fit model. This would increase the robustness of the method but would also significantly obstruct the physiological interpretation of the results, as the output volume would be nonconstant. The choice of speed regions is arbitrary. However, due to the number of free model parameters and their effect on the speed distribution of the model, it does not seem reasonable to have more than three regions.

In addition to the presented three-layered model with six free model parameters, the performance of a number of other models was also evaluated. The previously presented six-layered skin model^{10} with various number of free parameters and speed components was evaluated, as well as two-layered models and other variants of the three-layered model. However, the current model was chosen due to its high degree of generality—i.e., its ability to adjust to different microcirculatory conditions—and its robustness.

The ${\chi}^{2}$ -error is used for selecting the best model, but can also be used as an estimator of the reliability of the perfusion estimate. A large ${\chi}^{2}$ -error suggests that the three-layered model could not adapt to the measured spectra, and the estimated perfusion measures are then not reliable. This can be seen for example in Fig. 8, where ${\chi}^{2}=1.6$ and the deviation of the estimated total perfusion to the true perfusion is $-43\mathrm{\%}$ . We therefore suggest a ${\chi}^{2}$ threshold of 0.5. All of the estimated perfusion values from the example measurement in Sec. 4 resulted from models with ${\chi}^{2}<0.5$ (Fig. 10), and all adapted models to the random models in Sec. 3.2 except one had ${\chi}^{2}<0.5$ . For the one model with ${\chi}^{2}>0.5$ , the deviation of the total perfusion estimate was $-40\mathrm{\%}$ —i.e., the most underestimated case.

As seen in Fig. 11, the total quantitative perfusion measure increases more than the conventional perfusion measures as a response to the heat. The difference is especially large for the
$1.2\text{-}\mathrm{mm}$
-source–detector separation. This is due to the nonlinearity to the blood tissue fraction of the conventional perfusion estimates, which is foremost due to multiple Doppler shifts, but also due to a decreased measurement volume.^{5} A signal analysis that is supposed to correct for this nonlinearity has previously been presented.^{21} As that analysis was based on experiments in plastic flow models, there is no evidence that the correction of *in vivo* measurements has the expected effect. This type of nonlinearity is not present in the qLDF method.

The qLDF method requires a higher signal quality compared to cLDF. This can be accomplished by signal averaging or hardware improvements (higher power laser and/or more sensitive detector with less noise). Signal averaging limits its use when a high temporal resolution is required, as for example in spectral analysis of the perfusion signal.^{22}

The use of two different source–detector separations is of the utmost importance for the presented method, to avoid multiple solutions to the inverse problem. The robustness of the method would further improve if more measurement data were introduced into the method. On the other hand, it is likely that yet another source–detector separation would mostly give redundant information. One suggestion could be to use the detected dc signal as extra information, as that is also affected by the blood tissue fraction. Unfortunately, it is likely that the dc signal is even more affected by the melanin fraction in epidermis and other factors than by the blood tissue fraction, especially at
$780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, where the absorption of blood is low. Another suggestion is to add simultaneous diffuse reflectance spectroscopy (DRS)^{23, 24} measurements as an extra input to the model. This would be suitable for several reasons: DRS reveals a completely different type of information that would enhance the blood tissue fraction estimation; DRS is suitable for predicting the tissue optical properties that could be used for enhancing the model; and the calculations of the absorption effect are already included in the current qLDF methodology (Sec.
2.3.1).

## 5.2.

### Robustness Analysis

We performed a robustness analysis using simulations covering physiological variations of model geometry and tissue optical properties. For the random simulations presented in Sec. 3.2, it was shown that the robustness of the method is better than that of the conventional perfusion estimate, especially at the source–detector separation of $0.25\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , a separation that is used as standard in many commercially available LDF systems. For this analysis, the cLDF was normalized based on qLDF estimates in the linear region. It is shown that the conventional perfusion estimate from longer source–detector separations suffers from nonlinearities, but even when ignoring these nonlinearities, the robustness of the quantitative perfusion estimate is as least as good as the conventional perfusion estimate at $1.2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ .

It was shown (Table 2) that the quantitative perfusion estimate is relatively insensitive to variations in the optical properties of the static media. However, for other tissue types, the optical properties may differ considerably from those found in skin. For the method to work as accurately as possible for different tissue types, the model may need to be adapted to the actual optical properties, and also to the tissue geometry.

Since the absorption coefficient of blood is low at
$780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
and the blood tissue fraction is relatively low, a large relative change of the blood absorption coefficient has a negligible effect on the perfusion estimate (Table 3). The scattering properties have a higher impact, however, although it must be emphasized that the variations tested
$(\pm 40\mathrm{\%})$
are much larger than could reasonably be expected, and the variations of the anisotropy factor are also much larger than could normally be expected.^{19}

It is interesting to see that neither discrete capillaries, hairs, nor sweat glands introduce any significant errors in the perfusion estimate. Inhomogeneities in the form of larger vessels may have a higher impact, however, as demonstrated in Fig. 7. The largest errors occurred when a large vessel was located very close to the surface, which is not likely to occur in a real measurement situation. Furthermore, in a real measurement situation, it is not likely that the detected Doppler power spectrum is affected by only one discrete vessel. Nevertheless, large discrete blood vessels may add inhomogeneities that the three-layered model cannot handle. However, this is reflected by a large ${\chi}^{2}$ -error that can be used to discard such measurements. This mean of rejecting uncertain measurements is an important advantage compared to cLDF.

No validations have been performed using physical phantoms in this study. This is due to the extreme difficulties in constructing a controllable flow phantom that mimics the microcirculation in terms of geometry, optical properties, and flow to a sufficient degree. However, some of the basic principles of the method have previously been validated using a plastic flow phantom.^{8}

## 5.3.

### Physiological Implications

Generally, low blood flow speeds are found in small vessels and high speeds in larger vessels. This enables some differentiation between different types of blood vessels involved in the microcirculation in the probed volume, by evaluating the speed-resolved perfusion measure. However, in some situations, low speeds are also found in larger vessels—foremost in veins—a fact that complicates this interpretation.

For the presented example in Fig. 9, one can observe that the perfusion above $1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ is highly increased during the heat provocation, whereas the perfusion for low speeds remains unaffected or even decreases. This is expected since the heat regulation foremost involves relatively large vessels where the blood flow speed is high. This result supports the differentiating capabilities of qLDF.

The main advantages of the proposed qLDF method compared to conventional LDF are the output in physiologically relevant units ( $\mathrm{g}\phantom{\rule{0.3em}{0ex}}\mathrm{RBC}\u2215100\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ $\text{tissue}\times \mathrm{mm}\u2215\mathrm{s}$ ), the constant output volume, the differentiation in different speed regions, and the ${\chi}^{2}$ -indicator of uncertain estimates. Together, this has the potential to significantly facilitate the physiological interpretation of the results. A first clinical study that gives an indication of the clinical use of the method is presented in Ref. 25, where microvascular changes in diabetes type 2 are examined. Specifically, the speed-resolved measure has been used in that study to draw conclusions about which microvascular compartments are affected by the disease, conclusions that would be impossible to draw using cLDF.

## Acknowledgments

The study was financed by Perimed AB and Linköping University through the Center for Excellence NIMED-CBDP (Center for Biomedical Data Processing) and by VINNOVA and Perimed AB through the SamBIO research collaboration program between companies and academia within bioscience (VINNOVA D. No. 2008-00149).

## References

*Optical Diagnostics and Sensing VI*, G. L. Cote and A. V. Priezzhev, Eds., Proc. SPIE0277-786X 6094, 60940A (2006). 10.1117/12.659206 Google Scholar

*in vivo*laser Doppler flowmetry,” J. Biomed. Opt.1083-3668 13(1), 014015 (2008). 10.1117/1.2854691 Google Scholar