Lookup-table-based inverse model for mapping oxygen concentration of cutaneous microvessels using hyperspectral imaging

Wei Feng; Rui Shi; Chao Zhang; Tingting Yu; Dan Zhu

doi:10.1364/OE.25.003481

1. Introduction

The microvascular oxygen saturation (SO₂) 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 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₂ 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.

Fig. 1 The schematic of the experimental set-up.

Download Full Size | PDF

The hyperspectral raw data was normalized to obtain the reflectance value of the specimen ( $R_{E x p}^{λ}$ ) as follow:

R^{λ} = (I_{r a w}^{λ} - I_{b l a c k}^{λ}) / (I_{s t d}^{λ} - I_{b l a c k}^{λ}) .

R_{E x p}^{λ} = R^{λ} / R_{B a c k g r o u n d}^{λ} .

where

I_{r a w}^{λ}

is the raw data radiance value of a given pixel,

I_{s t d}^{λ}

and

I_{b l a c k}^{λ}

are the radiance values of the certified reflectance standard and the dark current for each spectral band, respectively.

R^{λ}

is the calculated reflectance value for each wavelength λ,

R_{B a c k g r o u n d}^{λ}

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

Fig. 2 The flow chart of LUT-NLF.

Download Full Size | PDF

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₂ 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_simulation) [37].

Fig. 3 Simplified tissue geometry.

Download Full Size | PDF

For the two-layer LUT establishment, the number of photons was 10⁶; 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 ( $μ_{a . e p i}^{λ} \in [\begin{matrix} 0 c m^{- 1}, & 100 c m^{- 1} \end{matrix}], μ_{a . d e r m}^{λ} \in [\begin{matrix} 0 c m^{- 1}, & 50 c m^{- 1} \end{matrix}]$ ), respectively, and the number of reduced scattering coefficients was 20 values for epidermis and dermis ( $μ_{s}^{'} (λ) \in [\begin{matrix} 2 c m^{- 1}, & 70 c m^{- 1} \end{matrix}]$ ). It is well known that $μ_{a . e p i}^{λ}$ , $μ_{a . d e r m}^{λ}$ and $μ_{s}^{'} (λ)$ 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 ( $μ_{a . e p i}^{λ}$ ) could be calculated as follows [38, 39]:

μ_{a . e p i}^{λ} = 6.60 \times 10^{11} λ^{- 3.33} • M + 7.84 \times 10^{8} λ^{- 3.255} • (1 - M) .

where M refer to the melanin volume fraction, which determined the epidermal absorption (µ_a,epi). In the dermis, the hemoglobin is main absorption contributor, which can be expressed as [38, 39]:

μ_{a . d e r m}^{λ} = B • (ε_{o x y}^{λ} • S O_{2} + ε_{d e o x y}^{λ} • (1 - S O_{2})) .

where

ε_{o x y}^{λ}

and

ε_{d e o x y}^{λ}

are wavelength-dependent molar extinction coefficient of oxy- and deoxy-hemoglobin [39], respectively. B refers to the volume fraction of blood. The reduced scattering coefficient is a wavelength-dependent function as [40]:

μ_{s}^{,} (λ) = μ_{s}^{,} (λ_{0}) • {(λ / λ_{0})}^{- γ} .

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

According to Zhong et al.’s previous work, the simulated reflectance (R_simulation) could be calibrated as follows [29]:

R_{s i m u l a t i o n}^{c a l i b r a t i o n} = 10^{k} • R_{s i m u l a t i o n} .

where the factor 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 $\vec{P} = 〈 \begin{matrix} B • S O_{2}, & B • (1 - S O_{2}), & γ, & k \end{matrix} 〉$ . The goal of the inverse problem was to estimate the $\vec{P}$ from the reflectance at each pixel of hyperspectral image ( $R_{E x p}^{λ_{i}}$ ). To estimate vector $\vec{P}$ , we solved this inverse problem by minimizing the sum of squared residuals δ expressed as follows:

δ = \sum_{i = 1}^{w} {[R_{L U T} (\vec{P}, λ_{i}) - R_{E x p}^{λ_{i}}]}^{2} .

where R_LUT can be obtained by interpolation of the LUT. Here, the function “lsqcurvefit” was used to solve this inverse nonlinear least-squares problem in MATLAB. And, the SO₂ could be calculated as follows:

S O_{2} = \vec{P} (1) / (\vec{P} (1) + \vec{P} (2)) = B • S O_{2} / (B • S O_{2} + B • (1 - S O_{2})) .

If the Pearson correlation coefficient between the measured and fitted values was larger than 0.9, the estimated SO₂ 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 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.3 MLR and NNLS methods for estimating SO₂

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

A (λ) = \log \frac{1}{R_{E x p}^{λ}} = B • (ε_{o x y}^{λ} • S O_{2} + ε_{d e o x y}^{λ} • (1 - S O_{2})) + ξ .

where A(λ) is the wavelength-dependent absorbance of blood from the reflectance

R_{E x p}^{λ}

; B and SO₂ are the volume fraction and oxygen saturation of blood, respectively;

ε_{o x y}^{λ}

and

ε_{d e o x y}^{λ}

are the wavelength-dependent molar extinction coefficient of oxy- and deoxy-hemoglobin, respectively [39]; ζ is a constant of scattering; λ 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]:

