## Abstract

Hyperspectral imaging combining with skin optical clearing technique provides a possible way to non-invasively monitor hemodynamics of cutaneous microvessels. In order to estimate microvascular blood oxygen saturation, in this work, a lookup-table-based inverse model was developed to extract the microvascular optical and physiological properties using hyperspectral analysis. This approach showed a higher fitting degree than currently existing hyperspectral analysis methods (*i.e.* multiple linear regression and non-negative least square fit) in estimating blood oxygen saturation. Hypoxic stimulation experiment showed that calculated results were in accordance with physiological changes, and the relative changes of estimated oxygen saturation indicated this method appeared to be more sensitive to blood oxygen fluctuation. And a simulated blood model was used for verification here, indicating this method also showed a good accuracy in determining oxygen saturation from the simulated spectra.

© 2017 Optical Society of America

## 1. Introduction

The microvascular oxygen saturation (*SO _{2}*) can indicate the rate of oxygen delivery and consumption by skin and underlying tissues. It is very important to monitor the

*SO*in clinical diagnosis of microcirculation and metabolism, or basic research of chronic diseases and tumor development [1–4]. Hyperspectral imaging (HSI) technique has shown a powerful tool for noninvasive disease diagnosis [5–7], including eye dark circles evaluation [8], diabetic foot assessment [9], and skin wounds characterization [10, 11]. But both the image depth and contrast suffer from the high scattering characteristic of turbid skin. Recently, the developed

_{2}*in vivo*skin optical clearing techniques evidently improved the imaging performance of various optical image techniques including optical coherence tomography [12–14], photoacoustic microscopy [15,16], laser speckle contrast imaging [17–21]. It will show a great potential for enhancing the hyperspectral imaging quality for cutaneous microvessels.

However, the hyperspectral analysis technique is not less than satisfactory. To extract the physiological properties from the hyperspectral data sets, some analysis methods have been developed. For instance, the multiple linear regression (MLR) [21–24] and non-negative least square method (NNLS) [25–27], were proposed to solve least-square problem. The former assumed that the variables are linearly independent [28], which implies that there is little correlation among of variables. Besides, the MLR may induce negative estimators. NNLS could solve a constrained version of the least squares problem. The algorithm implementation of NNLS is similar to MLR, but the estimators only are not allowed to be negative [25–27]. However, the both could not consider the multilayer structure and optical properties of tissue. The lookup-table-based (LUT) inverse model was proposed to realize reflectance spectroscopic analysis for solving inverse problems [29–32]. The multilayer structure, the absorption and scattering properties of tissue can be considered in LUT-based inverse model, thus, it may be more accurate for spectral analysis [29]. Therefore, we can establish different LUT according to different optical properties and structure of tissues. A similar inverse method based on a semi-empirical model also has been proposed for hyperspectral analysis [33, 34]. It is well known that Monte Carlo (MC) simulations, as a golden standard of calculation modeling, has been used to investigate those complex systems and processes, and validate new theoretical models or evaluate the newly-developed techniques [35]. Thus, the established LUT based on MC simulation would have better performance on extracting the physiological properties. Zhong *et. al.* established a LUT based on MC simulation to extract the optical properties and physiological properties of skin. Even though the multilayer structure of skin was considered in their work, the spatial distribution could not be obtained based on fiber optical spectrometer [29]. In contrast to the fiber probe-based spectroscopy technique, HSI cannot only access the physiological parameters of blood, *e.g.* oxygen saturation, and hemoglobin, but also their spatial distribution. Thus HSI permit us to more visually and easily obtain the characterization of skin [6].

In this study, we used *in vivo* skin optical clearing technique to improve HSI imaging quality. And we applied the LUT-based inverse model based on MC simulation for hyperspectral analysis to estimate distribution of blood oxygen saturation. Further, we compared differences between LUT-based inverse method and common hyperspectral analysis methods (*i.e.* MLR and NNLS) for extracting *SO _{2}* from hyperspectral data sets. And the simulated blood model was used to verify the accuracy and correlation of these methods.

## 2. Theory and methods

#### 2.1 Hyperspectral imaging

The HSI system mainly consists of a charge coupled device (CCD) camera (Pixelfly USB, PCO Company, Germany), a liquid crystal tunable filter (LCTF, CRi Varispec VIS, Perkin Elmer, USA), a stereo microscope (SZ61TR, Olympus, Japan) and a ring of LED light with polarizer as shown in Fig. 1. The bandwidth of the LCTF was 7 nm. A certified reflectance standard (SRS-99-020, Labsphere, USA) was used to acquire the illuminated calibration hyperspectral images. All the hyperspectral images were acquired from 500 nm to 620 nm with a wavelength-interval of 10 nm. The dark noise was measured by keeping the camera shutter closed, and there were no motion artefacts in the hypercube. The exposure times of CCD was 100 milliseconds, and the total acquisition time for a single measurement was about 1.3 seconds.

The hyperspectral raw data was normalized to obtain the reflectance value of the specimen (${R}_{Exp}^{\lambda}$) as follow:

*λ*, ${R}_{Background}^{\lambda}$represents reflectance value of the avascular region.

The *in vivo* skin clearing technique can be used to improve the performance of HSI, and allow to monitor the distribution of microvascular *SO _{2}* more clearly. The skin optical clearing agent (SOCA) (37°C) was applied on the dorsal region of interest (ROI) of

