## 1.

## Introduction

Diffuse reflectance spectroscopy (DRS) can be used for assessing tissue chromophores. Using white light in the visible to near infrared wavelength range and small source-detector separations will allow superficial sampling of tissue which is perfused mainly by the vessels in the microcirculation. Furthermore, in this wavelength range the absorption spectra of oxygenized and reduced hemoglobin show distinct features. This is advantageous for the determination of tissue oxygenation. With a model-based analysis of calibrated DRS, where knowledge on the absorption spectra of included tissue chromophores is assumed, quantitative measures of the chromophores can be attained.^{1}2.3.^{–}^{4} For accurate estimations, this kind of model must include all chromophores of importance in the tissue, or else included chromophores will compensate for omitted ones.^{1}^{,}^{5} Additionally, the fact that blood is confined to blood vessels and not homogeneously distributed must be taken into account.^{6}^{,}^{7}

For setups using small source-detector separations in the visible and near infrared wavelength range, where the absorption is on the same scale as the reduced scattering coefficient, diffusion theory fails to accurately describe light propagation. Instead, numerical simulations using Monte Carlo techniques provide a way to overcome these deficiencies. In addition, the Monte Carlo technique can be used to accurately model arbitrary complex structures.

It has previously been shown that a homogeneous model fails to describe the light propagation in skin at multiple source-detector separations in the visible wavelength range.^{5} This fact is a strong indication that a multilayered model should be used when analyzing DRS data from skin in order to accurately estimate the chromophore content. Both Wang et al.^{8} and Yudovsky et al.^{9} used Monte Carlo simulations of a two layered model in combination with a neural network analysis approach for estimating tissue optical properties (OP). Yudovsky et al. demonstrated that spatial frequency domain reflectance from a two-layered tissue model could be used for estimating dermis absorption and reduced scattering, while only epidermal optical thickness could be determined and this with limited accuracy. Wang et al. used spatially resolved diffuse reflectance and found that upper layer thickness improved upper layer OP estimation while lower layer OP accuracy deteriorated. These studies essentially determine OP using absolute calibrated reflectance at one wavelength but imply extension to spectroscopic use.

The approach proposed in this article differs in that it utilizes a subset of Monte Carlo simulated data for a limited number of tissue geometry and scattering parameters while adding the unique absorption characteristics of each chromophore directly in the inverse Monte Carlo algorithm. This is done by applying Beer-Lambert’s law on multiple path length-distributions in each layer in the forward calculation. Recordings were taken at two source-detector distances to make the method more sensitive to a layered structure of the tissue. The inverse algorithm is designed not to be sensitive to the average intensity detected at the two distances. This eliminates the need for an absolute calibration of the DRS system, which is difficult to perform accurately, and hence, eliminates the need for a well characterized calibration phantom that is stable over time.

The aim of this article is to present an inverse Monte Carlo method that uses a multilayered tissue model to estimate microcirculatory parameters in a diffuse reflectance spectroscopy setup. The accuracy of the method is evaluated using Monte Carlo simulations of tissue-like models containing discrete blood vessels. The performance of the method is compared to the performance of a similar method using a homogeneous model with absolute and relative intensity calibration. Examples of *in vivo* measurements are also given.

## 2.

## Methods, Models, and Measurements

## 2.1.

### Three Layered Bio-Optical Model

Photon propagation in tissue was modeled using Monte Carlo simulations of a three layered bio-optical model with nine free parameters. The design was based on previously presented models,^{10}^{,}^{11} and consisted of one epidermal layer with a variable thickness ${t}_{\mathrm{epi}}$ (parameter 1) and two dermal layers where the upper had a fix thickness of 0.5 mm and the lower an infinite thickness. The scattering coefficient [${\mathrm{mm}}^{-1}$] was equal for all layers and was modeled as

## (1)

$${\mu}_{\text{s}}=\alpha [(1-\gamma ){(\lambda /{\lambda}_{0})}^{-\beta}+\gamma {(\lambda /{\lambda}_{0})}^{-4}],$$^{12}as and the absorption coefficient [${\mathrm{mm}}^{-1}$] of the epidermis layer was modeled as

## (3)

$${\mu}_{\text{a},0}(\lambda )={f}_{\mathrm{mel}}\text{\hspace{0.17em}}{\mu}_{\text{a},\mathrm{mel}}(\lambda ).$$The dermis layers contained different tissue fractions of blood (${f}_{\text{blood},n}$ where $n$ is the layer number, parameters 6-7), where parameter 6 determines the average blood tissue fraction of the two layers and parameter 7 determines the relative ratio, ${f}_{\text{blood},\mathrm{rel}}$, between the layers. The same blood oxygen saturation was used in both layers ($s$ [-], parameter 8). The absorption spectra of fully oxygenated and deoxygenated blood^{13} (${\mu}_{\text{a},\mathrm{oxy}}$ and ${\mu}_{\text{a},\mathrm{deoxy}}$, respectively, hematocrit of 43%) are shown in Fig. 1 together with the absorption spectrum of melanin. The absorption coefficient of blood was given by

