Translator Disclaimer
1 September 2010 Model-based quantitative laser Doppler flowmetry in skin
Author Affiliations +
Laser Doppler flowmetry (LDF) can be used for assessing the microcirculatory perfusion. However, conventional LDF (cLDF) gives only a relative perfusion estimate for an unknown measurement volume, with no information about the blood flow speed distribution. To overcome these limitations, a model-based analysis method for quantitative LDF (qLDF) is proposed. The method uses inverse Monte Carlo technique with an adaptive three-layer skin model. By analyzing the optimal model where measured and simulated LDF spectra detected at two different source-detector separations match, the absolute microcirculatory perfusion for a specified speed region in a predefined volume is determined. qLDF displayed errors <12% when evaluated using simulations of physiologically relevant variations in the layer structure, in the optical properties of static tissue, and in blood absorption. Inhomogeneous models containing small blood vessels, hair, and sweat glands displayed errors <5%. Evaluation models containing single larger blood vessels displayed significant errors but could be dismissed by residual analysis. In vivo measurements using local heat provocation displayed a higher perfusion increase with qLDF than cLDF, due to nonlinear effects in the latter. The qLDF showed that the perfusion increase occurred due to an increased amount of red blood cells with a speed >1 mm/s.



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.2mm ), 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 gRBC100g tissue×mms , resolved in three different speed regions, below 1mms , 1 to 10mms , and above 10mms . The perfusion is given for a well-defined output volume of 3.0mm3 (half-sphere with radius 1.13mm ). The model is developed for skin relevant optical properties at 780nm 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.


Method Description


Data Collection and Calibration

A two-channel Periflux 5000 (Perimed AB, Järfälla, Sweden) modified for 50-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 780nm . The probe was custom built by Perimed AB, with one light-emitting fiber placed in the center, one light-collecting fiber placed 0.25mm (center-to-center distance) from the light-emitting fiber, and another 11 light-collecting fibers placed in a circle 1.2±0.1mm from the light-emitting fiber. All fibers were made of silica with a 0.125-mm diameter and 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 dc2 -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 100Hz .

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 0Hz 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 30s , in order to cancel out high-power, low-frequency fluctuations in the calibration spectra.


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.5mm thick, and the second dermal layer has an infinite thickness. The optical properties of epidermis and dermis at 780nm 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.400.43.5230.85

The blood in the model was assumed to have 80% hemoglobin oxygenation, 43% hematocrit, a hemoglobin value of 145g Hb/liter blood, and mean cell hemoglobin concentration (MCHC) of 345g Hb/liter RBC. Given those values, the absorption coefficient of blood at 780nm is 0.43mm1 (Ref. 11). The scattering coefficient μs was set to 222mm1 , and the phase function was set to a Gegenbauer kernel phase function with αGk=1.0 and gGk=0.948 , resulting in an anisotropy factor g of 0.991 and a reduced scattering coefficient μs=μs(1g)=2.0mm1 (Ref. 10).


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 .

Fig. 1

Flowchart showing the forward problem of calculating a Doppler power spectrum from a given set of skin model parameters.



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:

Eq. 1

where I0 is the detected intensity in the nonabsorbing medium, ppl is the normalized path-length distribution, dj is the path length (mm) for path-length bin j , D is the number of path-length bins, and μa is the absorption coefficient (mm1) .

For two or more layers, the same principle can be adopted, although it becomes more complex. Now, the path-length distribution ppl is separated into the distributions, ppl,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)2 distributions are thus required. In the first step, all these distributions are normalized to unity, i.e., jppl,m,n(dj)=1 , and then they are modified using the Beer-Lambert law:

Eq. 2

where μa,m is the absorption coefficient of layer m .

We also need another N parameter, I0,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:

Eq. 3

The total detected intensity is calculated as the sum of these intensities:

Eq. 4


The parameters ppl,m,n and I0,n are given from a Monte Carlo simulation. To account for epidermal thickness variations in the model, ppl,m,n and I0,n are interpolated from a series of Monte Carlo simulations covering a range of epidermal thicknesses.


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 ppl,m,n(d)=ppl,m,n(d)jppl,m,n(dj) , the blood tissue fraction in the layer ζblood,m , and the scattering coefficient of blood, μs,blood :