[\begin{matrix} A (λ_{1}) \\ A (λ_{2}) \\ . \\ . \\ A (λ_{w}) \end{matrix}] \approx [\begin{matrix} 1 & ε_{o x y}^{λ_{1}} & ε_{d e o x y}^{λ_{1}} \\ 1 & ε_{o x y}^{λ_{2}} & ε_{d e o x y}^{λ_{2}} \\ . & . & . \\ . & . & . \\ 1 & ε_{o x y}^{λ_{n}} & ε_{d e o x y}^{λ n} \end{matrix}] [\begin{matrix} ξ \\ B • S O_{2} \\ B • (1 - S O_{2}) \end{matrix}] .

where w is the number of wavelength. The matrix notation could be simplified as:

A \approx X • C .

We could obtain C by solving Eq. (11):

C = {(X^{T} X)}^{- 1} X^{T} A .

So, the SO₂ for each pixel could be estimated as follows:

S O_{2} = C (2) / (C (2) + C (3)) = B • S O_{2} / (B • S O_{2} + B • (1 - S O_{2})) .

Here, the estimator C may be negative, in this case, we took the absolute value of C for obtaining positive SO₂. 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:

\min \sum_{i = 1}^{w} {(C • X (λ_{i}) - A (λ_{i}))}^{2}, \begin{matrix} B \in [\begin{matrix} 0 & + \infty \end{matrix}), & S O_{2} \in [\begin{matrix} 0 & 1 \end{matrix}] . \end{matrix}

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₂ 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⁶, 1.4, and 0.9, respectively. And the absorption coefficient of blood model was determined by the SO₂ as follows:

μ_{a, b l o o d}^{λ} = μ_{o x y}^{λ} • S O_{2} + μ_{d e o x y}^{λ} • (1 - S O_{2}) .

where

μ_{o x y}^{λ}

and

μ_{d e o x y}^{λ}

are wavelength-dependent molar extinction coefficients of oxy- and deoxy-hemoglobin, respectively [39]; the SO₂ values were preset as 0.3, 0.5, 0.7 and 0.9; the wavelength λ ranged from 500nm to 620nm with a wavelength interval of 5 nm. The scattering coefficient was obtained through Eq. (5), λ₀ and γ were set to be 630nm and 1.4, respectively. The reflectance R for each SO₂ value could be obtained by the CONV program.

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.

Fig. 4 The LUTs established based on MC simulation. (a) One-layer LUT. (b) Two-layer LUT (R represents the simulated diffusion reflectance (R_simulation)) (µ_a,epi and µ_a,derm represent the absorption coefficients of epidermis and dermis).

