1. Introduction

The microcirculation is that part of the vascular tree that comprises blood vessels smaller than 100 μm in diameter, including arterioles, capillaries, and venules. One of its primary roles is to distribute oxygen throughout tissues. The main oxygen carrier within the circulatory system is red blood cells (RBCs), and the oxygen transport rate in various organs is controlled by microvascular geometry, hemodynamics, and RBC hemoglobin (Hb) oxygen saturation. Since diseases and injuries can lead to certain local modifications in microvasculature, and hence microvascular oxygen transport, local tissue may be hypoxic despite normal global Hb oxygenation. For example, sepsis leads to the termination of flow in a significant number of capillaries, such that the microcirculation fails to compensate for decreased functional capillary density [1]. This change elevates the likelihood of local organ failure. Furthermore, the steady state level of Hb oxygenation in the microcirculation represents the balance between cellular O₂ consumption and the supply of oxygenated blood. This, in turn, provides indirect access to the metabolic characteristics of local tissue, an important physiological indicator of tissue viability. Therefore, understanding and assessing Hb oxygenation in the microcirculation provides valuable insights into a given patient’s level of health. It is not surprising, then, that local Hb oxygenation is highly-sought information in medical research, as well as in clinical care [2–7].

Common in vitro Hb oxygenation sensors can be divided into two categories: electrochemical sensors and optical sensors. In the field of medicine, the optical sensors often are selected because of their non-destructive and/or non-invasive nature. The optical sensors utilize the oxygenation-dependent absorption properties of Hb to quantify its oxygenation. Depending upon its implementation, transmitted light at two or more wavelengths will be captured. When the assumptions of constant Hb concentration and optical pathlength are valid, a simple Beer-Lambert law may be used to convert light attenuation to the relative absorption properties of Hb and, hence, Hb oxygenation. A conventional CO-oximeter is typical of such an approach [8].

Accurately measuring Hb oxygenation from in vivo tissue is a difficult task because biological tissues are highly scattering. In addition, diffusely reflected light signals (i.e., diffuse reflectance spectroscopy) often are measured because of several limitations imposed by tissue geometry and measurement environment. These two conditions, collectively, complicate the process of assessing in vivo Hb oxygenation; scattering-induced light attenuation and pathlength alteration have to be considered by the photon migration model that converts diffuse reflectance signals to optical properties and then to Hb oxygenation. Several groups have modified Beer-Lambert law either empirically or theoretically to accommodate those scattering induced effects [9,10]. The more sophisticated models have been proposed, including the diffuse approximation model [11], the Monte Carlo (MC) based inverse model [12], the optical pathlength estimation model [13,14], and the empirical model [15,16]. Within these models, a set of coefficients is used to define a diffuse reflectance spectrum, and these coefficients are adjusted systematically until the model-predicted diffuse reflectance spectrum closely resembles the one actually measured. From the optimal coefficients, the absorption properties of the measured subject can be estimated and its Hb oxygenation determined.

In this article, we propose a new algorithm that estimates Hb oxygenation by analyzing the Hb oxygenation-induced spectral profile alterations in diffuse reflectance spectra within the visible wavelength region. The algorithm utilizes the differential wavelet transform (DWT) to enhance the tracking of the spectral profile alterations and to reduce and monitor noises in the spectra at the same time. The validation and evaluation using MC simulation and tissue phantoms concluded that the algorithm is applicable to various tissue types (i.e., scattering properties) and can be easily used in conjunction with a diverse range of probe geometries for real-time monitoring of Hb oxygenation without adjustment.

2. Methods and theory

The Hb oxygenation extraction algorithm proposed in this study was developed based on an important absorption characteristic of Hb: the profile of the extinction coefficient spectrum of Hb from 450 nm to 600 nm alters drastically between the oxygenated and the deoxygenated states. The extinction coefficient spectrum of Hb at a given oxygen saturation (SatO₂) level can be calculated using Eq. (1).

where λ is wavelength [nm], ${\epsilon }^{oxy}(\lambda )$is the extinction coefficient spectrum of oxy-Hb [liter/mole/cm], and ${\epsilon }^{deoxy}(\lambda )$is the extinction coefficient spectrum of deoxy-Hb [liter/mole/cm]. When comparing ${\epsilon }^{Sat{O}_2}(\lambda )$ at various SatO₂ levels, it clearly shows that SatO₂ level influences both the spectral feature types (e.g., peaks and valleys) as well as their locations in ${\epsilon }^{Sat{O}_2}(\lambda )$. Therefore SatO₂ information may be extracted from ${\epsilon }^{Sat{O}_2}(\lambda )$ by quantifying the location of a specific spectral feature in ${\epsilon }^{Sat{O}_2}(\lambda )$, which is the central theorem of the proposed Hb oxygenation extraction algorithm.