## (4)

$${\mu}_{\text{a},\text{blood}}(\lambda )=s{\mu}_{\text{a},\mathrm{oxy}}(\lambda )+(1-s){\mu}_{\text{a},\mathrm{deoxy}}(\lambda )$$## (5)

$${\mu}_{\text{a},n}(\lambda )={f}_{\text{blood},n}\text{\hspace{0.17em}}{c}_{\mathrm{vp}}(\lambda ){\mu}_{\text{a},\text{blood}}(\lambda ),$$## (6)

$${c}_{\mathrm{vp}}(\lambda )=\frac{1-\mathrm{exp}[-D{\mu}_{\text{a},\text{blood}}(\lambda )]}{{D}_{{\mu}_{\text{a},\text{blood}}}(\lambda )}.$$Mathematically, the nine parameters ${p}_{1}\dots {p}_{9}$ were used to calculate the skin model parameters as:

Here, ${t}_{\mathrm{epi}}$ is the epidermal thickness in mm, ${f}_{\mathrm{mel}}$ is the volume fraction of melanin, ${f}_{\text{blood},\mathrm{avg}}$ is the average tissue volume fraction of blood, ${f}_{\text{blood},\mathrm{rel}}$ is the relative difference in blood tissue fraction from the average in the two dermal layers (ranging between $-1$ and 1), ${f}_{\text{blood},1}$ and ${f}_{\text{blood},2}$ are the volume tissue fraction of blood in the two dermis layers, respectively, $s$ is the fraction of oxygenated blood, and $D$ is the average vessel diameter in mm.

## 2.1.1.

#### Base simulations

Monte Carlo simulations were run for all combinations of six levels of the epidermis thickness ([mm])

and 12 levels of the scattering coefficients ([${\mathrm{mm}}^{-1}$])In total 72 simulations were thus performed.

A Henyey-Greenstein phase function with $g=0.8$ was used in all layers of the model. All layers had a refractive index of 1.4, whereas the refractive index of the source, detector and surrounding media was 1.5. The numerical aperture of the source and detector fibers was set to 0.37. All photons propagating outside a 5 mm radius of the source were terminated to increase simulation speed. In each simulation five million photons were detected between 0.01 and 1.8 mm from the pencil beam light source. The position of detection, photon weight and the path length in each layer was stored in the simulations.

## 2.1.2.

#### Path length distributions