Download Full Size | PDF

For the one-layer LUT, it can be found that when the absorption coefficient (µ_a) is constant, the absorption intensity of blood model (log₁₀(1/R)) decreases with the increase of scattering coefficient (µ_s). In contrast, when scattering coefficient (µ_s) is constant, the absorption intensity of blood model (log₁₀(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₁₀(1/R)). If the optical parameters of one-layer LUT (derma or epidermis) is constant, the variation trend of the absorption intensity (log₁₀(1/R)) is same as that of the one-layer LUT. And for a constant µ_a,epi and µ_a,derm, the absorption intensity of two-layer skin model (log₁₀(1/R)) increase with the reduction of µ_s. For the LUT of different tissue model, the variation trend of absorption (log₁₀(1/R)) of model with absorption coefficient and scattering coefficient is similar, but the amplitude of variation of absorption (log₁₀(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).

Fig. 5 The skin imaging before and after treatment with SOCA. (a) The image at 540nm before treatment with SOCA. (b) The image at 540nm after treatment with SOCA. (c) The arteries and veins are extracted from hyperspectral images (The labeled (A-1 to A-5, V-1 to V-5) white and yellow rectangles represent the ROIs of artery and vein for SO₂ statistic analysis, respectively).

Download Full Size | PDF

3.3 The changes of microvascular SO₂ during hypoxic stimulation

Here we monitored microvascular oxygen changes during hypoxic stimulation experiment. Figure 6 shows typical SO₂ 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 SO₂ of cutaneous arteries and veins could return to normal almost.

Fig. 6 The spatiotemporal changes of SO₂ maps obtained based on LUT-NLF in cutaneous microvessels before hypoxia, during hypoxia and recovery. (The color bar represents SO₂).

Download Full Size | PDF

When using MLR to estimate SO₂, 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 SO₂ values form the same hyperspectral data sets based on the LUT-NLF, MLR (|Estimators|) and NNLS methods. Figure 7 shows the changes of 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 SO₂ of vascular pixels (A-1 and V-1) and relative 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 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 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 SO₂ (the first row in Fig. 7(a)) and venous SO₂ (the second row in Fig. 7(a)) estimated using the three methods were similar. However, the estimated SO₂ values determined by NNLS were higher than other two methods in Fig. 7(a). And, the SO₂ decreasing extent determined by NNLS was lower than other two methods during the hypoxic process.

Fig. 7 The time-lapse changes of microvascular oxyhemoglobin saturation in hypoxic stimulation. (a) The typical changes of SO₂ values in the arterious A-1 ROI (the first row) and venous V-1 ROI (the second row) (Mean ± Standard Deviation). (b) The relative changes of SO₂ in the arterious ROIs (form A-1 to A-5) and venous ROIs (form V-1 to V-5), and the curves and height of the colored shadows indicate mean and standard deviation, respectively. In (a) and (b) gray shade areas indicates the hypoxic stimulation period. Additionally, each column represents different methods (LUT-NLF, MLR (|Estimators|) and NNLS). (Hypoxic process was monitored at 190s, and process of recovery was monitored at 370s as shown on time axes in statistical graphs)

Download Full Size | PDF

Further, we performed the statistical analysis of relative change in SO₂ 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 SO₂ based on LUT-NLF appears to be more sensitive to the hypoxic stimulation.

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₂ 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 SO₂ maps estimated based on LUT-NLF and the other two common methods (MLR (|Estimators|) or NNLS).

Fig. 8 Fitting error and correlation of LUT-NLF, NNLS, MLR and MLR (|Estimators|). (a) The SO₂ maps obtained from the same hyperspectral data by LUT-NLF, MLR (|Estimators|) and NNLS (The color bar represents the SO₂). (b) Analysis on the same vascular pixels in (a). (i) The polar axis represents the MAE, and the different color marked circles represent different methods used. (ii) The proportion of pixels in different ranges of the Pearson correlation coefficients based on LUT-NLF, NNLS, MLR and MLR (|Estimators)|.

Download Full Size | PDF

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:

M A E = \frac{1}{W} \sum_{i = 1}^{w} | A_{E x p}^{λ_{i}} - \overset{\land}{A_{λ_{i}}} | .

Where

A_{E x p}^{λ_{i}}

and

\overset{\land}{A_{λ_{i}}}

represent the wavelength-dependent measured and fitted values, respectively; 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₂ 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₂, 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 SO₂ values 0.3, 0.5, 0.7, and 0.9. Furthermore, the estimated SO₂ values were obtained based on LUT-NLF, MLR and NNLS. And when using MLR for obtaining 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 SO₂ values and the estimated SO₂ values obtained by these methods. The estimated 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 SO₂ values obtained by LUT-NLF are closer to the pre-set SO₂ values than MLR and NNLS. Thus, the LUT-NLF has a good accuracy in determining SO₂ from this simulated spectra.

Fig. 9 The simulated blood model for verifying the analysis methods. (a) The simulated spectra of different SO₂ models (SO₂ was preset as 0.3, 0.5, 0.7, and 0.9). (b) The correlational analysis between pre-set SO₂ values and estimated SO₂ values.

Download Full Size | PDF

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₂ values for hyperspectral analysis.

When using LUT-NLF, MLR (|Estimators|) and NNLS to estimate SO₂ in the same ROIs, even though the changing tendency of SO₂ was similar during hypoxic stimulation, the decreasing degree of SO₂ were different, as shown in Fig. 7. The estimated 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 SO₂ was from 0 to 1) and made estimated SO₂ close to 1, especially the estimated SO₂ in arteries. Thus, it could reduce the sensitivity to arterial changes of 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 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.

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₂ values and estimated SO₂ values based on these methods, shown in Fig. 9(b). Thus, the estimated SO₂ values based on LUT-NLF, MLR and NNLS can effectively reflect the relative changes of SO₂. Additionally, Fig. 9(b) showed estimated 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 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 SO₂ in some way. For the further work, the absolute SO₂ value based on LUT-NLF will further carry on.

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₂ values based on MLR. The difference of spectra characteristics in vascular pixels may be due to the noise of CCD itself.

Fig. 10 Analysis on the data fitting processes of four vascular pixels. The estimated SO₂ of Pixel-1 and Pixel-2 are positive values. Conversely, the estimated SO₂ of Pixel-3 and Pixel-4 are negative values.

Download Full Size | PDF

Additionally, the estimated SO₂ 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 SO₂ values are too extreme to reflect sensitively changes in SO₂. Meanwhile, estimated SO₂ values of Pixel-3 and Pixel-4 based on LUT-NLF were obtained with a good precision. It can be found that estimated 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 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 SO₂ values for a large amount of pixels, as shown in Fig. 7.

Table 1. The estimated SO₂ of Pixel-1, Pixel-2, Pixel-3 and Pixel-4 from Fig. 10.

View Table

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₂ 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₂ 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 SO₂ based on LUT-NLF was largest among these methods, which appeared to indicate that estimated SO₂ based on LUT-NLF was more sensitive to the fluctuation of SO₂. In addition, the simulated blood model was used to verify, which indicated that estimated SO₂ values based on LUT-NLF were much closer to pre-set SO₂ values. That means that LUT-NLF has a good accuracy in determining 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.

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]

Methods	*SO₂*
Methods	Pixel-1	Pixel-2	Pixel-3	Pixel-4
LUT-NLF	0.76	0.72	2.78 × 10⁻⁷	4.31 × 10⁻⁵
MLR	0.65	0.69	−0.02	−0.15
MLR (\|Estimators\|)	0.65	0.69	0.02	0.11
NNLS	1	0.99	0	0

Lookup-table-based inverse model for mapping oxygen concentration of cutaneous microvessels using hyperspectral imaging

Abstract

1. Introduction