To facilitate the identification and localization of a specific spectral feature in ${\epsilon }^{SatO2}(\lambda )$, a DWT technique was employed in the Hb oxygenation extraction algorithm. DWT of ${\epsilon }^{SatO2}(\lambda )$ is denoted by $W_s^n{\epsilon }^{SatO2}(u)$ where u is a translation, s is a scale parameter, and n is a differential parameter. From a signal processing point of view, the wavelet transform resembles a windowed band-pass filter with a passing band, its center frequency and a window spread controlled by s. Moreover, the wavelet transform also can detect and measure variations in the input signal in the area surrounding u whose size is proportional to s. The detailed definition and explanation of DWT can be found in Appendix. In the algorithm, DWT removes those unwanted frequency components from the input signal (i.e., processed signal) and various n values were used in DWT to identify different spectral feature types in ${\epsilon }^{SatO2}(\lambda )$; zero-crossing points of $W_s^n{\epsilon }^{SatO2}(u)$, denoted by $Z_{j,s,n}^{\epsilon }$, were used to localize these features [17]. Specifically, $Z_{j,s,n=1}^{\epsilon }$ locates the j-th local maximum or minimum in filtered ${\epsilon }^{SatO2}(\lambda )$; $Z_{j,s,n=2}^{\epsilon }$ locates the j-th inflection point; $Z_{j,s,n=3}^{\epsilon }$ locates the j-th local maximum curvature.

Based on the central theorem mentioned above, the variations of $Z_{j,s,n}^{\epsilon }$ (i.e., movements of the spectral features in ${\epsilon }^{SatO2}(\lambda )$) should be strongly related to the SatO₂ level of Hb. This relationship can be formulated by a function $S_{j,s,n}$ that is written as:

In order to apply the relationship derived in Eq. (2) to diffuse reflectance spectra, preprocessing of diffuse reflectance spectra is necessary. According to the theory of photon migration, absorption events reduce the weight of a photon package during the migration process. Therefore, the remaining energy of a diffusely reflected photon package may be described as:

where r is the distance between the source and detector [cm], E ₀ is the initial energy of the photon package [J], r ₀ is the entrance point of the photon package (i.e., source location), $E(r,\lambda )$is the remaining energy of the photon package [J] when it remerges at r, ${\mu }_a(\lambda )$ is the absorption coefficient of the medium [1/cm], and $l(r,\lambda )$ is the migration pathlength of the photon package [cm]. Note that $l(r,\lambda )$ is a function of the reduced scattering properties of the medium ${\mu }_s{}^{\text{'}}(\lambda )$ [1/cm] as well as r [18]. A steady-state diffuse reflectance signal $Rd(r,\lambda )$ [W], therefore, is the sum of the all escaped photons at a given surface location per unit time. That is,

where x is the number of diffusely reflected photon packages remerged at r and it is a probability determined by ${\mu }_s{}^{\text{'}}(\lambda )$ and r, T is the time to deliver the entire incident photon packages, and $k(r,\lambda )$ is a constant calculated from $k(r,\lambda )=-\ln ({\displaystyle {\sum }_{i=1}^x{e}^{-{\mu }_a(\lambda )\cdot l(r,\lambda )}})/{\mu }_a(\lambda )$. Based on this theory, a natural logarithmic operation could be applied to $Rd(r,\lambda )$ to obtain the absorption characteristics (i.e., ${\mu }_a(\lambda )$) and the scattering characteristics (i.e., $k(r,\lambda )$). Assuming Hb is the only dominant chromophore of biological tissues between 450 nm and 600 nm and the spectral profile of $k(r,\lambda )$ is relatively monotonic, the spectral profile of $-\ln [Rd(r,\lambda )]$ between 450 nm and 600 nm should strongly resemble that of ${\epsilon }^{SatO2}(\lambda )$within the same spectral region.

To effectively detect the spectral profile alterations in $-\ln [Rd(r,\lambda )]$ induced by SatO₂ variations, DWT was again used on $-\ln [Rd(r,\lambda )]$. The transformed $-\ln [Rd(r,\lambda )]$ is denoted by $W_s^n\{-\ln [Rd(r,u)]\}$. The locations of the spectral features in $-\ln [Rd(r,\lambda )]$ are tracked by the zero-crossing points in $W_s^n\{-\ln [Rd(r,u)]\}$, which are denoted by $Z_{j,s,n}^{Rd}$. In the extraction algorithm, it is assumed that $Z_{j,s,n}^{\epsilon }$ acquired from the extinction spectra should be coincided with $Z_{j,s,n}^{Rd}$. Therefore, $S_{j,s,n}$ in Eq. (2) can be used as the lookup functions to estimate the SatO₂ level in $Rd(r,\lambda )$ using $Z_{j,s,n}^{Rd}$ as the inputs. That is,

3. Algorithm validation and evaluation using Monte Carlo simulation

To validate the proposed algorithm, an extensive diffuse reflectance spectra database was created using a MC simulation model of photon migration. In the simulation model, a semi-infinite medium with homogenous optical properties was used to mimic a biological tissue. In each simulation, a total of three million photons were injected into the medium at the origin (i.e., (r,z) = (0,0)). The energy of reflected photons was recorded as a function of the exit position r in order to create the diffuse reflectance spectra $R{d}_{sim}(r,\lambda )$.

Within the spectral region of 400 nm to 650 nm, Hb is the dominant chromophore in biological tissues and processes unique spectral profile characteristics. In contrast, other biological chromophores do not have distinct spectral features in the same wavelength region. In the MC simulation model, the only absorber used was Hb. The absorption property of the simulated medium hence was defined by the following equation:

where BVF is the blood volume fraction, $\overline{Sat{O}_2}$ is the user-defined SatO₂ level of Hb [%], ${\mu }_a^{oxy}(\lambda )$ is the absorption coefficient of oxygenated blood [1/cm], and ${\mu }_a^{deoxy}(\lambda )$is the absorption coefficient of deoxygenated blood [1/cm]. ${\mu }_a^{oxy}(\lambda )$and ${\mu }_a^{deoxy}(\lambda )$ were obtained from ${\epsilon }^{oxy}(\lambda )$and ${\epsilon }^{deoxy}(\lambda )$, respectively, with the whole blood Hb concentration of 150 gram/liter. The reduced scattering coefficients in biological tissue were approximated by the following equation [11,19]. That is,

where $w=\lambda \times {10}^3$. The ranges of all variables associated with the optical properties are summarized in Table 1 . Note that the ranges of the optical properties used to construct $R{d}_{sim}(r,\lambda )$ encompass the majority of biological tissue types and hemodynamic conditions. Furthermore, $R{d}_{sim}(r,\lambda )$ was recorded from r = 0.0075 cm to r = 0.9975 cm, which covers the source-detection separations commonly used in the in vivo diffuse reflectance spectroscopy studies.

Table 1. Variables used in constructing $R{d}_{sim}(r,\lambda )$ and their ranges used

View Table | View all tables in this article

To retrieve the SatO₂ information from$R{d}_{sim}(r,\lambda )$, the algorithm described in the previous section was applied to$R{d}_{sim}(r,\lambda )$. The performance of the algorithm was evaluated by calculating the absolute error between the actual and the estimated SatO₂ levels. That is,

where $\overline{Sat{O}_2}$ is the actual SatO₂ level of the Hb used in the MC simulations, as defined in Eq. (6), and $Sat{O}_2{}_{(j,s,n)}^{est}$ is the estimated SatO₂ level of Hb derived from the extraction algorithm using Eq. (5). $Sat{O}_2{}_{(j,s,n)}^{error}$ was calculated with extended combinations of (j,s,n). The $S_{j,s,n}$, as defined in Eq. (5), producing the lowest mean $Sat{O}_2{}_{(j,s,n)}^{error}$ over the entire range of the simulation parameters was identified as the ideal lookup function to retrieve the SatO₂ level of Hb from a diffuse reflectance spectrum using the proposed algorithm.

4. Algorithm validation and evaluation using tissue phantoms and in vivo experiments

To demonstrate its clinical utility, the proposed Hb oxygenation extraction algorithm was applied to the diffuse reflectance spectra acquired from tissue phantoms and in vivo tissue, using a fiber-optic spectroscopy system. The system utilized a tungsten-halogen lamp (LS-1, Ocean Optics, Dunedin, FL) as the excitation source and a spectrometer (USB 2000, Ocean Optics, Dunedin, FL) to record diffuse reflectance spectra. A fiber-optic probe was employed for excitation light delivery and reflected light collection. All the experimentally-acquired diffuse reflectance spectra were calibrated against a standard (FGS-20-02c, Avian Technologies, NH), so as to remove the spectral profile alterations induced by the spectroscopy system, and then re-sampled to reduce the spectral data interval to 2 nm.

The tissue phantoms used in the experiments consisted of a mixture of whole blood and scatterers in a saline solution. The blood samples were heparinized after their acquisition from healthy human subjects. The scatterers added to the solution were 1μm microspheres (Polysciences Inc., Warrington, PA). Two tissue phantoms were prepared for the experiments. The BVF of whole blood samples were 0.05 in Phantom 1 and 0.03 in Phantom 2. The concentrations of the scatters were 0.0159 particle/μm³ in Phantom 1 and 0.0108 particle/μm³ in Phantom 2. According to Mei theory estimates, the scattering properties resulting from the added scatterers were 31.7 1/cm at 500nm and 29.0 1/cm at 600nm for Phantom 1; and 21.5 1/cm at 500nm and 19.7 1/cm at 600nm for Phantom 2. The SatO₂ level of Hb in the phantom was first elevated by bubbling oxygen into the phantom. Then the phantom was deoxygenated by bubbling nitrogen into it. During the de-oxygenation process, the SatO₂ level of the Hb was estimated by applying the proposed algorithm to the diffuse reflectance signals acquired from the phantom surface using a fiber-optic diffuse reflectance spectroscopy system. In the probe, two identical fibers, each having a 200μm core-diameter were used, one for excitation and the other for detection, and the source-detector separation distance of the probe was 1652μm. Meanwhile, a small amount of the tissue phantom was extracted and placed into a separate apparatus to measure its collimated transmission, T_c(λ). The same measurement also was carried out with a reference sample containing the same concentration of scatterers, but no whole blood, to yield T_c,ref(λ). From the ratio of T_c(λ) to T_c,ref(λ), the absorption properties of the phantom sample ${\mu }_{a,phantom}(\lambda )$ were determined using modified Beer’s law. Sequentially, the SatO₂ level of Hb in the phantom sample (SatO_2,ref) was derived from ${\mu }_{a,phantom}(\lambda )$ using the Hb extinction coefficient fitting routine. Please note that the fitting routine was applied to ${\mu }_{a,phantom}(\lambda )$ between 650 nm to 850 nm, in order to avoid the effect of Hb packing in red blood cells [20].

In the first in vivo experiment, the algorithm was used to extract the local Hb oxygenation information at a fingertip under normal and ischemic conditions. The fiber-optic probe used in this experiment consisted of three 600 μm core-diameter optical fibers: one was used for excitation and the other two (fiber A and fiber B) for detection. The distances between the center of the excitation fiber and those of the detection fibers were 810 μm for fiber A and 1116 μm for fiber B. During the experiment, the probe was placed in direct contact with the fingertip, using a probe holder, to maintain stability of the contact position and pressure. Two spectrometers were used to acquire the diffuse reflectance spectra from the two detection fibers separately and simultaneously. A short-term local ischemic condition of the fingertip was created by tightening a rubber band placed at the base of the finger to reduce blood flow for 400 seconds. Diffuse reflectance spectra were acquired continuously at the rate of 0.5 Hz prior to the occlusion to establish a normal reference, during the ischemia and after releasing the rubber band.

In the second in vivo experiment, the algorithm was applied to diffuse reflectance spectra measured from the in vivo brains of six pediatric patients who were undergoing a craniotomy for brain tumor resection. The clinical study was conducted in accordance with a protocol approved by the Western Institutional Review Board and the Institutional Review Board of Florida International University, with consent from the patients and their parents. Methodological details of this in vivo study have already been published [21]. The fiber-optic probe used consisted of seven 300 μm optical fibers; one was used for excitation, one remained unused, and the remaining five were utilized for detection. The average source-detector separation distance was about 460 μm. The probe was held by the surgeon and placed in direct contact with the brain tissue. For each patient study, diffuse reflectance spectra were collected at a rate of 33 Hz for 12 seconds from normal cortex (i.e., outside the resection zone) and from the tumor (later confirmed by histology).

5. Spectral signal quality evaluation

MC simulation is a statistical-based model and requires a large number of repetitions (i.e., input photons) to reduce the noises. Diffuse reflectance spectra acquired by both in vivo experiments and MC simulation inevitably contain noises. For this reason, the noise level of $-\ln [Rd(r,\lambda )]$ and its impacts on the performance of the Hb oxygenation extraction algorithm were evaluated. To accomplish this goal, the signal to noise ratios (S/N) of $-\ln [Rd(r,\lambda )]$ had to be defined first. As mentioned previously, DWT is capable of filtering. With s = 1, DWT of $-\ln [Rd(r,\lambda )]$ isolates the highest frequency component in $-\ln [Rd(r,\lambda )]$, which is treated as the noise component. Therefore, depending on the s and n parameters used in DWT, S/N of a transformed $-\ln [Rd(r,\lambda )]$ with s > 1 can be conveniently estimated using the following equation:

6. Results

According to Table 1, the database of $R{d}_{sim}(r,\lambda )$ consists of 350240 spectra. During the evaluation and the validation, however, some spectra, especially those with large source-detector separations r, were considered as unusable because they contained zero data points (i.e., no signal). The average S/N of transformed $R{d}_{sim}(r,\lambda )$ at each r was used to evaluate this phenomenon. It was found that the number of unusable spectra increases and the S/N of transformed spectra decreases when r exceeds 0.3 cm. Figure 1 demonstrates two examples of the average S/N versus number of unusable spectra at various r with different s and n combinations. In order to avoid the errors induced by poor S/N over the majority optical property ranges defined in Table 1, the following evaluation and validation were carried out for r ≤ 0.3 cm only. This limitation results in a total of 103840 spectra, 9440 spectra per SatO₂ level, and 1760 at each r in the reduced $R{d}_{sim}(r,\lambda )$ database.

 

Fig. 1 The correlation between the number of unusable spectra and S/N(s,n) at various r. The long dash line represents the number of unusable spectra in $R{d}_{sim}(r,\lambda )$. The solid curve is the mean S/N of transformed$R{d}_{sim}(r,\lambda )$, and the error bars represent the standard deviations. Panel (a) represents the condition of s = 6 and n = 2 in DWT. Panel (b) represents the condition of s = 6 and n = 3 in DWT

Download Full Size | PDF

Search of the ideal $S_{j,s,n}$ for the proposed Hb oxygenation extraction algorithm was carried out using n = 1 to 3 and s = 2 to 9 in DWT of $R{d}_{sim}(r,\lambda )$ within the spectral region of 450 nm to 600nm. The results of the search show that, among all (j,s,n) combinations, (j,s,n) = (2,6,2) produces the minimal mean $Sat{O}_{2(}{}_{2,6,2)}^{error}$ defined by Eq. (8). This indicates that $S_{2,6,2}$is the ideal lookup function for the proposed algorithm. According to Eq. (2), $S_{2,6,2}$ was constructed using the location of the second inflection point identified by $Z_{2,6,2}^{\epsilon }$. Figure 2 (a) depicts how $Z_{2,6,2}^{\epsilon }$ shifts in according with the spectral profile alteration in ${\epsilon }^{SatO2}(\lambda )$ induced by SatO₂. This relationship was fitted nicely (r² = 0.99) by a rational decay function, $S_{2,6,2}$, as shown in Fig. 3 (a) . $Z_{2,6,2}^{Rd}$ retrieved from $R{d}_{sim}(r,\lambda )$ at each SatO₂ level defined in Table 1 were also plotted in Fig. 3 (a), which were found to follow the trend of $S_{2,6,2}$. This result verifies the validity as well as the accuracy of the proposed extraction algorithm.

 

Fig. 2 Panel (a) describes the effects of the SatO₂ induced spectral profile alterations in ${\epsilon }^{SatO2}(\lambda )$ on $Z_{2,6,2}^{\epsilon }$ for SatO₂ ≥ 40%. The top and the bottom groups of curves represent the filtered ${\epsilon }^{SatO2}(\lambda )$ and the DWT of ${\epsilon }^{SatO2}(\lambda )$, respectively, at various SatO₂ levels. The dot vertical lines show that the inflection points in the filtered ${\epsilon }^{SatO2}(\lambda )$ are accurately localized by $Z_{2,6,2}^{\epsilon }$. Panel (b) describes the similar capability of $Z_{4,6,3}^{\epsilon }$ for SatO₂ < 40%

Download Full Size | PDF

 

Fig. 3 In panel (a), the solid dots show $Z_{2,6,2}^{\epsilon }$ as a function of SatO₂ (≥ 40%). The solid line is the curve fit of all $Z_{2,6,2}^{\epsilon }$ using a rational decay function, which becomes the lookup function $S_{2,6,2}$. The empty squares represent the medians of $Z_{2,6,2}^{Rd}$ and the error bars represent the corresponding 10th and 90th percentiles. In panel (b), the solid dots are $Z_{4,6,3}^{\epsilon }$ at various levels of SatO₂ that are less than 40%. The solid line is the curve fit of all $Z_{4,6,3}^{Rd}$ using a 5th order polynomial function, which becomes the lookup function $S_{4,6,3}$. The empty squares represent the medians of $Z_{4,6,3}^{Rd}$ and the error bars represent the corresponding 10th and 90th percentiles. These numbers are calculated from the entire reduced $R{d}_{sim}(r,\lambda )$ database.

Download Full Size | PDF

It was also noted that in Fig. 2 (a) $Z_{2,6,2}^{\epsilon }$ vanishes when the SatO₂ level drops below 40%. Therefore a second lookup function is required for low SatO₂ (i.e., < 40%) in the proposed SatO₂ extraction algorithm. The results of the feature search and validation suggested that $S_{4,6,3}$ is the ideal lookup function to track low SatO₂ levels because it produced the minimal mean $Sat{O}_2{}_{(4,6,3)}^{error}$ among all combinations evaluated. Figure 2 (b) depicts the correlation between $Z_{4,6,3}^{\epsilon }$ and the SatO₂ induced spectral profile alterations in ${\epsilon }^{SatO2}(\lambda )$. Similarly, this relationship was fitted nicely (r² = 0.99) using a 5th degree polynomial function, $S_{4,6,3}$, as shown in Fig. 3 (b). $Z_{4,6,3}^{Rd}$ extracted from $R{d}_{sim}(r,\lambda )$ at each SatO₂ level were also plotted in Fig. 3 (b) which shows $S_{4,6,3}$ is a good lookup function for $Z_{4,6,3}^{Rd}$.

Upon identifying the two ideal lookup functions, their corresponding estimation errors, namely $Sat{O}_{2(}{}_{2,6,2)}^{error}$ and $Sat{O}_2{}_{(4,6,3)}^{error}$, were used to evaluate the performance of the proposed Hb oxygenation extraction algorithm. The overall estimation error of the algorithm is 1.23% ± 1.35% (mean ± standard deviation) for the entire reduced $R{d}_{sim}(r,\lambda )$ database.

Influences of individual variable defined in Table 1 on the performance of the algorithm also were individually evaluated. The analysis was carried out by calculating the estimation errors at all levels of one given variable, while the remaining variables were adjusted over the entire ranges.

The estimation errors at various SatO₂ levels are shown in Fig. 4 (a) . The results indicate that the estimation errors are the highest at the SatO₂ level of 30%. This level is the transition point at which one lookup function is replaced by the other one to interpret the SatO₂ level. The accuracy of the algorithm was also evaluated as a function of BVF, and the results are shown in Fig. 4 (b). In general, the SatO₂ estimation errors do not change significantly in accordance with the variations in BVF. The estimation error increases as the BVF level decreases; the highest estimation error (1.7%) is found at BVF = 0.01.

 

Fig. 4 (a) SatO₂ –dependence and (b) BVF-dependence of the SatO₂ estimation errors. The solid bars represent the average estimation errors; the error bars represent the one standard deviation.

Download Full Size | PDF

The SatO₂ estimation errors were evaluated in according with the parameters A and B, which are used to define the reduced scattering coefficient${\mu }_s^{\text{'}}(\lambda )$ in Eq. (7). The outcomes, as shown in Fig. 5 , again suggest that the estimation errors are constant for all ${\mu }_s^{\text{'}}(\lambda )$ levels.

 

Fig. 5 ${\mu }_s^{\text{'}}(\lambda )$-dependence of the SatO₂ estimation errors. The solid bars represent the average estimation errors; the error bars represent the one standard deviation.

Download Full Size | PDF

As mentioned previously, the decision of choosing the ideal lookup functions and the evaluation of the algorithm were determined using the reduced $R{d}_{sim}(r,\lambda )$ database (r < 0.3 cm). When applying this algorithm to the entire $R{d}_{sim}(r,\lambda )$ database (i.e., r = 0.0075 cm to r = 0.9975 cm), the overall estimation error of the algorithm increases to 2.39% ± 3.64%. Figure 6 illustrates the estimation errors as a function of r, derived from the entire $R{d}_{sim}(r,\lambda )$ database. As mentioned before, the number of usable spectra drops when r exceeds 0.3 cm. The estimation errors increase when r approaches the two extremes and the lowest estimation error was found between r = 0.06 cm to r = 0.15 cm.

 

Fig. 6 r-dependence of the SatO₂ estimation errors. The solid line represents the average estimation errors; the error bars represent the one standard deviation. The long dash line represents the number of usable spectra in $R{d}_{sim}(r,\lambda )$.

Download Full Size | PDF

Through the observations described above, it was found that the average estimation errors tend to increase when the mean S/N declines. To illustrate this negative correlation, the estimation errors were plotted against S/N(2,6,2) for SatO₂ ≥ 40% and S/N(4,6,3) for SatO₂ < 40% respectively in Fig. 7 . Figure 7 also reveals that the average estimation errors rise above the overall estimation error when S/N(2,6,2) for SatO₂ ≥ 40% and S/N(4,6,3) for SatO₂ < 40% descend below 700 and 70, respectively.

 

Fig. 7 The relationship between the estimation errors versus S/N(2,6,2) for SatO₂ ≥ 40% and S/N(4,6,3) for SatO₂ < 40%. The solid bars represent the average estimation errors within a specific range of S/N; the error bars represent the one standard deviation. The solid line represents the overall estimation error (2.39%).

Download Full Size | PDF

The results of tissue phantom experiments are plotted in Fig. 8(a) , which depicts very good agreement between the SatO₂ levels estimated by the proposed algorithm and the SatO_2,ref (reference) levels for both tissue phantoms employed. The Pearson correlation coefficient between the two variables was 0.99, with p < 0.05. Using SatO_2,ref as the error-free standard, the average estimation error of the algorithm was calculated as 2.4% ± 2.4% (n = 20) for Phantom 1 and 2.3% ± 2.5% (n = 22) for Phantom 2.

 

Fig. 8 (a) Results of the tissue phantom study. Each symbol represents the SatO₂ level estimated by the algorithm and the SatO_2,ref level determined by collimated transmission during the de-oxygenation process of two tissue phantoms. The line indicates unity. (b) Time histories of local SatO₂ variations at the fingertip during experimentally-induced ischemia. The SatO2 levels were derived from the diffuse reflectance spectra, recorded through fiber A and fiber B, using the proposed algorithm.

Download Full Size | PDF

For the in vivo fingertip ischemia experiment, the spectral data were taken from the tip of the middle finger of a volunteer (P. C.): the results are shown in Fig. 8 (b). The proximity of the two curves indicates that the SatO₂ levels extracted using the diffuse reflectance spectra from the two detection fibers are fairly similar. For the first 400 seconds (normal reference), the average SatO₂ level of the fingertip was 68.45% ± 2.55% from fiber A and 70.49% ± 2.55% from fiber B. The SatO₂ level started to decline right after the rubber band was tightened and eventually plateaued. During the steady-state ischemic condition (i.e., 600 sec to 800 sec), the average SatO₂ levels from fiber A and fiber B were 16.77% ± 2.90% and 17.01% ± 1.87%, respectively. The mean difference between fibers A and B was reduced during this period, relative to the normal reference level. The SatO₂ level returned to normal range, as expected, once the rubber band was released.

The demographic and histological information of the six brain tumor patients who participated in the clinical study are provided in Table 2 . Figure 9 represents the estimated SatO₂ level by the algorithm for each patient in normal cortex and the tumor area within the 12 second acquisition period. Figure 9 illustrates how the SatO₂ level may fluctuate during acquisition. However, in general, normal cortex exhibited relatively higher SatO₂ levels than tumor areas within the same patient. The average SatO₂ level of the six patients was 66.42% ± 8.86% within normal cortex and 50.89% ± 13.09% in brain tumor. The S/N of all the experimental spectra were evaluated using Eq. (9) and found above the threshold depicted in Fig. 7.

Table 2. Demographic data on the 6 patients who participated in the clinical study

View Table | View all tables in this article

 

Fig. 9 Representative time histories of the estimated SatO₂ levels of normal cortex (black) and tumor (grey) from all six patients.

Download Full Size | PDF

7. Discussion and conclusion

MC simulation for photon migration has become a standard for predicting photon migration in biological tissues [22]. One critical advantage of the MC model is that it can be used to evaluate a broad range of tissue optical characteristics efficiently. While constructing the $R{d}_{sim}(r,\lambda )$ database for the algorithm validation, the ranges of the optical properties selected encompass the majority of in vivo biomedical tissues. The outcome of validation shows that the algorithm provides satisfactory accuracy for r is less than 0.3 cm. Furthermore, the estimation performance of the algorithm is fairly constant over the entire optical property range, as indicated by low and even estimation errors. This demonstrates the applicability of the algorithm to extract the SatO₂ level of Hb from diffuse reflectance spectra measured for various tissue types in vivo.

The tissue phantom experiment was performed to further ascertain the validity of the algorithm. The algorithm exhibits adequate estimation accuracy versus SatO_2,ref in both tissue phantoms with different BVF and scattering properties. This finding further supports the validation and evaluation outcomes using simulated spectra generated by the MC simulation. As mentioned previously, a whole blood sample was used in the tissue phantom study and a packing effect of Hb in red blood cells (i.e., inhomogeneity in Hb distribution) was evident. According to Finlay et al. [20], the packing effect has a relatively small influence on the visible absorption peaks and can be ignored in the near-infrared region. SatO_2,ref was calculated from ${\mu }_{a,phantom}(\lambda )$ in the wavelength region of 650 to 850nm. The phantom study demonstrated that the SatO₂ levels of Hb estimated by our algorithm closely match SatO_2,ref levels. This suggests that the packing effect has little impact on our algorithm.

Through the derivation and validation processes, several advantages of the proposed algorithm were noticed. The operation of the algorithm involves conducting two transformations to the input (i.e., $Rd(r,\lambda )$) and applying the target spectral feature location (i.e., $Z_{j,s,n}^{Rd}$) to its corresponding lookup function $S_{j,s,n}$. This process is much simpler than the repeated fitting routine (i.e., inverse approach) required by many other models [11–16]. The computational resources required by the proposed algorithm are minimal, which makes it suitable for real-time applications. Furthermore, with an appropriate selection of the spectral profile features and hence the lookup functions, the algorithm is insensitive to the scattering properties and BVF of the biological media, as well as to the source-to-detector separation. In other words, the proposed Hb extraction algorithm with a set of established lookup functions can be immediately applied to the in vivo diffuse reflectance spectra acquired from various tissue types using optical probes with a wide range of source-detection geometry without adjustment.

The utilization of DWT in the proposed extraction algorithm also reduces the impact of spectral noise on the accuracy of the SatO₂ level assessment. As mentioned before, the DWT process is capable of signal filtering; unwanted noise in the signal is removed by a tunable linear band-pass filter. In addition, the noise level in signals can be assessed by DWT with a small scale, which provides a convenient way to access the quality of the signals (i.e., S/N) instantaneously and further ensure the accuracy of the estimation. In terms of S/N, we believe the S/N level of the signal, as defined by (9), should be higher than 700 and 70 for $S_{2,6,2}$ and $S_{4,6,3}$, respectively, in order to achieve acceptable algorithm accuracy. This conclusion is derived based upon the S/N-estimation error analysis shown in Fig. 7. It is found that the information regarding how signal quality affects the performance of other reported models is insufficient.

All other models citied in the Introduction section access SatO₂ level of Hb by relying on identifying the absorption properties of the tissue, which needs to take the excitation intensity of the light source into consideration. To adapt the experimental diffuse reflectance to different models, various intensity calibration procedures need to be employed for each method. Some methods calculate the attenuation/optical density spectrum by using a reference spectrum acquired from an intensity-dependent reference sample [9,10,13]. Other methods calibrate their diffuse reflectance spectrum intensities directly to the measurements acquired from one or multiple phantoms with known optical properties [12,15,16]. And still other models avoid this calibration procedure by introducing an intensity coupling/proportionality factor into their fitting routines [11,14]. The algorithm that we propose uses the spectral profile features in a diffuse reflectance spectrum to determine its corresponding SatO₂ level of Hb. Hence, knowing the absolute intensity of the light source is not needed. In other words, only simple calibration to the diffuse reflectance spectroscopy system is required to remove all possible spectral profile alterations induced by the instrumentation itself.

To ensure the accuracy of the extraction algorithm, the spectral profile characteristics of the input signals (i.e., $Rd(r,\lambda )$) must be well preserved. In other words, the success of applying the algorithm hinges upon the spectral resolution. During the development and evaluation processes, a 2 nm spectral resolution was used in both ${\epsilon }^{SatO2}(\lambda )$ and $R{d}_{sim}(r,\lambda )$. This spectral resolution is easily achievable with the spectrometers in today’s market. Further investigation is required to elucidate the impact of lower spectral resolution on the accuracy of the extraction algorithm.

The clinical relevance of the SatO₂ extraction algorithm also was confirmed through a couple of in vivo experiments. In the fingertip ischemia experiment, the SatO₂ level of the fingertip under normal conditions is about 70%. The algorithm immediately reported regional tissue hypoxia once the ischemia was induced and the recovery of SatO₂ level once the occlusion was cleared. These results are similar to the studies published previously [3,13]. As shown in Fig. 8, the time histories of the SatO₂ variations derived from the two source-detector separations are very similar, which demonstrates the insensitivity of the algorithm to the source-detector separation. The further reduction in the discrepancy between the two SatO₂ trends in Fig. 8 during the steady-state ischemic period may be attributed to the disappearance of the normal oxygen gradient across the epidermis under ischemic conditions [23].

It was noticed that, in some cases, the SatO₂ level of Hb in cortex fluctuated significantly during the recording period (see Fig. 9). Since the optical probe was held by surgeon's hand, we suspect that the hand movement is the primary cause of the SatO₂ fluctuation. In a previous study, excessive movement along the axial direction was found to sometimes induce large variations in probe contact pressure, thereby altering local hemodynamics [24]. The reported SatO₂ level in normal cortex exhibits certain degree of variability. This phenomenon was also observed in a large scale study which shows a wide distribution of cortical oxygenation [25]. Nevertheless, the SatO₂ levels detected from the brain tumor areas were consistently lower than those from the normal cortex, which verifies a well known physiological characteristic of brain tumors [26]. This trend may be used as an additional indicator for in vivo brain tumor demarcation using optical spectroscopy.

The extraction algorithm operates within the spectral region, 450 nm to 600 nm, where the Hb absorption spectrum is rich in spectral profile features. Unfortunately, this is also the region where the Hb absorption is very strong. This characteristic impedes the photons from penetrating deeply into tissue. According to the results generated using the MC simulations, the penetration depth of light within this region usually falls within one or two millimeters, depending on the tissue type and the Hb concentration in the tissue. In addition, the high absorption characteristics of Hb in this spectral region also limits the source-detector separation to preserve the S/N of the signals. These two factors, together, limit the volume of tissue investigated by the diffuse reflectance spectroscopy system employing the extraction algorithm for SatO₂ assessment. Under this circumstance, the suitable instrumentation partners for the algorithm would be systems with fiber-optic based setups that acquire diffuse reflectance from tissues using surface probes, catheters or needles. Such a combination will be ideal for studies of in vivo microcirculation for which focal information is required.