When analyzing the base simulations, the pencil beam source was smeared out to mimic the size of the real fiber core by convoluting the detection positions with a circular area of 200 *μ*m diameter.^{15} Two ring detectors were used, one between 0.3 and 0.5 mm (from the source), and one between 1.1 and 1.3 mm. The photons were re-weighted to reflect the detected intensity in a fiber rather than the ring detectors used in the simulations.^{16} Separate path length distributions were created for photons detected by these two detectors. For each detector, photons were grouped according to the deepest layer that they had reached. For each group, distributions were created over the path lengths that these photons had been propagated in each of the three layers. This resulted in a set of six distributions per simulation and detector: the distribution of path lengths in layer one for the photons that had only been propagating in layer one; the distribution in layer one for photons that had only been propagating in layer one and two; the distribution in layer two for the same photons; and so on. The resulting distributions (for one source detector separation) is denoted ${p}_{\mathrm{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, $m\le n$. An example of the six path length distributions for the short source-detector separation (base simulation with ${t}_{\mathrm{epi}}=0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ and ${\mu}_{\text{s}}=8\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$) is shown in Fig. 2.

## 2.1.3.

#### Forward problem

The forward problem of calculating a model spectrum from the three layered model was based on interpolation of the path length distributions from the base simulations using Beer-Lambert’s law to add the absorption for each path length. We have previously described the important steps of this multi-path-length multilayered Beer-Lambert algorithm.^{11}

Initially, a set of path length distributions ${p}_{\mathrm{pl},m,n}$($d$) were interpolated linearly in two dimensions from the 72 base simulations, depending on the epidermis thickness of the model and the scattering coefficient of the wavelength of interest. From these distributions, the total intensity ${I}_{0,n}$ for photons that had reached layer $n$, was calculated as

where ${I}_{0,n}$ is independent of $m\le n$. The path length distributions for varying $n$ were then normalized to unityThe path length distributions were modified for all path lengths $d$ by adding the absorption using Beer-Lamberts law

## (22)

$${p}_{\mathrm{pl},m,n}^{\prime \prime}(d)={p}_{\mathrm{pl},m,n}^{\prime}(d)\mathrm{exp}(-d{\mu}_{\mathrm{a},m}),$$## (23)

$${I}_{n}={I}_{0,n}\underset{m=0}{\overset{n}{\mathrm{\Pi}}}\sum _{j}{p}_{\mathrm{pl},m,n}^{\prime \prime}({d}_{j}),$$^{11}With this assumption, a Monte Carlo-based technique that is sufficiently fast to be used in a real time inverse engineering approach with a multilayered model is introduced ($<1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ms}$ for calculating a single 32-point spectra at the two source-detector separations with our implementation using an ordinary workstation).

The flow chart in Fig. 3 shows an overview of the forward problem.

## 2.2.

### Inverse Problem

We utilized the Levenberg-Marquardt method^{17} to find an optimal fit between measured (or simulated) spectra and the modeled spectra resulting from the forward calculation of the three layered analysis model. The nine model parameters were optimized for a least square fit at the two source detector separations 0.4 and 1.2 mm simultaneously. Before comparing the spectra, they were normalized by the average intensity of the short and long source-detector separations, i.e.,

## (24)

$${I}_{s}^{\prime}(\lambda )=\frac{{I}_{s}(\lambda )}{\left(\langle {I}_{0.4}{\rangle}_{\lambda}+\langle {I}_{1.2}{\rangle}_{\lambda}\right)/2},$$The choice of merit function, as described in Sec. 2.2.1, is crucial for model convergence. To further improve the convergence properties of the inverse problem a canonical correlation analysis was performed to find linear combinations of the nine model parameters with maximal correlation to linear combinations of the $32\times 2$ spectral points. The Levenberg-Marquardt method was then applied on the resulting linear combinations of the model parameters instead of the original parameters. This is further described in Sec. 2.2.2. Furthermore, a multiple starting point strategy was utilized where it was assumed that the global minimum was found when the same minimum point was found from five different start positions.^{11}

## 2.2.1.

#### Merit function

The merit function to minimize is given by

## (25)

$${\chi}^{2}=\sum _{s}\sum _{\lambda}{\left[\frac{{I}_{s,\mathrm{meas}}^{\prime}(\lambda )-{I}_{s,\text{model}}^{\prime}(\lambda )}{{\sigma}_{s}(\lambda )}\right]}^{2}+{\left[\frac{\phi ({I}_{0.4,\mathrm{meas}}^{\prime})-\phi ({I}_{0.4,\text{model}}^{\prime})}{2}\right]}^{2}+{\left[\frac{\phi ({I}_{1.2,\mathrm{meas}}^{\prime})-\phi ({I}_{1.2,\text{model}}^{\prime})}{0.5}\right]}^{2},$$## (26)

$$\phi (I)=\frac{\langle I{\rangle}_{\lambda \in [535,540],[575,580]}}{\langle I{\rangle}_{\lambda \in [550,570]}}.$$Adding the function $\phi (I)$ to the ${\chi}^{2}$-error adds importance to the shape of the spectrum at the double-peak interval of the oxygenated hemoglobin spectrum, which improves the oxygen saturation estimation of the method. To increase the convergence properties and to emphasize important aspects of the spectra, the weighting factor ${\sigma}_{s}$ vector should be chosen carefully. For example, by setting ${\sigma}_{s}$ to a constant, the optimization is based on the *absolute* difference between the measurement and the model, whereas the optimization is based on the *relative* difference if ${\sigma}_{s}={I}_{s,\mathrm{meas}}$. In the current study, the weighting factor was calculated as

## (27)

$${\sigma}_{s}(\lambda )=\frac{{\zeta}_{s}(\lambda )}{({\langle {\zeta}_{0.4}(\lambda )\rangle}_{\lambda}+\langle {\zeta}_{1.2}(\lambda ){\rangle}_{\lambda})/2},$$## (28)

$${\zeta}_{s}(\lambda )=\frac{w(\lambda )}{v(S)}\left[\langle {I}_{s,\mathrm{meas}}{\rangle}_{\lambda \in ]500,600[}+{I}_{s,\mathrm{meas}}(\lambda )\right],$$## (29)

$$w(\lambda )=\{\begin{array}{cc}1& \text{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\lambda \notin ]500,600[\\ \frac{((\lambda -550)/50{)}^{4}+1}{2}& \text{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\lambda \in [500,600]\end{array}.$$By calculating $\sigma $ in this way the ${\chi}^{2}$-error will be influenced both by the absolute and the relative difference where the importance of the absorption band of hemoglobin between 500 mm and 600 nm is amplified [$w(\lambda )$], and the contribution to the ${\chi}^{2}$-error is higher for the longer source-detector separation [$v(S)$].

## 2.2.2.

#### Parameter de-correlation

Although the nine model parameters have been chosen in a way that they do not apparently change the spectrum in a similar way, there is still a high degree of correlation between the parameters in the way they affect the spectra at the two source-detector separations. That is not ideal for the convergence properties of the nonlinear optimization problem. To minimize this correlation, a canonical correlation analysis (CCA) was performed. 5000 sets of the model parameters ($\mathbf{X}$) were randomly chosen between some relevant limits, and the spectra at the two source-detector separations ($\mathbf{Y}$) were then calculated for each set of the parameters. The CCA was then performed for $\mathbf{X}$ and $\mathbf{Y}$, which resulted in linear combinations ${\mathbf{a}}_{k}$ of parameters in $\mathbf{X}$ which were maximally correlated to linear combinations ${\mathbf{b}}_{k}$ of the intensities in $\mathbf{Y}$, at the same time as ${\mathbf{a}}_{k}\mathbf{X}$ was uncorrelated with ${\mathbf{a}}_{l}\mathbf{X}$ for all $k\ne l$, and ${\mathbf{b}}_{k}\mathbf{Y}$ was uncorrelated with ${\mathbf{b}}_{l}\mathbf{Y}$ for all $k\ne l$. In other words, nine new parameters were created, which were linear combinations of the original nine parameters. These new parameters affected the spectra in a maximally uncorrelated manner.

## 2.2.3.

#### Output

We chose to present the average RBC tissue fraction within a 1.13 mm radius ($3.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ half sphere) from the center of the light emitting fiber, as the RBC tissue fraction was not homogeneously distributed in this analysis model. This volume reflects the main sampling volumes of the two source-detector separations rather well.^{11}^{,}^{18} In addition, the estimated oxygen saturation is presented and evaluated.

## 2.3.

### Homogeneous Analysis Model

For comparison, a homogeneous semi-infinite model containing a turbid medium that scatters light according to a Henyey Greenstein phase function with an anisotropy $g=0.8$, was evaluated. The geometry and the numerical aperture of the illuminating and detecting fibers, and the refractive indexes of the model, were all identical to those in the multilayer model. By using a white Monte Carlo approach^{19} the effect of changing the scattering and absorption coefficients was evaluated by post-processing of a single simulation consisting of half a million backscattered photons. In this way a dataset that describes how the absorption and scattering affect the detected intensity was constructed.

The forward problem of calculating a simulated spectrum was solved using linear interpolation on this dataset, assuming a scattering and absorption model equal to that in the multilayer model [i.e., Eqs. (1) through (17)]. To solve the inverse problem, where model spectra are fitted to measured spectra (i.e., validation simulation spectra), a trust-region-reflective optimization algorithm was employed. Optimization using this algorithm was done in either one (relative calibration) or two steps (absolute calibration). First, a shape-sensitive merit function that discarded any overall magnitude difference was used. This function was defined as

## (30)

$${\chi}_{s}^{2}=\sum _{\lambda}{(\frac{{I}_{s,\text{model}(\lambda )}}{{I}_{s,\mathrm{meas}(\lambda )}}/\langle \frac{{I}_{s,\text{model}}(\lambda )}{{I}_{s,\mathrm{meas}(\lambda )}}{\rangle}_{\lambda}-1)}^{2}.$$In the second step the optimization parameters ($\alpha $, $\beta $, $\gamma $, ${f}_{\mathrm{mel}}$, ${f}_{\text{blood}}$, $s$ and $D$) where fine-tuned by including also the overall magnitude difference. This merit function was defined as

## (31)

$${\chi}_{s}^{2}=\sum _{\lambda}{[\frac{{I}_{s,\text{model}}(\lambda )}{{I}_{s,\mathrm{meas}}(\lambda )}-1]}^{2}.$$The last step demands that measured and modeled spectra are absolute calibrated in magnitude.

## 2.4.

### Validation Simulations

In order to test if a three layered model can be used to describe the light propagation in geometrically much more complex tissue, a set of 50 random validation models, all containing more than 100 individual blood vessels, were Monte Carlo simulated. By fitting spectra from the three-layered model to the spectra resulting from these simulations, the accuracy of the method could be evaluated. The details on how these models were randomly created can be found in the Appendix. These models consisted of an epidermis layer with a thickness varying between 0.045 and 0.47 mm with a melanin tissue fraction varying between 0.28% and 9.2%.

A semi-infinite nonabsorbing dermis layer was located below the epidermis layer. The dermis layer contained between 100 and 1750 infinitely long discrete blood vessels located with a random orientation parallel to the surface layer. The center positions of these blood vessels where randomly placed within a 3 mm radius of the light source. The blood vessels contained blood with a wavelength dependent absorption coefficient as given by Fig. 1, and the vessel diameters ranged between 6 and 160 *μ*m. The blood vessel diameter generally increased with depth in the model. The oxygen saturation level of the blood in the vessels ranged between 0% and 100%.

The scattering coefficient was equal for both layers and modeled according to Eq. (1), with resulting values ranging between 1.1 and $17\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. A Henyey-Greenstein phase function with $g=0.85$ in the epidermis layer and 0.7 in the dermal layer was used. The reduced scattering coefficient of the blood was set to $2.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for all wavelengths, and a Gegenbauer kernel phase function with ${g}_{\mathrm{Gk}}=0.948$ and ${\alpha}_{\mathrm{Gk}}=1.0$ (resulting in $g=0.991$) was used.^{10}

The refractive index of both layers and the blood vessels was set to 1.4, whereas the refractive index of the source, detectors and surrounding media was set to 1.5. The numerical aperture of both source and detector was 0.37, the source had a diameter of $200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$ and the detectors were set to ring detectors 0.3 to 0.5 and 1.1 to 1.3 mm from the center of the source, respectively. For each of the 32 wavelengths, 150,000 photons were detected in total for the two detectors. A cross sectional view of one of the models is shown in Fig. 4.

In order to calculate the accuracy of the estimated RBC tissue fraction and oxygen saturation of the fitted three layered analysis model, we calculated the RBC tissue fraction and oxygen saturation within a $3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ half sphere in the 50 simulated models. This calculation was done by randomly probing 15 million homogeneously distributed points within the $3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ volume, and calculating the average RBC tissue fraction and oxygen saturation of these points. The $3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ half sphere is marked in Fig. 4.

## 2.5.

### Testing Parameter Importance

The spectral fitting of various analysis models including different sets of the parameters ${p}_{1}-{p}_{9}$, were statistically compared by evaluating the ${\chi}^{2}$-statistics [Eq. (25)]. The full analysis model (${\chi}_{\text{all}}^{2}$ including parameters ${p}_{1}-{p}_{9}$), was compared to the analysis model with a reduced number of parameters: indefinitely thin epidermis thickness (${\chi}_{{t}_{\mathrm{epi}=0}}^{2}$), constant scattering parameter $\alpha =2$ (${\chi}_{\alpha =2}^{2}$), removing $\beta $ by setting scattering parameter $\gamma =1$ (${\chi}_{\gamma =1}^{2}$), removing scattering parameter $\gamma $ (${\chi}_{\gamma =0}^{2}$), no melanin in the epidermis (${\chi}_{{f}_{\mathrm{mel}}=0}^{2}$), constant average RBC tissue fraction (${\chi}_{{f}_{\text{blood},\mathrm{avg}}=0.4}^{2}$), same RBC tissue fraction in both dermis layers (${\chi}_{{f}_{\text{blood},\mathrm{rel}}=0}^{2}$), constant blood oxygen saturation (${\chi}_{s=0.5}^{2}$), and no compensation for the vessel packaging effect (${\chi}_{D=0}^{2}$). Note that the model with an indefinitely thin epidermis thickness still contains melanin [see Eq. (11)]. The ${\chi}^{2}$-error for the full analysis model versus any of the reduced analysis models was tested using the test variable

## (32)

$$v=\frac{{\chi}_{\mathrm{red}}^{2}-{\chi}_{\text{all}}^{2}}{{n}_{{\chi}_{\text{all}}^{2}}-{n}_{{\chi}_{\mathrm{red}}^{2}}}/\frac{{\chi}_{\text{all}}^{2}}{N-{n}_{{x}_{\text{all}}^{2}}-1},$$^{20}

## 2.6.

### In Vivo Measurement

One measurement was performed on the volar side of the forearm of a healthy volunteer (male, 31 years old, fair skin). A blood pressure cuff was placed on the upper arm and inflated to a pressure of 170 mm Hg after 1 min of the measurement, completely occluding the blood flow. The pressure was held for 3 min and the measurement lasted for another 3 min.

A custom made optical fiber probe was used to deliver light to the tissue and to transmit the backscattered light to the detectors. The illuminating fiber was placed in the middle of the probe tip, and two light collecting fibers were placed at 0.4 and 1.2 mm, respectively, away from the center of the illuminating fiber. The fibers were made of fused silica, had a diameter of $200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$ and a numerical aperture of 0.37. The illuminating fiber was connected to a broadband white light source (Avalight-HAL-S, Avantes BV, The Netherlands), and the light collecting fibers were connected to two spectroscopes (AvaSpec-ULS2048L, Avantes BV, The Netherlands). The measurement probe was placed in a plastic probe holder (PH 08; Perimed AB, Järfälla, Sweden) which was fixated to the skin using double-sided adhesive rings (PF 105-1; Perimed AB). This assured a minimum contact pressure between the probe tip and the examined skin.

Before analyzing the recorded spectra, the spectra corresponding to each source-detector separation, were processed in three calibration steps:

1. A dark reference spectrum was subtracted.

2. The spectrum was normalized with a white reference spectrum.

3. The spectrum was normalized by the mean intensity at 612 to 617 nm from a measurement (calibrated according to steps 1 and 2 above) where the face of the probe was evenly illuminated.

Steps 2 and 3 are both relative calibration steps that enable inter-wavelength and inter-channel intensity comparisons. No additional calibration between modeled and measured intensity is necessary as this difference is removed when solving the inverse problem by normalizing with their average intensity [Eq. (24)].

## 3.

## Results

## 3.1.

### Forward Calculation

We simulated a three layered model with RBC tissue fraction of 0.32% and blood oxygen saturation of 80%. The detected intensities at each wavelength were compared with the forward calculated spectra of the three layered analysis model with the exact same properties. The comparison is shown in Fig. 5. In this particular example, the relative RMS deviation between the simulated and forward calculated spectra was 1.0%.

## 3.2.

### Necessary Parameters

The comparison of analysis models [Eq. (32)] was performed on the 50 simulated models. The results indicated that the full analysis model had a significantly smaller ${\chi}^{2}$-error than the reduced analysis models in ${n}_{\text{significant}}$ of 50 models according to Table 1. The effect of the reduced analysis models on the absolute RMS deviation between the true and estimated RBC tissue fraction and blood oxygen saturation was also evaluated on the 50 simulated models. The RMS deviations for all models are also found in Table 1.

## Table 1

Performance of various analysis models.

χall2 | χtepi=02 | χα=22 | χγ=12 | χγ=02 | χfmel=02 | χfblood,avg=0.42 | χfblood,rel=02 | χs=0.52 | χD=02 | |
---|---|---|---|---|---|---|---|---|---|---|

nsignificant | — | 37 | 50 | 50 | 49 | 47 | 45 | 28 | 41 | 39 |

RMS RBCa | 0.087 | 1.1 | 3.8 | 8.1 | 0.12 | 4.7 | 0.22 | 0.10 | 0.10 | 0.14 |

RMS oxyb | 4.2 | 9.9 | 19 | 40 | 13 | 38 | 14 | 4.0 | 11 | 3.5 |

## a

Number of simulated models that had a significantly lower χ2-error with the full analysis model than the reduced model.

## b

Absolute RMS deviation between estimated and true RBC tissue fraction or blood oxygen saturation, in percentage units.

The spectral fit was significantly better for at least 28 of the 50 simulated models for the full analysis model compared to any of the reduced analysis models. In addition, the RMS deviation between the true and estimated RBC tissue fraction and blood oxygen saturation was generally much lower for the full than for the reduced analysis models, and only marginally higher in two cases (blood oxygen saturation for ${\chi}_{{f}_{\text{blood},\mathrm{rel}}=0}^{2}$ and ${\chi}_{D=0}^{2}$). Hence, the full analysis model was used in the further analysis.

## 3.3.

### Accuracy of RBC Tissue Fraction and Oxygen Saturation

## 3.3.1.

#### Three layered analysis model

The plots in Fig. 6 show the estimated RBC tissue fraction and blood oxygen saturation, respectively, compared to the values calculated from the $3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ half sphere in the simulated models. The absolute RMS deviation between the true and estimated RBC tissue fraction was 0.087 percentage units and the absolute RMS deviation of the oxygen saturation was 4.2 percentage units (see also Table 1).

Half of the simulated models had an RBC tissue fraction $<0.3\%$. For these models, the absolute RMS deviation was 0.055 and 5.3 percentage units for the RBC tissue fraction and blood oxygen saturation, respectively. The *relative* RMS deviation for the RBC tissue fraction was 24% for these models. For the models with an RBC tissue fraction $>0.3\%$, the absolute RMS deviation was 0.11 and 2.9 percentage units for the RBC tissue fraction and blood oxygen saturation, respectively, and the relative RMS deviation for the RBC tissue fraction was 19%.

An example of the model fit is given in Fig. 7. The chosen spectra had the median ${\chi}^{2}$-deviation of all 50 models.

## 3.3.2.

#### Homogeneous analysis model

The 50 simulated models were also analyzed using the homogeneous analysis model described in Sec. 2.3. The data was analyzed in three different ways: 1) a single detector using absolute calibration; 2) a single detector using relative calibration; and 3) two detectors using relative calibration. When the absolute intensity of a single spectrum was considered [1, merit function Eqs. (30) and (31)], the absolute RMS deviation of the RBC tissue fraction was 0.20 and 0.15 percentage units, respectively, analyzed for the source-detector separations of 0.4 and 1.2 mm. When the single spectra were relative calibrated [2, by using Eq. (30) in the merit function], corresponding numbers were 0.20 and 0.14 percentage units. These results should be compared to an RMS deviation of 0.087 percentage units for the three layered analysis model. Corresponding numbers for the RMS deviation of the blood oxygen saturation were 11 and 5.8 percentage units for absolute calibrated spectra and 12 and 5.3 percentage units for relative calibrated spectra, compared to 4.2 percentage units for the three-layer model. Paired t-tests indicated that the average deviation between estimated and true RBC tissue fraction and blood oxygen saturation was significantly ($p<0.01$) lower for the three layered analysis model than for the homogeneous analysis model at both source-detector separations. The estimated RBC tissue fraction and blood oxygen saturation from the homogeneous analysis model with absolute calibration are shown in Fig. 8.

When using the homogeneous model and considering relative calibrated spectra at both separations (3), the spectral fit was poor (average spectral RMS deviation of 9.7%, compared to 1.3% for the three layered model), with obvious incorrect estimations of especially the blood oxygen saturation (RMS deviation of 30 percentage units).

## 3.4.

### In Vivo Measurement

The RBC tissue fraction and blood oxygen saturation during the occlusion measurement are shown in Fig. 9. A slight increase in RBC tissue fraction can be noticed during the first minute of the occlusion, likely due to redistribution of blood to the venous side. As the pressure is released, a short peak is seen in the RBC tissue fraction, before it approaches the baseline. The blood oxygen saturation is slowly decreased during the occlusion, and a massive and quick increase is seen when the pressure is released, before it slowly approaches the baseline. Two examples of the spectral fit, at 0.5 and 3 min, are shown in Fig. 10.

## 4.

## Discussion

There are two key findings in this study. First, a three layered multi parameter tissue model is needed to explain spectra recorded at two separate source-detector distances in realistic simulated tissue models. The light propagation cannot be explained using a homogeneous tissue model. Secondly, when detecting the light at two distances from the source, a relative calibration routine, which only determines the relative intensity gain between the detectors, is sufficient.

The proposed method has an important advantage in the fact that it is not dependent on an absolute intensity calibration. An accurate absolute calibration is complicated to perform as it is sensitive to changes in the calibration standards, and thus error prone. Furthermore, when using an absolute calibration the inverse algorithm may become sensitive to measurement imperfections that are not included in the model but affect the two distances similarly, such as an air gap between the probe and the tissue and to intensity variations in the light source.^{3} The need for use of an index matching liquid is thus reduced, which facilitates the measurement procedure.

Simulated tissue-like models were used to determine the accuracy of the method, since controlled optical phantoms cannot be constructed with enough geometrically complexity. We have previously proven that it is possible to simulate homogeneous optical phantoms using a very similar measurement setup with high accuracy.^{21} Wang et al. managed to construct a two-layered semi-infinite tissue phantom^{8} with controlled OP, but still their tissue model was too simplified to mimic skin. The aim of the validation in this study was to test if a three layered model can describe the light propagation in a geometrically more complex, tissue-like, structure, and therefore optical phantoms were ruled out. *In vivo* measurements are improper for initial method validation since the exact expected values of the RBC tissue fraction and oxygen saturation cannot be determined, and the accuracy can thus not be evaluated. However, *in vivo* studies where the method is used are important to show the applicability of the method, and such studies are ongoing.

The three layered analysis model consisted of nine variable parameters. It was concluded that all nine parameters contributed to a significantly improved model fit in at least 28 of the 50 simulated tissue-like models. Based on the data presented in Table 1, it may be questioned if the parameters controlling the relative amount of blood in the two dermal layers and the average vessel diameter are necessary (significantly better fit in 28 and 39 of 50 spectra, but comparable RMS deviation in RBC tissue fraction and oxygen saturation). It should be noted, however, that in most of the 50 simulated models the blood tissue fraction does generally not differ much between superficial and deeper parts of the dermis, and the average vessel diameters are generally small compared to what can be expected in real tissue. In models where this was not the case, both the spectral fit and the estimation of the RBC tissue fraction were much worse when any of these parameters were excluded. The data in Table 1 also indicate the importance of not excluding any significant parameter; as doing so will cause the other parameters to compensate for the missing one. By omitting significant parameters the estimations of the RBC tissue fraction and oxygen saturation will be poor (see for example the column ${\chi}_{\gamma =0}^{2}$ in Table 1). The importance of not excluding any significant parameter has also been demonstrated in other studies.^{1}^{,}^{5}

In measurements on real tissue, the set of model parameters may need to be revised to better reflect the actual content and structure of the tissue under investigation and to improve the spectral fit. In epidermis and subcutaneous fat, beta-carotene has a nonnegligible absorption up to 500 nm.^{22} In blood, small amounts of bilirubin and met-hemoglobin are present and these concentrations may be high enough to have a significant influence on the diffuse reflected spectra in the visible wavelength range $<500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ for bilirubin and between about 600 and 650 nm for met-hemoglobin) under some physiological conditions.^{13}^{,}^{22}^{,}^{23} Water and fat are present in high amounts in tissue, but their absorption is negligible in the wavelength range used in this study (450 to 850 nm).^{24}^{,}^{25} In muscle tissue, oxygenized and reduced myoglobin and cytochrome (c and aa3) both have significant effect on the diffusely reflected spectra. Myoglobin has absorption spectra rather similar to hemoglobin and it may be difficult to distinguish them. The cytochromes, which have characteristic absorption bands in the wavelength interval 500 to 650 nm, take a central part of the oxygen transport in the mitochondria in the cell, and monitoring their oxidation status is physiologically very interesting.^{1}^{,}^{26} The same strategy as adopted in this study to determine model parameters that has a significant effect on the spectral fit can be adopted in clinical studies to determine which additional chromophores should be included when analyzing various types of tissue, and if any of the nine parameters included here should be removed. Studying Fig. 10, it may be seen that the spectral fit is relatively poor below 500 nm for the long fiber separation, a fact indicating that the model might benefit from the inclusion of bilirubin.

