Inverse Monte Carlo in a multilayered tissue model: merging diffuse reflectance spectroscopy and laser Doppler flowmetry

Ingemar Fredriksson; Oleg Burdakov; Marcus Larsson; Tomas Strömberg

doi:10.1117/1.JBO.18.12.127004

18 December 2013 Inverse Monte Carlo in a multilayered tissue model: merging diffuse reflectance spectroscopy and laser Doppler flowmetry

Ingemar Fredriksson, Oleg Burdakov, Marcus Larsson, Tomas Strömberg

Author Affiliations +

Journal of Biomedical Optics, Vol. 18, Issue 12, 127004 (December 2013). https://doi.org/10.1117/1.JBO.18.12.127004

Abstract

The tissue fraction of red blood cells (RBCs) and their oxygenation and speed-resolved perfusion are estimated in absolute units by combining diffuse reflectance spectroscopy (DRS) and laser Doppler flowmetry (LDF). The DRS spectra (450 to 850 nm) are assessed at two source–detector separations (0.4 and 1.2 mm), allowing for a relative calibration routine, whereas LDF spectra are assessed at 1.2 mm in the same fiber-optic probe. Data are analyzed using nonlinear optimization in an inverse Monte Carlo technique by applying an adaptive multilayered tissue model based on geometrical, scattering, and absorbing properties, as well as RBC flow-speed information. Simulations of 250 tissue-like models including up to 2000 individual blood vessels were used to evaluate the method. The absolute root mean square (RMS) deviation between estimated and true oxygenation was 4.1 percentage units, whereas the relative RMS deviations for the RBC tissue fraction and perfusion were 19% and 23%, respectively. Examples of in vivo measurements on forearm and foot during common provocations are presented. The method offers several advantages such as simultaneous quantification of RBC tissue fraction and oxygenation and perfusion from the same, predictable, sampling volume. The perfusion estimate is speed resolved, absolute (% RBC×mm/s ), and more accurate due to the combination with DRS.

1. Introduction

The microcirculation involves the smallest vessels in the tissue: the arterioles, the capillaries, and the venules. These vessels are arranged into a refined microvascular network with a primary function of delivering oxygen to the surrounding tissue. This function depends strongly on the oxygen gradients in the microcirculation, a property that is governed by the balance between blood flow and tissue oxygen demand.¹ In tissues such as the brain, with a high-metabolic oxygen demand that is balanced by a high-blood flow, the oxygen transport to tissue mainly occurs in the capillaries. However, in resting muscle, with a lower oxygen demand and a low-blood flow, the oxygen transport occurs already in the arterioles. The above-mentioned study clearly indicates that both oxygenation and blood flow are important parameters to assess when evaluating the microcirculation. This article describes a method for measuring blood flow and oxygenation parameters in the microcirculation by integrating diffuse reflectance spectroscopy (DRS) and laser Doppler flowmetry (LDF).

In the visible to near-infrared wavelength range, the absorption spectra of oxygenized and reduced hemoglobin show distinct characteristics, which are advantageous for the determination of tissue oxygenation using DRS. However, the large dynamic range in the absorption coefficient, as displayed by the tissue chromophores in this wavelength range, has proven difficult to model accurately in the presence of a tissue-scattering coefficient in the same magnitude. For setups using small source–detector separations, at about 1 mm and below, diffusion theory fails to accurately describe light propagation in tissue.²^,³ Instead, numerical simulations using Monte Carlo techniques provide a way to overcome these deficiencies.

Quantitative measures of the amount of tissue chromophores can be attained by analyzing calibrated DRS data using an inverse Monte Carlo technique, assuming that the included tissue chromophores are known.⁴ To obtain accurate estimations, all tissue chromophores of importance need to be included in the model, or else included chromophores will falsely compensate for excluded ones.⁵^,⁶ The effect of confining blood to vessels rather than assuming a homogeneous distribution must also be taken into account.⁷^–⁹ Still, it has been shown that using inverse Monte Carlo with a single-layer model fails to describe how light propagates through skin tissue and how it is detected at multiple source–detector separations.⁶ This result strongly indicates that DRS data from skin should be analyzed using a multilayered model, in order to accurately estimate the chromophore content.

Inverse analysis techniques using Monte Carlo simulations of two-layered models for estimating the tissue optical properties (OPs) have previously been used by Wang et al.¹⁰ and Yudovsky.¹¹ By analyzing the spatial frequency domain reflectance, Yudovsky showed that it is possible to estimate both dermis absorption and reduced scattering, while epidermal optical thickness could only be determined with limited accuracy. Wang et al. showed that it is difficult to accurately estimate the OP of both layers when analyzing the spatially resolved diffuse reflectance. These studies essentially determine OP using absolute calibrated reflectance at one wavelength, but imply extension to spectroscopic use.

When light is scattered by moving red blood cells (RBCs), a small frequency shift will occur. The effect of this frequency shift can be studied with an LDF instrument in order to indirectly estimate the microcirculatory blood flow. Commonly, an estimate is formed by calculating the first moment of the detected Doppler power spectrum.¹²^,¹³ This measure increases linearly with the RBC speed and nonlinearly with the RBC concentration.¹⁴ It is also highly affected by the tissue OPs and the structure of the tissue and the microvasculature. Consequently, it is given in nonabsolute units, a fact that complicates comparisons between and within individuals. Neither can this measure be used to differentiate between flow speeds. Previous results indicate that these shortcomings can be overcome by using a model-based data analysis approach based on presimulated Monte Carlo data.¹⁵^,¹⁶

The combination of LDF and DRS in the same fiber-optic probe has been presented by other groups.¹⁷^,¹⁸ These approaches involved conventional LDF and a modified Beer–Lambert law¹⁷ or Kubelka–Munk algorithm¹⁸^,¹⁹ for DRS. These simplified light transport algorithms have problems compensating for the effect of tissue OPs, reduced scattering and melanin absorption in epidermis, and source–detector separation.¹⁷ The total hemoglobin concentration during occlusion release was, for example, highly affected by source–detector separation in Ref. 17.

The two studies by Wang et al. and Yudovsky mentioned above indicate that there is a need for a more sophisticated light transport model when merging DRS and LDF. The approach of joining DRS and LDF, as proposed in this article, was based on previously presented algorithms, where either of the two modalities has been used in combination with an adaptive tissue model.¹⁵^,²⁰ Each of these algorithms uses a subset of presimulated Monte Carlo data for a limited number of tissue geometry and scattering parameters, while both Doppler and absorption effects are directly added in the inverse Monte Carlo algorithm. By using recordings at two source–detector separations, the sensitivity to the layered structure of the tissue is increased. Furthermore, the use of two source–detector separations enables a relative calibration, where only the relative intensity difference between the two DRS channels has to be calibrated. This eliminates the need for an absolute calibration of the DRS system, which is difficult to perform accurately in a clinical setting.

The aim of this article was to present an inverse Monte Carlo method that uses a multilayered tissue model to estimate the microcirculatory parameters in a setup based on a fiber-optic probe that integrates the DRS and LDF modalities. The tissue model parameters are found by comparing measured DRS spectra at two different source–detector separations (0.4 and 1.2 mm) and measured LDF spectra at one separation (1.2 mm) with calculated DRS and LDF spectra from the tissue model. When the optimal fit is found, the output data are given from the model. We present the tissue fraction of RBCs in mass percentage (%) [equal to (g $RBC / 100 g$ tissue)], the hemoglobin oxygen saturation in (%), and the perfusion in the absolute unit (% $RBC \times mm / s$ ), resolved in three different speed regions: below $1 mm / s$ , 1 to $10 mm / s$ , and above $10 mm / s$ . A significant part of the work has been spent on finding a fast and robust way of solving the inverse nonlinear Monte Carlo optimization problem.

The accuracy of the method is evaluated using Monte Carlo simulations of well-known tissue-like models containing discrete blood vessels. The performance of the method is compared with conventional perfusion estimates. Examples of in vivo measurements are also given.

2. Method

2.1.

Adaptive Skin Model

The skin was modeled in the wavelength range 450 to 850 nm as a three-layer structure consisting of a bloodless epidermis layer and two dermis layers containing different amounts of blood with variable speed distributions and equal blood oxygen saturations. The thickness of the epidermis layer was variable, whereas the upper dermis layer had a fixed thickness of 0.5 mm and the lower had an infinite thickness. All layers had the same wavelength-dependent scattering coefficient. The model also contained an average vessel diameter, twice as large in the deeper dermis layer as in the upper, which was used to compensate for the vessel-packaging effect.⁷^,⁸ The model has previously been described in Ref. 15, 20, and 21, and a brief summary of the parameters is given below.

The thickness of the bloodless epidermis layer was given by one parameter.

The wavelength-dependent scattering coefficient, $μ_{s} (λ)$ , was described by three parameters: $α$ , $β$ , and $γ$ according to

Eq. (1)

μ_{s} (λ) = α [(1 - γ) {(λ / λ_{0})}^{- β} + γ {(λ / λ_{0})}^{- 4}],

where

λ

is the wavelength in nm and

λ_{0} = 600 nm

. A fixed Henyey–Greenstein scattering phase function with

g = 0.8

was used.

The melanin content of the epidermis layer was given by one parameter, $f_{mel}$ , and the absorption coefficient of the epidermis layer was modeled as

Eq. (2)

μ_{a, 0} (λ) = f_{mel} μ_{a, mel} (λ),

where²²

Eq. (3)

