2.1.

Experimental Setup

The schematic of the experimental setup is shown in Fig. 1. The light source is a halogen lamp (Ocean Optics LS-1) that is stabilized (variations less than 1 after 40 min of warming up). Other light sources can be selected by an optical switch. The optical switch has one output connected to the illuminating fiber of the optical probe and another one to a reference channel that is connected directly to the spectrometer. Temporal variations of the source spectrum are corrected by the measurement made on the reference channel.

Figure 1

General scheme of the acquisition setup.

The optical probe is made of eleven optical fibers, one for illumination and ten for detection [see Fig. 2(e)]. The fibers usedare multimode silica fibers (FVP-200 PF from Thorlabs Inc.) with a numerical aperture of 0.22 and a core diameter of 200 μm. They are coated with a 20-μm polyimide sheath. The ten detection fibers are placed at various distances (noted as ρ_i, i=1–10) from the illumination fiber, ranging approximately from 0.3 to 1.35 mm with a step of approximately 0.1 mm. They are held in place in the probe by two laser micromachined alignment pieces. This procedure leads to a high accuracy in fiber positioning (approximately 15 μm), which allows the construction of several identical probes. All the components have been chosen to with stand routine disinfections and/or sterilization without damage.

Figure 2

Diagram of the algorithm for the automatic spectral determination of μ_a, μ_s ^′, and γ. (a) First, a set of reflectance curves simulated by the Monte Carlo method is produced for a wide range of discrete values for each optical coefficient (i.e., μ_a, μ_s ^′, and γ). (b) Second, the simulated curves are interpolated by cubic B-splines to provide a reflectance curve (c) for any value of μ_a, μ_s ^′, and γ. Finally, for each wavelength of the spectrum, a fit (d) of the measured reflectance curve (e) is performed onto the interpolated set of simulated curves, allowing for the determination of μ_a, μ_s ^′, and γ. The fitting algorithm used is Levenberg-Marquardt.

The ten detecting fibers are connected to an imaging spectrometer (Jobin-Yvon CP 200) with a Peltier-cooled CCD camera (Hamamatsu C7041 and S7031-1008). A mechanical shutter sets the integrating time of the CCD camera (usually about 600 ms). The image acquired on the camera is digitized by a 12-bit analog/digital card (National Instrument PCI-MIO-16E-4). A single computer is used to control the instrumentation and to process the data. A typical measurement takes 4 s, including a first acquisition while the tissue is illuminated by the source fiber and a second acquisition made with the source fiber switched off. This allows the subtraction of the dark current noise and the background light.

In order to obtain absolute reflectance spectra, a two-step calibration procedure is performed. First, the effect of the spectral responses of the source, fibers, grating, and detector is corrected by performing a measurement on a spectrally flat reflectance standard. Typically, the probe is inserted inside an integrating sphere (made of Spectralon SRS-99, Labsphere Inc.). Second, the effective source intensity is obtained by performing a series of measurements on a solid turbid siloxane phantom whose optical properties (μ_a, μ_s ^′, and γ) were previously determined at 675 nm by combining frequency-domain and spatially resolved measurements.⁷ The comparison of the phantom measurement with the expected absolute reflectance from a Monte Carlo simulation allows the calibration of the absolute intensity of the experimental reflectance. Note that this calibration factor is wavelength independent, since all the spectral features of the instrumentation have been corrected in the first step. This calibration routine is performed every time the probe is reconnected to the setup, before or after the measurement session. Repeated measurements on the calibration phantom over a six-month period showed a standard deviation of approximately 4, which was mostly due to the variability of the coupling between the probe and the phantom. Once the absolute reflectance spectra are calculated, the optical coefficients μ_a(λ), μ_s ^′(λ) and γ(λ) are determined by the algorithm presented in Sec. 2.4.

2.2.

Optical Model at Short Source–Detector Distances: The γ Parameter

Consider a tissue illuminated by a fiber with a unit of incident power. The spatially resolved reflectance R(ρ) is defined as the power of the backscattered light per unit of area detected by a second fiber at the surface of the tissue at a distance ρ from the source. R(ρ) depends on the optical properties of the sample, i.e., the absorption coefficient μ_a, the scattering coefficient μ_s, and the phase function p(cos θ) (where θ is the scattering angle), the refractive index n, and the numerical aperture NA of the detecting fibers.

The phase function can be expanded into a series of Legendre polynomials P_n(cos θ):