Eq. 5

where Po(s,μ) denotes the probability of obtaining s number shifts from a Poisson distribution with parameter μ=djζblood,mμs,blood (Ref. 12).

For a known speed distribution in each layer, a single-shifted optical Doppler spectrum G1(β) 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 G1(β) , β denotes the frequency. This spectrum is calculated as a histogram from 2.5million Doppler-shifted photons, where each Doppler shift fD is stochastically generated using

Eq. 6

where v=|v| is the speed of the RBCs; q is the difference between the wave vectors ki and ks of the incident and scattered light wave, respectively; λ is the wavelength in the medium; θ is the scattering angle; and φ is the angle between q and v . In this stochastic generation of Doppler shifts, v is randomly chosen from the speed distribution (see Sec. 2.3.3), θ is randomly chosen from the scattering phase function of the RBCs, and cosφ is randomly chosen from a uniform distribution between 1 and 1.

The optical Doppler spectrum for any s number of shifts can be calculated using the cross-correlation:

Eq. 7

G0(β) , i.e., the non-Doppler-shifted spectrum, equals the Dirac delta function. Since Gs(β) can be assumed to be symmetric around β=0 (equal amount of positive and negative Doppler shifts), the cross-correlation can be replaced by a convolution. By Fourier transforming (F{}) the optical Doppler spectrum, we can thus write:

Eq. 8

We now extend the notation to Ym,n(γ) , where m and n are used in analogy with Eq. 2. The Fourier transformed multiple Doppler-shifted spectrum for all photons in layer m that have been propagating down to layer n can then be calculated as

Eq. 9

The next step is to summarize the spectra from all layers by cross-correlating them, i.e., calculating the product in the Fourier domain:

Eq. 10

where In=InmIm . Last, the Doppler power spectrum is calculated as

Eq. 11

where the square is introduced because the Doppler power spectrum equals the autocorrelation of the optical Doppler spectrum13.

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.

Fig. 2

Example of the forward-calculated spectra at two different source–detector separations ( 0.25mm lower, 1.2mm upper) from a three-layered model (light thinline), compared to Monte Carlo simulated spectra of the same setup (jagged blackline). The spectral difference was calculated to χ2=0.12 (see Sec. 2.4.2).



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.1mms is zero, and the slope at 200mms 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.1mms and the slope of 50 at 200mms . Those four slopes and smooth transitions between them constitute the speed distribution.

Fig. 3

Illustration of the speed distribution used; the thin lines indicate the slope or derivative at various speeds.



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.2mm 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.5s 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 .

Fig. 4

Flowchart showing how the inverse problem is solved ensuring global convergence.



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.


Difference between calculated and measured spectra

Both calculated and measured spectra are given in 211 points evenly distributed between 0 and 25kHz . Before calculating the difference between the spectra, they are log-transformed and reduced in length according to

Eq. 12

where ⟨…⟩ denotes the average, and fi denotes the frequency of frequency bin i . The indices jn are calculated as

Eq. 13

where i16kHz=1311 is the last frequency bin below 16kHz , and […] denotes the rounding to the closest integer. The limit at 16kHz is chosen due to the effect of the low-pass filtering of the time-varying part of the detector current in the measurement system. Although the filter characteristics are compensated for, the compensation becomes uncertain above 16kHz , where the signal is reduced more than 5dB by the filter. The indices jn are denser at low frequencies, where the dynamic of the Doppler power spectrum is higher, which emphasizes the impact of lower frequencies when comparing spectra. The averaging in Eq. 12 also prevents the calculated difference between the spectra from being influenced by noise.

A χ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 χ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 χ2 -error is expressed as


where the subscript “red” indicates the reduced spectra; “calc” and “meas” indicate the calculated and measured spectra, respectively; and ρ denotes the center–center source–detector separation (0.25 or 1.2mm ). The spectral noise generally increases with frequency and becomes unacceptably high at high frequencies relative to the signal. Therefore, the maximal index N is chosen as the first index where the difference between two adjacent frequency bins is higher than one decade, i.e., |log10[Pmeas,ρ(fN)]log10[Pmeas,ρ(fN+1)]|> 1 , where |…| denotes the absolute value.


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 χ2 -error function. Therefore, as the χ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 χ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.


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 3mm3 (radius 1.13mm ); 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 1mms , 1 to 10mms , and above 10mms . 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 gRBC100g tissue, the speed in mm/s, and the perfusion in gRBC100g tissue×mms . 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 145g Hb/liter blood, and MCHC of 345g Hb/liter RBC. A tissue density of 1.04kgliter is assumed.15

Fig. 5

Illustration of the output volume (half-sphere). The derived quantities are calculated from the model properties within this volume. The banana shapes represent the main probed volumes for the short (ρ=0.25) and long (ρ=1.2) source–detector separations. The emitting fiber is labeled e .



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.


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% and 3.4% for the total perfusion. For the separate perfusion in the three speed regions, the deviation was between 8.1% 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%.


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.3mm , 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.

Fig. 6

True total perfusion (thick line, dots) from the 50 simulated random models versus the estimated perfusion using the new method [(a), crosses] and versus the conventional perfusion estimates (b) at source–detector separations 0.25mm (crosses) and 1.2mm (plus signs).


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 1mms , 16 and 21% in the interval 1 to 10mms , and 34 and 41% for speeds above 10mms . 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% , for the highest speed region 3% , and for the total perfusion 5% .


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% for the lowest speed region, 1.1% 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.7mm1 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.3mm1 for μs in epidermis, 0.064 to 0.136mm1 for μa in dermis, and 10.4 to 15.6mm1 for μ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
Epidermis μa 2.3 0.6 9.5 8.1 1.6 0.3 7.98.9
Epidermis μs 2.4 0.9 8.5 4.1 2.3 0.1 7.38.3
Dermis μa 7.0 8.4 8.8 4.2 4.2
Dermis μs 2.7 0.9 8.6 6.7 2.6


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.58mm1 (+35%) , whereas totally oxygenized hemoglobin decreases it from 0.43 to 0.39mm1 (9%) . 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 ( ±40% , 0.26 to 0.6mm1 ) absorption coefficient of blood was evaluated. In analogy, the effect of a changed ( ±40% , 133 to 311mm1 ) 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 1g rather than on g , 1g was varied ±40% , resulting in g from 0.9874 to 0.9946. For the applied Gegenbauer kernel phase function using αGk=1.0, this interval corresponds to a gGk 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
Absorption coefficient 2.4 1.2 8.8 7.2 1.9 0.9 7.47.9
Scattering coefficient 32 34 16 5.8 16 12 0.2 16
Anisotropy factor 27 22 16 1.0 11 9.10.214


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μm were distributed with a density corresponding to approximately 20 capillary loops per mm2 .

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-mm3 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-mm3 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-μm vessel with mean speed of 3mms (parabolic flow profile) was successively moved downward in the six-layered skin model, between 0.15 and 1.7mm 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.25mm . 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 χ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.1mm ) are shown in Fig. 8 .

Fig. 7

Estimated (squares) and true (circles) perfusion are plotted in (a), and the relative difference between them (circles) and the χ2 -error of the model fit (squares) in (b), as a function of the depth of the 50-μm vessel. The outer radius of the output volume (1.13mm) is marked on the depth scale.


Fig. 8

Spectral fit between simulation (jagged blackline) and adapted model (smooth lightline) when the 50-μm vessel is located at the depth 0.5mm [(a), χ2=1.6 ] and 1.1mm [(b), χ2=0.22 ].



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-mm -long cylinder, with a diameter of 170μm , tilting away from the closest detector fiber.20 The hair shaft was modeled as a 0.7-mm -long cylinder, with a diameter of 40μm , that extended the hair root and penetrated the epidermis surface about 0.3mm away from the illuminating fiber. The optical properties of the hair were approximated to μa=2.0mm1 , μs=13mm1 , g=0.85, and n=1.55 .

The sweat gland was modeled as a 20-μm -diam, 0.1-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 μa=0mm1 , μs=2mm1 , 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.


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 15min . After the first minute, the measurement site was heated to 38°C using the thermostatic probe holder PF 450 (Perimed AB), and after an additional 5min , the skin was heated to 44°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.

Fig. 9

Absolute RBC tissue fraction (a) and perfusion (b) in the three speed intervals <1mms (thin blackline), 1 to 10mms (medium grayline), and > 10mms (thick blackline).