μ_{a, mel} (λ) = 6.6 \times 10^{10} λ^{- 10 / 3} .

The blood was assumed to have a hematocrit of 43%, a hemoglobin concentration of $145 g / l$ blood, and a mean cell hemoglobin concentration of $345 g / l$ RBC. The wavelength-dependent absorption coefficients for oxygenated and deoxygenated blood ( $μ_{a, oxy}$ and $μ_{a, deoxy}$ ) were derived from Ziljstra et al.²³ A scattering coefficient of blood of $222 {mm}^{- 1}$ at the laser wavelength (780 nm) was used, and at the same wavelength, a Gegenbauer kernel phase function with parameters $g_{Gk} = 0.948$ and $α_{Gk} = 1.0$ , resulting in an anisotropy factor of 0.991, was used.²¹ The oxygen saturation ( $s$ ) was given by one parameter and was assumed to be equal in the two dermal layers. The absorption coefficient of blood in layer $n$ was calculated as

Eq. (4)

μ_{a, blood, n} (λ) = [s μ_{a, oxy} (λ) + (1 - s) μ_{a, deoxy} (λ)] c_{vp, n},

where

c_{vp, n}

is a factor compensating for the so called vessel-packaging effect. This factor is given by⁷^,⁸^,²⁴

Eq. (5)

c_{vp, n} (λ) = \frac{1 - \exp [- D_{n} μ_{a, blood} (λ)]}{D_{n} μ_{a, blood} (λ)},

where

D_{n}

is the average vessel diameter in layer

n

(variable parameter twice as high for the deeper dermis layer). It can be realized that

c_{vp, n} (λ) \to 1

when

D_{n} \to 0

or

μ_{a, blood} (λ) \to 0

. The absorption coefficient of the two dermal layers was given by (

n \in [1, 2]

)

Eq. (6)

μ_{a, n} (λ) = f_{blood, n} μ_{a, blood, n} (λ),

where

f_{blood, n}

was the volume fraction in the two dermis layers given by two parameters.

The blood in the dermis layers had a speed distribution that was given by 10 parameters. All these 10 parameters had a rectangular speed distribution between 0 and twice the mean speed of each parameter, in order to resemble the parabolic flow profile that can be expected in the blood vessels.

In total, the model was controlled by 1 (epidermal thickness) $+ 3$ (scattering) $+ 1$ (melanin fraction) $+ 2$ (blood tissue fraction) $+ 1$ (oxygen saturation) $+ 1$ (mean vessel diameter) $+ 10$ (speed distribution), i.e., 19 variable parameters. As will be seen in Sec. 2.3, the first nine parameters are determined by fitting DRS spectra and the 10 speed parameters by fitting the LDF spectra.

2.2.

Forward Problem

The forward problem of calculating DRS²⁰ and LDF spectra¹⁵ from the model has previously been described. A summary is given in the following subsections. The vessel-packaging effect for LDF is given extra attention, since that has not been described before but yet has a high impact on the resulting spectra.

2.2.1.

Generation of path length distributions

The forward calculations are based on a limited number of base simulations with various levels of the scattering coefficient of the three layers ( $μ_{s} = 2^{x}, x \in [0, 0.5, \dots, 6.5] {mm}^{- 1}$ , 14 levels) and various thicknesses of the epidermis layer ( $t_{epi} \in {[0.05, 0.1, \dots, 0.7]}^{2} mm$ , 14 levels), i.e., in total 196 simulations.

The path lengths in each of the three layers were stored for each detected photon in the base simulations. These were used to generate path length distributions that were used in the forward calculations. Path length distributions were generated for the epidermis layer for photons that had only been propagating in the epidermis layer, for the epidermis and first dermis layers for photons that had only been propagating in those two layers, and for all three layers for photons that had been propagating in all three layers. In total this gives six path length distributions per simulation. The notation used for these distributions throughout this article is $p_{pl, m, n} (d)$ , where $d$ is the path length, $m$ is the layer number, and $n$ is the number of the deepest layer that the photon reached, i.e., $m \leq n$ (see Ref. 20 for more details and for an example). These simulations, as well as the evaluation simulations presented in Sec. 2.5, were performed with the previously presented and validated in-house Monte Carlo software.²⁵

2.2.2.

DRS forward calculation

The DRS spectra are calculated by first interpolating the base simulations based on the scattering coefficient and the epidermis thickness of each wavelength of the current model. Linear interpolation was used in the dimension of the scattering coefficient, whereas the thickness dimension was interpolated using a cosinus-shaped kernel covering the three closest thickness levels. For the latter dimension, this kind of smoothing was aimed at preventing the appearance of local minima in the inverse solution. The absorption effect was then added in each layer by applying Beer–Lambert’s law for each path length in the path length distributions. Equations (7)–(9) describe this mathematically.

For each wavelength, the total intensity $I_{0, n}$ for photons that had reached layer $n$ was calculated as

Eq. (7)

I_{0, n} = \sum_{d} p_{pl, m, n} (d),

where

I_{0, n}

is independent of

m \leq n

. The path length distributions were then normalized to unity