The diffusion approximation generally holds for μ_s ^′ρ>5 and μ_s ^′>100μ_a. For shorter source–detector separations, typically when 0.5<μ_s ^′ρ<5, we showed¹⁷ that the second moment of the phase function g₂ must also be taken into account. This case corresponds to the measurement performed with the probe described in Sec. 2.1.

Using the similarity relations introduced by Wyman et al.;,²¹ it is easy to show that besides the refractive index, only three parameters are needed to accurately describe the light propagation at such short source–detection distances¹⁷: μ_a, μ_s ^′, and γ where

2.3.

Monte Carlo Simulation

We used Monte Carlo simulations as our basis for the model of light propagation. As already mentioned, standard diffusion theory is inadequate for the small source–detector separation we used in experiments¹⁷ The Monte Carlo simulations allow us to generate reflectance curves noted as R _sim (ρ_i,μ_a,μ_s ^′,γ) [see Fig. 2(a)]. The code we used has been extensively tested and allows the use of any phase function given in discretized form,^{23 24} We used a modified form of the Henyey-Greenstein²⁵ phase function p _HG (θ,g _HG ), noted as p _MHG (θ,g _HG ,ξ) (see Ref. 17):

In order to decrease the Monte Carlo noise and limit the number of photons (one million photons have been simulated for each reflectance curve), the numerical aperture was set at a higher value than in the experiments (0.279 instead of 0.157). This use of the higher numerical aperture was chosen as a tradeoff between simulation time, noise in the simulation, and accuracy degradation. Several tests showed that after renormalization of the intensity, this mismatch introduced a typical error smaller than 2 and a maximum error of 4. Note that in both the simulation and the experiments, the numerical aperture corresponds to a narrow solid angle (a half angle of 16° instead of 9°) and the angular distribution is close to a Lambertian distribution for such small-angle ranges.

Approximately forty discrete values were chosen for μ_a and μ_s ^′, ranging respectively from 0.003 to 10 mm⁻¹ and from 0.5 to 10 mm⁻¹; twenty values were chosen for γ, ranging from 1.0 to 2.9. This resulted in a four-dimensional matrix of simulated reflectance curves, noted as R _sim (ρ_i,μ_a,μ_s ^′,γ), each coefficient taking only discrete values. Using these parameters, 32,000 reflectance curves have been computed, each for a specific triplet of optical coefficients (μ_a,μ_s ^′,γ). Since the path lengths and exit positions of each simulated photon were stored, a full Monte Carlo simulation was required only for each different value of γ (thus twenty simulations only). Assuming the tissue is homogeneous, scaling relationships¹⁷ allow us to derive any reflectance curve R _sim (ρ_i,μ_a,μ_s ^′,γ) from a single simulation R _sim (ρ_i,μ_a=0,μ_s ^′=1,γ), γ being kept constant.

2.4.

Algorithm for the Determination of μ_a, μ_s ^′, and γ

We denote ρ_i, i=1–10, the distances between the ten detecting fibers and the illuminating fiber, and R _mes (ρ_i) as the measured reflectance at these distances for a given wavelength. The inverse problem corresponds to determining the optical coefficient triplet (μ_a,μ_s ^′,γ) that is the closest to the reflectance curve R _mes (ρ_i).

The basic idea of the algorithm, as shown in Fig. 2, is to fit the measured reflectance curve R _mes (ρ_i) onto a set of previously simulated reflectance curves, forming a matrix noted as R _sim (ρ_i,μ_a,μ_s ^′,γ). An algorithm based on a lookup table can be performed as a first step. However, in such a case, the number of simulations performed naturally limits the accuracy of the result. Since the reflectance curves are continuous in all their parameters, this problem could be overcome by interpolating the basic set of simulations. The result of the interpolation is thus a function R¯ _sim (ρ_i,μ_a,μ_s ^′,γ), which is a numerical function defined for any value of the optical coefficient triplets (μ_a,μ_s ^′,γ) within their definition ranges and for each discrete distance ρ_i, i=1–10.

The matrix R _sim (ρ_i,μ_a,μ_s ^′,γ) was interpolated with cubic B-splines.²⁶ Scho¨nberg et al.;²⁷ have shown that any discrete function R(x), (x∈ℤ) may be interpolated by a unique combination c(k) of B-splines βⁿ(x):

Cubic B-splines β³(x) are defined like pieces of polynomials of degree 3, smoothly connected:

If the function R(x) is known on x values that are not integers—as is the case for the reflectance matrix R _sim (ρ_i,μ_a,μ_s ^′,γ) —it must be made uniform. This step is performed by interpolating the function with nonuniform splines, resampling it on the integers, and finally interpolating it again with uniform B-splines. This process, fully described in Ref. 28, has been applied on R _sim (ρ_i,μ_a,μ_s ^′,γ).

In terms of performance, this interpolation method is very efficient because the initial matrix R _sim (ρ_i,μ_a,μ_s ^′,γ) is economically replaced by a matrix of coefficients c(ρ_i,k,l,m) of the same dimension, and any value of R¯ _sim (ρ_i,μ_a,μ_s ^′,γ) is calculated by the summation of 3×4 B-splines weighted by their coefficients c(ρ_i,k,l,m). * The interpolated function R¯ _sim (ρ_i,μ_a,μ_s ^′,γ) could then be properly used to fit any measured reflectance curve in order to determine its optical coefficients. To perform the fit, the Levenberg-Marquardt algorithm was chosen because of its efficiency and its robustness. Appropriately, this algorithm minimizes the quantity noted as χ² and defined as

As required by the algorithm, the Jacobean of the R¯ _sim (ρ_i,μ_a,μ_s ^′,γ) function can be calculated economically, as mentioned before, by the combination of square B-splines. For the lowest wavelength, a first set of parameters (μ_a,μ_s ^′,γ) is determined using a lookup table based on the discrete set of simulations R _sim (ρ_i,μ_a,μ_s ^′,γ). These values are then used as initial guesses for the Levenberg-Marquardt algorithm. For the next wavelength, the initial values are taken directly from the previous wavelengths (typically from an average of the four closest wavelengths). Nevertheless, we found the converged values to be insensitive to the initial guesses, which is related to the fact that only single minima were found in the χ² function, as discussed in the next section. The determination of the optical coefficients for a whole spectrum, composed of 510 wavelengths (480 to 950 nm) takes about 50 s (on a 400-MHz Pentium personal computer).

3.1.

Algorithm Robustness to Noise

Various assessments of the algorithm are reported in this section. First, in order for the algorithm to work properly, its solution must be unique. This is guaranteed if only a single minimum is present in the χ² function. An inspection of the χ² function over the whole range of values of μ_a, μ_s ^′, and γ for which the algorithm is defined (i.e., μ_a from 0.003 to 10 mm⁻¹, μ_s ^′ from 0.5 to 10 mm⁻¹, and γ from 1.0 to 2.9) revealed only a single minimum for each set of parameters. Second, the robustness of the algorithm was studied by assessing the impact of noise on the accuracy of the optical coefficient determination. For this, noise (uniformly distributed) was added to simulated experimental reflectance curves. The relative errors between the result of the fit and the expected optical properties were estimated and are shown in Fig. 3 for noise of 2, 4, and 8 of the reflectance curve. Such noise figures are realistic for clinical tissue measurements, where tissue heterogeneity plays a significant role as noise in the reflectance curve. For each noise value, both the maximal and average errors of determination are plotted. Typically, the mean error on μ_s ^′ and γ is lower than 3, 5, and 10 for a noise of 2, 4, and 8, respectively. The maximum error on μ_a is lower than 6, 15, and 30 for a noise of 2, 4, and 8, respectively. The determination of μ_a is the most sensitive to noise because the sensitivity of the measurement to the absorption is not optimal owing to the relatively short path length at a short source–detector separation. Also, the error increases significantly at low absorption, especially for absorptions lower than 0.05 mm⁻¹. Specifically, μ_a values smaller than 0.02 mm⁻¹ cannot be reasonably estimated accurately with such a small source–detector separation. Therefore this method is best adapted to wavelengths in the visible (400 to 700 nm) or wavelengths higher than 900 nm, where μ_a is typically higher than 0.02 mm⁻¹ in tissues. Overall, these results demonstrate that μ_a, μ_s ^′, and γ can be determined simultaneously from the reflectance curve at short source–detector separations, which is a remarkable result.

Figure 3

Accuracy of the algorithm to determine μ_a, μ_s ^′, and γ from noisy reflectance curves. The relative amplitudes of uniformly distributed noise added to the simulated reflectance curves were 2, 4, and 8, and both the maximal (thin lines) and average (bold lines) inaccuracies are plotted.

3.2.

Experiments With Phantoms

