1.

Introduction

Spectroscopy of biological tissues in near-infrared (NIR) wavelengths are based on diffused photon density waves, optical absorption, and scattering and is extended for functional NIR (fNIR) imaging. Photons in the NIR wavelength range diffuse through tissue by scattering movement.¹ Since the individual photon undergoes a random walk within the scattering media (i.e., each photon travels in any spherical direction, resulting in a combined photon density wave); these photons are being scattered and a favorable path² is allocated at each modulation frequency. Various fNIR system topologies are reported based on the continuous wave (CW),³ time domain (TD),⁴ and frequency domain (FD).⁵ The CW is based on unmodulated photons, and the optical absorption coefficient⁶ can be extracted⁵ only by monitoring only the photon amplitude. TD systems depend on tunable mode-locked lasers to generate pulses with different time samples,⁴ from which both absorption and scattering coefficients⁷ are efficiently extracted by monitoring amplitude changes as a function of time due to propagation delay of photons in a scattered media.⁴ FD is based on either narrow-band (e.g., single modulating frequency I/Q system) or broadband frequency-swept diffused photon waves,³ and both absorption and scattering coefficients^8–10 are extracted by monitoring the change in the amplitude and IP of modulated photons. Note that short optical pulse and broadband frequency modulation are related through Fourier transformation.

Photon dynamics has been modeled using the transport theory of particles and extending the diffusion equation (DE) of photons to explain the photon frequency modulation or photon pulse modulation.^11–13 Extraction of the optical parameters of absorption (${\mu }_a$) and reduced scattering (${\mu }_s^{\prime }$) is an inverse problem of predicting amplitude and phase [i.e., insertion loss (IL) and insertion phase (IP) in frequency-modulated systems at each frequency] using the DE of NIR photons. A closed-form analytical expression is only reported for a homogenous medium¹⁴ and the DE solution of a multilayer and inhomogeneous medium is not very challenging. Alternatively, a time-consuming Monte Carlo (MC) photon counting could be used to extract the optical parameters; however, this technique is reported only for two-layer, turbid media.¹⁵

Our goal is to develop, as depicted in Fig. 1, an fNIR medical imaging system of a multilayer human head with multioptical transmitters and receivers that are strategically located on the head to ultimately monitor any brain physiological activity. The limitation of current modeling and extraction to either homogenous medium using DE, or two layers using MC, introduces an important drawback in addressing blind optical parameter extraction in inhomogeneous and multilayer tissue. The inverse problem calculation of the optical parameters, particularly, is very challenging for inhomogeneous and multilayer media [e.g., the human head with layers of skin, skull, cerebrospinal fluid (CSF), and inhomogeneous brain matter] and an alternative numerical solution is required to meet this challenge. For the first time, this paper proposes and reports on the accuracy of using finite element method (FEM)-based modeling for the blind extraction of four phantoms resembling brain matter. The challenge of accurate optical parameter extraction of inhomogeneous and multilayer brain material is addressed using FEM numerical calculations that are demonstrated to be reliable, accurate, and robust. As the accuracy of blind parameter extraction is demonstrated using FEM modeling, FEM-based modeling with appropriate boundary conditions can be extended to multilayer phantoms that resemble a multilayer human head. Successful blind extractions are required before development of accurate fNIR imaging that could potentially lead to predictions of traumatic brain injury (TBI). A technically viable, helmet-mounted fNIR imaging system is envisioned that can potentially be applied for TBI assessment through the extraction of optical parameters of the multilayer structure of the head and inhomogeneous brain.

Fig. 1

Photon frequency modulated 3-D sectional image of a human head and multioptical transmitters and receivers.

The objective of this work is to validate the use of a commercially available, FEM-based tool using a multiphysics package called COMSOL (COMSOL, Inc., Burlington, Massachusetts) to solve the inverse problem and to extract absorption and reduced scattering coefficients of various homogenous phantoms by fitting predicted results to broadband frequency modulated (from 30 to 1000 MHz) measurements. This paper is structured as follows. It begins by discussing frequency domain NIR spectroscopy and the need of FEM for the extraction of optical phantom parameters as background material to justify the utility of the presented procedures. First, a review of DE-based analytical modeling is presented as the basis of MC numerical simulations. Fundamental equations for solving the inverse problem using FEM are also reviewed. In Sec. 2, experimental results and both DE- and FEM-based signal processing are presented. More specifically, this paper employs the reported¹⁶ differential measurement procedures and novel signal processing techniques for accurate extraction of optical phantom parameters using DE. In Sec. 3, blind extraction of optical parameters for four solid phantoms that closely resemble different regions of brain are demonstrated using DE, and then these parameters are applied to numerical simulations in FEM to demonstrate accuracy and fast predictions of photon dynamics. Discussions and concluding remarks are presented in Sec. 4, emphasizing that this powerful simulation procedure is applicable to custom-designed TBI assessment using broadband frequency measurement of IP and IL.

2.

Modeling and Optical Parameters Extraction Techniques

2.1.

Transmission Mode Procedure

In this section, procedures for extracting the optical parameters are introduced. Blind extraction is performed for four different homogenous phantoms, and results were compared with the manufacturer’s optical properties reported at 630 nm. Table 1 summarizes the optical properties of the four phantoms provided by the phantom manufacturer at 630 nm using an alternative technique, which were blindly extracted using the procedure outlined later in this article in Sec. 3. For each phantom, the extraction of optical parameters is accomplished by minimizing error between the predicted and measured IL and IP of modulated photons. Cylindrically shaped solid phantoms with a diameter of 2.5 cm and height of 9 cm (see Fig. 2) are employed for the experimental measurements and blind optical parameter extractions.

Table 1

Manufacturer-provided optical parameters for phantoms at 630 nm.

	${\mu }_a({\mathrm{cm}}^{-1})$	${\mu }_s^{\prime }({\mathrm{cm}}^{-1})$
Phantom 1	0.45	19.9
Phantom 2	0.67	19.9
Phantom 3	1.08	19.9
Phantom 4	3.01	19.9

Fig. 2

Experimental setup for the characterization of a solid phantom using fNIR spectroscopy: network analyzer as an accurate measurement of IL and IP, solid phantom, hybrid optical Rx, hybrid optical Tx, and mechanical manipulator for accurate positioning of optical Tx and Rx for fixed separations.

In transmission mode, the photons are forward scattered, as opposed to reflection mode, where they are backscattered. The reflection mode technique requires a very sensitive receiver, especially when the absorption coefficient is high in the homogenous media. As a result, a transmission mode technique is adopted here because there is a higher signal-to-noise ratio for a similar receiver separation. The overall measurement setup with custom-designed optical Tx and Rx is depicted in Fig. 2. An automatic network analyzer (Anritsu MS4623B; Test Equipment Connection, Lake Mary, Florida) was employed to measure the forward IP and IL of modulated photons, where the rf source used for modulating optical Tx is compared to a very sensitive received rf. The customized optical Tx used the rf switch from Hittite (HMC245QS16) to drive three high-powered vertical cavity surface-emitting lasers (VCSELs) (at 670, 795, and 850 nm). A printed board was designed and fabricated using a commercial substrate (FR4) to accommodate all surface mount technology (SMT) components. Details of this optical Tx was discussed in Ref. 16, and it includes a tri-wavelength source from Vixar (Plymouth, MN; Module V3WLM) with pin designations for three wavelengths: 680, 795, and 850 nm. High-powered lasers are employed with each VCSEL’s threshold currents: 8.5 mA for the 680 nm, 9.5 mA for the 795 nm, and 2.8 mA for the 850 nm. Even though in other publications,^26,27 we have described a custom-designed optical Rx with a lower-operating voltage, better light collection efficiency, and better sensitivity than the available commercial avalanche photodiode (APD); however, throughout this set of measurements, we have used optical Rx based on APD module C5658 (Hamamatsu Phonics, Tokyo, Japan) along with a built-in trans-impedance amplifier (TIA). This commercial module provides high-sensitivity optical measurements with a spectral response in the range of 400 to 1000 nm. The built-in TIA provides 34 dB of power amplification gain up to 1 GHz. The overall optical sensitivity is $-40\mathrm{dBm}$ at 800 nm, with a compression dynamic range of about 70 dB for all three NIR wavelengths of interest. We previously showed in Ref. 16 the full capability of our built-in system. More specifically, it demonstrated a novel extraction procedure using the analytical DE method for homogenous phantom resembling white brain matter for the fNIR three wavelengths of interest.

2.2.

Signal Processing for DE Using the Two Separations Measurement Method

The solution of DE is presented in terms of IL and IP for a given phantom optical parameter values. The optical parameters ${\mu }_a$ and ${\mu }_s^{\prime }$ were extracted by best fitting of the predicted DE results to the measured experimental results. For greater accuracy, a curve-fitting procedure is proposed by approximating the IL and IP as a function of modulating frequency (i.e., $f$) to the complex $a\sqrt{\mathrm{freq}}+b$-form; the least-square curve algorithm [e.g., lsqcurvefit from MATLAB (Mathworks, Natick, Massachusetts)] is used to minimize differences between measured and calculated values, based on the Levenberg-Macquart technique. In particular, the measured IL and IP of modulated photons at a 680-nm wavelength were curve fitted, and then the least-square fit was employed to extract the optical parameters at 680 nm. Since the rf electrical loss was the square of the optical loss, a conversion factor of 2 dB electrical loss was made for every 1 dB optical loss. Noise sources from the optical transmitter and receiver were individually identified and their predicted values are combined to quantify the level of measurement noise floor. Systemic errors at the interface of the optical Tx and optical Rx with the solid phantom and any artifact associated with frequency response deviation of a broadband laser driver and tarnsimpedance amplifier boards were removed for more accurate fitting and extraction.

A two-separation measurement method¹⁶ was employed to remove any interface artifacts and accurately measure IL and IP. More specifically, as shown in Fig. 3, two separations of 0.5 and 0.8 cm were considered for transmission mode measurement and calculation of incremental IL and IP as

Eq. (7)

$$\mathrm{\Delta }\mathrm{IL}(\mathrm{dB})={\mathrm{IL}}_{d2=0.8\mathrm{cm}}-{\mathrm{IL}}_{d1=0.5\mathrm{cm}}.$$

Eq. (8)

$$\mathrm{\Delta }\mathrm{IP}(\mathrm{Deg})={\mathrm{IP}}_{d2=0.8\mathrm{cm}}-{\mathrm{IP}}_{d1=0.5\mathrm{cm}}.$$

Fig. 3

Structure of transmission mode measurement and simulation. (a) The mesh structure used for a 2-D modeling of DPDW in a solid phantom. (b) A two-separation method setup for eliminating common mode noise and artifacts.

2.3.

FEM Simulation and Best-Fitting Using COMSOL

The accuracy of finite element simulation results depends on the resolution of FEM mesh and associated points in space where the diffused photon density wave (DPDW) solution is to be computed using the partial difference equation of the Helmholtz equation. A higher-resolution meshing structure is required to accurately simulate the photon dynamics, which naturally results in longer computation times or even causes the computation system to crash; therefore, a compromise in accuracy and mesh resolution has to be made. The scope of simulation and time required depend on the mesh resolution and size of structure being modeled. Here, we will focus on a two-dimensional (2-D) structure [as opposed to three-dimensional (3-D)] to perform simulation in less time. As confidence is established here, the work will be extended later to the 3-D model. Our efforts is focused on random (as opposed to uniform) meshing to achieve great accuracy in a shorter simulation time, and they also extend to 3-D modeling of large phantom geometries.

The optical parameter extraction is based on subtracting the common source of error from the two consecutive distance measurements using Eqs. (7) and (8). The challenge of the FEM is the required computation resources for reaching an accurate solution of diffused photon density wave dynamics for given homogenous and inhomogeneous media. To meet the typical computation time of 45 s using an Intel Core i3 processor PC with a reasonable error of less than 5%, mesh element sizes of 80 and 1000 nm are considered in regions around optical Tx and Rx and everywhere else, respectively, as depicted in Fig. 3. These meshing element sizes are also recommended when human tissue is approximated.²⁸ In COMSOL, the boundary condition is set either through Dirichlet or Neumann boundary conditions, where they express the fluence rate $u$ at desired boundaries. Dirichlet and Neumann boundary conditions are chosen appropriately for air-dielectric (or dielectric-dielectric) and radiation conditions in the region surrounding the phantom. The mathematical representation is expressed as

Eq. (9)

$$n·[{(c▿u)}_1]-[{(c▿u)}_2]+q·u=g$$

Eq. (10)

$$h·u=r,$$

where $h,r,q$, and $g$ are the boundary coefficients for the phantom modeling. For the diffused photon flux, $h=1$ and $r=0$. For the matching boundary of air dielectric and dielectric and dielectric interfaces $q=0$ and $g=0$, it is important to mention that extrapolated boundary condition take into account the phantom-air boundary as a radiation boundary by setting the fluence rate to zero using Eq. (10) without any impact on the error achieved.

The simulation results are shown in Figs. 4 and 5 for both IL and IP at frequencies of 100 MHz and 1000 MHz for optical source wavelength of 680 nm from the Vixar VCSEL. Full broadband frequency modulation from 30 to 1000 MHz has been analyzed using COMSOL, and the next section presents results for the wavelength of 680 nm. Note that the phantom manufacturer had provided the optical parameters of ${\mu }_s^{\prime }=19{\mathrm{cm}}^{-1}$ and ${\mu }_a=0.45{\mathrm{cm}}^{-1}$ only at the wavelength of 630 nm; and therefore, our accuracy comparison is best made to the nearest VCSEL wavelength of 680 nm.

Fig. 4

Insertion gain simulation results for 100 and 1000 MHz. The scale is from 0 to 30 dB, depending on the loss in phantom.

Fig. 5

Insertion phase simulation results for 100 and 1000 MHz. The scale is from 0 deg to 75 deg, depending on the phase shift.

3.

Blind Extraction of Optical Parameters and Comparison of Results

This section addresses the experimental and blind extraction results by simulating the analytical solution using DE and numerical predictions using FEM for broadband frequency modulation of diffused photon density waves (DPDWs).

Comparison of analytical, numerical, and experimental results of IL and IP are depicted in Figs. 6 and 7 at 680 nm for all four solid phantoms covering a frequency range of 30 to 1000 MHz. The least-square-method error fitting of DE using MATLAB and FEM results using COMSOL is used, along with differential measurements, in order to perform blind extraction of the phantom optical parameters. This IL and IP loss provides information related to absorption and scattering coefficients of the four phantoms. Initial measurements were done to give a sense of blind extraction, through monitoring the DC level in the media to perform extraction of absorption coefficient using the Beer-Lambert law ${\mu }_a=\frac{A}{l}$, where $A$ is the DC transmission ratio and $l$ is the path length of the photons, as shown in Table 2. This result is obtained by comparing the DC level obtained from measured new phantoms to that of the known optical parameter phantom. The absorption level depends on the amount of India ink²⁹ used in each phantom. Therefore, the darker the phantom, the higher IL is expected.

Fig. 6

(a) Phantom 1 IL, (b) phantom 1 IP, (c) phantom 2 IL, and (d) phantom 2 IP.

Fig. 7

(a) Phantom 3 IL, (b) phantom 3 IP, (c) phantom 4 IL, and (d) phantom 4 IP.

Table 2

Beer-Lambert extraction.

The results, shown in Table 3, render errors between the manufacturer data at 630 nm with respect to the extracted values using broadband and narrow-band measurements of IL and IP at 680 nm. Error analyses of the extracted results are presented based on the differences between the known data at 630 nm from the manufacturer and the best-fit values of absorption and scattering coefficients at 680 nm using the measurement bandwidths. Since all experimental results were collected at the fNIR optical Tx wavelength of 680 nm, the reduced scattering coefficient of a higher error could be due to this wavelength mismatch, not necessarily due to measurement or extraction procedure errors. Note that in all the phantoms, an IP shift of approximately 19 deg is observed in the optical Rx positions of $d_1=0.5\mathrm{cm}$ and $d_2=0.8\mathrm{cm}$, respectively (see Fig. 5). Observation of this constant phase indicates that all four phantoms have similar scattering parameters, ${\mu }_s^{\prime }$. These phantoms represent a homogenous solution of scattering and absorption that leads to unique IL and IP behavior in the frequency domain.

Table 3

Optical parameter extraction for Phantoms 1, 2, 3, and 4, with broadband and narrow-band extraction with error analysis.

Measurements of IL over broadband (depicted in Figs. 6 and 7) show that phantom 1 (with a low optical absorption parameter) has a low IL compared to phantom 4 (a high-absorption phantom). It is important to note that phantom 4 has the greatest amount of ink, which leads to high absorption and high IL. The measured IL response depicts a poor performance at around 600 MHz for phantom 4. This could be due to hardware limitations of a higher thermal noise and poor optical receiver responsivity at frequencies above 600 MHz and potential bandwidth limitation of APD. In other words, phantoms with ${\mu }_a>3.01{\mathrm{cm}}^{-1}$ would require a more sensitive receiver than the current APD-based optical Rx module.^26,27 As rendered in Table 4, the extracted optical parameters using DE in MATLAB depict measurements over a broadband frequency modulation with a small amount of error.

Table 4

Optical parameter extraction summary for all four phantoms at 630 nm.

4.

Discussion and Conclusions

FEM-based numerical modeling of brain phantoms is reported with important significance associated with our ability to analyze and predict photon scattering. The blind extraction of optical parameters is compared to the extraction results using DE for the homogenous medium, where a closed-form analytical solution with infinite spatial resolution is available. A similar extracted optical parameter value using DE- and FEM-based modeling with limited spatial resolution of meshing in homogenous solid phantoms placed in air justifies the extension of FEM modeling to inhomogeneous and multilayer phantoms. Our novel signal processing of the raw signal was done by curve fitting to enhance the free space measurements and a two-separation subtraction to eliminate the common sources of error. The combination of curve fitting and the two-separation subtraction technique was applied to the experimental measurements and both the DE analytical and FEM numerical simulations. The advantage of both techniques is that they extract the accurate optical parameters from unknown phantoms, a process that is called blind extraction.

Accurate extraction of the optical parameters depends on the sensitivity of the modeling method used to explain the photon behavior in the media. Broadband and narrow-band optical parameters extraction (shown in Table 3) is for the homogenous media (phantoms with a mix of absorption and scattering materials). The analysis and the experimental results show that extraction based on broadband frequency modulation would give an accurate assessment of the optical parameters of the homogenous media. Results for Phantom 1 show a lowest error of 2% and a highest error of 8% for absorption and scattering coefficients, respectively. The extraction conducted and reported in this paper at NIR wavelength of 680 nm is primarily meant to demonstrate compatibility of DE- and FEM-based modeling rather than parameter extraction of all NIR wavelengths of interest; note that parameter extraction for NIR tri-wavelengths of interest has been reported for a solid brain phantom, resembling white brain matter.¹⁶

The novel aspect of this paper is the comparison of the extracted results with the manufactured data for the purpose of validating the numerical simulation and modeling using FEM (i.e., COMSOL). Table 4 summarizes the extracted optical parameters using analytical-based DE over broadband frequency measurements of 30 to 1000 MHz. The optical parameters extracted using DE are then used through FEM using COMSOL to simulate IL and IP for different narrow-band and broadband frequency modulation. Then error analysis between analytical-based DE and numerical FEM are computed based on standard error analysis between two curves in different narrow-band and broadband frequency ranges, as summarized in Table 5. This small error for broadband measurements provides a high level of confidence in FEM numerical simulation, even for coarse spatial resolution, compared to infinite resolution of the DE solution. The low percentage of errors show that FEM-based modeling in 2-D with our selected meshing has a strong utility as a replacement for DE-based analytical modeling.

Table 5

Error percentage between FEM using COMSOL and analytical DE using MATLAB.

The advantage of the FEM numerical extraction is shown by the computational time required for simulating the modulated NIR photons, in which $T_{\mathrm{simulations}}<45$ s. The confidence established in accurate and robust 2-D simulation can now be extended to 3-D homogenous and inhomogeneous and multilayer phantoms. Our ability to demonstrate a 3-D multilayer numerical simulation and its blind extraction is being reported elsewhere.