An alternative to the multi-path-length multilayered Beer-Lambert algorithm presented in this paper is to use a five dimensional look-up-table (LUT), where the five dimensions are constituted by epidermis thickness, scattering coefficient and absorption coefficient in the three layers, respectively. The LUT could be constructed from the same base simulations as used here, adding the effect of the absorption directly on each photon in a post processing step performed once (doing that for each forward calculation would be too slow). Using a LUT could even be faster than the proposed algorithm, but the absorption range in all three layers has to be restricted and the resolution in the absorption dimensions has to be fine enough over the whole range for accurate results.

One main focus during the development of the proposed multilayered inverse Monte Carlo technique has been to optimize the convergence properties. Equation (1) as well as Eqs. (7) through (17) are carefully chosen to ensure that the parameters affect the spectral shape in a somewhat linear manner over a wide range of parameter values, for example with the exponential in Eqs. (8) and (13), and the squares in Eqs. (11) and (12). Even if the parameters should affect the spectra in a perfectly linear manner, large correlations between how different parameters affect the spectra may exist. Such correlations have a negative influence on the convergence properties and should be avoided—such as by performing the canonical correlation analysis described in Sec. 2.2.2. Note that this type of de-correlation is not very efficient if the parameters affect the spectra in a highly nonlinear fashion.