p_{pl, m, n}^{'} (d) = p_{pl, m, n} (d) / I_{0, n}

. The path length distributions were modified for all path lengths

d

by adding the absorption using Beer–Lambert’s law

Eq. (8)

p_{pl, m, n}^{''} (d) = p_{pl, m, n}^{'} (d) \exp (- d μ_{a, m}),

where

μ_{a, m}

is the absorption coefficient of layer

m

[Eqs. (2) and (6) for the epidermis layer and the dermis layers, respectively].

The detected intensity for all photons that had propagated down to layer $n$ and back was then calculated by multiplying $I_{0, n}$ with the total absorption effect from all layers the photon had propagated through as

Eq. (9)

I_{n} = I_{0, n} \prod_{m = 1}^{n} \sum_{j} p_{pl, m, n}^{''} (d_{j}),

and the total detected intensity

I

for each wavelength was then simply calculated as the sum of the

I_{n}

’s

Eq. (10)

I = \sum_{n} I_{n} .

An example of the accuracy of this forward calculation can be found elsewhere.²⁰

2.2.3.

LDF forward calculation

The forward calculation of LDF spectra starts with the same steps as for the DRS forward calculation, but only with the laser wavelength (780 nm in this study). In Ref. 15, it is described how single-shifted optical Doppler spectra are calculated and summarized using shift distributions that are based on the path length distributions. Here, taking the vessel-packaging effect into account, the optical Doppler spectra resulting from the passing through a single vessel are calculated and summarized using a distribution of the number of vessels the light passes through. This is a small but important difference from Ref. 15 and essentially compensates for a vessel-packaging effect in LDF. A mathematical description of the forward calculation of LDF spectra is given in Eqs. (11)–(17).

Consider a single-shifted optical Doppler spectrum originating from a certain evenly distributed speed distribution between 0 and $2 v mm / s$ , $y_{1, v} (f)$ . The calculation of a single-shifted optical Doppler spectrum for a single speed is described in the Appendix. The optical Doppler spectrum for any number of shifts is given by cross-correlating the single-shifted optical Doppler spectrum as

Eq. (11)

y_{s, v} (f) = y_{s - 1, v} (f) ⋆ y_{1, v} (f),

where

⋆

denotes the cross-correlation function.

The non-Doppler shifted spectrum, $y_{0, v} (f)$ , equals the Dirac’s delta function. Since $y_{s, v} (f)$ can be assumed to be symmetric around $f = 0$ (equal amount of positive and negative Doppler shifts), the cross-correlation is equal to a convolution. Fourier transforming ( $F {\cdot}$ ) the optical Doppler spectrum results in $Y_{v} (γ) \equiv F {y_{1, v} (f)}$ and thus $F {y_{s, v} (f)} = Y_{v}^{s} (γ)$ [i.e., $Y_{v} (γ)$ to the power of $s$ ], which is a more convenient form.

Then consider a vessel with diameter $D$ with a parabolic flow profile with mean speed $v$ . When light passes perpendicularly through the center axis of that vessel, the light will be Doppler shifted a number of times following a Poisson distribution with expectation value ${\hat{n}}_{shifts} = D μ_{s, blood}$ (neglecting a slightly increased path length due to scattering). The light will thus on average be Doppler shifted ${\hat{n}}_{shifts}$ times when passing the vessel, under the assumption that the average path length through the vessel is $D$ (some will propagate through the periphery of the vessel and thus a shorter way, and some will propagate through the vessel with a rather oblique angle and thus a longer way). The Fourier-transformed optical Doppler spectrum from a single vessel is then given by

Eq. (12)

Y_{vessel, v} (γ) = \sum_{s} Y_{v}^{s} (γ) p_{Po} (s; {\hat{n}}_{shifts}),

where

p_{Po} (s; {\hat{n}}_{shifts})

denotes the probability of

s

shifts given by the Poisson distribution with

{\hat{n}}_{shifts}

. Summarizing these Fourier-transformed optical Doppler spectra over the whole distribution of vessels with various flow speeds in the layer results in

Eq. (13)

Y_{vessel} (γ) = \sum_{j} p_{v, j} Y_{vessel, v_{j}} (γ),

where

p_{v, j}

denotes the fraction of blood moving with speed

v_{j}

.

For each path length $d_{j}$ in a certain layer $m$ , the distribution of the number of vessels that are intersected can be assumed to follow a Poisson distribution with parameter $μ_{vessels, m} = d_{j} f_{blood, m} / D_{m}$ . Thus, the distribution of intersected vessels can be calculated from the normalized absorption-modified path length distribution $p_{pl, m, n}^{'''} (d) = p_{pl, m, n}^{''} (d) / \sum_{j} p_{pl, m, n}^{''} (d_{j})$ (normalized since the distribution of intersected vessels is not dependent on the total light intensity of $p_{pl, m, n}^{''}$ ) by

Eq. (14)

p_{vessels, m, n} (s) = \sum_{j} p_{pl, m, n}^{'''} (d_{j}) p_{Po} (s; μ_{vessels, m}) .

The Fourier-transformed optical Doppler-shifted spectrum for light that has passed several vessels for all photons in layer $m$ that have been propagating down to layer $n$ can then be calculated as

Eq. (15)

Y_{ms, m, n} (γ) = \sum_{v = 0}^{V} p_{vessels, m, n} (v) Y_{vessel, m, n}^{v} (γ),

where

V

is the maximum number of vessels that the light has passed through.

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

Eq. (16)

Y_{ml} (γ) = \sum_{n = 1}^{N} I_{n}^{'} \prod_{m = 1}^{n} Y_{ms, m, n} (γ),

where

I_{n}^{'} = I_{n} / \sum_{m} I_{m}

, and

N

is the index of the deepest dermis layer (i.e., 2). Finally, the Doppler power spectrum

y_{D} (f)

is calculated as

Eq. (17)

y_{D} (f) = F^{- 1} {Y_{ml}^{2} (γ)} .

2.3.

Inverse Problem

The inverse problem is to find a model and its parameters that fit measured DRS spectra at two source–detector separations (0.4 and 1.2 mm) and LDF spectra at one source–detector separation (1.2 mm). The use of two source–detector separations is needed in the DRS case in order to avoid multiple solutions to a single measurement. In the case of DRS, the use of two source–detector separations also makes the calibration routine more simple and reliable in practice, as only the relative intensity difference between the two separations then has to be considered when solving the inverse problem, not the absolute levels. In fact, the importance of an accurate relative intensity calibration between the two separations can be reduced by introducing an adaptive amplification factor $q$ for one of the separations, too (see Sec. 2.3.1). The negative impact on the solution that the introduction of $q$ may lead to is limited by penalizing solutions which requires a $q$ other than a unity.

2.3.1.

Fitting DRS spectra

The inverse problem of finding a set of parameters that generates forward-calculated DRS spectra as similar to the measured spectra as possible is solved by using the trust region reflective algorithm (Matlab R2012b, MathWorks Inc., Natick, Massachusetts).

All parameters except the speed parameters, i.e., nine parameters, are determined by solving the inverse problem of fitting DRS spectra. The spectral interval used in this study is 450 to 850 nm, and 32 wavelengths within this interval are chosen. The wavelengths are chosen more densely where the absorption spectra of hemoglobin display a highly dynamic shape, such as around the absorption peaks at 520 to 600 nm. The resulting optimization problem is given by

Eq. (18)

\min_{x} F_{obj, DRS} (x) {subject to c}_{i} (x) \leq 0, i \in I,

where

F_{obj, DRS} (x)

is the objective function of the nine parameters

x

, and

I

is a set of inequality constraints determined by lower and upper bounds for each of the nine parameters. Furthermore, the objective function is given by

Eq. (19)

F_{obj, DRS} (x) = ‖ \begin{matrix} F_{intensity} (x) \\ F_{shape} (x) \\ F_{q} (x) \\ F_{parameters} (x) \end{matrix} ‖_{2}^{2},

where

{‖ \cdot ‖}_{2}

denotes the Euclidian norm. The first part,

F_{intensity} (x)

, represents a weighted intensity difference between the calculated spectra from the model with parameters

x

and the measurement at both source–detector separations with emphasis for wavelengths at the hemoglobin absorption peaks between about 520 and 600 nm. The second part,

F_{shape} (x)

, represents some shape properties within the same wavelength interval, e.g., the relative amplitude difference between the absorption peaks and valleys. The third part,

F_{q} (x)

, represents the deviation from unity in the parameter

q

[see Eq. (20)]. The last part,

F_{parameters} (x)

, serves for penalizing unwanted behavior in the model parameters

x

directly including a very uneven distribution of blood between the upper and lower dermis layers and a too small vessel diameter (

< 15 μ m

).

The spectral fitting is performed in three steps. First, the epidermal thickness, the three scattering parameters, the melanin fraction, and the average blood tissue fraction are fitted (preliminary) to the spectra for wavelengths above 700 nm, in order to quickly and robustly find appropriate starting points for the second and third optimization steps. The rationale for this is that hemoglobin absorption is relatively weak for wavelengths above 700 nm and therefore has a rather low impact on the spectra here. Thus, the details about the hemoglobin absorption (difference in fraction between layers, oxygen saturation, and vessel diameter) are unnecessary to fit in this step. This is done for several random starting points of these six parameters, and the best candidates (those with best fit above 700 nm but with essentially different optimal values of the six parameters) are chosen for the second step.

In the second step, the average blood tissue fraction, the relative blood tissue fraction, the oxygen saturation, and the average vessel diameter are fitted. The epidermal thickness, the three scattering parameters, and the melanin fraction are locked to the values found in the first step. The average blood tissue fraction value has a starting point that is the result from the first step, and random starting points are chosen for the relative blood tissue fraction, the oxygen saturation, and the average vessel diameter.

In the third step, all nine parameters are fitted to the whole spectrum with starting points from the best candidates in the second step. The point with the smallest error is then regarded as the global minimum when this point is repeatedly obtained. The convergence properties are improved when performing the optimization in three steps like this.

Before calculating the objective function, the measured $[I_{meas, ρ} (λ)]$ and modeled $[I_{model, ρ} (λ)]$ spectra at the two source–detector separations ( $ρ \in [0.4, 1.2]$ ) are normalized according to

Eq. (20)

I_{meas, 0.4} (λ) \leftarrow q I_{meas, 0.4} (λ), where q = \frac{{〈 I_{model, 0.4} 〉}_{λ > 700} / {〈 I_{model, 1.2} 〉}_{λ > 700}}{〈 I_{meas, 0.4} 〉_{λ > 700} / 〈 I_{meas, 1.2} 〉_{λ > 700}},

where

〈 \cdot 〉

denotes the average, and then

Eq. (21)

I_{m, ρ} (λ) \leftarrow \frac{I_{m, ρ} (λ)}{{〈 {〈 I_{m, ρ_{0}} (λ) 〉}_{λ > 700} 〉}_{ρ_{0} \in [0.4, 1.2]}}, m \in [model, meas], ρ \in [0.4, 1.2] .

As mentioned above, the normalization in Eq. (21) allows for avoiding an absolute calibration of the spectrometers, and the introduction of the factor $q$ in Eq. (20) releases the necessity of an exact relative intensity calibration between the two source–detector separations. The factor $q$ is penalized when deviating from unity in order to avoid multiple solutions.

2.3.2.

Fitting LDF spectra

The 10 speed parameters, $p_{v, j}$ ( $j \in [0, \dots, 9]$ ) are fitted to the LDF spectra. The other nine model parameters have, in this step, fixed values obtained from fitting the DRS spectra. Each of the 10 speed parameters consists of speeds between 0 and $2 v$ , where $v$ is the average speed of each speed parameter evenly distributed in the logarithmic plane between 0.2 and $75 mm / s$ .

The resulting optimization problem is given by

Eq. (22)

\min_{p_{v}} F_{obj, LDF} (p_{v}) subject to p_{v, j} \geq 0, \sum_{j} p_{v, j} \leq 1, j \in [0, 9],

where

p_{v} = {[p_{v, 0}, \dots, p_{v, 9}]}^{T}

, and the objective function has the form

Eq. (23)

F_{obj, LDF} (p_{v}) = {‖ F_{diff} (p_{v}) ‖}_{2}^{2} .

The inner part, $F_{diff} (p_{v})$ , of the objective function is given by $F_{diff} (p_{v}) = {W [F_{M 1, model} (p_{v}) - F_{M 1, meas}]}$ , where $W$ is a diagonal weight matrix. Here, $F_{M 1, model} (p_{v})$ is given by

Eq. (24)

F_{M 1, model} (p_{v}) = [\begin{matrix} \int_{f_{0}}^{f_{1}} f y_{D, model} (f, p_{v}) d f \\ ⋮ \\ \int_{f_{i_{\max} - 1}}^{f_{i_{\max}}} f y_{D, model} (f, p_{v}) d f \\ \int_{f_{0}}^{f_{i_{\max}}} f y_{D, model} (f, p_{v}) d f \end{matrix}],

where

y_{D, model} (f, p_{v})

is the Doppler power spectrum calculated from the model and

Eq. (25)

f_{i} = {(f_{ϵ} + i Δ f)}^{2},

where

f_{ϵ}

is a low frequency above 0, and

Δ

is a frequency step that depends on the sampling frequency and the number of points used in the FFT. In analogy,

F_{M 1, meas}

is given by replacing

y_{D, model} (f, p_{v})

with

y_{D, meas} (f)

in Eq. (24). Note that the last element of the vectors

F_{M 1, meas}

and

F_{M 1, model}

equals the conventional perfusion estimate, i.e., the first-order moment of the Doppler power spectrum. The other elements represent the first-order moment in certain frequency intervals.

We observed that $f_{M 1, model} (p_{v})$ changes approximately linearly with the 10 speed parameters $p_{v, j}$ , i.e.,

Eq. (26)

F_{M 1, model} (p_{v}) + J (p_{v}) δ = F_{M 1, model} (p_{v} + δ) + ε,

where

J (p_{v})

denotes the Jacobian matrix, i.e., the derivative of each element in

F_{M 1, model} (p_{v})

with respect to

p_{v}

for each element in

p_{v}

. Here,

δ

is a change in

p_{v}

, and

ε

is the deviation from the linear approximation. Based on this, the optimization problem is solved by iterating the following steps:

1. Choose an initial $p_{v}$ ;
2. Calculate $F_{M 1, model} (p_{v})$ ;
3. Approximate the Jacobian at $p_{v}$ using finite differences;
4. Use constrained regression for finding $δ$ , which approximately solves the optimization problem;
5. Calculate $F_{M 1, model} (p_{v} + δ)$ ;
6. Calculate $ε$ from Eq. (26);
7. If $ε$ is small enough, STOP;
8. Else, use Broyden’s update to estimate the Jacobian at $p_{v} + δ$ ;
9. Go to step 4.

This iterative procedure is based on a special case of Broyden’s quasi-Newton method with a full-step size (i.e., no line search) in each iteration.²⁶ The full-step size is possible due to the close to linear behavior of $F_{M 1, model} (p_{v})$ . In Broyden’s method, the Jacobian is successively approximated over the step $δ$ with use of the updating formula

Eq. (27)

\hat{J} (p_{v} + δ) = \hat{J} (p_{v}) + \frac{[F_{M 1, model} (p_{v} + δ) - F_{M 1, model} (p_{v}) - J (p_{v}) δ] δ^{T}}{δ δ^{T}} .

In the proposed algorithm, the time-consuming part is to calculate the spectra from the model, which is done once in step 2, 10 times in step 3, and once in step 5. Since only step 5 is iterated among the expensive steps, this optimization procedure is known to be very fast, which is also the case in practice as it converges to the global minimum on average in about four iterations from a randomly generated initial $p_{v}$ .

2.4.

Model Output

When the optimal model is identified, some output parameters are directly obtained from that model. The quantities that will be considered in detail in this study are the RBC tissue fraction, the RBC oxygen saturation, and the perfusion (RBC tissue fraction times speed) in three different speed regions (0 to $1 mm / s$ , 1 to $10 mm / s$ , and above $10 mm / s$ ). The oxygen saturation is given directly by one model parameter, whereas the RBC tissue fraction and perfusion may be different in the two dermal layers and must therefore be averaged depending on the actual sampling volume.

The sampling volume used to calculate the output parameters is generated by randomly choosing at least 25,000 points among the trajectories of each detected photon in the Monte Carlo simulations of the fitted three-layered model. This is done for all wavelengths and source–detector separations included in the spectral fit and results in a point cloud. Based on the number of points located in each layer, the desired quantities are calculated. The same point cloud is applied on the evaluation simulations described in the next section in order to calculate the “true” parameters in those models. In this way, the accuracy of the output parameters from the fitted model can be determined.

2.5.

Evaluation Simulations

In order to evaluate the accuracy of the proposed method, 250 models of complex geometry, including an epidermis layer and a dermis layer containing 150 to 2000 individual blood vessels of various orientation, diameter (6 to 400 μm), oxygen saturation (0% to 100%), and speed distribution (0 to $200 mm / s$ ), were generated using a number of probability distributions. The models were generated similarly to the models described in Ref. 20. The average and standard deviation of the diameter varied between the 250 models (mean 12 to 50 mm, standard deviation 5 to 43 mm), that of the oxygen saturation (29% to $66 % \pm 15 %$ to 32%), and that of the speed distribution (0.005 to $7.4 \pm 0.07$ to $20 mm / s$ ). The epidermis thickness ranged from 7 to 220 μm and the melanin fraction from 0.1% to 17%. The reduced scattering coefficient was different in the epidermis and in the dermis layers and followed the shape given by Eq. (1). In the epidermis layer, it ranged from 1.2 to $11 {mm}^{- 1}$ with a difference between 450 and 850 nm ranging from 1.6 to $5.5 {mm}^{- 1}$ . Corresponding numbers for the dermis layer was 0.9 to $10 {mm}^{- 1}$ and 1.8 to $6.8 {mm}^{- 1}$ , respectively.

Photon transport in these models was then Monte Carlo simulated with OPs valid at 33 selected wavelengths (32 for DRS and 1 wavelength for LDF) from 450 to 850 nm at the two source–detector separations 0.4 and 1.2 mm. The DRS and LDF spectra were generated from the detected photons at the two source–detector separations for DRS and one source–detector separation for LDF, and the three-layered model was fitted to the resulting spectra. In this way, the accuracy could be evaluated for tissue-like models.

For the 32 DRS simulations and for each model and each source–detector separation, enough photons were simulated in order to achieve a signal-to-noise ratio (SNR) of at least 100, where SNR was approximated by

Eq. (28)

SNR = \frac{〈 W 〉}{std (W)} \sqrt{n}, W_{j} = \sum_{i = \frac{j N}{n}}^{i = \frac{(j + 1) N}{n} - 1} w_{i}, j \in [0, \dots, n - 1],

where

w_{i}

denotes the weight of detected photon

I

, and

n

was set to 10. For the LDF simulations, at least 180,000 photons were detected.

2.6.

In Vivo Measurements

Measurements were performed on a healthy male with Caucasian skin of age 37. Two different measurements were conducted on the volar side of the lower forearm: a systolic and a venous occlusion, and one measurement was performed on the volar surface of the foot: a heating provocation where the tissue was heated to 44°C locally. The two occlusions lasted for 5 min, followed by a 5-min reperfusion phase, and the heating provocation for 25 min, all during which the subject was sitting comfortable with arms and legs fully rested. All provocations were preceded by a 5-min baseline recording. This study was approved by the regional ethics committee of Linköping (D. no M83-09).

The measurements were performed using a custom-made fiber-optic probe with two emitting fibers and nine detecting fibers: two detecting fibers connected to one spectroscope each (AvaSpec-ULS2048L, Avantes BV, The Netherlands) and six detecting fibers connected to a single detector in a modified Periflux 5000 system (Perimed AB, Järfälla, Sweden). The emitting fibers consisted of one fiber delivering light from a broadband white-light source (Avalight-HAL-S, Avantes BV, Apeldoorn, The Netherlands) and the other delivering light from a laser light source (780 nm). The multiple-detecting fibers for the DRS system were placed at a distance of 0.4 and 1.2 mm from the white-light illuminating fiber, and the six detecting fibers for the LDF system were placed at a distance of 1.2 mm to the laser light source. The fibers used in the probe were made of fused silica with a radius of 100 μm and a numerical aperture of 0.37. A notch filter to attenuate wavelengths of $790 \pm 20 nm$ was added between the probe and the two spectroscopes to ensure minimal influence from the laser light on the DRS spectra. The DRS spectra were analyzed in the 450- to 850-nm wavelength interval. During the measurements, the probe was placed in a thermostatic probe holder (PF 450, Perimed AB, Järfälla, Sweden) to ensure a stable skin temperature and a good tissue contact.

The spectrometers were calibrated in three steps as previously described.²⁰ First, a dark reference spectrum was subtracted. Then, the spectrum was normalized to a white reference spectrum, and finally, a relative calibration between the two detector channels was performed by normalizing the spectrum with the average intensity of a calibration measurement, where the fibers at both source–detector separations were evenly illuminated.

The laser Doppler units were calibrated by measuring noise levels at various DC levels, subtracting the estimated noise level for the DC level at the time of the measurement, and normalizing with the frequency characteristics of the noise that can be assumed to be white. Finally, the spectra were normalized with the first-order moment of a spectrum originating from a measurement in motility standard. This calibration routine is described in detail elsewhere²⁷ and results in intensity and frequency-calibrated spectra that can be compared with the spectra calculated from the model.

2.7.

DRS Analysis for Comparison

For comparison, an existing state-of-the-art DRS analysis algorithm based on a modified Beer–Lambert’s law expression, according to Bargo et al.,²⁸ was implemented. In that method, an absolute calibration was performed using a grid of phantoms with varying $μ_{a}$ and $μ_{s}^{'}$ covering the range of expected values. For each level of $μ_{s}^{'}$ , three parameters ( $C_{1}$ , $L_{1}$ , and $C_{2}$ ) were fitted to all $μ_{a}$ levels in the intensity grid $I (μ_{a}, μ_{s}^{'})$ in the following modified Beer–Lambert expression

Eq. (29)

I (μ_{a}, μ_{s}^{'}) = C_{1} (μ_{s}^{'}) \exp [- μ_{a} L_{1} (μ_{s}^{'})] + C_{2} (μ_{s}^{'}) .

The parameters $C_{1}$ , $L_{1}$ , and $C_{2}$ were then fitted to polynomials of orders 4, 15, and 15, respectively. For more details, see the original paper.²⁸

Bargo’s method was essentially used to calculate the forward problem in an alternative way. The original method was, however, slightly adapted to better reflect the setting in the evaluation simulations using the same source–detector separations (one at a time), same wavelengths, same scattering model [Eq. (1)] and absorbers, and the same model for describing the vessel packaging [Eq. (5)]. The model parameters describing the epidermis thickness and the relative difference in blood tissue fraction between the two layers were excluded, as the original Bargo method is based on a homogeneous medium assumption. Thus, seven free parameters remained to fit to the measured spectra in this alternative method. For comparison, the same fitting strategy, as outlined in Sec. 2.3.1, was used.

Before comparing the two methods, the Bargo algorithm was calibrated using a range of Monte Carlo-simulated homogeneous phantoms covering the range of OPs found in the 250 evaluation simulations. When calculating the true RBC tissue fraction from the simulated models, the sampling volume for the current source–detector separation was considered, in contrast to our proposed method where the average sampling volume from both source–detector separations was considered.

3. Results

3.1.

Evaluation Simulations

The average relative root mean square (RMS) error in the fit of the DRS spectra for the 250 evaluation simulations was 0.9% for the 0.4-mm fiber separation and 2.0% for the 1.2-mm fiber separation. The average relative RMS error in the fit for LDF spectra for frequencies [Eq. (25)], where the power exceeded $10^{- 5} {Hz}^{- 1}$ , was 4.5%. These RMS errors indicate rather good spectral fits for both DRS and LDF for the 250 simulated models. An example for the simulated model with the median DRS fit [Eq. (18)] is shown in Fig. 1 and with the median LDF fit in Fig. 2.

Fig. 1

Simulated (thick black) and fitted (thin yellow) diffuse reflectance spectroscopy (DRS) spectra at 0.4 (a) and 1.2 (b) mm source–detector separations. The example shown is the one of the 250 models with the median DRS fit.

Fig. 2

Simulated (thick black) and fitted (thin yellow) laser Doppler flowmetry (LDF) spectra at source–detector separation 1.2 mm. The example shown is the one of the 250 models with the median LDF fit.

The RMS deviation between estimated and true RBC tissue fraction ${[〈 {({RBC}_{est} - {RBC}_{true})}^{2} 〉]}^{1 / 2}$ in the actual sampling volume in the evaluation simulations was 0.099 percentage units, which should be compared with an average RBC tissue fraction in the 250 simulated models of 0.56%. When only considering fractions below 0.5%, the absolute RMS deviation was 0.058 percentage units. The relative RMS deviation between the estimated and true RBC tissue fraction ${[〈 {({RBC}_{est} / {RBC}_{true} - 1)}^{2} 〉]}^{1 / 2}$ was 19% (16% for fraction above 0.5%). The absolute RMS deviation of the oxygen saturation was 4.1 percentage units (4.9 and 2.7 percentage units for RBC tissue fractions below and above 0.5%, respectively), which should be compared with an average oxygen saturation of 45% in the 250 simulated models. Scatter plots of true and estimated RBC tissue fraction and oxygen saturation are shown in Fig. 3.

Fig. 3

Scatter plots of estimated versus true red blood cell (RBC) tissue fraction (a) and oxygen saturation (b, open circles—RBC tissue fractions below 0.5%, filled circles—RBC tissue fractions above 0.5%).

The relative RMS deviation between estimated and true total perfusion was 23%. By normalizing the conventional perfusion estimate with the mean of the true perfusion for perfusion values below 0.1% $RBC \times mm / s$ , the accuracy of the normalized conventional perfusion can also be calculated. This deviation was 49%. For the perfusion for the three individual speed regions (0 to 1, 1 to 10, and above $10 mm / s$ , respectively), the relative RMS deviation in relation to the total perfusion ${〈 [{({perf}_{v, est} - {perf}_{v, true}) / {perf}_{tot, true}]}^{2} 〉}^{1 / 2}$ was 7.7%, 17%, and 18%, respectively. Scatter plots of true and estimated perfusions are shown in Fig. 4.

Fig. 4

Scatter plots of true total perfusion versus estimated and normalized conventional perfusions (a, black dots—model-based estimation, red crosses—normalized conventional perfusion) and perfusions for the individual speed regions (b).

The estimated average vessel diameter is a quantity which could potentially be clinically very interesting. This quantity also affects how accurately the other quantities, especially the RBC tissue fraction and perfusion, are estimated. The absolute RMS deviation of the average vessel diameter was 10 μm (8.3 μm for RBC tissue fractions above 0.5%), which should be compared with an average value of 27 μm in all the 250 models. There was a strong correlation between the accuracy in estimated vessel diameter and accuracy in estimated RBC tissue fraction and total perfusion (correlation coefficients of 0.83 and 0.75, respectively). Only a weak correlation was seen to the accuracy of the estimated RBC oxygen saturation. Omitting the vessel-packaging effect, i.e., setting that model parameter to 0, leads to a worsened spectral fit of DRS spectra and a general underestimation of RBC tissue fraction and perfusion of 34% and 37%, respectively, for the 250 models.

The correlation between estimated and true values for a number of estimated model parameters and their RMS deviations are given in Table 1.

Table 1

Accuracy of estimated model parameters for all evaluation models and those with a RBC tissue fraction (t.f.) >0.5%.

				All 250 models		RBC t.f. $> 0.5 %$ ( $n = 137$ )
	Correlation coefficient	Average	Range	Absolute RMS	Relative RMS (%)	Absolute RMS	Relative RMS (%)
RBC tissue fraction (%)	0.98	0.56	0.06 to 3.2	0.099	19	0.13	16
Oxygen saturation (%)	0.94	45	21 to 75	4.1	9.6	2.7	6.7
Total perfusion (% $\times mm / s$ )	0.96	0.61	$8 \times 10^{- 4}$ to 6.7	0.25	23	0.35	20
Perf. $< 1.0 mm / s$ (% $\times mm / s$ )a	0.93	0.088	$8 \times 10^{- 4}$ to 0.5	0.037	7.7	0.047	6.9
Perf. 1 to $10 mm / s$ (% $\times mm / s$ )a	0.90	0.27	0 to 1.6	0.15	17	0.19	16
Perf. $> 10 mm / s$ (% $\times mm / s$ )a	0.96	0.25	0 to 5.5	0.16	18	0.22	13
Vessel diameter (μm)	0.60	27	12 to 81	10	40	8.3	31
Average ( ${mm}^{- 1}$ )	0.89	3.3	2.0 to 5.1	0.31	11	0.32	11
Relative RBC between layers (-)	0.72	0.0	$- 1.0$ to 1.0	0.33	b	0.29	b
Melanin content (% $\times mm$ )	0.96	0.29	$4 \times 10^{- 3}$ to 1.8	0.11	b	0.13	b
Epidermis thickness (μm)	0.32	79	7.1 to 218	81	b	62	b

^aRelative values given relative to total perfusion.

^bNot relevant because values close to zero cause large relative differences.

The relaxation factor $q$ [Eq. (20)] was $0.998 \pm 0.043$ ( $mean \pm standard deviation$ ) for the 250 simulated models. For six of the models, $q$ was below 0.9. In Eq. (21), the normalization factor was $4.4 % \pm 7.6 %$ lower for the simulated spectra than for the fitted spectra. The simulated models were also fitted in an absolute manner, i.e., by omitting Eq. (21) and/or by omitting the relaxation factor $q$ . Interestingly enough, the accuracy of the estimations of the oxygen saturation, the RBC tissue fraction, as well as the total perfusion decreased when omitting any or both of these equations. The decrease was significant in all cases (paired $t$ -test, $p < 0.02$ ), except for the RBC tissue fraction, when only omitting Eq. (20), i.e., the relaxation parameter $q$ .

3.2.

Evaluation Simulations with Alternative DRS Analysis

The Bargo method for analyzing DRS spectra, outlined in Sec. 2.7, was applied on the 250 evaluation simulations, with one source–detector separation at a time. As with our proposed method, a good spectral fit (average RMS deviation of 5.7% and 3.9% for 0.4 and 1.2 mm source–detector separation, respectively) was achieved. The relative RMS deviations between true and estimated RBC tissue fraction were 123% and 52% for the source–detector separations of 0.4 and 1.2 mm, respectively. This can be compared with 19% for our proposed method. The absolute RMS deviations for the oxygen saturation were 18 and 6.0 percentage units, compared with 4.1% for our method. Corresponding numbers for the average vessel diameter were 40 and 25 μm, compared with 10 μm for our method. Our proposed method was significantly more accurate in all cases except for the oxygen saturation at 1.2 mm source–detector separation ( $p < 0.001$ ; Wilcoxon’s matched pairs test).

3.3.

In Vivo Measurements

The RBC tissue fraction, oxygen saturation, and perfusion during the heat provocation are plotted as functions of time in Fig. 5. The heating started at 5 min and lasted throughout the measurement. Examples of the spectral fit at 5 and 25 min, respectively, are shown in Fig. 6. Note especially the correlation between the oxygen saturation and the total perfusion and the somewhat different characteristics of the three speed-resolved perfusion estimates.

Fig. 5

Time-resolved RBC tissue fraction and oxygen saturation (a) and speed-resolved perfusion (b) during a heat provocation on skin (foot).

Fig. 6

Example of spectral fit during heat provocation for DRS spectra at 0.4 mm source–detector separation (a), DRS spectra at 1.2 mm source–detector separation (b), and LDF spectra (c) at 5 (before the start of provocation) and 25 (near the end of provocation) min. Black thick curves correspond to measured spectra, and yellow thin curves correspond to fitted (modeled) spectra.

The results of the venous occlusion are shown in Fig. 7. The occlusion started at 5 min and lasted for another 5 min. Note the reduction in perfusion although the RBC tissue fraction is doubled during the provocation.

Fig. 7

Time-resolved RBC tissue fraction and oxygen saturation (a) and speed-resolved perfusion (b) in forearm skin during a venous occlusion.

The results of the systolic occlusion are shown in Fig. 8. The occlusion started at 5 min and lasted for another 5 min. Note that the oxygen saturation settles at 0% toward the end of the occlusion (the model allows for negative saturations) and the much more rapid drop in perfusion component for speeds above $10 mm / s$ compared with the component 1 to $10 mm / s$ in the reperfusion phase. The increase in RBC tissue fraction during the occlusion phase is likely a result from a regional redistribution of blood from the arterial (deeper) vessels to the venous (more superficial) vessels.

Fig. 8

Time-resolved RBC tissue fraction and oxygen saturation (a) and speed-resolved perfusion (b) in forearm skin during a systolic occlusion.

4. Discussion

The most important difference between the proposed model-based method and other attempts to combine DRS and LDF¹⁷^,¹⁸ is that the proposed method integrates not only the hardware (probe), but also the analysis of the measured spectra. In our analysis, the tissue model parameters fitted to the DRS spectra are used when fitting LDF spectra, which make the LDF analysis more accurate as scattering and vessel-packaging effects can be accounted for. Thus, the quantitative LDF measures from the method proposed in this article are more accurate than in the method we have previously presented, where we made use of LDF spectra at two source–detector separations but without the support from DRS measurements.¹⁵ It is possible to use two source–detector separations for the LDF measurements. This enables the use of a model parameter differentiating the flow speeds in the two dermal layers. In contrast to the previous attempts,¹⁵ two LDF source–detector separations are not needed in the combined method in order to avoid multiple solutions to the inverse problem. Therefore, only one separation was used in this article.

Compared with methods which are not based on an adaptive tissue model, our proposed method has several advantages. It enables determination of absolute quantities of RBC tissue fraction and oxygen saturation, by accounting for scattering and geometrical effects such as layered structures and vessels. It also enables quantitative speed-resolved estimations of the RBC flow by feeding the LDF analysis algorithm with the output from the DRS analysis. Other advantages include the possibility to give estimates on the uncertainty of the estimated quantities by analyzing the effect on the objective function by a small change in the model parameters or to reject the measurements where an acceptable model fit cannot be found. An estimate of the sampling volume can also be presented.²⁹

The importance of the vessel-packaging effect for accurate estimations of the output parameters was outlined in Sec. 3.1, which states that ignoring this effect in general causes about 35% underestimation of the RBC tissue fraction and perfusion. The reason why no correlation was found between the accuracy of the vessel diameter and oxygen saturation estimations is probably due to a too small range of oxygen saturation. Generally, the oxygen saturation level was between 25% and 60% in the 250 models. For oxygenation levels lower than 15%, the results from ignoring the vessel-packaging effect indicated a general overestimation of about 5 percentage units, and a general underestimation of about 10 percentage units was indicated for saturation levels above 70%. The latter is well in line with previous results.⁵

While the vessel-packaging effect and its importance have previously been thoroughly examined for DRS,⁵^,⁷^–⁹^,³⁰ it has previously not been accurately described for LDF. In DRS, the vessel-packaging effect is an effect of inhomogeneous light absorption, but for LDF, the main reason for the vessel-packaging effect is that consecutive Doppler shifts are caused by RBC’s moving with similar speed (similar in relation to the speed distribution in the whole tissue). Furthermore, larger vessels imply that less light will be Doppler shifted for a constant total RBC tissue fraction, but the light that is Doppler shifted will be more multiple shifted, and this also affects the Doppler power spectrum. An even better approximation of the vessel diameter would be beneficial for the proposed method.

The adaptive tissue model contains in total 19 free parameters. We have previously shown that all nine parameters that are fitted to the DRS spectra are needed in order for a good spectral fit at both source–detector separations, i.e., that all nine parameters contribute to a significant decrease of the objective function.²⁰ It is possible that choosing different geometrical properties could improve the spectral fit further. For example, allowing a variable thickness of the upper dermis layer, instead of the epidermis, could enhance the results, as indicated by the rather poor correlation between estimated and true epidermis thickness (Table 1).

The 10 speed parameters describe a speed distribution, where a high degree of cross-talk between neighboring parameters can be expected. A reduction in the number of speed parameters could probably reduce over-fitting issues for the 10 parameters. However, this should not be regarded as crucial, since the output perfusion is condensed into only three different flow speed regions.

A possible expansion of the model is to allow for an adaptive exponent in the absorption coefficient for melanin [Eq. (3)]. Reported values on this exponent are inconsistent, and it will depend on the fraction of eumelanin in relation to pheomelanin.³¹ Changing this value from $- 3 / 10$ to $- 2$ or $- 5$ in the three-layered adaptive model increased the value of the objective function at the solution with about 75% (not presented in Sec.3), which indicates that it may be a good idea to let this parameter be variable in the adaptive model when used on measured spectra. The importance of this should be evaluated in a larger in vivo study.

Water is often included as a chromophore³¹^,³² in tissue models similar to the one presented here, because water can be a prominent absorber for wavelengths above 750 nm. We have chosen not to include water, since it has a negligible influence on the backscattered intensity in the wavelength interval (450 to 850 nm) and source–detector separations (0.4 and 1.2 mm) used. Similarly, various “yellow” chromophores such as bilirubin, methemoglobin, and beta-carotene may also be included as they can have a significant effect on the backscattered intensity below about 500 nm.³¹ We have chosen not to include such chromophores, but instead reduce the influence on the objective function for overestimations of the backscattered light intensity for wavelengths below 500 nm. The effect of this relaxation in the objective function can be noticed in Figs. 1(a) and 6(b).

In the spectral interval 450 to 650 nm, both oxygenated and deoxygenated hemoglobin have a strong and unique characteristic footprint. The interval 650 to 850 nm, where tissue absorption is minimal, strengthens the possibility of determining the level of scattering throughout the complete spectral range. These are important reasons for choosing the wavelength interval 450 to 850 nm in the DRS part of this method. Furthermore, since we primarily want to study the microcirculation, it is an advantage to use the visible wavelength region compared with the near-infrared region, in order to reduce the sampling volume. The reduced sampling volume is also the reason for choosing the rather small, but well separated, source–detector separations of 0.4 and 1.2 mm.

The design of the objective functions and the order of search for the optimal fit (three steps in global search for DRS) are important in order to achieve good convergence properties when solving the inverse problem. A large part of the development of the method has focused on this. In the three-step search strategy, approximate values for the parameters that are mainly affecting the spectra above 700 nm are found relatively quickly in the first step. In the second step, those parameters are locked using their approximate values from step 1, while the other parameters are fitted to the entire spectral range. By locking these parameters in step 2, potential raising conditions, caused by starting points located far from their optimum, are avoided. In the third step, all parameters, which are now likely to be located relatively close to their optimal values, are freely fitted a last time. In the first and the second steps, multiple random starting points are used for the fitting parameters until the same best optimal point is found several times. This is necessary in order to be sure that the global optimum is found, but is also the reason that a global search is rather slow. If all nine parameters were to be found simultaneously in a single step with multiple starting points for all parameters, many more random starting points would be needed (one to two orders of magnitude more), which would increase the time needed to find the global minimum.

Broyden’s method was used for updating the Jacobian in each iteration. This has a high impact on the time needed for solving the inverse problem compared with using partial derivatives which requires calculating the forward problem several times. For Broyden’s method to work smoothly, prescaling of the model parameters has to be performed, so that a change in each parameter affects the objective function to a similar extent.²⁶ For LDF, the prescaling was done by normalizing each speed parameter with its average speed. For DRS, empirical scaling factors were used.

The use of Eq. (21) leads to an important simplification of the calibration routine, as no absolute calibration is required, whereas Eq. (20) limits the effect of inaccuracies in the relative calibration. This is important since especially an absolute calibration and also an exact relative calibration are difficult to perform in practice. The absolute levels in real measurements may also be largely influenced by the contact between the probe and the skin, which would influence the results without these relaxations.

The accuracy of the estimated quantities is, in general, improved by applying a calibration relaxation according to Eqs. (20) and (21), even for the simulated evaluation simulations that are absolute calibrated by nature. The fact that the accuracy and, obviously, the goodness of fit are improved by introducing Eqs. (20) and (21) indicates that they, to some extent, compensate for inhomogeneities that cannot accurately be accounted for by the three-layered model. Although deviations from unity in the factor $q$ are penalized by the objective function, it settled at values below 0.9 for 6 of the 250 simulations. Common for these cases is a high concentration of blood vessels located directly in front of the light-emitting fiber or in a layer that is much thinner than 0.5 mm (the thickness of the upper dermis layer in the three-layered adaptive model). Modifications to the adaptive three-layered model, such as introducing an adaptive parameter that controls the thickness of the upper dermal layer, may avoid this behavior in $q$ . It should also be noted that the introduction of a relaxation parameter can result in an impaired estimation of the scattering properties in contradiction to the general improvement seen in other parameters, especially when $q$ deviates too much from unity.

For DRS, a variety of methods have been presented, most of them capable of estimating relative or absolute changes in RBC tissue fraction (or equivalent) and/or the absolute oxygen saturation.⁵^,⁶^,¹⁷^,²⁸^,³³^,³⁴ We have performed a comparison to one of these methods,²⁸ which is in many senses a state-of-the-art DRS method based on an absolute calibrated modified Beer–Lambert’s law model. For the longer source–detector separation (1.2 mm), the accuracy in RBC oxygen saturation was approximately the same as for our method. The accuracy in RBC tissue fraction and average vessel diameter was much worse in that method compared with ours. This indicates that for estimating the oxygen saturation only, a relatively simple method may be good enough, but in order to accurately estimate absolute RBC tissue fractions and many other parameters, a more complex method, as our proposed method, is necessary. Our method also has the advantage in that it can be combined with LDF and that it does not require a cumbersome absolute calibration which is needed for many other methods.

For LDF, on the other hand, all commercial products use the first moment of the Doppler power spectrum to calculate the perfusion index. By normalizing that measure, we were able to compare our results with the conventional measure [Fig. 4(a)]. The known nonlinearity in the conventional measure can clearly be seen, but except from this, the “noise” to the true values is similar for our model-based method and the conventional method. Nevertheless, other previously mentioned advantages with the model-based approach including the quantitative and speed-resolved nature of the estimated perfusion are important.

In this study, simulations of complex models containing discrete blood vessels were used to evaluate the accuracy of the proposed method. This evaluation indicates certain accuracy, but it is not a hard proof of that the method gives accurate estimates for in vivo measurements. As a gold standard does not exist for measurements in the microcirculation, proving the accuracy of the method is cumbersome. We have previously shown that this type of model-based approach works well for measurements in bio-optical phantoms both for LDF³⁵^,³⁶ and DRS.³⁷ Those phantom measurements constitute the first link in the chain of proof for the model-based method. The second link is the accuracy shown by the evaluation simulations in Sec. 3.1 in this article. The third link is merely the fact that a single model that fits measured spectra at two source–detector separations for DRS and one source–detector separation for LDF well is possible to find, as shown in Sec. 3.2. Some of the estimated parameters in the fitted model can also be validated by comparing them with the results obtained using other bio-optical measurement techniques,³⁸^,³⁹ which can further strengthen the trustworthiness of the proposed method. Finally, expected behavior of the output parameters during in vivo measurements serves as an indicator that the method works accurately. An example of such expected behavior can be seen in Fig. 8(a), where the oxygen saturation settles at the zero level at the end of the occlusion.

Methods based on adaptive tissue models are more complex than comparable conventional methods based on direct calculations on the measured spectra. The iterative fitting process is also time consuming and may require recorded spectra with less noise and more accurate calibration processes. The latter issue is solved in the presented method by the relaxations in Eqs. (20) and (21) for DRS and by the previously presented calibration process²⁷ for LDF. The time resolution issues, both for acquiring high-quality spectra and for solving the inverse problem, may be a limiting factor in some specialized high-dynamic measurement situations. Utilizing smart algorithms and hyper-parallel computing on graphic computing units have reduced the time to solve the inverse problem of both DRS and LDF on a standard laptop (NVIDIA Quadro NVS 160M with eight cores) to the order of 100 ms, reducing this limitation.

5. Conclusions

We have presented an adaptive model-based method for a combined analysis of DRS and LDF data. The accuracy of the method was evaluated on complex, tissue-like models containing a high amount of individual blood vessels. The absolute RMS deviation between estimated and true oxygenation was 4.1 percentage units, whereas the relative RMS deviations for the RBC tissue fraction and perfusion were 19% and 23%, respectively. Advantages compared with stand-alone conventional DRS and LDF methods include simultaneous RBC tissue fraction, oxygen saturation, and speed-resolved perfusion from the same sampling volume given in absolute quantitative units.

Appendices

Appendix:

Analytical Calculation of Single-Shifted Optical Doppler Spectrum

The Doppler frequency shift $f$ that results when light is scattered an angle $θ$ by an RBC moving with velocity $v$ can be expressed by

Eq. (30)

f = v \cdot q = \frac{2 n v}{λ} \sin (θ / 2) \cos φ = \frac{2 n v}{λ} \sqrt{\frac{1 - \cos θ}{2}} \cos φ = [\begin{matrix} μ = \cos θ \\ ξ = \cos φ \end{matrix}] = \frac{2 n v}{λ} \sqrt{\frac{1 - μ}{2}} ξ,

where

v = | v |

,

q

is the difference between the wave vectors

k_{i}

and

k_{s}

of the incident and scattered light waves, respectively,

λ

is the laser wavelength,

n

is the refractive index of the medium,

θ

is the scattering angle, and

φ

is the angle between

v

and

q

. The single-shifted optical Doppler spectrum

p_{o, 1} (f)

for a given set of

v

,

λ

, and

n

can thus be calculated by considering the probability distributions of the angles

θ

and

φ

. Specifically, the probability to be calculated for each

f

is the probability of

Eq. (31)

\sqrt{\frac{1 - μ}{2} ξ} = \frac{f λ}{\underset{A}{\underset{︸}{2 n v}}} \Rightarrow \sqrt{\frac{1 - μ}{2}} = x and ξ = \frac{A}{x},

for all values of

x

between

A

and 1, i.e., for

μ \in [- 1, 1 - 2 A^{2}]

. The probability of

μ

is given by the scattering phase function of blood (Gegenbauer kernel with

g_{Gk} = 0.948

and

α_{Gk} = 1

is used), and the probability of

ξ

can be assumed to be rectangularly distributed between 0 and 1 for positive frequency shifts. The probability density function for

ξ = A / x

is given by

Eq. (32)

p_{ξ} (x) = \frac{A}{x} .

The probability density function for $μ$ is given by the scattering phase function (Gegenbauer kernel) which with $α_{G k} = 1$ is

p_{μ} (μ) = \frac{K}{{(1 + g_{Gk}^{2} - 2 g_{Gk} μ)}^{2}}, where K = \frac{g_{Gk} {(1 - g_{Gk}^{2})}^{2}}{2 [{(1 + g_{Gk})}^{2} - {(1 - g_{Gk})}^{2}]} .

The probability density function $p_{o, 1} (f)$ is thus calculated by

Eq. (33)

p_{o, 1} = \frac{1}{f} \int_{- 1}^{μ_{0}} p_{μ} (μ) p_{ξ} (\frac{\sqrt{1 - μ}}{2}) d μ = \dots = [B = \frac{{(1 - g_{Gk})}^{2}}{4 g_{Gk}}] = \frac{K λ}{8 n v g_{Gk}^{2}} \int_{A}^{1} \frac{1}{{(B + x^{2})}^{2}} d x = \dots = \frac{K λ}{\underset{a}{\underset{︸}{2 n v g_{Gk} {(1 - g_{Gk})}^{2}}}} (\underset{b}{\underset{︸}{\frac{1}{(1 + B)} + \frac{1}{\sqrt{B}} \arctan \frac{1}{\sqrt{B}}}} - \frac{x}{\underset{c (f)}{\underset{︸}{(A^{2} + B)}}} - \underset{d (f)}{\underset{︸}{\frac{1}{\sqrt{B}} \arctan \frac{A}{\sqrt{B}}}}) .

In practice, the probability density function is rather wanted for a frequency bin $f_{i}$ to $f_{i + 1}$

Eq. (34)

p_{o, 1} (f_{i}) = a \int_{f_{i}}^{f_{i + 1}} [b - c (f) - d (f)] d f = \dots = \frac{K λ}{2 n v g_{Gk} {(1 - g_{Gk})}^{2}} {[f (\frac{1}{(1 + B)} + \frac{1}{\sqrt{B}} \arctan \frac{1}{\sqrt{B}} - \frac{1}{\sqrt{B}} \arctan \frac{λ f}{2 n v \sqrt{B}})]}_{f_{i}}^{f_{i + 1}} .

Note that this expression is only valid for a Gegenbauer kernel phase function with $α_{Gk} = 1$ . Note also that the expression is only valid for positive frequency shifts, but the shape is identical (mirrored) for negative frequency shifts.

Acknowledgments

This study was financed by VINNOVA and Perimed AB through the SamBIO research collaboration program between companies and academia within bioscience (VINNOVA D.no. 2008-00149) and the Research&Grow program (VINNOVA D.no. 2011-03074), and also by NovaMedTech, supported by the European Union Regional Development Fund. The authors also want to acknowledge PhD student Hanna Karlsson for valuable assistance during the in vivo measurements.

References

1.

A. G. TsaiP. C. JohnsonM. Intaglietta, “Oxygen gradients in the microcirculation,” Phys. Rev., 83 (3), 933 –963 (2003). PHREA7 0031-9333 Google Scholar

2.

V. VenugopalanJ. S. YouB. J. Tromberg, “Radiative transport in the diffusion approximation: an extension for highly absorbing media and small source-detector separations,” Phys. Rev. E, 58 (2), 2395 –2407 (1998). http://dx.doi.org/10.1103/PhysRevE.58.2395 PLEEE8 1063-651X Google Scholar

3.

J. L. HollmannL. V. Wang, “Multiple-source optical diffusion approximation for a multilayer scattering medium,” Appl. Opt., 46 (23), 6004 –6009 (2007). http://dx.doi.org/10.1364/AO.46.006004 APOPAI 0003-6935 Google Scholar

4.

G. M. PalmerN. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: theory and validation on synthetic phantoms: erratum,” Appl. Opt., 46 (27), 6847 –6847 (2007). http://dx.doi.org/10.1364/AO.46.006847 APOPAI 0003-6935 Google Scholar

5.

T. Lindberghet al., “Improved model for myocardial diffuse reflectance spectra by including mitochondrial cytochrome aa3, methemoglobin, and inhomogenously distributed RBC,” J. Biophotonics, 4 (4), 268 –276 (2011). http://dx.doi.org/10.1002/jbio.201000048 JBOIBX 1864-063X Google Scholar

6.

H. Karlssonet al., “Can a one-layer optical skin model including melanin and inhomogeneously distributed blood explain spatially resolved diffuse reflectance spectra?,” Proc. SPIE, 7896 78962Y (2011). http://dx.doi.org/10.1117/12.873134 PSISDG 0277-786X Google Scholar

7.

L. Svaasandet al., “Therapeutic response during pulsed laser treatment of port-wine stains: dependence on vessel diameter and depth in dermis,” Laser. Med. Sci., 10 (4), 235 –243 (1995). http://dx.doi.org/10.1007/BF02133615 LMSCEZ 1435-604X Google Scholar

8.

R. L. P. van VeenW. VerkruysseH. J. C. M. Sterenborg, “Diffuse-reflectance spectroscopy from 500 to 1060 nm by correction for inhomogeneously distributed absorbers,” Opt. Lett., 27 (4), 246 –248 (2002). http://dx.doi.org/10.1364/OL.27.000246 OPLEDP 0146-9592 Google Scholar

9.

I. FredrikssonM. LarssonT. Strömberg, “Accuracy of vessel diameter estimated from a vessel packaging compensation in diffuse reflectance spectroscopy,” Proc. SPIE, 8087 80871M (2011). http://dx.doi.org/10.1117/12.889684 PSISDG 0277-786X Google Scholar

10.

Q. WangK. ShastriT. J. Pfefer, “Experimental and theoretical evaluation of a fiber-optic approach for optical property measurement in layered epithelial tissue,” Appl. Opt., 49 (28), 5309 –5320 (2010). http://dx.doi.org/10.1364/AO.49.005309 APOPAI 0003-6935 Google Scholar

11.

D. Yudovsky, “Spatial frequency domain spectroscopy of two layer media,” J. Biomed. Opt., 16 (10), 107005 (2011). http://dx.doi.org/10.1117/1.3640814 JBOPFO 1083-3668 Google Scholar

12.

G. E. NilssonT. TenlandP. Å. Öberg, “Evaluation of a laser Doppler flowmeter for measurement of tissue blood flow,” IEEE Trans. Biomed. Eng., BME-27 (10), 597 –604 (1980). http://dx.doi.org/10.1109/TBME.1980.326582 IEBEAX 0018-9294 Google Scholar

13.

R. BonnerR. Nossal, “Model for laser Doppler measurements of blood flow in tissue,” Appl. Opt., 20 (12), 2097 –2107 (1981). http://dx.doi.org/10.1364/AO.20.002097 APOPAI 0003-6935 Google Scholar

14.

I. FredrikssonM. LarssonT. Strömberg, “Laser Doppler flowmetry,” Microcirculation Imaging, Wiley-Blackwell, Weinheim, Germany (2012). Google Scholar

15.

I. FredrikssonM. LarssonT. Strömberg, “Model-based quantitative laser Doppler flowmetry in skin,” J. Biomed. Opt., 15 (5), 057002 (2010). http://dx.doi.org/10.1117/1.3484746 JBOPFO 1083-3668 Google Scholar

16.

I. Fredrikssonet al., “Reduced arteriovenous 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). http://dx.doi.org/10.2337/db10-0080 DIAEAZ 0012-1797 Google Scholar

17.

H. LiuM. Kohl-BareisX. Huang, “Design of a tissue oxygenation monitor and verification on human skin,” Proc. SPIE, 8087 80871Y (2011). http://dx.doi.org/10.1117/12.889197 Google Scholar

18.

S. Beckertet al., “The impact of the micro-lightguide O2C for the quantification of tissue ischemia in diabetic foot ulcers,” Diabetes Care, 27 (12), 2863 –2867 (2004). http://dx.doi.org/10.2337/diacare.27.12.2863 DICAD2 0149-5992 Google Scholar

19.

K. H. Franket al., “The Erlangen micro-lightguide spectrophotometer EMPHO I,” Phys. Med. Biol., 34 (12), 1883 –1900 (1989). http://dx.doi.org/10.1088/0031-9155/34/12/011 PHMBA7 0031-9155 Google Scholar

20.

I. FredrikssonM. LarssonT. Strömberg, “Inverse Monte Carlo method in a multilayered tissue model for diffuse reflectance spectroscopy,” J. Biomed. Opt., 17 (4), 047004 (2012). http://dx.doi.org/10.1117/1.JBO.17.4.047004 JBOPFO 1083-3668 Google Scholar

21.

I. FredrikssonM. LarssonT. 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). http://dx.doi.org/10.1117/1.2854691 JBOPFO 1083-3668 Google Scholar

22.

S. L. Jacques, Skin Optics, (1998) http://omlc.ogi.edu/news/jan98/skinoptics.html Google Scholar

23.

W. G. ZijlstraA. BuursmaO. 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

24.

I. FredrikssonM. LarssonT. Strömberg, “Accuracy of vessel diameter estimated from a vessel packaging compensation in diffuse reflectance spectroscopy,” Proc. SPIE, 8087 80871M (2011). http://dx.doi.org/10.1117/12.889684 PSISDG 0277-786X Google Scholar

25.

I. Fredriksson, “Monte Carlo simulations,” Quantitative Laser Doppler Flowmetry, Linköping University, (2009). Google Scholar

26.

J. A. NocedalS. J. Wright, Numerical Optimization., 2nd ed.Springer, New York, NY (2006). Google Scholar

27.

I. Fredrikssonet al., “Improved calibration procedure for laser Doppler perfusion monitors,” Proc. SPIE, 7906 790602 (2011). http://dx.doi.org/10.1117/12.871938 PSISDG 0277-786X Google Scholar

28.

P. R. Bargoet al., “In vivo determination of optical properties of normal and tumor tissue with white light reflectance and an empirical light transport model during endoscopy,” J. Biomed. Opt., 10 (3), 034018 (2005). http://dx.doi.org/10.1117/1.1921907 JBOPFO 1083-3668 Google Scholar

29.

I. FredrikssonM. LarssonT. Strömberg, “Measurement depth and volume in laser Doppler flowmetry,” Microvasc. Res., 78 (1), 4 –13 (2009). http://dx.doi.org/10.1016/j.mvr.2009.02.008 MIVRA6 0026-2862 Google Scholar

30.

S. L. Jacques, “Optical assessment of cutaneous blood volume depends on the vessel size distribution: a computer simulation study,” J. Biophotonics, 3 (1–2), 75 –81 (2010). http://dx.doi.org/10.1002/jbio.200900085 JBOIBX 1864-063X Google Scholar

31.

S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol., 58 (11), R37 –R61 (2013). http://dx.doi.org/10.1088/0031-9155/58/11/R37 PHMBA7 0031-9155 Google Scholar

32.

I. V. MeglinskyS. J. Matcher, “Modelling the sampling volume for skin blood oxygenation measurements,” Med. Biol. Eng. Comput., 39 (1), 44 –50 (2001). http://dx.doi.org/10.1007/BF02345265 MBECDY 0140-0118 Google Scholar

33.

J. W. Featheret al., “A portable scanning reflectance spectrophotometer using visible wavelengths for the rapid measurement of skin pigments,” Phys. Med. Biol., 34 (7), 807 –820 (1989). http://dx.doi.org/10.1088/0031-9155/34/7/002 PHMBA7 0031-9155 Google Scholar

34.

C. E. Thornet al., “Is mean blood saturation (SmbO2) a useful marker of tissue oxygenation?,” Am. J. Physiol. Heart Circ. Physiol., 296 (5), H1289 –H1295 (2009). http://dx.doi.org/10.1152/ajpheart.01192.2008 0363-6135 Google Scholar

35.

M. LarssonT. Strömberg, “Toward a velocity-resolved microvascular blood flow measure by decomposition of the laser Doppler spectrum,” J. Biomed. Opt., 11 (1), 014024 (2006). http://dx.doi.org/10.1117/1.2166378 JBIEDR 0141-5425 Google Scholar

36.

I. FredrikssonM. LarssonT. Strömberg, “Absolute flow velocity components in laser Doppler flowmetry,” Proc. SPIE, 6094 60940A (2006). http://dx.doi.org/10.1117/12.659206 PSISDG 0277-786X Google Scholar

37.

H. Karlssonet al., “Inverse Monte Carlo for estimation of scattering and absorption in liquid optical phantoms,” Opt. Express, 20 (11), 12233 –12246 (2012). http://dx.doi.org/10.1364/OE.20.012233 OPEXFF 1094-4087 Google Scholar

38.

R. B. Saageret al., “Method for depth-resolved quantitation of optical properties in layered media using spatially modulated quantitative spectroscopy,” J. Biomed. Opt., 16 (7), 077002 (2011). http://dx.doi.org/10.1117/1.3597621 Google Scholar

39.

E. JonathanJ. EnfieldM. J. Leahy, “Correlation mapping method for generating microcirculation morphology from optical coherence tomography (OCT) intensity images,” J. Biophotonics, 4 (9), 583 –587 (2011). http://dx.doi.org/10.1002/jbio.201000103 JBOIBX 1864-063X Google Scholar

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Ingemar Fredriksson, Oleg Burdakov, Marcus Larsson, and Tomas Strömberg "Inverse Monte Carlo in a multilayered tissue model: merging diffuse reflectance spectroscopy and laser Doppler flowmetry," Journal of Biomedical Optics 18(12), 127004 (18 December 2013). https://doi.org/10.1117/1.JBO.18.12.127004

Published: 18 December 2013

Access the abstract

JOURNAL ARTICLE
15 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY