1.

Introduction

Near-infrared (NIR) (650 to 1000 nm) optical properties of tissue offer insight into biological processes that occur in normal and pathological conditions. Diffuse reflectance spectroscopy is a technique that provides noninvasively quantitative information about the absorption and reduced scattering coefficients (${\mu }_a$ and ${\mu }_s^{\prime }$, respectively) of thick and multiple scattering tissues. It is used to quantify the concentration of important chromophores, such as water, hemoglobin, and fat, and to understand the composition and organization of tissue.¹

In vivo macroscopic optical properties ${\mu }_a$ and ${\mu }_s^{\prime }$ can be measured by performing a series of reflectance measurements using one of three major approaches.² All of them are based on the diffusion approximation of the radiative transport theory of light propagation, but experimentally, they are different in that the source is either a short pulse of light, sinusoidally modulated light, or a constant-intensity light. In the first approach, often referred to as the time domain (TD) method, as the pulse propagates through tissue, scatterers and absorbers attenuate, and temporally broaden the pulse. Information about the pulse attenuation and the temporal broadening is sufficient to extract ${\mu }_a$ and ${\mu }_s^{\prime }$ from measurements made at a single source-detector distance.^3,4 In the second approach, referred to as the frequency-domain (FD) method, sinusoidally intensity modulated light in the MHz range is used to probe tissue optical properties.^5–8 FD photon migration (FDPM) techniques measure the phase and amplitude of the light coming out of the tissue at one or more points on the surface, depending on how many source modulation frequencies are used. In the third approach, referred to as the steady-state (SS) method, tissue is excited by light of constant intensity, and measurements of light intensity at multiple distances from the source are performed.^9,10

SS and FD reflectance measurements can be used together to obtain broadband wavelength coverage with increased penetration depth.² In a typical FDPM instrument, most of the wavelength coverage is provided by a white-light SS measurement, whereas the FD data are acquired at a few selected wavelengths using laser diodes. Absorption and reduced scattering coefficients derived from the FD data are used to calibrate the intensity of the SS measurements and to estimate ${\mu }_s^{\prime }$ at all wavelengths in the spectral region of interest. Finally, ${\mu }_a$ can be determined by comparing the corrected SS reflectance values with the diffusion theory predictions for the wavelength range used.⁸

The number of wavelength points used to obtain the FD data basically limits this approach. We developed a method to improve the acquisition of optical parameters of the examined tissues, based on digital parallel acquisition in FD. With our system, both FD and SS measurements are performed using a pulsed super-continuum white laser alone. Moreover, the pulsed white laser allows a continuous scan of frequencies in the spectral medical window (650 to 1000 nm), thus providing additional wavelength information. At each wavelength, we extract what we refer to as the tissue phasor, a mathematical representation of phase shift and demodulation of the output light with respect to the source light. The estimated absorption and scattering coefficients are obtained by fitting phase shift and demodulation to the standard diffusion equation with semi-infinite geometry boundary conditions.

2.

Background and Theory

2.2.

Calibration and Data Analysis

The FDPM approach uses both FD and SS measurements in order to obtain the optical parameters over a large range of wavelengths. From the combination of these parameters, quantitative broadband ${\mu }_a$ spectra can be obtained. Each FD measurement contains instrumental artifacts (which can generate phase delay and demodulation). At each modulation frequency, the measured phase is the sum of the sample and the instrument phase, while the measured amplitude is the product of the sample and instrument amplitude. In order to estimate the instrument phase and modulation, a calibration is performed through the acquisition of FD data from a sample whose optical parameters are known a priori (from a set of two distance measurements). The instrumental response can be obtained by looking at the differences between the measured and predicted phase and the modulation amplitude from a sample of known parameters.

After calibration, data depend only on two unknown quantities, ${\mu }_a$ and ${\mu }_s^{\prime }$. Their values are obtained from the best fit of the data to Eqs. (3), (4) and/or (5), using a nonlinear least-squares method.¹⁶ Once the values of ${\mu }_s^{\prime }$ at the selected wavelengths measured in FD have been calculated, they are fitted to a power function of the form ${\mu }_s^{\prime }=A·{\lambda }^{-B}$ (where A and B are constants calculated through the fit process), as this function well describes the scattering wavelength dependence over the spectral medical window,¹⁷ providing a good estimates of ${\mu }_s^{\prime }$ at all wavelengths in the necessary range.

Single-distance SS reflectance measurements are usually made using a halogen lamp (a white laser in our case), and the light is analyzed by a spectrometer. The spectrum of the light source is also measured separately by using an integrating sphere or a standard reflectance disk (99% reflectance; Edmund Optics, Barrington, NJ). Finally, the relative reflectance is obtained by dividing the sample spectrum by the source spectrum. The instrumental factors, due to cable phase delays and demodulation, are independent of the wavelength; thus, in order to recover the actual ${\mu }_a$ spectrum, theoretically we could estimate them at only one wavelength. This is readily done at any of the FD wavelengths, as ${\mu }_a$ and ${\mu }_s^{\prime }$ are both known and absolute reflectance is a function of only these variables. Using all the FD data, we can calculate the value of the instrumental factor ($k$ factor) that scales the measured SS reflectance values to match the predicted reflectance as closely as possible (in the least-squares sense). With ${\mu }_s^{\prime }$ calculated and the reflectance correctly scaled, we can compute ${\mu }_a$ by comparing the calibrated SS reflectance measurements and the theoretical reflectance predicted by Eq. (2) for a given trial value of ${\mu }_a$ (Fig. 1).

Fig. 1

Steps for the determination of broadband ${\mu }_a$ spectrum. (a) Measured SS reflectance spectrum; (b) amplitude calibration of the measured SS reflectance spectrum through the determination of the instrumental factor $k$; (c) SS reflectance spectrum scaled to fit the discrete FD reflectance values; (d) least-squares power-law fit to the ${\mu }_s^{\prime }$ values measured by the FD; (e) actual broadband ${\mu }_a$ spectrum.

3.

Materials and Methods

In a typical FDPM instrument, the wavelength range is provided by a white-light SS spectral measurement and the FD data are acquired using several laser diodes that are sequentially modulated. Typically, this is done through a network analyzer or a RF circuitry.⁸ In our approach, we minimize the hardware complexity by introducing a single pulsed white laser as a source of excitation (Optical Supercontinuum Source SC450, Fianium Ltd, Southampton, UK) that generates a train of pulses with a repetition frequency of approximately 20 MHz in a wavelength range from 450 nm and to greater than 2500 nm, giving a total output power of more than 2 W ($\approx 1\mathrm{mW}/\mathrm{nm}$). The pulse train in the TD is represented by a number of harmonics in the FD, whose range depends inversely on the pulse width. The pulse width of the white laser is 3 ps, corresponding to a harmonic content of several hundred gigahertz. The detector has a rise time of 0.57 ns, corresponding to a bandwidth given by $BW\approx \frac{0.34}{t_r}\approx 600\mathrm{MHz}$.¹⁸ The FLIMbox card used for the measurements shown in this paper is characterized by the frequency response shown in Fig. 2. A faster card is also available, allowing greater signal harmonic content. Tissue is the important factor in the determination of the spectral information available. In particular, the intensity of the signal at higher modulation frequencies is much less than at lower frequencies (the modulation decreases), as predicted by Eqs. (3) and (4). In Fig. 2, it has been considered a tissue with typical physical parameters.¹⁵

Fig. 2

Frequency response of the system components: FLIMbox, detector, laser source, and tissue. Harmonic intensities are normalized with respect to the DC value (at $f=0$). Data are processed for values above the threshold (dashed line).

From Fig. 2, it is evident that although the limitation factor of the system in terms of amount of harmonics available is the FLIMbox card, this does not affect the potential information present in the light diffusing through the tissue and acquired by the detector, as tissue does not respond at high frequencies. As a consequence of the train of narrow pulses, we have simultaneous access to the range of modulation frequencies needed for the FD measurements.

The selection of a specific wavelength of the white laser is achieved through a motor-controlled monochromator. The pulsed white light is directed via the built-in optical fiber into the optical system, where an H10 monochromator (ISA Instruments S.A., Longjumeau, France) sequentially selects the wavelengths from 680 to 850 nm. The monochromator has an exit slit of 1 mm with a dispersion of $8\mathrm{nm}/\mathrm{mm}$. From the monochromator, the light is coupled into an optical fiber bundle to be delivered to the sample. The measured power at the sample is about 0.3 mW at 750 nm, in an area of about $10{\mathrm{mm}}^2$. Diffusely reflected light is collected and processed at two detecting channels simultaneously. Each channel consists of a photomultiplier (PMT) detector (Models R7400U-20 and H10721-20; Hamamatsu Corporation, Bridgewater, NJ) that is included in a custom-made, handheld probe. Moreover two high-pass glass filters (Model NT66-045; Edmund Optics, Barrington, NJ) have been used for the detection to filter out photons below a certain wavelength ($\approx 620\mathrm{nm}$). The probe also allows multiple source detector distances (Fig. 3) by placing the source fiber in the different available positions. The PMT modules have a gain curve depending on the voltage control that goes from ${10}^4$ to $3\times {10}^6$. The modules have a sensitivity of $78\mathrm{mA}/\mathrm{W}$ at 630 nm, a sensing element with an 8-mm diameter, and very low dark noise. The detector sensitivity and the high brightness of the source allow us to obtain satisfactory signals for distances greater than those commonly seen in conventional FDPM instruments. The signals from the PMTs are converted in “photon pulses” using a constant fraction discriminator and directed to the FLIMbox (Laboratory for Fluorescence Dynamics, University of California Irvine, Irvine, CA) that is synchronized with the white laser frequency of 20 MHz. The FLIMbox performs the digital FD acquisition based on the heterodyning principle using a field programmable gate array (FPGA).⁸ It replaces the RF mixing of analog signals with a digital parallel mixer scheme in the FPGA (rather than using analog mixers, which are noisy and have a low conversion efficiency) by using multiple sampling windows within each excitation period, as shown in Fig. 4. Because the mixing is performed digitally, the signal can be processed by an arbitrary number of sampling windows without any loss of signal (photocurrent pulses) at a 100% duty cycle. The PMT can be operated at full power during the entire acquisition. The photon counting operation mode further improves the signal-to-noise (S/N) ratio of the detection system. Measurements have been obtained with an acquisition of about ${10}^6\mathrm{photons}/\text{s}$. We maintain the input laser power at the sample constant using a feedback circuit. The changes in intensity at the detector depend on the sample absorption and scattering. Since data collection terminates when a given number of photons is collected, the only consequence of the transmissions at different wavelengths is to modify the integration time in our acquisition protocol. This means that all data points are collected with the same statistics and the same input power. Figure 5(a) shows a picture of the actual A320 FLIMbox module (ISS Inc., Champaign IL). Figure 5(b) shows the schematics of the FLIMbox module, which includes the preamplifiers, constant fraction discriminators for the two independent input channels, the FPGA chip (Spartan 3E; Xilinx Inc., San Jose, CA), and the connections for synchronization with the laser frequency.

Fig. 3

(a) Schematic of the hand-held probe used for the measurements, (b) showing the two PMTs and the available source-detector distances. (c) Picture of the actual probe with the relative distances employed during the measurements presented in this paper.

Fig. 4

Illustration of the digital heterodyning principle, with exaggerated heterodyning such that the cross-correlation frequency (fcc) is $5/8$ of the sampling frequency (fs) instead of fcc $5\mathrm{fs}/256$, which was used in our implementation. Arriving photons (dots) are assigned to one of four sampling windows according to their arrival time. In the real case, the sampling windows slide through the entire period of the emission response due to the slight difference in frequencies, for a total of 256 steps.

Fig. 5

(a) Picture of the A320 module housing the actual FPGA where the FLIMbox circuit is implemented. (b) Schematic of the modules inside the FLIMbox unit, showing the preamplifiers, constant fraction discriminators, the FPGA card, and the synchronization for the external laser.

The software SimFCS (Laboratory for Fluorescence Dynamics, University of California Irvine, Irvine, CA) calculates the phase shift and demodulation (relative to the reference channel). An HR4000 spectrometer (Ocean Optics, Dunedin, FL) with a sensitivity range from 680 to 1100 nm (not shown in the schematic) is used to compare the scatter-corrected absorption obtained from FDPM measurements with the SS signal.

A schematic of the parallel digital FD spectroscopy system is shown in Fig. 6. The instrument measures the phase and amplitude of the photon density waves at different modulation frequencies over a range of optical wavelengths.

Fig. 6

Schematic of the electronics and optical components of the digital FD parallel spectroscopy instrument. The flow of the signals is shown: the optical path (dashed line) goes from the white laser to the phantom and the electrical paths (solid lines) from the detectors to the FLIMbox to the computer.

4.

Results and Discussion

Tissue phantoms with known optical properties were used to determine the FDPM capability to quantify scattering. Scatter-corrected absorption spectra from FDPM were compared to those taken with the spectrometer. Reflectance measurements were performed with semi-infinite geometry placing the probe in contact with the surface of the phantom. The source fiber and detectors were at two relative distances, 3.5 and 5.5 cm. Measured optical properties were compared to the expected ${\mu }_a$ and ${\mu }_s^{\prime }$ parameters to determine the accuracy in characterizing optical properties. Figure 7 shows the actual phase and modulation data obtained at 750 nm, one of the 170 wavelengths available (from 680 to 850 nm, each).

Fig. 7

Phase shift and demodulation at one wavelength obtained from a phantom with known optical properties.

To determine the uncertainty in the measurements, we can calculate the standard deviation of the measured phase and modulation for each harmonic. Considering a Poissonian error, the standard deviations for the phase and for the modulation at each frequency are proportional to the inverse of the square root of the number of counts $N$ and can be expressed using the following mathematical model:¹⁹

A nonlinear least-squares approach is used to fit the model functions [described by Eqs. (3) to (5)] to the calibrated data. In principle, phase and modulation data, fit with their corresponding functions, provide a good estimation of optical properties.⁸

We noticed a correlation between the recovered optical coefficients, and we assumed that such behavior could be due to the model function pair used for the fitting. Consequently, we performed a chi-square ($\chi 2$) statistical analysis in order to find both the model function pair and the number of modulation frequencies necessary to best fit the acquired data. We used the following $\chi 2$ function for simultaneous fitting of different pairs: ${({\chi }_{\text{tot}})}^2={({\chi }_{f1})}^2+{({\chi }_{f1})}^2.$ In the present study, we compared six fitting model function pairs for robustness and accuracy in extracting optical properties: fitting phase ($\mathrm{\Phi }$) and modulation ($M$), fitting $\mathrm{\Phi }$ and $Y=(1-M)/\mathrm{\Phi }$, fitting $\mathrm{\Phi }$ and amplitude of the detected signal ($ac$), fitting $M$ and $ac$, fitting $Y$ and $ac$, and fitting the slope of $\mathrm{\Phi }$ and the slope of $ac$.

We computed the surfaces of the $\chi 2$ for each model function pair and evaluated the best pair that generates a parabolic shape with a well-defined minimum. Simulations were performed by selecting broad sets of combinations of ${\mu }_a$ and ${\mu }_s^{\prime }$. Combinations of Eqs. (3)–(5) were solved as photon density wave model function pairs in semi-infinite media. Noise with a normal distribution was added to the analytical solution of Eqs. (3)–(5) at each frequency point. Furthermore, we studied how the number of modulation frequencies affects the fits. We did not notice improvements by considering frequencies larger than 260 to 300 MHz; therefore, all simulations presented in this paper have a frequency range from 20 to 260 MHz. Also, the simulations show that about 12 to 15 frequencies equally spaced in this range are sufficient to determine the optical parameters and that further increase in the number of frequencies does not result in smaller errors.

The simulated surfaces are obtained by using experimentally derived uncertainties. Chi-squared surfaces for $\mathrm{\Phi }$ and $M$ [Fig. 8(a)] and $\mathrm{\Phi }$ and $Y$ [Fig. 8(b)], as well as for $M$ and $ac$, and $Y$ and $ac$ fits flatten with increasing absorption and measurement uncertainty, making it difficult for fitting algorithms to find a true global minimum. In contrast, the chi-squared surface for the fit using $\mathrm{\Phi }$ and $ac$ and their slope is parabolic with a well-defined minimum. Advantages of the fit using $\mathrm{\Phi }$ and $ac$, as well as their slopes as model functions, can be appreciated by visualizing the chi-squared surfaces, shown in Fig. 8(c) and 8(d), respectively.

Fig. 8

Chi-squared surfaces of fits performed using different pair functions: (a) $\varphi$ and $M$, (b) $\varphi$ and $Y$, (c) $\varphi$ and $ac$, and (d) slope $\varphi$ and slope $ac$.

Figure 9(a) shows the results of an FDPM fit from data acquired on a phantom with known optical properties, using $\mathrm{\Phi }$ and $M$ as the model function pair. The recovered ${\mu }_a$ and ${\mu }_s^{\prime }$ appear to be strongly correlated. Using $\mathrm{\Phi }$ and $ac$ Fig. 9(b) for the fit, the correlations eventually disappears, confirming the results obtained with the chi-square analysis.

Fig. 9

Comparison of FDPM fits with two different model function pairs: (a) $\varphi$ and $M$, (b) $\varphi$ and $ac$.

$\mathrm{\Phi }$ and $ac$ were used as model functions to fit FDPM data. Figure 10 shows the FDPM measurements of ${\mu }_a$ and ${\mu }_s^{\prime }$, the fit to obtain the scattering factor using the power law,⁷ and the comparison between the corrected spectrum taken with the spectrometer and the absorption spectrum obtained through the digital system. The spectrum taken with the spectrometer has been corrected through the correction factor ($k$ factor) that takes into account the calibration process described previously. Figure 9(b) shows that the large number of wavelengths used for the FDPM acquisition allows us to obtain a scattering curve approaching the one predicted by the power law approximation.

Fig. 10

(a) Absorption and scattering obtained from FDPM measurements. (b) Scattering from FDPM measurements overlapped by the scattering fit (power law) and scattering-corrected absorption. (c) Spectrometer spectrum, correction factor and corrected spectrum. (d) Comparison between spectrometer corrected spectrum and FDPM spectrum.

Multiple phantom measurements were performed to determine the repeatability of the FDPM to quantify ${\mu }_a$ and ${\mu }_s^{\prime }$. Overall, the results clearly show the capability of the new digital parallel acquisition system to perform FDPM measurements.

In this paper, we have discussed the factors that govern the accuracy and sensitivity of the FDPM-measured optical properties. The number of photons collected at each wavelength determines the noise of a single determination. Since data are collected until a given number of photon is reached, the noise is constant and predictable, according to Eqs. (6) and (7). The light propagation model determines the absolute accuracy of the recovered parameters. The type of data used in the fitting algorithms to recover the optical properties from measured data determines the correlations among parameters, as shown in the discussion of Fig. 8. We performed a data fit of the FD response after acquiring the entire spectrum to obtain the optical parameters. This fitting is almost instantaneous (less than 1 s for all wavelengths); therefore, it does not limit the throughput of our approach.

The digital mixer has no noise and no losses. Figure 11 shows the absence of noise (the uncertainty is smaller than the size of the symbols) in the digital approach. In the case of the digital system, the noise is predictable and it depends only on the number of photons detected, which determines the minimal noise attainable. The analog system has other sources of noise that are difficult to evaluate because they are specific to the implementation and to the instrument. For example, two “identical” analog systems could have quite different levels of noise. Two FLIMboxes measuring the same signal will have the same amount of noise. The FLIMbox has no “dead” time, per second; the dead time is due to the constant fraction discriminator that in our system is about 7 ns. The FLIMbox input circuit operates internally at about 320 MHz on two independent channels. There is no specific sampling rate. For this specific study, we divided the laser repetition period (50 ns) into 256 phase bins (of about 200 ps each) to measure the phase of the photon with respect to the laser period.

Fig. 11

Phase (dots) and modulation (squares) measurements showing the absence of noise. The uncertainty is smaller than the size of the symbols.

In conclusion, our system collects FD data at 170 wavelengths and allows accurate separation between absorption and scattering and the recovery of highly resolved spectra. The system performs parallel acquisition, which decreases the acquisition time required for standard serial systems that need to scan through all the modulation frequencies. Furthermore, the digital modulation removes analog noise (Fig. 11), and the analog mixer does not create radio frequency interference or emission. The brightness of the single laser source and the large detection area, allow a deep penetration with extremely low power on the sample. Finally, the detection scheme allows very sensitive acquisition due to the low S/N ratio.

The ability of the system to penetrate deep into the tissues makes it a good candidate for use in breast cancer detection and monitoring. In previous works,²⁰ we showed the possibility of noninvasively predicting chemotherapy response prior to treatment based on biomarkers obtained from tumor spatial heterogeneities of spectral features measured using a conventional DOS instrument. With the method presented in this paper, we are obtaining absorption spectra with a high number of FDPM wavelength points; thus, we expect to increase the accuracy of the obtained spectra. The more detailed spectra should improve the characterization of breast tumor spatial heterogeneities, and thus the prediction of chemotherapy outcomes.

There is ample room for improvement and future development of the system. First, in order to obtain multiple spatial points at the same time, an array of PMTs could be introduced to scan an area instead of a single point. Furthermore, the monochromator could be removed by employing an acousto-optic filter, thus removing the mechanical component and allowing even faster acquisition.