Fig. 10

Spectral difference χ2 over the measurement.


Fig. 11

Absolute total RBC perfusion (thin blackline) versus conventional perfusion estimates at the two source–detector separations 0.25mm (medium grayline) and 1.2mm (thick blackline).


Fig. 12

Spectral fit between the measured spectra (jagged blackline) and the calculated spectra (smooth lightline) at the two source–detector separations 0.25mm (lower spectra) and 1.2mm (upper spectra) at 30s [(a), χ2=0.12 ] and 14min [(b), χ2=0.16 ].


In Fig. 9, the change of the RBC tissue fraction and perfusion in the three speed regions below 1mms , 1 to 10mms , and above 10mms during the heat provocation is shown. The spectral difference χ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 30s ) and at the end (at 14min ) of the measurement.




Methodological Considerations

The presented qLDF measures are separated into three speed regions; below 1mms , 1 to 10mms , and above 10mms . These measures are given for a half-sphere with a volume of 3mm3 . Additional volumes of 1, 2, 4, and 5mm3 were also tested. However, since the 3-mm3 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 780nm using the 0.25- and 1.2-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 model10 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 χ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 χ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 χ2=1.6 and the deviation of the estimated total perfusion to the true perfusion is 43% . We therefore suggest a χ2 threshold of 0.5. All of the estimated perfusion values from the example measurement in Sec. 4 resulted from models with χ2<0.5 (Fig. 10), and all adapted models to the random models in Sec. 3.2 except one had χ2<0.5 . For the one model with χ2> 0.5 , the deviation of the total perfusion estimate was 40% —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-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 780nm , 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).


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.25mm , 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.2mm .

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 780nm 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 (±40%) 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 χ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


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 1mms 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 ( gRBC100g tissue×mms ), the constant output volume, the differentiation in different speed regions, and the χ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.


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).



G. E. Nilsson, T. Tenland, and P. Å. Öberg, “Evaluation of a laser Doppler flowmeter for measurement of tissue blood flow,” IEEE Trans. Biomed. Eng., BME–27 (10), 597 –604 (1980). 0018-9294 Google Scholar


R. Bonner and R. Nossal, “Model for laser Doppler measurements of blood flow in tissue,” Appl. Opt., 20 (12), 2097 –2107 (1981). 0003-6935 Google Scholar


M. Larsson, W. Steenbergen, and T. Strömberg, “Influence of optical properties and fiber separation on laser Doppler flowmetry,” J. Biomed. Opt., 7 (2), 236 –243 (2002). 1083-3668 Google Scholar


I. Fredriksson, “Quantitative laser Doppler flowmetry,” Department of Biomedical Engineering, Linköping University, (2009). Google Scholar


I. Fredriksson, M. Larsson, and T. Strömberg, “Measurement depth and volume in laser Doppler flowmetry,” Microvasc. Res., 78 (1), 4 –13 (2009). 0026-2862 Google Scholar


M. Larsson and T. Strömberg, “Toward a velocity-resolved microvascular blood flow measure by decomposition of the laser Doppler spectrum,” J. Biomed. Opt., 11 (1), 014024 (2006). 1083-3668 Google Scholar


A. Liebert, N. Zolek, and R. Maniewski, “Decomposition of a laser-Doppler spectrum for estimation of speed distribution of particles moving in an optically turbid medium: Monte Carlo validation study,” Phys. Med. Biol., 51 (22), 5737 –5751 (2006). 0031-9155 Google Scholar


I. Fredriksson, M. Larsson, and T. Strömberg, “Absolute flow velocity components in laser Doppler flowmetry,” Proc. SPIE, 6094 60940A (2006). 0277-786X Google Scholar


S. Wojtkiewicz, A. Liebert, H. Rix, N. Zolek, and R. Maniewski, “Laser-Doppler spectrum decomposition applied for the estimation of speed distribution of particles moving in a multiple scattering medium,” Phys. Med. Biol., 54 (3), 679 –697 (2009). 0031-9155 Google Scholar


I. Fredriksson, M. Larsson, and T. Strömberg, “Optical microcirculatory skin model: assessed by Monte Carlo simulations paired with in vivo laser Doppler flowmetry,” J. Biomed. Opt., 13 (1), 014015 (2008). 1083-3668 Google Scholar