In order to complete the validation of the method, we performed a series of measurements on tissuelike phantoms made of suspended and dyed microspheres. Ge´le´bart et al.;²⁹ have suggested that tissue phase functions could be advantageously simulated with a mixture of spheres with different diameters. Following this approach, we mixed polystyrene microspheres of five different diameters: 1, 0.5, 0.2, 0.1, and 0.05 μm. An identical concentration of 1.3 [g/liter] was set for each diameter, corresponding to a sphere density decreasing as d⁻³, where d is the sphere’s diameter. The scattering spectra μ_s ^′(λ) and γ(λ) were calculated by Mie theory,³⁰ assuming a refractive index of 1.60 for polystyrene and 1.33 for water.

The absorption properties of the solution were obtained by adding Nigrosin to the microsphere suspensions with a concentration of 59.09±1.65 μg/ml. The Nigrosin absorption spectrum μ_a(λ) was determined by spectrophotometry. Figure 4 shows a comparison of the predicted spectra of optical coefficients μ_a(λ), μ_s ^′(λ), and γ(λ) and the corresponding spectra measured experimentally with the probe.

Figure 4

Comparison of the predicted spectra of optical coefficients (dashed curve) and those measured (solid curve) on the phantoms. The bold line is the average and the thin solid ones are the standard deviations.

A preliminary analysis allowed us to make two adjustments in the experimental procedure that eliminated two sources of systematic errors. First, we found that the use of the two closest fibers (ρ≈0.3) did not fit well with the simulations and was causing a systematic error. This effect was most likely due to a mismatch in higher moments of the phase function between our theoretical and experimental models. Our theoretical model assumes that only the first two moments of the phase function affect the reflectance. We showed in Ref. 17 that such an assumption is valid for an optical distance ρμ_s ^′>0.5. Shorter distances can be used only if higher moments are matched between the simulation and experiment. This is indeed generally not the case between the average phase function of microsphere mixtures made of only five sizes of spheres and the modified Henyey-Greenstein phase function used in the simulation. Thus, these two fibers were removed from the analysis of the phantom measurements. However, the first two fibers have been kept when analyzing tissue measurements, since in this case the high-order moments of the tissue phase function are expected to match significantly better with the ones of the modified Henyey-Greenstein phase function. Moreover, adding these two fibers improved the robustness of the fit. Second, a systematic difference of 5 was found in the reflectance intensity between experiments and simulations, which was mainly due to the refractive index difference between the phantom (n=1.33) and the refractive index assumed in the simulations (n=1.40, corresponding to tissue).

Figure 4 shows the comparison between the predicted and experimental μ_a(λ), μ_s ^′(λ), and γ(λ) values. The agreement on μ_s ^′(λ) and γ(λ) is excellent, with an error typically less than 5. A good agreement is also found for μ_a(λ) with an error less than 10. These results confirmed the algorithm test described in the previous section.

3.3.

In Vivo Measurements of Human Gastrointestinal Tracts

In this section we present measurements performed in vivo on the human stomach, in the antrum (lower part of the stomach) and fundus (upper part) during gastroscopic examination. These results represent an example of the clinical feasibility of the method. Other studies are currently being performed on the colon and the uterus. One of the medical interests of the method lies in the possibility of instantly discriminating gastritis from normal cases in vivo during gastroscopy. This necessitates the determination of normal values for the mucosa in the antrum and in the fundus, which are presented here. In a companion paper,³¹ we will report further studies on the stomach, showing that this technique is sensitive to different stages of gastritis and to subtypes as well.³²

A set of 35 patients, 21 females and 14 males, aged from 23 to 87 with an average age of 50 were examined. For each patient, four sites were usually selected to be measured, two in the antrum and two in the fundus. The optical probe was inserted into the working channel of the gastroscope at the end of the examination and gently applied on the tissue. The pressure could not be monitored during acquisition, but other controlled measurements on mouse skin showed that pressure variations did not significantly alter the results.

The endoscope illumination was turned off after the probe placement because of its high intensity. During this time, four to six measurements were taken in approximately 20 s. At the end, the endoscope illumination was turned on again and measurements were accepted and averaged if the probe stayed at the same location. The accept rate was approximately 80. Immediately after that, a biopsy was taken on the measurement site, which could be precisely located by a light mark on the epithelium. The data were analyzed after the clinical session.

Figure 5 shows the average optical properties, μ_a(λ), μ_s ^′(λ), and γ(λ), of twenty-two measurements performed in the antrum and sixty measurements performed in the fundus. All sites were found normal after biopsy. First, as expected, μ_a(λ) clearly shows the presence of hemoglobin, characterized by absorption peaks at 541 and 576 nm. The absorption decreases significantly in the near-infrared range. The slight increase for wavelengths larger than 850 nm is due to water and hemoglobin. The concentration of total oxyhemoglobin (HbO ₂) and deoxyhemoglobin (Hb) were obtained by a least-square fit of μ_a(λ). The total hemoglobin concentration (i.e., [HbO ₂]+[Hb]) was 20.0±2.6 μM in the antrum, with a saturation (sO ₂) around 48. In the fundus, a similar total hemoglobin concentration was found, 21.0±2.6 μM, whereas the saturation drops to 21. This surprising very low saturation could nevertheless be explained by a high concentration of CO ₂ in the mucosa. Further investigations are needed to validate such results, since no other data on the oxygenation and saturation of hemoglobin in the fundus mucosa have been published yet. The μ_s ^′(λ) spectra measured in both the fundus and antrum are monotonically decreasing, corresponding approximately to a power function, as expected for tissue.^{9 29 33} The scattering spectrum μ_s ^′(λ) is found to be significantly lower in the antrum than in the fundus over nearly the entire spectrum.

Figure 5

Average (bold lines) and standard deviations (thin lines) of optical coefficients measured for normal antrum (solid curve) and normal fundus (dashed curve).

The γ(λ) spectra are approximately constant. This is consistent in the case of a large distribution of scatterer sizes, as found in tissue. As a matter of fact, similar findings were displayed in the experiment and simulation using the microsphere suspensions. Computations with Mie scatterers showed that the γ(λ) value is related to the scatterer size distribution.^{22 32} Therefore the absolute γ value might be important in differentiating tissue types.

For the antrum and the fundus, the γ values are 1.9±0.1 and 2.0±0.1, respectively. Assuming a Henyey-Greenstein phase function, these values would correspond to a g of 0.9 and higher. Such values are basically in agreement with g values commonly accepted for tissue. However, the γ values are found in some cases to be larger than 2, which could not be produced by a simple Henyey-Greenstein phase function. This seems to confirm previous results showing that a Henyey-Greenstein function alone does not approximate the tissue phase function adequately at large angles, but needs to be completed by either a second Henyey-Greenstein function with a negative g value or a cos ² θ function.^{17 20}

The differences between the fundus and antrum found in the scattering properties μ_s ^′(λ) and γ(λ) can be correlated with the histological analysis of the biopsy samples. Figure 6 shows typical histological images of tissue resected after measurements. These images clearly show a difference in the thickness of the foveolar and glandular regions between the antrum and the fundus. These differences are likely to be the source of the differences observed in the μ_s ^′(λ) and γ(λ) coefficients. We recently showed³⁴ correlations between spectra of optical coefficients measured with this setup and histological analysis of a mouse skin alteration model.

Figure 6

Histological slices of antrum (a) and fundus (b). Whereas the layered structure is clearly visible in the fundus with the glandular region (G) clearly distinct from the foveolar (F), the structure of the epithelium is thinner and less contrasted in the antrum. Moreover, the mucosa is two times thicker in the fundus than in the antrum.

Globally, these results demonstrate that the scattering, μ_s ^′(λ) and γ(λ), and absorption, μ_a(λ), properties of the tissue can be separated quantitatively. Only a fairly small coupling effect with the μ_a(λ) spectrum can be observed in both μ_s ^′(λ) and γ(λ) spectra. Such a coupling is unavoidable considering that actual tissue is not homogeneous, which is the opposite of the theoretical model. As already mentioned, the layer structure is obvious in Fig. 6.

Two approaches are possible to further diminish such a coupling and improve accuracy. First, the μ_a(λ), μ_s ^′(λ), and γ(λ) spectral shapes could be constrained. The μ_a(λ) spectrum can be assumed to result from known chromophores such as Hb, HbO ₂, water and fat. As discussed earlier, in tissue, μ_s ^′(λ) approximately follows a power law and γ(λ) is approximately a linear function. Incorporating such conditions in the inverse problem would significantly decrease coupling effects. Second, a layer model could be used to more closely model the tissue structure. Spatially resolved measurements seem appropriate to determine some of the properties of a two-layer geometry.¹⁹