*Babl/c*mice (28 ± 2g, 8 weeks old) for 5 minutes. The skin was effectively optical transparent so that the dermal vessels were visible. Here, the SOCA was based on our previous study, and the details can refer to the literature [17–19].

After skin was optically transparent, mouse was kept to inhale the air mixed with 1.5% isoflurane prior to the hypoxic phase, and the hyperspectral data were recorded for 3 minutes in the normal state. Then, the inhalation was switched to a mixed gas (10% oxygen + 88.5% nitrogen) together with 1.5% isoflurane to create a hypoxic stimulation, and this phase lasted for 3 minutes followed by a 4 minutes recovery period in which the inhalation was switched back to the former. The hyperspectral data was recorded every 10 seconds. Figure 1 showed the schematic of the experimental set-up.

#### 2.2 LUT-based inverse model for hyperspectral analysis

The LUT of tissue model includes the absorption, scattering and reflectance properties. MC simulation has been used for analyzing diffuse reflectance spectra as well [36]. And the LUT established by MC simulation can be used in least-squares sense to solve the nonlinear fitting (NLF) problems. Here, the LUT-NLF method was used for hyperspectral analysis. Figure 2 shows the flow chart of hyperspectral analysis for each pixel based on LUT-NLF.

### 2.2.1 Establishment of LUT for tissue model

Figure 3 shows the schematic of two-layer tissue model and one-layer tissue model used in this work. Both tissue models include absorption and scattering. For the two-layer tissue model, the absorption is from the melanin in epidermal layer and blood vessel in dermal layer, respectively. The one-layer tissue model neglect completely the epidermis, only the blood layer is considered. We established both one-layer and two-layer LUT to verify the simulated blood model and estimate the *SO _{2}* of cutaneous microvessels, respectively. Here, GPU-assisted MC simulation was used for simulating light distribution in skin, and Wang

*et al.*’s CONV program was used to obtain the diffuse reflectance (

*R*) [37].