When dealing with algorithms based on model fitting, one faces a possibility of rejecting models with a poor spectral fit to the measured spectra. One can also evaluate how sensitive the modeled spectra are to variations in the model parameters. If a parameter has a low impact on the model spectra, it is likely to be less accurately estimated. In both cases, one cannot assume that the fitted model reflects the actual quantities accurately. This is an advantage compared to methods that only calculate some type of index in the measured spectrum.

Other model parameters than the RBC tissue fraction and blood oxygen saturation may be clinically interesting in various clinical situations. Both the average vessel diameter and scattering parameters, obtained with related methods, has been used in studies involving various types of cancer.^{27}28.29.^{–}^{30} The epidermal thickness, volume fraction of melanin and depth distribution of blood (${f}_{\text{blood},\mathrm{rel}}$) may also be interesting, but before using any of these parameters as a clinical parameter, their accuracy should be evaluated. We have chosen not to evaluate the accuracy of these model parameters due to the fact that the method has been designed to estimate only the RBC tissue fraction and oxygen saturation accurately, especially when choosing the ${\chi}^{2}$-merit function (Sec. 2.2.1).

In conclusion, the three-layered analysis model is able to describe diffusely reflected spectra at two close but well separated source-detector separations simultaneously and is able to estimate the RBC tissue fraction and blood oxygen saturation more accurately than a homogeneous analysis model. An RMS deviation of 0.087 and 4.2 percentage units in the estimated RBC tissue fraction and blood oxygen saturation, respectively, was achieved on 50 simulated tissue-like models containing discrete blood vessels. No absolute calibration routine is needed for the method, and it is fast enough for real time analysis of recorded data.