W. G. Zijlstra, A. Buursma, and O. W. van Assendelft, Visible and Near Infrared Absorption Spectra of Human and Animal Haemoglobin: Determination and Application, VSP Books, Utrecht, Boston, Köln, Tokyo (2000). Google Scholar


A. Serov, W. Steenbergen, and F. F. M. de Mul, “Prediction of the photodetector signal generated by Doppler-induced speckle fluctuations: theory and some validations,” J. Opt. Soc. Am., 18 (3), 622 –630 (2001). 0030-3941 Google Scholar


A. T. Forrester, “Photoelectric mixing as a spectroscopic tool,” J. Opt. Soc. Am., 51 (3), 253 –259 (1961). 0030-3941 Google Scholar


M. Hazen, “Search strategies for global optimization,” Department of Electrical Engineering, University of Washington, (2008). Google Scholar


A. R. C. Lee and K. Tojo, “An experimental approach to study the binding properties of vitamin E (alpha-tocopherol) during hairless mouse skin permeation,” Chem. Pharmaceut. Bull., 49 (6), 659 –663 (2001). Google Scholar


S. L. Jacques, “Skin optics,” Oregon Medical Laser Center News, (1998) Google Scholar


A. J. Thody, E. M. Higgins, K. Wakamatsu, S. Ito, S. A. Burchill, and J. M. Marks, “Pheomelanin as well as eumelanin is present in human epidermis,” J. Invest. Dermatol., 97 (2), 340 –344 (1991). 0022-202X Google Scholar


E. Salomatina, B. Jiang, J. Novak, and A. N. Yaroslavsky, “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt., 11 (6), 064026 (2006). 1083-3668 Google Scholar


A. M. K. Enejder, J. Swartling, P. Aruna, and S. Andersson-Engels, “Influence of cell shape and aggregate formation on the optical properties of flowing whole blood,” Appl. Opt., 42 (7), 1384 –1394 (2003). 0003-6935 Google Scholar


N. Otberg, H. Richter, H. Schaefer, U. Blume-Peytavi, W. Sterry, and J. Lademann, “Variations of hair follicle size and distribution in different body sites,” J. Invest. Dermatol., 122 (1), 14 –19 (2004). 0022-202X Google Scholar


G. E. Nilsson, “Signal processor for laser Doppler tissue flowmeters,” Med. Biol. Eng. Comput., 22 (4), 343 –348 (1984). 0140-0118 Google Scholar


U. Hoffmann, U. K. Franzeck, M. Geiger, A. Yanar, and A. Bollinger, “Variability of different patterns of skin oscillatory flux in healthy controls and patients with peripheral arterial occlusive disease,” Int. J. Microcirc.: Clin. Exp., 12 (3), 255 –273 (1993). 0167-6865 Google Scholar


G. M. Palmer and N. Ramanujam, “Monte Carlo–based inverse model for calculating tissue optical properties. Part I: theory and validation on synthetic phantoms,” Appl. Opt., 45 (5), 1062 –1071 (2006). 0003-6935 Google Scholar


G. M. Palmer, C. Zhu, T. M. Breslin, F. Xu, K. W. Gilchrist, and N. Ramanujan, “Monte Carlo–based inverse model for calculating tissue optical properties. Part II: application to breast cancer diagnosis,” Appl. Opt., 45 (5), 1072 –1078 (2006). 0003-6935 Google Scholar


I. Fredriksson, M. Larsson, F. H. Nyström, T. Länne, C. J. Östgren, and T. Strömberg, “Reduced arterio-venous shunting capacity after local heating and redistribution of baseline skin blood flow in type 2 diabetes assessed with velocity-resolved quantitative laser Doppler flowmetry,” Diabetes, 59 (7), 1578 –1584 (2010). 0012-1797 Google Scholar
©(2010) Society of Photo-Optical Instrumentation Engineers (SPIE)
Ingemar Fredriksson, Marcus Larsson, and Tomas Strömberg "Model-based quantitative laser Doppler flowmetry in skin," Journal of Biomedical Optics 15(5), 057002 (1 September 2010).
Published: 1 September 2010

Back to Top