_{simulation}For the two-layer LUT establishment, the number of photons was 10^{6}; the refractive index of the skin was 1.4; the anisotropy factor was 0.9; and the mice epidermal layer thickness was set to 20 μm (Male *Balb/c* mice, 8 weeks old) [12]. In addition, the number of absorption coefficients of epidermis and dermis was 40 and 30 values chosen from following set (${\mu}_{a.epi}^{\lambda}\in [\begin{array}{cc}0c{m}^{-1},& 100c{m}^{-1}\end{array}],{\mu}_{a.derm}^{\lambda}\in [\begin{array}{cc}0c{m}^{-1},& 50c{m}^{-1}\end{array}]$), respectively, and the number of reduced scattering coefficients was 20 values for epidermis and dermis (${\mu}_{s}^{\text{'}}(\lambda )\in [\begin{array}{cc}2c{m}^{-1},& 70c{m}^{-1}\end{array}]$). It is well known that${\mu}_{a.epi}^{\lambda}$, ${\mu}_{a.derm}^{\lambda}$ and ${\mu}_{s}^{\text{'}}(\lambda )$ are wavelength-dependent. The values used in MC simulation cover the range of absorption coefficient and scattering coefficient of skin in the range of visible-NIR wavelength [29] [39].

The main absorption in epidermis results from melanin packaged in melanosomes, which is stronger in the short wavelength. The absorption coefficient of epidermis (${\mu}_{a.epi}^{\lambda}$) could be calculated as follows [38, 39]:

*M*refer to the melanin volume fraction, which determined the epidermal absorption (

*µ*). In the dermis, the hemoglobin is main absorption contributor, which can be expressed as [38, 39]:

_{a,epi}*B*refers to the volume fraction of blood. The reduced scattering coefficient is a wavelength-dependent function as [40]:

*λ*is equal to 630 nm, and γ is the scattering power [41].

_{0}According to Zhong *et al.*’s previous work, the simulated reflectance (*R _{simulation}*) could be calibrated as follows [29]:

*k*is one of the estimated parameters.

For the one-layer LUT, only the absorption of blood was under consideration in MC simulation. Thus, the absorption of blood model (*μ _{a}*) was equal to the Eq. (4), and the scattering and other simulated parameters were the same as the two-layer LUT.

### 2.2.2 Extracting physiological parameters based on LUT-NLF

The four estimated parameters could be expressed by a property vector$\overrightarrow{P}=\u3008\begin{array}{cccc}B\u2022S{O}_{2},& B\u2022(1-S{O}_{2}),& \gamma ,& k\end{array}\u3009$. The goal of the inverse problem was to estimate the $\overrightarrow{P}$ from the reflectance at each pixel of hyperspectral image (${R}_{Exp}^{{\lambda}_{i}}$). To estimate vector$\overrightarrow{P}$, we solved this inverse problem by minimizing the sum of squared residuals *δ* expressed as follows:

*R*can be obtained by interpolation of the LUT. Here, the function “

_{LUT}*lsqcurvefit*” was used to solve this inverse nonlinear least-squares problem in MATLAB. And, the

*SO*could be calculated as follows:

_{2}If the Pearson correlation coefficient between the measured and fitted values was larger than 0.9, the estimated *SO _{2}* would be accepted, otherwise, the

*SO*should be rejected instead. To reduce the computing time-consuming, the microvascular distribution could be extracted from hyperspectral images at 540 nm by gray threshold segmentation in advance because there is strong absorption of vessels at 540 nm. The gray threshold was chosen by artificial discrimination. Additionally, the multithread parallels were employed to accelerate the LUT-NLF method in MATLAB. In this work, the time-consumption of the entire calculation was about 10 minutes for obtaining one

_{2}*SO*map. The information of the computational hardware is as follows: Intel(R) Xeon(R) E5-2620 CPU at 2. 0 GHz (two CPU), 32GB of physical RAM, windows 7. The time-consumption of calculation would be decreased if high performance computer device was used or vascular pixels is decreased.

_{2}#### 2.3 MLR and NNLS methods for estimating SO_{2}

The *SO _{2}* assessment by HSI was based upon the modified Beer Lambert law (MBLL), which can be expressed as follows [42]:

*A(λ)*is the wavelength-dependent absorbance of blood from the reflectance${R}_{Exp}^{\lambda}$;

*B*and

*SO*are the volume fraction and oxygen saturation of blood, respectively; ${\epsilon}_{oxy}^{\lambda}$ and ${\epsilon}_{deoxy}^{\lambda}$are the wavelength-dependent molar extinction coefficient of oxy- and deoxy-hemoglobin, respectively [39]; ζ is a constant of scattering;

_{2}*λ*represents the wavelength.

To solve the MBLL equation, many previous works mainly adopted the MLR and NNLS. As a matrix approach for MLR analysis was expressed as follows [28]:

*w*is the number of wavelength. The matrix notation could be simplified as:We could obtain

*C*by solving Eq. (11):So, the

*SO*for each pixel could be estimated as follows:

_{2}Here, the estimator *C* may be negative, in this case, we took the absolute value of *C* for obtaining positive *SO _{2}*. However, the NNLS was used in previous studies to avoid negative estimators by limiting the range of

*C*. Thus, we could solve Eq. (11) by NNLS as follows:

Both MLR and NNLS could be realized by matrix calculation with the help of MATLAB, and we can use “*nnls*” function in MATLAB to implement NNLS method directly. The estimated *SO _{2}* value for each pixel was accepted or rejected based on the Pearson correlation coefficient between the measured and fitted values. If the correlation coefficient was larger than 0.9, the estimated values would be accepted, otherwise, the values would be rejected instead.

#### 2.4 Simulated blood model

In order to verify the performance of LUT-NLF, we established simulated blood model based on Beer Lambert Law by means of MC simulation. In MC simulation, the number of photons, the refractive index of blood, and the anisotropy factor were set to be 10^{6}, 1.4, and 0.9, respectively. And the absorption coefficient of blood model was determined by the *SO _{2}* as follows:

*SO*values were preset as 0.3, 0.5, 0.7 and 0.9; the wavelength

_{2}*λ*ranged from 500nm to 620nm with a wavelength interval of 5 nm. The scattering coefficient was obtained through Eq. (5),

*λ*and

_{0}*γ*were set to be 630nm and 1.4, respectively. The reflectance

*R*for each

*SO*value could be obtained by the CONV program.

_{2}## 3. Results

#### 3.1 The LUT for hyperspectral analysis

MC simulation has been widely employed to simulate the light distribution in turbid media [43–46]. Figure 4 shows that established one-layer LUT and two-layer LUT based on MC simulation, which covers the range of absorption coefficient and scattering coefficient of skin in the range of visible-NIR wavelength.

For the one-layer LUT, it can be found that when the absorption coefficient (*µ _{a}*) is constant, the absorption intensity of blood model (log

_{10}(1/

*R*)) decreases with the increase of scattering coefficient (

*µ*). In contrast, when scattering coefficient (

_{s}*µ*) is constant, the absorption intensity of blood model (log

_{s}_{10}(1/

*R*)) increases with the increase of absorption coefficient (

*µ*).

_{a}For the two-layer LUT, the color-bar shows range of the absorption intensity of two-layer skin model (log_{10}(1/*R*)). If the optical parameters of one-layer LUT (derma or epidermis) is constant, the variation trend of the absorption intensity (log_{10}(1/*R*)) is same as that of the one-layer LUT. And for a constant *µ _{a,epi}* and

*µ*the absorption intensity of two-layer skin model (log

_{a,derm},_{10}(1/

*R*)) increase with the reduction of

*µ*. For the LUT of different tissue model, the variation trend of absorption (log

_{s}_{10}(1/R)) of model with absorption coefficient and scattering coefficient is similar, but the amplitude of variation of absorption (log

_{10}(1/R)) is different.

#### 3.2 Hyperspectral imaging combining with skin optical clearing technique

The skin optical clearing technique has been widely used to improve the imaging depth and contrast. Here, we employed HSI to non-invasively monitor the microvascular oxygen saturation with the assistance of skin optical clearing technique. Figures 5(a) and 5(b) showed the images at 540 nm before and after treatment with SOCA, respectively. It can be found that the image contrast can be obviously improved after SOCA treatment. The extinction coefficients of oxy- and deoxy-hemoglobin at 540 nm are strong, thus both arteries and veins at 540 nm can be clearly visualized after SOCA treatment. The extinction coefficient of deoxy-hemoglobin is much larger than that of oxy-hemoglobin at 620 nm [39], which help to distinguish arteries and veins. Therefore, based on hyperspectral images at 540 nm and 620 nm, the arteries and veins could be extracted as shown in Fig. 5(c).

#### 3.3 The changes of microvascular SO_{2} during hypoxic stimulation

Here we monitored microvascular oxygen changes during hypoxic stimulation experiment. Figure 6 shows typical *SO _{2}* maps obtained based on LUT-NLF with two-layer LUT. There were visible decreases in microvascular

*SO*of cutaneous arteries and veins during the hypoxic stimulation. After hypoxic stimulation, the

_{2}*SO*of cutaneous arteries and veins could return to normal almost.

_{2}When using MLR to estimate *SO _{2}*, the estimators for some pixels could be negative, thus under such conditions we took absolute values of estimators for obtaining positive

*SO*values. We also compared

_{2}*SO*values form the same hyperspectral data sets based on the LUT-NLF, MLR (|Estimators|) and NNLS methods. Figure 7 shows the changes of

_{2}*SO*values estimated by LUT-NLF at different time points, which correspond with the maps shown in Fig. 6. The statistical results of Fig. 7(a) and 7(b) are from the

_{2}*SO*of vascular pixels (A-1 and V-1) and relative

_{2}*SO*changes of five ROIs (A-1 to A-5 and V-1 to V-5) indicated in Fig. 5, respectively. The time axe of Fig. 7 corresponds to the time point of

_{2}*SO*maps in first column of Fig. 6, and at 190s hypoxic process, and the recovery were monitored at 190 s and 370 s, respectively. The gray shade areas indicate the hypoxic stimulation period in Fig. 7. Figure 7(a) demonstrates the typical changes of microvascular

_{2}*SO*during hypoxic stimulation in the ROIs of an artery and a vein labeled with the white and yellow rectangles (A-1 and V-1), respectively, as indicated in Fig. 5(c). The time-lapse variations of arterial

_{2}*SO*(the first row in Fig. 7(a)) and venous

_{2}*SO*(the second row in Fig. 7(a)) estimated using the three methods were similar. However, the estimated

_{2}*SO*values determined by NNLS were higher than other two methods in Fig. 7(a). And, the

_{2}*SO*decreasing extent determined by NNLS was lower than other two methods during the hypoxic process.

_{2}Further, we performed the statistical analysis of relative change in *SO _{2}* of arteries and veins for five ROIs indicated in Fig. 7(b). It can be found that venous responses to hypoxia are stronger than those of arteries according to the trends of relative changes in Fig. 7(b) (

*p*> 0.05). In addition, during the hypoxic stimulation, the microvascular relative decreasing extent of

*SO*based on LUT-NLF, MLR (|Estimators|) and NNLS methods decreased sequentially. Therefore, the estimated

_{2}*SO*based on LUT-NLF appears to be more sensitive to the hypoxic stimulation.

_{2}#### 3.4 Fitting error and correlation of LUT-NLF

To evaluate the fitting error and the correlation of LUT-NLF, NNLS, MLR and MLR (|Estimators|), one of hyperspectral data sets of hypoxic experiment was used to analyze. Figure 8(a) shows *SO _{2}* maps obtained based on LUT-NLF, MLR (|Estimators|) and NNLS. The correlation coefficient of

*SO*maps between LUT-NLF and MLR (|Estimators|) was ~0.9. And the correlation coefficient of that between LUT-NLF and NNLS was ~0.9 as well. It indicated that there was a strong correlation between

_{2}*SO*maps estimated based on LUT-NLF and the other two common methods (MLR (|Estimators|) or NNLS).

_{2}Moreover, we analyzed entire vascular pixels based on LUT-NLF, NNLS, MLR and MLR (|Estimators|) during the fitting process, as shown in Fig. 8(a). And we compared their mean absolute error (MAE) and the Pearson correlation coefficient. The MAE of each pixel can be expressed as:

*W*is the number of fitting points.

Figure. 8(b-i) shows that all the MAE based on LUT-NLF, MLR and NNLS seem to be small. To obtain positive *SO _{2}* values based on MLR, we took the absolute values of estimators. Unfortunately, it caused the MAE increase (red circles in Fig. 8(b-i)). Since we only kept the pixels with Pearson correlation coefficient greater than 0.9, Fig. 8(b-ii) shows the proportion of pixels in different ranges of Pearson correlation coefficients from 0.9 to 1 (with a step size of 0.01) in the fitting processes with LUT-NLF, NNLS, MLR and MLR (|Estimators|). And the proportion of pixels in ranges of Pearson correlation coefficients from 0.97 to 1 based on LUT-NLF, NNLS, MLR and MLR (|Estimators|) were 54.2%, 23.0%, 34.1% and 34.0%, respectively. Hence, it indicated that LUT-NLF method had a higher fitting degree during the data fitting process.

#### 3.5 Validation of the performance of LUT-NLF

To verify the performance of LUT-NLF, MLR and NNLS in estimating *SO _{2}*, we employed MC simulation of blood model based on Beer Lambert Law. This simulated blood model was previously used to verify in Ye Yang

*et al.*’ work [47]. The pure blood model, as one-layer model, was determined by pre-set

*SO*, thus, the one-layer LUT established previously was used in LUT-NLF. Figure 9(a) shows the simulated spectra with different pre-set

_{2}*SO*values 0.3, 0.5, 0.7, and 0.9. Furthermore, the estimated

_{2}*SO*values were obtained based on LUT-NLF, MLR and NNLS. And when using MLR for obtaining

_{2}*SO*from the simulated spectra, no negative estimators values appeared, thus we need not take the absolute values of estimators. Figure 9(b) reveals there are good linear correlations between the pre-set

_{2}*SO*values and the estimated

_{2}*SO*values obtained by these methods. The estimated

_{2}*SO*based on MLR and NNLS were equal. The fitting line based on LUT-NLF (black line) was much closer to the standard line (blue dotted line) than those fitting lines based on MLR and NNLS (red dashed line and green line, respectively), this means that the estimated

_{2}*SO*values obtained by LUT-NLF are closer to the pre-set

_{2}*SO*values than MLR and NNLS. Thus, the LUT-NLF has a good accuracy in determining

_{2}*SO*from this simulated spectra.

_{2}## 4. Discussion

Currently, the major methods for estimating the physiological properties in hyperspectral analysis field are MLR [21–24] and NNLS [25–27]. NNLS is a bounded-variable least-square method, which limits the range of estimators, but MLR doesn’t limit it. Besides, MLR is unlike the simple linear regression, the inferences are made about the degree of interaction or correlation between each of the independent variables [28].

Recently, Zhong *et al.* has used the two-layer LUT to extract the skin properties from the fiber-based spectroscopic diffusion reflectance [29]. In this work, we further developed this LUT-based inverse method for hyperspectral analysis of the microvascular oxyhemoglobin saturation. This LUT-based inverse model considers the multilayer structure of skin, and combines the absorption and scattering properties based on MC simulation. Therefore, LUT-NLF method could obtain a more accurate estimated *SO _{2}* values for hyperspectral analysis.

When using LUT-NLF, MLR (|Estimators|) and NNLS to estimate *SO _{2}* in the same ROIs, even though the changing tendency of

*SO*was similar during hypoxic stimulation, the decreasing degree of

_{2}*SO*were different, as shown in Fig. 7. The estimated

_{2}*SO*values based on NNLS were generally higher than that based on MLR (|Estimators|) and LUT-NLF, as shown in Fig. 7(a), that may be due to NNLS solved the least-squares problem by limiting the range of estimators (the

_{2}*SO*was from 0 to 1) and made estimated

_{2}*SO*close to 1, especially the estimated

_{2}*SO*in arteries. Thus, it could reduce the sensitivity to arterial changes of

_{2}*SO*estimated based on NNLS. There are many vascular pixels of statistics (fromA-1 and V-1) in Fig. 7(a), and the relative change

_{2}*SO*of statistics is from five areas (A1 to A5 or V1 to V5) in Fig. 7(b). Thus, the standard deviation of Fig. 7(a) is higher than that of Fig. 7(b).And, the standard deviation of Fig. 7 may be related to size of statistical areas and physiological change.

_{2}The MC simulation, as a classical method for modeling, has been widely used for analyzing the diffuse reflectance spectra [29, 36]. Here we employed the simulated blood model based on Beer Lambert law via MC simulation to verify the performance of LUT-NLF (one-layer LUT), MLR and NNLS. There was a good linear correlation between the pre-set *SO _{2}* values and estimated

*SO*values based on these methods, shown in Fig. 9(b). Thus, the estimated

_{2}*SO*values based on LUT-NLF, MLR and NNLS can effectively reflect the relative changes of

_{2}*SO*. Additionally, Fig. 9(b) showed estimated

_{2}*SO*values obtained by MLR and NNLS were same, probably because both algorithms obtained the same estimators (${({X}^{T}X)}^{-1}{X}^{T}A$) by computing the pseudoinverse, therefore the solutions may be equal in some cases. Here, we verified the fitting accuracy and some differences of these methods for calculating the same fitting data rather than caring about the absolute

_{2}*SO*value. In order to analyze the fitting processes based on these approaches, the same stimulated blood spectra were used. Therefore, at least for the simulated spectra, the LUT-NLF has a good accuracy in determining

_{2}*SO*in some way. For the further work, the absolute

_{2}*SO*value based on LUT-NLF will further carry on.

_{2}In the data fitting process based on MLR, the estimators could turn negative. Thus, we selected four typical vascular pixels to analyze based on the positive and negative values of estimators. For Pixel-1 and Pixel-2 the estimated *C(2)* and *C(3)* were positive except for the residual *C(1)* in Eq. (13). However, the estimated *C(2)* and *C(3)* at Pixel-3 and Pixel-4 were both negative. The Pixel-1 and Pixel-2 in Fig. 10 could explain the reason why the MAE of fitting process based on MLR (|Estimators|) was larger when estimators of some pixels were changed into its absolute values. There were different spectra characteristics between Pixel-1, Pixel-2 and Pixel-3, Pixel-4 in Fig. 10, which may lead to negative estimated *SO _{2}* values based on MLR. The difference of spectra characteristics in vascular pixels may be due to the noise of CCD itself.

Additionally, the estimated *SO _{2}* values of Pixel-1, Pixel-2, Pixel-3 and Pixel-4 in Fig. 10 were shown in Table 1. The results indicate that the estimated

*SO*values of Pixel-1 and Pixel-2 based on NNLS are 1 and 0.99, respectively, and that of Pixel-3 and Pixel-4 based on NNLS are 0, which is similar with results of Fig. 7. The estimated

_{2}*SO*values are too extreme to reflect sensitively changes in

_{2}*SO*. Meanwhile, estimated

_{2}*SO*values of Pixel-3 and Pixel-4 based on LUT-NLF were obtained with a good precision. It can be found that estimated

_{2}*SO*values derived from negative estimators based on MLR may be close to 0. The proportion of the number of microvascular pixels at which negative estimated

_{2}*SO*values based on MLR to all the microvascular pixels was less. Hence, the absolute estimators based on MLR (|Estimators|) causes less influence on the statistical analysis of estimated

_{2}*SO*values for a large amount of pixels, as shown in Fig. 7.

_{2}When the two-layer LUT was established here, we considered the impact of skin optical clearing agent on two-layer skin model by MC simulation. Therefore, the refractive index of the epidermis and dermis were both set as 1.4 in MC simulation because of the refractive index matching after applying skin optical clearing agent. And in this LUT-based inverse model, we assumed homogeneous tissue and epidermal thickness, which could limit accuracy of the simulated reflectance. However, compared with the common hyperspectral analysis methods, LUT-NLF method shows a better performance for estimating *SO _{2}* values. And this method can be applied to tissues other than skin. We can establish different LUT by considering different optical properties and tissue structure for hyperspectral analysis. Additionally, a similar inverse method was used to improve to predict accuracy of the diabetic foot ulcer region for HSI [33]. Thus, LUT-based inverse model for hyperspectral analysis may have potential applications for non-invasive clinical diagnoses.

## 5. Conclusion

In this study, we used HSI to non-invasively monitor microvascular *SO _{2}* changes under the assistance of

*in vivo*skin optical clearing technique. To estimate

*SO*values from experimental data, we used LUT-NLF based on a two-layer LUT for hyperspectral analysis. Further, this method was compared with the other two common methods (MLR and NNLS) under the hypoxic stimulation experiment. The lower MAE and higher Pearson correlation coefficient indicted that LUT-NLF method had a higher fitting degree. For the hypoxic stimulation, the decreased degree of relative changes in

_{2}*SO*based on LUT-NLF was largest among these methods, which appeared to indicate that estimated

_{2}*SO*based on LUT-NLF was more sensitive to the fluctuation of

_{2}*SO*. In addition, the simulated blood model was used to verify, which indicated that estimated

_{2}*SO*values based on LUT-NLF were much closer to pre-set

_{2}*SO*values. That means that LUT-NLF has a good accuracy in determining

_{2}*SO*from the simulated spectra. And the LUT-based inverse model for hyperspectral analysis can take into consideration of the complex tissue structure and optical properties, which is expected to be applied in clinical tissue constituent monitoring and diseases diagnoses.

_{2}## Funding

This study was supported by the National Natural Science Foundation of China (NSFC) (31571002), the Science Fund for Creative Research Group of China (61421064) and the seed project of Wuhan National Laboratory for Optoelectronics.

## Acknowledgements

We would like to thank Valery V. Tuchin (Saratov National Research State University) for their helpful discussion. Especially, I would like to express my appreciation to Prof. Lin Z. Li (University of Pennsylvania) for helping us modify the manuscript.

## References and links

**1. **J. A. Lee, R. T. Kozikowski, and B. S. Sorg, “In vivo microscopy of microvessel oxygenation and network connections,” Microvasc. Res. **98**, 29–39 (2015). [CrossRef] [PubMed]

**2. **S. Schöneborn, C. Vollmer, F. Barthel, A. Herminghaus, J. Schulz, I. Bauer, C. Beck, and O. Picker, “Vasopressin V1A receptors mediate the stabilization of intestinal mucosal oxygenation during hypercapnia in septic rats,” Microvasc. Res. **106**, 24–30 (2016). [CrossRef] [PubMed]

**3. **B. S. Sorg, B. J. Moeller, O. Donovan, Y. Cao, and M. W. Dewhirst, “Hyperspectral imaging of hemoglobin saturation in tumor microvasculature and tumor hypoxia development,” J. Biomed. Opt. **10**(4), 44004 (2005). [CrossRef] [PubMed]

**4. **M. W. Dewhirst, “Concepts of oxygen transport at the microcirculatory level,” Semin. Radiat. Oncol. **8**(3), 143–150 (1998). [CrossRef] [PubMed]

**5. **S. Miclos, S. V. Parasca, M. A. Calin, D. Savastru, and D. Manea, “Algorithm for mapping cutaneous tissue oxygen concentration using hyperspectral imaging,” Biomed. Opt. Express **6**(9), 3420–3430 (2015). [CrossRef] [PubMed]

**6. **G. Lu and B. Fei, “Medical hyperspectral imaging: a review,” J. Biomed. Opt. **19**(1), 10901 (2014). [CrossRef] [PubMed]

**7. **M. A. Calin, S. V. Parasca, D. Savastru, and D. Manea, “Hyperspectral Imaging in the Medical Field: Present and Future,” Appl. Spectrosc. Rev. **49**(6), 435–447 (2013). [CrossRef]

**8. **K. Kikuchi, Y. Masuda, and T. Hirao, “Imaging of hemoglobin oxygen saturation ratio in the face by spectral camera and its application to evaluate dark circles,” Skin Res. Technol. **19**(4), 499–507 (2013). [PubMed]

**9. **D. Yudovsky, A. Nouvong, and L. Pilon, “Hyperspectral imaging in diabetic foot wound care,” J. Diabetes Sci. Technol. **4**(5), 1099–1113 (2010). [CrossRef] [PubMed]

**10. **A. Basiri, M. Nabili, S. Mathews, A. Libin, S. Groah, H. J. Noordmans, and J. C. Ramella-Roman, “Use of a multi-spectral camera in the characterization of skin wounds,” Opt. Express **18**(4), 3244–3257 (2010). [CrossRef] [PubMed]

**11. **S. A. Shah, N. Bachrach, S. J. Spear, D. S. Letbetter, R. A. Stone, R. Dhir, J. W. Prichard, H. G. Brown, and W. A. LaFramboise, “Cutaneous wound analysis using hyperspectral imaging,” Biotechniques **34**(2), 408–413 (2003). [PubMed]

**12. **L. Guo, R. Shi, C. Zhang, D. Zhu, Z. Ding, and P. Li, “Optical coherence tomography angiography offers comprehensive evaluation of skin optical clearing in vivo by quantifying optical properties and blood flow imaging simultaneously,” J. Biomed. Opt. **21**(8), 081202 (2016). [CrossRef] [PubMed]

**13. **K. V. Larin, M. G. Ghosn, A. N. Bashkatov, E. A. Genina, N. A. Trunina, and V. V. Tuchin, “Optical Clearing for OCT Image Enhancement and In-Depth Monitoring of Molecular Diffusion,” IEEE J. Sel. Top. Quantum Electron. **18**(3), 1244–1259 (2012). [CrossRef]

**14. **X. Wen, S. L. Jacques, V. V. Tuchin, and D. Zhu, “Enhanced optical clearing of skin in vivo and optical coherence tomography in-depth imaging,” J. Biomed. Opt. **17**(6), 066022 (2012). [CrossRef] [PubMed]

**15. **Y. Zhou, J. Yao, and L. V. Wang, “Optical clearing-aided photoacoustic microscopy with enhanced resolution and imaging depth,” Opt. Lett. **38**(14), 2592–2595 (2013). [CrossRef] [PubMed]

**16. **Y. Liu, X. Yang, D. Zhu, R. Shi, and Q. Luo, “Optical clearing agents improve photoacoustic imaging in the optical diffusive regime,” Opt. Lett. **38**(20), 4236–4239 (2013). [CrossRef] [PubMed]

**17. **D. Zhu, J. Wang, Z. Zhi, X. Wen, and Q. Luo, “Imaging dermal blood flow through the intact rat skin with an optical clearing method,” J. Biomed. Opt. **15**(2), 026008 (2010). [CrossRef] [PubMed]

**18. **J. Wang, R. Shi, and D. Zhu, “Switchable skin window induced by optical clearing method for dermal blood flow imaging,” J. Biomed. Opt. **18**(6), 061209 (2013). [CrossRef] [PubMed]

**19. **R. Shi, M. Chen, V. V. Tuchin, and D. Zhu, “Accessing to arteriovenous blood flow dynamics response using combined laser speckle contrast imaging and skin optical clearing,” Biomed. Opt. Express **6**(6), 1977–1989 (2015). [CrossRef] [PubMed]

**20. **J. Wang, N. Ma, R. Shi, Y. Zhang, T. Yu, and D. Zhu, “Sugar-induced skin optical clearing: from molecular dynamics simulation to experimental demonstration,” IEEE J. Sel. Top. Quantum Electron. **20**(2), 1–7 (2014). [CrossRef]

**21. **J. Wang, Y. Zhang, P. Li, Q. Luo, and D. Zhu, “Review: tissue optical clearing window for blood flow monitoring,” IEEE J. Sel. Top. Quantum Electron. **20**(2), 1–12 (2014). [CrossRef]

**22. **H. Ding and C. R. Chang, “Comparison of Photometric Stereo and Spectral Analysis for Visualization and Assessment of Burn Injury from Hyperspectral Imaging,” in *Computational Intelligence and Virtual Environments for Measurement Systems and Applications (CIVEMSA)* (IEEE, 2015), pp. 1–6.

**23. **K. J. Zuzak, M. D. Schaeberle, E. N. Lewis, and I. W. Levin, “Visible reflectance hyperspectral imaging: characterization of a noninvasive, in vivo system for determining tissue perfusion,” Anal. Chem. **74**(9), 2021–2028 (2002). [CrossRef] [PubMed]

**24. **I. Nishidate, N. Tanaka, T. Kawase, T. Maeda, T. Yuasa, Y. Aizu, T. Yuasa, and K. Niizeki, “Noninvasive imaging of human skin hemodynamics using a digital red-green-blue camera,” J. Biomed. Opt. **16**(8), 086012 (2011). [CrossRef] [PubMed]

**25. **D. W. Allen, J.-P. Bouchard, J. Wang, P. Ghassemi, A. Melchiorri, J. Ramella-Roman, S. A. Mathews, J. Coburn, B. Sorg, Y. Chen, and J. Pfefer, “3D printed biomimetic vascular phantoms for assessment of hyperspectral imaging systems,” Proc. SPIE **9325**, 932508 (2015). [CrossRef]

**26. **E. M. Hillman, A. Devor, M. B. Bouchard, A. K. Dunn, G. W. Krauss, J. Skoch, B. J. Bacskai, A. M. Dale, and D. A. Boas, “Depth-resolved optical imaging and microscopy of vascular compartment dynamics during somatosensory stimulation,” Neuroimage **35**(1), 89–104 (2007). [CrossRef] [PubMed]

**27. **G. M. Palmer, A. N. Fontanella, S. Shan, G. Hanna, G. Zhang, C. L. Fraser, and M. W. Dewhirst, “In vivo optical molecular imaging and analysis in mice using dorsal window chamber models applied to hypoxia, vasculature and fluorescent reporters,” Nat. Protoc. **6**(9), 1355–1366 (2011). [CrossRef] [PubMed]

**28. **S. H. Brown, “Multiple Linear Regression Analysis: A Matrix Approach with MATLAB,” Alabama J. Math. **34**, 1–3 (2009).

**29. **X. Zhong, X. Wen, and D. Zhu, “Lookup-table-based inverse model for human skin reflectance spectroscopy: two-layered Monte Carlo simulations and experiments,” Opt. Express **22**(2), 1852–1864 (2014). [CrossRef] [PubMed]

**30. **N. Rajaram, T. H. Nguyen, and J. W. Tunnell, “Lookup table-based inverse model for determining optical properties of turbid media,” J. Biomed. Opt. **13**(5), 050501 (2008). [CrossRef] [PubMed]

**31. **L. Lim, B. Nichols, N. Rajaram, and J. W. Tunnell, “Probe pressure effects on human skin diffuse reflectance and fluorescence spectroscopy measurements,” J. Biomed. Opt. **16**(1), 011012 (2011). [CrossRef] [PubMed]

**32. **I. Fredriksson, M. Larsson, and T. Strömberg, “Inverse Monte Carlo method in a multilayered tissue model for diffuse reflectance spectroscopy,” J. Biomed. Opt. **17**(4), 047004 (2012). [CrossRef] [PubMed]

**33. **D. Yudovsky, A. Nouvong, K. Schomacker, and L. Pilon, “Monitoring temporal development and healing of diabetic foot ulceration using hyperspectral imaging,” J. Biophotonics **4**(7-8), 565–576 (2011). [CrossRef] [PubMed]

**34. **D. Yudovsky and L. Pilon, “Simple and accurate expressions for diffuse reflectance of semi-infinite and two-layer absorbing and scattering media,” Appl. Opt. **48**(35), 6670–6683 (2009). [CrossRef] [PubMed]

**35. **Z. Qian, S. Victor, Y. Gu, C. Giller, and H. Liu, “Look-Ahead Distance of a fiber probe used to assist neurosurgery: Phantom and Monte Carlo study,” Opt. Express **11**(16), 1844–1855 (2003). [CrossRef] [PubMed]

**36. **G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt. **45**(5), 1062–1071 (2006). [CrossRef] [PubMed]

**37. **L. Wang, S. L. Jacques, and L. Zheng, “CONV--convolution for responses to a finite diameter photon beam incident on multi-layered tissues,” Comput. Methods Programs Biomed. **54**(3), 141–150 (1997). [CrossRef] [PubMed]

**38. **D. Yudovsky and L. Pilon, “Rapid and accurate estimation of blood saturation, melanin content, and epidermis thickness from spectral diffuse reflectance,” Appl. Opt. **49**(10), 1707–1719 (2010). [CrossRef] [PubMed]

**39. **S. A. Prahl, “Optical absorption of hemoglobin,” http://omlc.org/spectra/hemoglobin/summary.html.

**40. **R. T. Zaman, N. Rajaram, B. S. Nichols, H. G. Rylander 3rd, T. Wang, J. W. Tunnell, and A. J. Welch, “Changes in morphology and optical properties of sclera and choroidal layers due to hyperosmotic agent,” J. Biomed. Opt. **16**(7), 077008 (2011). [CrossRef] [PubMed]

**41. **A. Cerussi, N. Shah, D. Hsiang, A. Durkin, J. Butler, and B. J. Tromberg, “In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broadband diffuse optical spectroscopy,” J. Biomed. Opt. **11**(4), 044005 (2006). [CrossRef] [PubMed]

**42. **R. D. Shonat, E. S. Wachman, W. Niu, A. P. Koretsky, and D. L. Farkas, “Near-simultaneous hemoglobin saturation and oxygen tension maps in mouse brain using an AOTF microscope,” Biophys. J. **73**(3), 1223–1231 (1997). [CrossRef] [PubMed]

**43. **C. Jiang, H. He, P. Li, and Q. Luo, “Graphics processing unit cluster accelerated monte carlo simulation of photon transport in multi-layered tissues,” J. Innov. Opt. Health Sci. **5**(5), 476–482 (2012).

**44. **C. Zhu and Q. Liu, “Validity of the semi-infinite tumor model in diffuse reflectance spectroscopy for epithelial cancer diagnosis: a Monte Carlo study,” Opt. Express **19**(18), 17799–17812 (2011). [CrossRef] [PubMed]

**45. **B. Luo and S. He, “An improved Monte Carlo diffusion hybrid model for light reflectance by turbid media,” Opt. Express **15**(10), 5905–5918 (2007). [CrossRef] [PubMed]

**46. **L. Wang, S. L. Jacques, and L. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. **47**(2), 131–146 (1995). [CrossRef] [PubMed]

**47. **Y. Yang, O. Soyemi, P. J. Scott, M. R. Landry, S. M. Lee, L. Stroud, and B. R. Soller, “Quantitative measurement of muscle oxygen saturation without influence from skin and fat using continuous-wave near infrared spectroscopy,” Opt. Express **15**(21), 13715–13730 (2007). [CrossRef] [PubMed]