## Appendices

## Appendix:

### Creation of Validation Simulations

This appendix gives the details about how the validation models were randomly created.

The vessel diameter was given from a gamma distribution with parameters $\alpha $ and $\beta $. The probability density function of a gamma distribution $\mathrm{\Gamma}(\alpha ,\beta )$) is given by

## (33)

$$f(x;\alpha ,\beta )=\frac{{\beta}^{\alpha}}{\mathrm{\Gamma}(\alpha )}{x}^{\alpha -1}\text{\hspace{0.17em}}\mathrm{exp}(-\beta x),$$The approximate RBC tissue fraction $C$ was given from $\xi $ where $\xi $ was given from the gamma distribution $\mathrm{\Gamma}(1.5,3)$. The number $N$ of vessels were calculated from

where $r=3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, ${f}_{\mathrm{HCT}}$ is the hematocrit (set to 0.43) and is the expectation value of the vessel cross sectional area at $z=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$.A lower limit ${s}_{0}$ of the blood oxygen saturation was randomly chosen between 0 and 50% (rectangular distribution) for each model. For each vessel, the oxygen saturation level was calculated from ${\xi}^{2}(1-{s}_{0})+{s}_{0}$, where $\xi $ was $R({s}_{0},100)$.

The epidermis thickness was set to ${\xi}^{2}$ where $\xi $ was normally distributed with expectation value of 0.4 and standard deviation of 0.1 ($N(0.4,0.1)$). The melanin tissue fraction in % was given by $\xi $ where $\xi $ was $\mathrm{\Gamma}(1.75,0.5)$. The parameters $\alpha $, $\beta $, $\gamma $ giving the scattering coefficient were all given by ${\xi}^{2}$ where $\xi $ was $N(0.9,0.15)$, $N(0.8,0.2)$ and $N(0.25,0.05)$, respectively.

## Acknowledgments

The 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), through the Research and grow programme (VINNOVA D. no. 2011-03074), and also by NovaMedTech, supported by the European Union Regional Development Fund.

## References

*μ*m wavelength region,” Appl. Opt. 12(3), 555–563 (1973).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.12.000555 Google Scholar