Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Derivation of chromophore concentrations in turbid media using the rate of change of optical attenuation with respect to wavelength

Open Access Open Access

Abstract

A spectroscopic method is investigated for estimating the concentrations of absorbers in highly scattering media using measurements of the rate of change of optical attenuation with wavelength. Such measurements are independent of absolute intensity and thus may be significantly less influenced by changes in coupling which often cause artefacts in medical applications of near-infrared spectroscopy. The method has been explored using a combination of stochastic (Monte Carlo) and analytical (diffusion-based) models and experiments on samples of turbid fluids. Results suggest that the method is highly tolerant of changes in the measurement geometry. The accuracy of the derived concentrations of absorbers can be strongly influenced by the wavelength dependence of scattering, and an ad-hoc, empirically-derived correction for this dependency has been investigated and implemented with some success.

Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

1. Introduction

Near-infrared spectroscopy (NIRS) is a medical diagnostic method which usually involves estimating changes in the concentrations of oxy- and deoxy-hemoglobin in biological tissue from measurements of changes in the diffuse reflectance of near-infrared light. So-called functional NIRS is now widely used as a means of observing hemodynamic responses within the cerebral cortex to sensory stimuli or cognitive tasks [1,2]. Furthermore, measurements at multiple locations on the scalp enable images to be generated revealing the spatial distribution of blood volume and oxygenation changes in the cerebral cortex, using a technique known as diffuse optical tomography (DOT) [3].

Because scattering is the dominant interaction between near-infrared photons and most biological tissues, obtaining quantitative measurements of chromophore concentrations is challenging. A common approach to quantity tissue oxygenation, known as spatially resolved spectroscopy, involves measuring the diffuse reflectance as a function of distance from the source [4]. Assuming the wavelength dependence of scatter in the tissues is known, the wavelength dependence of absorption can be derived, which enables the absolute ratio of concentrations of oxy- and deoxy-hemoglobin to be determined. Experimental studies have suggested that although scatter varies considerably, its wavelength dependence is remarkably similar for all tissues [5].

Conversion of measurements of changes in diffuse reflectance into changes in chromophore concentrations requires an estimate of the mean optical pathlength of light between a source and detector on the tissue surface, which can be measured directly using time- or frequency-domain instruments [1]. Alternatively, Matcher et al. [6] demonstrated a simple method for measuring the pathlength utilizing the spectral absorption features of water. This involves obtaining the second derivative of the measured attenuation spectrum with respect to wavelength, which removes both the arbitrary constant offset (due to unknown coupling) and any linearly wavelength-dependent attenuation. The resulting spectrum contains only features corresponding to the second derivative of the absorption spectra of any chromophores in the tissue, including water, whose concentration is assumed for a given tissue. Observations that derivatives of spectra minimize the influence of scattering and nonspecific absorption while enhancing sensitivity to absorption features of chromophores of interest were reported earlier by Ferrari et al. [7].

The recent advent of sophisticated machine learning methods is enabling new techniques for determining optical properties of tissues, by analyzing spectra using algorithms trained on measurements on well-characterized ex vivo tissue samples as well as artificially-generated datasets generated using appropriate models [8,9].

Unfortunately, measurements of diffuse reflectance are highly sensitive to the uncertain and variable coupling of light into and out of the skin surface. This sensitivity largely inhibits the routine derivation of absolute concentrations of chromophores in vivo, and NIRS and DOT measurements are commonly contaminated by artefacts due to large and random changes in coupling due to relative motion between the tissue and the optical probe. As a consequence, much effort has been devoted to methods for identifying and/or eliminating such artefacts in recorded data [1012]. Meanwhile, some investigators have explored alternative types of optical measurements which might be less sensitive to coupling variability. Specifically, Dehghani et al. [13,14] and Pucci et al. [15] have noted the insensitivity of spectral derivatives to unknown coupling coefficients for both spectroscopy and imaging.

A new technique was recently proposed, known as wavelength modulated near-infrared spectroscopy (WM-NIRS) [16], which measures the spectral derivative directly. NIRS and DOT techniques commonly measure diffuse reflectance by modulating the source intensity at a fixed frequency (typically a few kHz) and use a lock-in amplifier (or Fourier analysis) to isolate the signal amplitude at the source frequency. WM-NIRS involves modulating the source wavelength (by a few nm, for example) at a fixed frequency. Variation in the tissue absorption over the wavelength range will result in a small modulation in detected intensity which can also be measured using lock-in or Fourier methods. It was originally hypothesized that this approach could facilitate estimates of absolute concentrations of chromophores if unknown surface coupling losses and geometrical factors are insensitive to small changes in wavelength. The potential of WM-NIRS has now been explored using a combination of stochastic and analytical models and experiments on samples of turbid fluids. In the following sections, the underlying theory is briefly presented, and the predictions of the theory are then investigated using an analytical model based on a Greens function solution to the diffusion equation and a separate Monte Carlo model of light propagation in turbid media. The results achieved with both models are then compared with experimental measurements using a spectrometer and a glass cell containing a mixture of inks and an aqueous suspension of polystyrene microspheres with precisely characterized absorbing and scattering properties.

2. Theory

The implementation of NIRS as a medical tool is founded on the modified Beer-Lambert law [17], which determines that for a single source fiber and a single detector fiber placed on the surface of the interrogated tissue, separated by a distance d, the detected intensity is given by:

$$I = k\; {I_0}{e^{ - {\mu _a}\beta d - G}}$$
where µa is the absorption coefficient of the tissue, β is the differential pathlength factor (DPF), G is an unknown scattering-dependent geometric factor, I0 is the intensity of the source, and k is the (usually unknown and highly variable) coupling efficiency. The value of µa depends on the concentrations cn of the constituent chromophores as follows:
$${\mu _a} = \mathop \sum \nolimits_n {\varepsilon _n}{c_n}$$
where n is the number of chromophores and εn are their molar extinction coefficients. Because k and G are unknown, NIRS typically relies on measuring changes in intensity in response to changes in tissue optical properties under conditions where it is assumed that k and G remain constant. Using Eq. (1) the attenuation A of the light between the source and detector can be expressed as follows:
$$A = \; - \ln \left( {\frac{I}{{{I_0}}}} \right) = {\mu _a}\beta d + G - \textrm{ln}(k )$$

Differentiating this equation with respect to wavelength λ yields:

$$\frac{{\partial A}}{{\partial \lambda }} = \frac{1}{{{I_0}}}\frac{{\partial {I_0}}}{{\partial \lambda }} - \frac{1}{I}\frac{{\partial I}}{{\partial \lambda }} = d\left( {\beta \frac{{\partial {\mu_a}}}{{\partial \lambda }} + \frac{{\partial \beta }}{{\partial \lambda }}{\mu_a}} \right) + \frac{{\partial G}}{{\partial \lambda }} - \frac{1}{k}\frac{{\partial k}}{{\partial \lambda }}$$

It is immediately evident that the rate of change of attenuation with respect to wavelength is a relative measure, independent of the absolute values of source intensity I0 and measured intensity I. As noted by Pucci et al. [15], assuming k is independent of wavelength, the use of the first derivative of spectra suppresses the influence of uncertain coupling. By employing this assumption Eq. (4) reduces to:

$$\frac{{\partial A}}{{\partial \lambda }} = d\beta \frac{{\partial {\mu _a}}}{{\partial \lambda }} + d\frac{{\partial \beta }}{{\partial \lambda }}{\mu _a} + \frac{{\partial G}}{{\partial \lambda }}.$$

It was previously proposed [16] that, if G can also be assumed to be independent (or only weakly dependent) on wavelength, then it should be possible to derive absolute concentrations of chromophores from measurements of ∂A/∂λ using prior knowledge of the derivatives of the extinction coefficients and the DPF. The primary goal of the modelling described in the following sections was to determine the significance of the term ∂G/∂λ compared to other terms in Eq. (5) under two scattering conditions: a) identical scattering properties at all wavelengths (unrealistic for biomedical tissues), and b) scattering whose wavelength dependence follows a simple power relationship. The results have important consequences regarding the potential of WM-NIRS to determine absolute values of chromophore concentrations.

3. Analytical model of photon migration for a turbid slab

3.1 Model geometry and optical properties

The diffusion approximation to the radiation transfer equation has been widely used to describe the transport of light through highly scattering biological tissues. The approximation typically holds well providing scattering interactions dominate over absorption events, and that measurements are made over distances much greater than the average photon pathlength between consecutive scatters. Scattering in homogenous media is normally represented by two parameters: the scatter coefficient µs (the number of scatters per unit length) and the anisotropy parameter g (the average cosine of the scattering angle). In the diffusion equation, these are combined into a single parameter known as the transport scattering coefficient (also known as the reduced scatter coefficient) µs′ = µs(1 – g). Contini et al. [18] present a comprehensive set of solutions to the diffusion equation for both a continuous and time-dependent (i.e. pulsed) source and for different measurement geometries. The continuous source solutions were used to estimate the wavelength-dependent attenuation A(λ) and differential pathlength factor β(λ) for a source and detector separated by a distance d = 20 mm on the same surface of a slab of thickness 25 mm, with other dimensions extending to infinity. The internal and external refractive indices of the slab were 1.4 and 1.0 respectively. The wavelength-dependent absorption coefficient of the slab was given by:

$${\mu _a} = {\epsilon _{Hb{O_2}}}{c_{Hb{O_2}}} + {\epsilon _{Hb}}{c_{Hb}} + {\mu _{a,Water}}W + {\mu _{a,Lipid}}L + {\mu _b},$$
where εHbO2 and cHbO2 are the molar extinction coefficient and molar concentration of oxyhemoglobin, εHb and cHb are the equivalent parameters for deoxyhemoglobin, µa,Water and µa,Lipid are the absorption coefficients of water and lipid, W and L are the fractions of water and lipid, and µb is a wavelength-independent background term. For the model employed here, the following values were chosen for the constant parameters (and which are roughly typical of what might be expected for the adult human brain [19]): cHbO = 56 µmol, cHb = 24 µmol, W = 0.8, L = 0.116, µb= 0.012 mm-1. Estimates of the four coefficients as a function of wavelength (in the range 740 nm – 860 nm) are displayed in Fig. 1 using data extracted from the online database compiled by Prof. Scott Prahl of the Oregon Institute of Technology [20]. Also shown are the first derivatives of these spectra with respect to wavelength, obtained by fitting a tenth-order polynomial to each spectrum, and then finding the gradient of the polynomial.

 figure: Fig. 1.

Fig. 1. a) Extinction coefficients of oxy- and deoxy-hemoglobin; b) Absorption coefficients of water and lipid; c) The first derivatives of the specific absorption coefficients of oxy- and deoxy-hemoglobin with respect to wavelength; d) The first derivatives of the absorption coefficients of water and lipid with respect to wavelength.

Download Full Size | PDF

Two scattering models were explored: a) µs′ = 0.5 mm-1 and b) µs′ = 0.5 (λ/800)p mm-1, where λ is expressed in nanometers. The wavelength dependence is in accordance with empirical models published elsewhere [21,22] and is normalized to yield a value of µs′ = 0.5 mm-1 at λ = 800 nm. Unless otherwise stated, a power p = 1.2 has been used.

Values of the absorption and transport scatter coefficient were estimated at 1 nm intervals in the wavelength range from λ = 740 nm to λ = 860 nm. Equation (45) of Contini et al. [18] yielded the diffuse reflectance R(λ) from which the attenuation is obtained using A(λ) = –ln [R(λ)]. Meanwhile Eq. (47) yielded the mean optical pathlength L(λ) from which the DPF is obtained using β(λ) = L(λ)/d.

The values of attenuation A(λ) and DPF β(λ) obtained for both constant and wavelength-dependent scattering are shown in Fig. 2. The geometric term G(λ) was then calculated in each case using Eq. (3) and k = 1. Finally, the first derivative of the attenuation with respect to wavelength ∂A/∂λ for both scattering models was computed using a 20th-order polynomial fit to the values of A(λ).

 figure: Fig. 2.

Fig. 2. Estimation of a) attenuation and b) differential pathlength factor for diffuse reflectance measurements on an analytical slab model with brain-like absorption and for both constant (black) and wavelength-dependent (blue) scattering.

Download Full Size | PDF

3.2 Results for constant scatter

Figure 3(a) shows a plot of the four terms in Eq. (5) for the constant scattering model.

 figure: Fig. 3.

Fig. 3. a) Estimation of the four terms in Eq. (5) for diffuse reflectance measurements on an analytical slab model with brain-like absorption and constant scatter; b) The sum of the terms in Eq. (7), and the difference between the terms in Eq. (8).

Download Full Size | PDF

The first observation to be made is that the term ∂G/∂λ is not negligible compared to the other terms, contrary to what has been proposed previously [16]. Instead, it is apparent that

$$\frac{{\partial G}}{{\partial \lambda }} + d\frac{{\partial \beta }}{{\partial \lambda }}{\mu _a} \approx 0$$
at all wavelengths, and consequently
$$\frac{{\partial A}}{{\partial \lambda }} \approx d\beta \frac{{\partial {\mu _a}}}{{\partial \lambda }}$$
at all wavelengths. This is illustrated in Fig. 3(b) which plots the sum of the terms in Eq. (7), and the difference between the terms in Eq. (8), which are approximately zero (and essentially identical as expected from Eq. (5)). So, for wavelength-independent scatter, Eq. (5) can be simplified as:
$$\frac{{\partial A}}{{\partial \lambda }} \approx d\beta \mathop \sum \nolimits_n {c_n}\frac{{\partial {\varepsilon _n}}}{{\partial \lambda }},$$
which suggests that measurements of ∂A/∂dλ and a reliable estimate of the DPF are sufficient to derive approximate values of the chromophore concentrations. This was tested by first extracting values of A(λ) and β(λ) at each wavelength, and computing the first derivatives of εHbO2, εHb, μa,Water, and μa,Lipid with respect to wavelength. Then, using only values of ∂A/∂λ and β(λ) at nine discrete wavelengths (λ = 760, 770, 780, …, 840 nm), values of cHbO2, cHb, W and L in Eq. (6) were derived by finding the minimum norm least-squares solution to Eq. (9) expressed as a linear matrix equation. The derived and true values of the parameters are shown in Table 1. Also shown is the oxygen saturation, equal to the percentage of hemoglobin in the oxygenated form. Values of the derived parameters are also shown when β(λ) is replaced with its mean value over the 760–840 nm range, <β> = 2.96.

Tables Icon

Table 1. True and estimated values of the parameters in Eq. (6), obtained by solving Eq. (9) using values of ∂A/∂λ at nine discrete wavelengths and either β(λ) or <β> derived from an analytical slab model with brain-like absorption and constant (wavelength-independent) scatter. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

All the parameters are derived within ten percent of their true values when either the wavelength dependent DPF or the mean value of the DPF is employed.

3.3 Results for wavelength-dependent scatter

Figure 4(a) shows a plot of the four terms in Eq. (5) for the wavelength-dependent scattering model, and Fig. 4(b) shows the corresponding sum of terms in Eq. (7), and the difference between the terms in Eq. (8). Comparison with the result for constant scatter (Fig. 3) reveals some curious similarities and differences. Figure 4(a) now exhibits a vertical offset between terms $\frac{{\partial A}}{{\partial \lambda }}$ and $d\beta \frac{{\partial {\mu _a}}}{{\partial \lambda }}$, and reveals that the sum of terms in Eq. (7) is no longer approximately zero. Thus Eq. (9) no longer holds and must be modified as follows:

$$\frac{{\partial A}}{{\partial \lambda }} = d\beta \mathop \sum \nolimits_n {c_n}\frac{{\partial {\varepsilon _n}}}{{\partial \lambda }} - r({\lambda ,p} )$$
where r(λ,p) represents a residual function which has already been shown to be (almost) negligible when scatter is constant (i.e. when scattering power p = 0).

 figure: Fig. 4.

Fig. 4. a) Estimation of the four terms in Eq. (5) for reflectance measurements on an analytical slab model with brain-like absorption and wavelength-dependent scatter. b) The sum of the terms in Eq. (7), and the difference between the terms in Eq. (8).

Download Full Size | PDF

Thus the apparent offset can be expressed as:

$$f(\lambda )= r({\lambda ,p} )- r({\lambda ,0} )$$

The form of offset f (λ) was investigated for different values of the scattering power, as shown in Fig. 5(a) for p = 0.8, 1.0, 1.2, 1.4, and 1.6. At each wavelength an almost linear dependence on p was observed. It was also observed that the wavelength dependence of the offset f (λ) was very similar to that of the transport scatter coefficient, and therefore the form of the quotient f(λ)/µs′(λ) was also investigated, as shown in Fig. 5(b) for the same values of p.

 figure: Fig. 5.

Fig. 5. a) The offset f(λ) derived for different values of the scatter power p. b) The offset f(λ) divided by the transport scatter coefficient µs′ for different values of the scatter power p.

Download Full Size | PDF

The apparent near-independence of this quotient on wavelength suggests that the offset can be represented by a constant multiplied by the transport scatter coefficient. Thus:

$$\frac{{\partial A}}{{\partial \lambda }} \approx d\beta \mathop \sum \nolimits_n {c_n}\frac{{\partial {\varepsilon _n}}}{{\partial \lambda }} - \alpha \mu _s^{\prime} = d\beta \mathop \sum \nolimits_n {c_n}\frac{{\partial {\varepsilon _n}}}{{\partial \lambda }} - \alpha {\lambda ^{ - p}}$$

Thus deriving the chromophore parameters from measurements of ∂A/∂λ requires fitting to two additional parameters (α and p), or one if a value of p is assumed. This was attempted using the same discrete values of ∂A/∂λ and β(λ) by expressing Eq. (12) as a matrix equation and finding the minimum norm least-squares solution. To minimize any dependency on the initial guess, fitting was repeated 100 times using initial values sampled randomly. The results presented in Table 2 represent the mean values. The third and fourth columns in Table 2 show the derived parameters when the known function β(λ) is used and when it is replaced with its mean value, respectively. The majority of the values are within ten percent of their true values, and using the average DPF does not have a major detrimental effect. Although the lipid fraction has a comparatively large percentage error when compared to the other parameters, the error in absolute terms is small. Furthermore, the contribution of lipid to the total absorption over the 760–840 nm range is only about 0.25%.

Tables Icon

Table 2. True and estimated values of the parameters in Eq. (6), obtained by solving Eq. (12) using values of ∂A/∂λ at nine discrete wavelengths and either β(λ) or <β> derived from an analytical slab model with brain-like absorption and wavelength-dependent scatter. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

 figure: Fig. 6.

Fig. 6. a) Derived values of the concentrations of oxy- and deoxy-hemoglobin for different oxygen saturations. The continuous lines represent the true values. b) Derived values of the fractions of water and lipid, and the estimated oxygenation saturation. The continuous lines represent the true values.

Download Full Size | PDF

The process was repeated for other values of oxygen saturation between 60% and 95%, and the derived parameter estimates are displayed in Fig. 6. The error bars represent the standard deviation of the derived values when the initial values employed in the fitting algorithm were varied randomly.

Figure 6 reveals a tendency for both cHbO2 and W to be overestimated at lower saturations and underestimated at higher saturations. The parameter cHb is also overestimated at lower saturations, while the derived saturation is consistently slightly underestimated. It is possible that the accuracy could be improved in future by implementing an alternative means of compensating for the wavelength dependence of scatter.

4. Monte Carlo model of photon migration for a turbid slab

4.1 Model geometry and optical properties

The analytical model presented in the previous section suggests that it may be possible to derive estimates of chromophore concentrations independently of the geometry factor. However, although the slab geometry had a finite thickness, the top and bottom planes of the slab were infinite. In order to explore the validity of the approach for finite geometries a Monte Carlo model was developed using MCX, a voxel-based open-source light transport simulator funded by the US National Institute of Health [23]. The chosen geometry is illustrated in Fig. 7.

 figure: Fig. 7.

Fig. 7. The slab geometry of a Monte Carlo model of photon migration. The source and detector are aligned along the y axis and separated by d = 20 mm, and the midpoint was varied between the center of the slab at (0,0,0) and a point near the edge at (49,0,0).

Download Full Size | PDF

The geometry and tissue-like optical properties were identical to those of the analytical slab model with the following exceptions:

  • • lateral boundaries exist at x = ± 50 mm and y = ± 50 mm;
  • • the scattering anisotropy parameter g = 0.95;
  • • the scatter coefficient µs = 10 mm-1 for constant scatter, and µs = 10 (λ/800)−1.2 mm-1 for wavelength-dependent scatter, where λ is expressed in nanometers;
  • • the detector has a circular collection area of radius of 2 mm;
  • • the distance between the center of the top surface of the slab (0,0,0) and the midpoint between the source and detector (Δx,0,0) was varied between Δx = 0 and Δx = 49 mm.

Using 108 injected photons, values of attenuation A(λ) and DPF β(λ) were estimated at eleven discrete wavelengths (λ = 750, 760, 770, …, 850 nm). The first derivative of the attenuation with respect to wavelength ∂A/∂λ was then computed using an eighth-order polynomial fit to the values of A(λ).

4.2 Results for constant scatter

Figure 8 shows plots of the four terms in Eq. (5) at nine discrete wavelengths (λ = 760, 770, 780, …, 840 nm) for the Monte Carlo slab model with constant scattering and the source-detector midpoint a) at the center of the slab, and b) very near the lateral edge of the slab.

 figure: Fig. 8.

Fig. 8. Estimation of the four terms in Eq. (5) for diffuse reflectance measurements on a Monte Carlo slab model with brain-like absorption and constant scatter, for the source-detector midpoint a) at the center of the slab (Δx = 0 mm), and b) very near the lateral edge of the slab (Δx = 49 mm).

Download Full Size | PDF

It is evident that the approximate relationships expressed by Eqs. (7) and (8) appear to hold at least as strongly for this model as they did for the analytical model. Furthermore, and perhaps most surprisingly, they hold just as strongly when the source-detector pair are located very close to the lateral edge of the slab. In other words, inevitable changes in the ∂G/∂λ term are apparently compensated for by corresponding changes in the term containing ∂β/∂λ.

Using the values of ∂A/∂λ and β(λ) at the same nine discrete wavelengths, values of cHbO2, cHb, W and L were derived by finding the minimum norm least-squares solution to Eq. (9). The derived and true values of the parameters are shown in Table 3. The correspondence between the derived and true values is remarkable, even when the source-detector pair is located adjacent to the edge of the slab. Generally, the results are superior in accuracy to the corresponding results for the analytical model shown in Table 1, despite the fact that Monte Carlo estimates of attenuation and DPF inevitably include noise due to the finite number of detected photons.

Tables Icon

Table 3. True and estimated values of the parameters in Eq. (6), obtained by solving Eq. (9) using values of ∂A/∂λ at nine discrete wavelengths and either β(λ) or <β> derived from a Monte Carlo slab model with brain-like absorption and constant (wavelength-independent) scatter. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

Note that the extinction coefficient gradients were estimated using the same sampling interval (10 nm) as the modelled values of attenuation and DPF. The errors in the derived values are significantly worse if, instead, the gradients are estimated at the same interval (1 nm) as used for the analytical model. The critical importance of spectral sampling is discussed below in section 6.

4.3 Results for wavelength-dependent scatter

Figure 9 shows the four terms in Eq. (5) at the nine discrete wavelengths for the Monte Carlo slab model with wavelength-dependent scattering and the source-detector midpoint a) at the center of the slab, and b) very near the lateral edge of the slab.

 figure: Fig. 9.

Fig. 9. Estimation of the four terms in Eq. (5) for diffuse reflectance measurements on a Monte Carlo slab model with brain-like absorption and wavelength-dependent scatter, for the source-detector midpoint a) at the center of the slab (Δx = 0 mm), and b) very near the lateral edge of the slab (Δx = 49 mm).

Download Full Size | PDF

As for the results achieved using the analytical model, Fig. 9 exhibits offsets between the terms $\frac{{\partial A}}{{\partial \lambda }}$ and $d\beta \frac{{\partial {\mu _a}}}{{\partial \lambda }}$, and indicates that the sum of terms in Eq. (7) is no longer approximately zero. The mean vertical offset is found to be approximately 30% larger for Δx = 49 mm. For both source-detector locations, an unsuccessful attempt was made to recover the chromophore parameters from the nine discrete values of ∂A/∂λ and β(λ) using Eq. (12). The combination of noise (due the finite number of photons) and the 10 nm discrete sampling of attenuation (leading to poor estimates of ∂A/∂λ) did not enable all six parameters (cHbO2, cHb, W, L, α and p) to be derived simultaneously with reasonable accuracy.

5. Experimental investigation

To explore whether the spectroscopic method presented above can be applied successfully to experimental data, a simple experiment was performed on a sample of turbid fluid consisting of an aqueous mixture of red and blue inks and polystyrene microspheres. The experimental set-up is illustrated in Fig. 10. The fluid was contained in a rectangular glass cuvette cell with internal dimensions of 40 × 45 × 10 mm, with a wall thickness of approximately 1.25 mm. The red and blue inks (“Deep Red” and “Ultramarine” drawing inks, Winsor & Newton, UK) are fully translucent in water. The microspheres (Bangs Laboratories, USA) have a specified diameter range of 0.95–1.05 µm (mean diameter = 1 µm) and are supplied in a 10% solution (100 mg per mL of water). Transmission spectra across the 10 mm thickness of the cell were acquired using a broadband tungsten-halogen light source (HL-2000-HP, Ocean Optics) and a spectrometer (HR2000+, Ocean Optics, USA). The output power of the source was about 8.8 mW, and a spot of approximately 3 mm diameter was illuminated at the center of one face of the cell. Transmitted light was delivered to the spectrometer via a 600 µm diameter optical fiber from a position immediately opposite the point of illumination as shown. The spectrometer provides a wavelength sampling interval of about 0.45 nm.

 figure: Fig. 10.

Fig. 10. Experimental set up for measurement of transmitted light across a glass cuvette cell containing inks and polystyrene microspheres in water.

Download Full Size | PDF

First, the absorbing characteristics of both inks were determined using known concentrations in water, and using a pure water measurement as a reference. The absorption coefficients of the (undiluted) inks, µa,red and µa,blue, are shown in Fig. 11(a). The scattering characteristics of the microspheres were estimated using a Mie scattering model based on the known wavelength-dependent refractive indices of polystyrene and water [24]. The cell was then filled with water containing a precisely measured amount of both inks and a quantity of microspheres calculated to provide a transport scatter coefficient as shown in Fig. 11(b). The combination of inks within the cell provided an overall absorption coefficient as illustrated in Fig. 11(c). Next, the differential pathlength factor β(λ) at nineteen discrete wavelengths (λ = 480, 490, 500, …, 660 nm) was estimated using a Monte Carlo model of the cell geometry and optical properties, and the values are shown in Fig. 11(d).

 figure: Fig. 11.

Fig. 11. a) The measured absorption coefficients of the undiluted red and blue inks; b) The transport scatter coefficient of the fluid predicted from Mie theory; c) The absorption coefficient of the fluid due to both inks; d) The differential pathlength factor for transmission measurements across the fluid predicted from a series of Monte Carlo simulations at 19 discrete wavelengths.

Download Full Size | PDF

The transmitted intensity across the fluid-filled cell as a function of wavelength I(λ) was recorded. For experimental measurements it is necessary to compensate for the spectral characteristics of the source and the detector. In this case this was easily achieved by recording the intensity IEC(λ) transmitted through an empty cell. The relative attenuation of light through the fluid-filled cell and the empty cell are defined as A(λ) = –ln [I(λ)] and AEC(λ) = –ln [IEC(λ)] respectively. It is important to note that for the analysis presented below, the absolute values of I(λ) and IEC(λ) are irrelevant, except that obviously higher intensities yields results with lower noise.

The objective of the experiment was to determine the absolute quantities of ink within the fluid using measurements of ∂A/∂λ and an estimate of the DPF. Given that the transport scatter coefficient µs′(λ) is known to be wavelength dependent, Eq. (12) was adapted as follows:

$$\frac{{\partial A}}{{\partial \lambda }} - \frac{{\partial {A_{EC}}}}{{\partial \lambda }} \approx d\beta \left( {R\frac{{\partial {\mu_{a,red}}}}{{\partial \lambda }} + B\frac{{\partial {\mu_{a,blue}}}}{{\partial \lambda }} + W\frac{{\partial {\mu_{a,water}}}}{{\partial \lambda }}} \right) + \alpha {\lambda ^{ - p}}$$
where R and B are the fractions of red and blue ink in the fluid and W is the fraction of water. The first derivatives of the A, AEC, µa,red, µa,blue and µa,water with respect to wavelength were computed using an 20th-order polynomial fit to each set of values sampled at 0.45 nm intervals over a wavelength range from 480 nm to 660 nm. A twelfth-order polynomial was also fitted to the DPF values obtained from the Monte Carlo model to generate an estimate of β(λ) for the same wavelength range and interval.

Figure 12 shows plots of the terms $\frac{{\partial A}}{{\partial \lambda }} - \frac{{\partial {A_{EC}}}}{{\partial \lambda }}$ and $d\beta \left( {R\frac{{\partial {\mu_{a,red}}}}{{\partial \lambda }} + B\frac{{\partial {\mu_{a,blue}}}}{{\partial \lambda }} + W\frac{{\partial {\mu_{a,water}}}}{{\partial \lambda }}} \right)$ on either side of Eq. (13), generated using the true values of fractions R, B, and W. Unlike the modelled data for wavelength-dependent scatter shown in Figs. 4(a) and 9, an offset between the terms is almost imperceptible, and thus there is unlikely to be a significant contribution from the scattering term in Eq. (13).

 figure: Fig. 12.

Fig. 12. Plots of terms on the left and right hand side of Eq. (13) generated from experimental measurements of transmission across a cell containing red and blue inks, water, and polystyrene microspheres.

Download Full Size | PDF

Using 221 discrete values of ∂A/∂λ and β(λ) over a wavelength range from 520 nm to 620 nm, values of R, B, W, α and p were derived by finding the minimum norm least-squares solution to Eq. (13). The results are summarized in Table 4 which also shows the ratio of R and B. Values of the derived parameters are also shown when β(λ) is replaced with its mean value over the 520–620 nm range, <β> = 3.469, and when the contribution of scatter is excluded (i.e. when α = 0). Note that water absorption and its wavelength derivative are extremely small compared to those of the inks at the visible wavelengths employed for this experiment. It was therefore necessary to restrict the water fraction W to values between 0 and 1, although imposing this restriction had minimal impact on the derived estimates of the ink fractions.

Tables Icon

Table 4. True and estimated values of the fractions of red (R) and blue (B) ink in the sample of turbid fluid, obtained by solving Eq. (13) using values of ∂A/∂λ at 221 discrete wavelengths in the range 520-620 nm, and either β(λ) or <β> derived from a Monte Carlo model. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

The estimated experimental error in the true values of R and B is about 1 percent. The accuracy of the derived values is remarkably good, even when the DPF is replaced with its mean value. As predicted from the comparison of terms in Fig. 12, including the scatter term in Eq. (13) evidently produces very little benefit.

To verify that the method is tolerant of changes in measurement geometry, the whole experiment was repeated using a much narrower cuvette cell, of dimensions 40 × 10 × 10 mm. The only difference in the analysis was that a new Monte Carlo model was created to derive β(λ) for the smaller cell, which inevitably had a smaller mean value of 2.999 over the 520–620 nm range. The derived values of R and B were no less accurate, with errors of less than five percent. Again, the inclusion of the ad hoc scatter correction had little influence on the observed accuracy.

Currently we can offer no definitive explanation for the apparent lack of offset in the experimental data shown in Fig. 12, despite the wavelength dependent scatter. The maximum absolute value of ∂µa/∂λ for the fluid was roughly an order of magnitude greater than the corresponding value for the modelled slabs, which may be a contributing factor.

6. Discussion and conclusions

It is important to appreciate that all of the above analysis is based on modelled and experimental values of a measurable quantity ∂A/∂λ that is independent of absolute intensity, and the derived parameter estimates are of absolute concentrations of absorbers. This is potentially a major advance on conventional NIRS that relies on measured changes in intensity to generate values of changes in concentrations. Furthermore, in principle measurements of ∂A/∂λ may be significantly less influenced by changes in coupling, thus reducing if not eliminating the motion artefact which contaminates much NIRS and DOT data acquired on human subjects. The assumption that the measurements of ∂A/∂λ are independent of coupling was indirectly verified for the turbid fluid experiments, which took no account of the inevitable and variable losses occurring at both surfaces of the cell. As previously proposed [16], rapid measurements of ∂A/∂λ can be made using a wavelength-modulated source.

The most surprising result is that the derivation of absorber concentrations from measurements of ∂A/∂λ appears to be highly robust to changes in measurement geometry providing that a sufficiently accurate estimate (or measurement) of the DPF is available, i.e. no independent assessment of the scattering-dependent geometric parameter G is required. As mentioned above, the DPF can be measured at discrete wavelengths directly using time- or frequency-domain instruments [1,17], or by using the water absorption spectrum [6]. However, the results presented here suggest that replacing the wavelength dependent DPF with its mean value does not hugely degrade the accuracy of the results. Furthermore, the derived values scale with the DPF, and the variation of DPF in the human head, for example, is relatively small [25].

The least satisfactory aspect of the proposed method is the ad-hoc, empirically-derived correction for the wavelength dependency of scattering. While the analytical and Monte Carlo models both demonstrated that measurements of ∂A/∂λ are influenced by this dependency, the experimental results revealed that that there are circumstances where this influence can be considered negligible. A more detailed theoretical investigation of the origin of the apparent offset between terms is warranted.

Yeganeh et al. [26] have demonstrated that crosstalk between derived parameters can be reduced by employing the second and well as the first derivative of spectra. While the first derivative eliminates the contribution of wavelength-independent effects (such as coupling), the second derivative eliminates effects which have a linear dependent on wavelength (such as εHbO2 over the 740 nm – 860 nm wavelength range).

As highlighted in a previous publication [16], not only is the measurement of ∂A/∂λ dependent on the spectral resolution of the system, but so are the estimates of the chromophore extinction coefficients and their derivatives. If the two resolutions are significantly different, then the fitting process is unlikely to produce realistic values of the concentrations. This problem was carefully avoided for the analytical and Monte Carlo models by using exactly the same chromophore data to create the estimates of attenuation and to solve for the concentrations. Likewise, the experiments used the same system to measure the absorption coefficients of the inks and to measure the transmitted intensity across the fluid-filled cell. When estimating spectral gradients using polynomial fitting, for example, the same spectral sampling interval and the same degree of polynomial should be used. However, if different sampling intervals or different order polynomials are used, the derived concentrations may be significantly less accurate. So, for any practical application of the method, the derivatives of the chromophore extinction spectra are required at the same spectral resolution as the system. Potentially, smoothing could be employed to achieve a common spectral resolution. The sampling intervals for the experimental measurements and the diffusion model were 0.45 nm and 1 nm respectively, but because of the high computational burden, the Monte Carlo model used a sampling interval of 10 nm. This greater interval inevitably led to greater uncertainty in the values of ∂A/∂λ and ultimately prevented realistic values being derived when the scatter was wavelength dependent.

A second sampling issue concerns the number of discrete values of ∂A/∂λ used to estimate the chromophore concentrations to ensure uniqueness. The number of samples should exceed the number of fit parameters (between four and six). This number was nine for both the analytical and Monte Carlo models, and was 221 for the experimental measurements. To ensure a unique solution, it is critically important that the discrete wavelengths at which ∂A/∂λ is measured should correspond to regions of the spectrum where first derivatives of the chromophore extinction coefficients are most significant and are distinct from each other. This and other practical aspects of the proposed method need to be explored further.

Funding

Philippine Department of Science and Technology - Science and Education Institute (DOST-SEI).

Acknowledgments

The authors would like to thank Dr. Rob Cooper for helpful discussions.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” NeuroImage 63(2), 921–935 (2012). [CrossRef]  

2. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. M. Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” NeuroImage 85(1), 6–27 (2014). [CrossRef]  

3. Y. Hoshi and Y. Yamada, “Overview of diffuse optical tomography and its clinical applications,” J. Biomed. Opt 21(9), 091312 (2016). [CrossRef]  

4. S. Suzuki, S. Takasaki, T. Ozaki, and Y. Kobayashi, “A tissue oxygenation monitor using NIR spatially resolved spectroscopy,” Proc. SPIE 3597, 582–592 (1999). [CrossRef]  

5. S. J. Matcher, M. Cope, and D. T. Delpy, “In vivo measurements of the wavelength dependence of tissue-scattering coefficients between 760 and 900 nm measured with time-resolved spectroscopy,” Appl. Opt. 36(1), 386–396 (1997). [CrossRef]  

6. S. J. Matcher, M. Cope, and D. T. Delpy, “Use of the water absorption spectrum to quantify tissue chromophore concentration changes in near-infrared spectroscopy,” Phys. Med. Biol. 39(1), 177–196 (1994). [CrossRef]  

7. M. Ferrari, D. A. Wilson, D. F. Hanley, J. F. Hartmann, M. C. Rogers, and R. J. Traystman, “Noninvasive determination of hemoglobin saturation in dogs by derivative near-infrared spectroscopy,” Am. J. Physiol. Heart Circ. 256(5), H1493–H1499 (1989). [CrossRef]  

8. L. Fernandes, S. Carvalho, I. Carneiro, R. Henrique, V. Tuchin, H. Oliveira, and L. Oliveira, “Diffuse reflectance and machine learning techniques to differentiate colorectal cancer ex vivo,” Chaos 31(5), 053118 (2021). [CrossRef]  

9. A. R. Botelho, H. F. Silva, I. S. Martins, I. C. Carneiro, S. D. Carvalho, R. M. Henrique, V. V. Tuchin, and L. M. Oliveira, “Fast calculation of spectral optical properties and pigment content detection in human normal and pathological kidney,” Spectrochim. Acta, Part A 286(5), 122002 (2023). [CrossRef]  

10. L. Pollonini, H. Bortfeld, and J. S. Oghalai, “PHOEBE: a method for real time mapping of optodes-scalp coupling in functional near-infrared spectroscopy,” Biomed. Opt. Express 7(12), 5104–5119 (2016). [CrossRef]  

11. I. Tachtsidis and F. Scholkmann, “False positives and false negatives in functional near-infrared spectroscopy: issues, challenges, and the way forward,” Neurophotonics 3(3), 031405 (2016). [CrossRef]  

12. A. von Lühmann, Z. Boukouvalas, K. R. Müller, and T. Adalı, “A new blind source separation framework for signal analysis and artifact rejection in functional near-infrared spectroscopy,” Neuroimage 200, 72–88 (2019). [CrossRef]  

13. H. Dehghani, F. Leblond, B. W. Pogue, and F. Chauchard, “Application of spectral derivative data in visible and near-infrared spectroscopy,” Phys. Med. Biol. 55(12), 3381–3399 (2010). [CrossRef]  

14. H. Dehghani, F. Leblond, B. W. Pogue, and F. Chauchard, “Application of spectral derivative data in spectral near infrared tomography,” Proc. SPIE 7896, 78960I (2011). [CrossRef]  

15. O. Pucci, V. Toronov, and K. St. Lawrence, “Measurement of the optical properties of a two-layer model of the human head using broadband near-infrared spectroscopy,” Appl. Opt. 49(32), 6324–6332 (2010). [CrossRef]  

16. J. C. Hebden, “Exploring the feasibility of wavelength modulated near-infrared spectroscopy,” J. Biomed. Opt. 25(11), 110501 (2020). [CrossRef]  

17. D. T. Delpy, M. Cope, P. van der Zee, S. Arridge, S. Wray, and J. Wyatt, “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol. 33(12), 1433–1442 (1988). [CrossRef]  

18. D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36(19), 4587–4599 (1997). [CrossRef]  

19. T. Correia, A. Gibson, and J. C. Hebden, “Identification of the optimal wavelengths for optical topography: a photon measurement density function analysis,” J. Biomed. Opt. 15(5), 056002 (2010). [CrossRef]  

20. https://omlc.org/spectra

21. R. B. Saager, A. Quach, R. A. Rowland, M. L. Baldado, and A. J. Durkin, “Low-cost tissue simulating phantoms with adjustable wavelength-dependent scattering properties in the visible and infrared ranges,” J. Biomed. Opt 21(6), 067001 (2016). [CrossRef]  

22. J. R. Mourant, J. P. Freyer, A. H. Hielscher, A. A. Eick, D. Shen, and T. M. Johnson, “Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics,” Appl. Opt. 37(16), 3586 (1998). [CrossRef]  

23. http://mcx.space/

24. C. F. Bohren and D. R. Huffmann, Absorption and scattering of light by small particles, 1st ed. (Wiley1998).

25. F. Scholkmann and M. Wolf, “General equation for the differential pathlength factor of the frontal human head depending on wavelength and age,” J. Biomed. Opt 18(10), 105004 (2013). [CrossRef]  

26. H. Z. Yeganeh, V. Toronov, J. T. Elliott, M. Diop, T.-Y. Lee, and K. St. Lawrence, “Broadband continuous-wave technique to measure baseline values and changes in the tissue chromophore concentrations,” Biomed. Opt. Express 3(11), 2761–2770 (2012). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (12)

Fig. 1.
Fig. 1. a) Extinction coefficients of oxy- and deoxy-hemoglobin; b) Absorption coefficients of water and lipid; c) The first derivatives of the specific absorption coefficients of oxy- and deoxy-hemoglobin with respect to wavelength; d) The first derivatives of the absorption coefficients of water and lipid with respect to wavelength.
Fig. 2.
Fig. 2. Estimation of a) attenuation and b) differential pathlength factor for diffuse reflectance measurements on an analytical slab model with brain-like absorption and for both constant (black) and wavelength-dependent (blue) scattering.
Fig. 3.
Fig. 3. a) Estimation of the four terms in Eq. (5) for diffuse reflectance measurements on an analytical slab model with brain-like absorption and constant scatter; b) The sum of the terms in Eq. (7), and the difference between the terms in Eq. (8).
Fig. 4.
Fig. 4. a) Estimation of the four terms in Eq. (5) for reflectance measurements on an analytical slab model with brain-like absorption and wavelength-dependent scatter. b) The sum of the terms in Eq. (7), and the difference between the terms in Eq. (8).
Fig. 5.
Fig. 5. a) The offset f(λ) derived for different values of the scatter power p. b) The offset f(λ) divided by the transport scatter coefficient µs′ for different values of the scatter power p.
Fig. 6.
Fig. 6. a) Derived values of the concentrations of oxy- and deoxy-hemoglobin for different oxygen saturations. The continuous lines represent the true values. b) Derived values of the fractions of water and lipid, and the estimated oxygenation saturation. The continuous lines represent the true values.
Fig. 7.
Fig. 7. The slab geometry of a Monte Carlo model of photon migration. The source and detector are aligned along the y axis and separated by d = 20 mm, and the midpoint was varied between the center of the slab at (0,0,0) and a point near the edge at (49,0,0).
Fig. 8.
Fig. 8. Estimation of the four terms in Eq. (5) for diffuse reflectance measurements on a Monte Carlo slab model with brain-like absorption and constant scatter, for the source-detector midpoint a) at the center of the slab (Δx = 0 mm), and b) very near the lateral edge of the slab (Δx = 49 mm).
Fig. 9.
Fig. 9. Estimation of the four terms in Eq. (5) for diffuse reflectance measurements on a Monte Carlo slab model with brain-like absorption and wavelength-dependent scatter, for the source-detector midpoint a) at the center of the slab (Δx = 0 mm), and b) very near the lateral edge of the slab (Δx = 49 mm).
Fig. 10.
Fig. 10. Experimental set up for measurement of transmitted light across a glass cuvette cell containing inks and polystyrene microspheres in water.
Fig. 11.
Fig. 11. a) The measured absorption coefficients of the undiluted red and blue inks; b) The transport scatter coefficient of the fluid predicted from Mie theory; c) The absorption coefficient of the fluid due to both inks; d) The differential pathlength factor for transmission measurements across the fluid predicted from a series of Monte Carlo simulations at 19 discrete wavelengths.
Fig. 12.
Fig. 12. Plots of terms on the left and right hand side of Eq. (13) generated from experimental measurements of transmission across a cell containing red and blue inks, water, and polystyrene microspheres.

Tables (4)

Tables Icon

Table 1. True and estimated values of the parameters in Eq. (6), obtained by solving Eq. (9) using values of ∂A/∂λ at nine discrete wavelengths and either β(λ) or <β> derived from an analytical slab model with brain-like absorption and constant (wavelength-independent) scatter. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

Tables Icon

Table 2. True and estimated values of the parameters in Eq. (6), obtained by solving Eq. (12) using values of ∂A/∂λ at nine discrete wavelengths and either β(λ) or <β> derived from an analytical slab model with brain-like absorption and wavelength-dependent scatter. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

Tables Icon

Table 3. True and estimated values of the parameters in Eq. (6), obtained by solving Eq. (9) using values of ∂A/∂λ at nine discrete wavelengths and either β(λ) or <β> derived from a Monte Carlo slab model with brain-like absorption and constant (wavelength-independent) scatter. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

Tables Icon

Table 4. True and estimated values of the fractions of red (R) and blue (B) ink in the sample of turbid fluid, obtained by solving Eq. (13) using values of ∂A/∂λ at 221 discrete wavelengths in the range 520-620 nm, and either β(λ) or <β> derived from a Monte Carlo model. Numbers in parenthesis are the percentage differences between true and derived values (positive for overestimates, negative for underestimates).

Equations (13)

Equations on this page are rendered with MathJax. Learn more.

I = k I 0 e μ a β d G
μ a = n ε n c n
A = ln ( I I 0 ) = μ a β d + G ln ( k )
A λ = 1 I 0 I 0 λ 1 I I λ = d ( β μ a λ + β λ μ a ) + G λ 1 k k λ
A λ = d β μ a λ + d β λ μ a + G λ .
μ a = ϵ H b O 2 c H b O 2 + ϵ H b c H b + μ a , W a t e r W + μ a , L i p i d L + μ b ,
G λ + d β λ μ a 0
A λ d β μ a λ
A λ d β n c n ε n λ ,
A λ = d β n c n ε n λ r ( λ , p )
f ( λ ) = r ( λ , p ) r ( λ , 0 )
A λ d β n c n ε n λ α μ s = d β n c n ε n λ α λ p
A λ A E C λ d β ( R μ a , r e d λ + B μ a , b l u e λ + W μ a , w a t e r λ ) + α λ p
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.