15 January 2013 Accurate optical parameter extraction procedure for broadband near-infrared spectroscopy of brain matter
Author Affiliations +
J. of Biomedical Optics, 18(1), 017008 (2013). doi:10.1117/1.JBO.18.1.017008
Abstract
Modeling behavior of broadband (30 to 1000 MHz) frequency modulated near-infrared (NIR) photons through a phantom is the basis for accurate extraction of optical absorption and scattering parameters of biological turbid media. Photon dynamics in a phantom are predicted using both analytical and numerical simulation and are related to the measured insertion loss (IL) and insertion phase (IP) for a given geometry based on phantom optical parameters. Accuracy of the extracted optical parameters using finite element method (FEM) simulation is compared to baseline analytical calculations from the diffusion equation (DE) for homogenous brain phantoms. NIR spectroscopy is performed using custom-designed, broadband, free-space optical transmitter (Tx) and receiver (Rx) modules that are developed for photon migration at wavelengths of 680, 780, and 820 nm. Differential detection between two optical Rx locations separated by 0.3 cm is employed to eliminate systemic artifacts associated with interfaces of the optical Tx and Rx with the phantoms. Optical parameter extraction is achieved for four solid phantom samples using the least-square-error method in MATLAB (for DE) and COMSOL (for FEM) simulation by fitting data to measured results over broadband and narrowband frequency modulation. Confidence in numerical modeling of the photonic behavior using FEM has been established here by comparing the transmission mode’s experimental results with the predictions made by DE and FEM for known commercial solid brain phantoms.
Sultan, Najafizadeh, Gandjbakhche, Pourrezaei, and Daryoush: Accurate optical parameter extraction procedure for broadband near-infrared spectroscopy of brain matter

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.1 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 path2 is allocated at each modulation frequency. Various fNIR system topologies are reported based on the continuous wave (CW),3 time domain (TD),4 and frequency domain (FD).5 The CW is based on unmodulated photons, and the optical absorption coefficient6 can be extracted5 only by monitoring only the photon amplitude. TD systems depend on tunable mode-locked lasers to generate pulses with different time samples,4 from which both absorption and scattering coefficients7 are efficiently extracted by monitoring amplitude changes as a function of time due to propagation delay of photons in a scattered media.4 FD is based on either narrow-band (e.g., single modulating frequency I/Q system) or broadband frequency-swept diffused photon waves,3 and both absorption and scattering coefficients89.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.1112.13 Extraction of the optical parameters of absorption (μa) and reduced scattering (μs) 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 medium14 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.15

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.

JBO_18_1_017008_f001.png

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 reported16 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.

1.1.

Inverse Problem Solving Using DE and MC

The simplified DE is used to model frequency-modulated photon migration uniform homogenous media by extending the Boltzman transport theory.17 It has been demonstrated that the DE could predict a behavior as accurately as the MC technique for modulation frequencies up to 1 GHz.18 The extracted optical parameters are used to describe the brain activity at a certain location due to oxygenated and deoxygenated hemoglobin concentrations.3,6,89.10 The DE has been used in our analysis for broadband frequency modulation analysis and extraction of the optical parameters. The derived DE states:

(1)

1cϕ(r,t)tD2ϕ(r,t)(μa)ϕ(r,t)=S(r,t),
where ϕ(r,t) is the fluence rate (w/mm2) for reflectance or transmittance mode, and c is the speed of light in the media. The photon density wave (PDW) is expressed for a sinusoidal point source modulated at an angular frequency of ω=2πf in a semi-infinite medium as:11,19

(2)

ϕ(r,t)=ϕDC(r,t)+ϕAC(r,t)ej(krωt)

(3)

ϕ(r,t)=Adc4πDe(rδ)r+Aac4πDe(kimagr)re[j(krealrωt)]

(4)

Kreal=32μaμs{[1+(ωcμa)2]12+1}12

(5)

Kimag=32μaμs{[1+(ωcμa)2]121}12,
where Adc and Aac are, respectively, the DC (unmodulated) and AC (modulated) components of the source; δ is the penetration depth for unmodulated light; and kreal and kimag are the real and imaginary components of the complex wave number that describes spatial evolution for the PDW. Parameter D is the diffusion coefficient D=13·μs, and r is the separation distance between transmitter and receiver.

The DE is attractive because of its simplicity and because it has a closed-form solution in a homogenous medium. Since the DE assumes that the radiance is weakly anisotropic,20 then an anisotropic factor is introduced to provide a better explanation for the photon migration. Even though the anisotropic factor, g, was studied in different homogenous biological media,21 the complexity of the inverse problem solving for inhomogeneous media is time consuming, and accuracy of the extracted optical parameters is still under investigation.9,10,22 Therefore, a variety of approximations, such as modeling the photon behavior through MC simulation, has been reported.23 The basic premise of MC simulation is that complex interactions of particles and biological matter can be treated as a stochastic process, with simulated random movement samples from probability density functions. Even though the MC simulation of two layers15 shows a good result compared to two-layer DE,24 it requires a great amount of computational time.19

1.2.

Novel Inverse Problem Solving Using FEM Simulation

Finite element analysis is a numerical method used to analyze the partial differential equation solutions in discrete points with appropriate material boundaries. This FEM-based numerical technique could be applied for the first time to a homogenous, single-layer medium and then extended to inhomogeneous multilayer media with appropriate boundary conditions. The FEM numerical solutions are based on the discretized DE and are applied to multilayer problems where a closed-form solution does not exist, and MC is computationally time consuming even when applied to two-layer media.25 Commercial multiphysics-based FEM software, such as COMSOL, is used as a robust numerical tool to predict the photon behavior in highly scattering media and to solve the inverse problem. The Helmholtz equation is used for FEM numerical analysis25 in COMSOL as

(6)

[Dϕ(r,t)]+μaϕ(r,t)=S(r,t),
where ϕ(r,t) is the photon flux,D=13·μs , and in a frequency-modulated NIR system, S(r,t)=So(1+mejωt), where m is the amplitude modulation index at a time harmonic modulation frequency of ω.

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.

μa(cm−1)μs′(cm−1)
Phantom 10.4519.9
Phantom 20.6719.9
Phantom 31.0819.9
Phantom 43.0119.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.

JBO_18_1_017008_f002.png

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 40dBm 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 μa and μs 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 afreq+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 method16 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

(7)

ΔIL(dB)=ILd2=0.8cmILd1=0.5cm.

(8)

ΔIP(Deg)=IPd2=0.8cmIPd1=0.5cm.

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.

JBO_18_1_017008_f003.png

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.28 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

(9)

n·[(cu)1][(cu)2]+q·u=g

(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 μs=19cm1 and μa=0.45cm1 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.

JBO_18_1_017008_f004.png

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.

JBO_18_1_017008_f005.png

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 μa=Al, 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 ink29 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.

JBO_18_1_017008_f006.png

Fig. 7

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

JBO_18_1_017008_f007.png

Table 2

Beer-Lambert extraction.

μa(cm−1)
Phantom 10.44
Phantom 20.66
Phantom 31.05
Phantom 42.96

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 d1=0.5cm and d2=0.8cm, respectively (see Fig. 5). Observation of this constant phase indicates that all four phantoms have similar scattering parameters, μs. 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.

Phantom 1Phantom 2
Extraction using DE @680 nmError %Extraction using DE @680 nmError %
μa(cm−1)μs′(cm−1)μa(cm−1)μs′(cm−1)μa(cm−1)μs′(cm−1)μa(cm−1)μs′(cm−1)
30–1000 MHz0.4418.4280.6518.4837
30–300 MHz1.347.30>100631.1419.79701
100–500 MHz0.915.78>100210.6420.0051
500–1000 MHz0.3215.6029220.4413.773431
30–500 MHz>1001.66>10092>1000>100100
100–1000 MHz0.6419.444220.719.2853
Phantom 3Phantom 4
Extraction using DE @680 nmError %Extraction using DE @680 nmError %
μa(cm−1)μs′(cm−1)μa(cm−1)μs′(cm−1)μa(cm−1)μs′(cm−1)μa(cm−1)μs′(cm−1)
30–1000 MHz1.0218.73662.9318.7736
30–300 MHz5.16.59>10067>1000.19>10099
100–500 MHz1.920.467631.452.2053>100
500–1000 MHz1.1329.735502.85314.94525
30–500 MHz>1000.16>10099>1000.19>10099
100–1000 MHz2.1419.659810.3413.628932

The bold values are the most accurate values when compared to the manufacture data.

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 μa>3.01cm1 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.

Manufacturer’s data at 630 nmDE using MATLAB at 680 nm
μa(cm−1)μs′(cm−1)μa(cm−1)μs′(cm−1)
Phantom 10.4519.80.4418.4
Phantom 20.67019.80.6518.48
Phantom 31.0819.81.0218.73
Phantom 43.0119.82.9318.77

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.16

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.

Phantom 1Phantom 2Phantom 2Phantom 3
Error %Error %Error %Error %
ILIPILIPILIPILIP
30 to 1000 MHz2.681.232.811.081.781.412.711.17
30 to 300 MHz2.951.183.712.012.712.622.832.59
100 to 500 MHz3.11.312.932.053.661.853.671.84
500 to 1000 MHz3.091.513.821.722.782.713.862.13
30 to 500 MHz2.961.382.981.632.811.683.12.26
100 to 1000 MHz2.891.142.871.191.911.562.781.21

The advantage of the FEM numerical extraction is shown by the computational time required for simulating the modulated NIR photons, in which Tsimulations<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.

Acknowledgments

This research is partly supported by the Center for Neuroscience and Regenerative Medicine (CNRM), and the Intramural Research Program (IRP) of Eunice Shriver National Institute of Child Health and Human Development (NICHD) of the National Institute of Health. We are grateful to Professor Jessica Rammella-Roman, at the Catholic University of America, for her support in providing the phantoms.

References

1. T. J. FarrellM. S. PattersonB. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992).MPHYA60094-2405 http://dx.doi.org/10.1118/1.596777 Google Scholar

2. L. T. Perelmanet al., “Time-dependent photon migration using path integrals,” Phys. Rev. E 51(6), 6134–6141 (1995).PLEEE81063-651X http://dx.doi.org/10.1103/PhysRevE.51.6134 Google Scholar

3. K. Izzetogluet al., “Functional optical brain imaging using near-infrared during cognitive tasks,” Int. J. Hum. Comput. Interact. 17(2), 211–227 (2004).1071-5819 http://dx.doi.org/10.1207/s15327590ijhc1702_6 Google Scholar

4. D. Grosenicket al., “Time-domain optical mammography: initial clinical results on detection and characterization of breast tumors,” Appl. Opt. 42(16), 3170–3186 (2003).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.42.003170 Google Scholar

5. T. H. Phamet al., “Broad bandwidth frequency domain instrument for quantitative tissue optical spectroscopy,” Rev. Sci. Instrum. 71(6), 2500–2513 (2000).RSINAK0034-6748 http://dx.doi.org/10.1063/1.1150665 Google Scholar

6. H. Ayazet al., “Optical brain monitoring for operator training and mental workload assessment,” Neuroimage 59(1), 36–47 (2011) 10.1016/j.neuroimage.2011.06.023.NEIMEF1053-8119 Google Scholar

7. J. C. Hebdenet al., “Three-dimensional time-resolved optical tomography of a conical breast phantom,” Appl. Opt. 40(19), 3278–3287 (2001).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.40.003278 Google Scholar

8. J. Choiet al., “Noninvasive determination of the optical properties of adult brain: near-infrared spectroscopy approach,” J. Biomed. Opt. 9(1), 221–229 (2004).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1628242 Google Scholar

9. F. Bevilacquaet al., “In vivo local determination of tissue optical properties: applications to human brain,” Appl Opt. 38(22), 4939–4950 (1999).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.38.004939 Google Scholar

10. J. Zhaoet al., “In vivo determination of the optical properties of infant brain using frequency-domain near-infrared spectroscopy,” J. Biomed. Opt. 10(2), 024028 (2005).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1891345 Google Scholar

11. M. S. PattersonB. ChangeB. C. Wilson, “Time-resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties.” Appl. Opt. 28(12), 2331–2336 (1989).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.28.002331 Google Scholar

12. F. Bevilacquaet al., “Broadband absorption spectroscopy in turbid media by combined frequency-domain and steady state methods,” Appl. Opt. 39(34), 6498–6507 (2000).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.39.006498 Google Scholar

13. A. SassaroliF. MartelliS. Fantini, “General perturbative approach to the diffusion equation in the presence of absorbing defects: frequency-domain and time-domain results,” Proc. SPIE 7896, 78961N (2011).PSISDG0277-786X http://dx.doi.org/10.1117/12.876943 Google Scholar

14. M. S. Pattersonet al., “Diffusion equation representation of photon migration in tissue,” in IEEE MTT-S Int. Microwave Symposium Digest, Vol. 2, pp. 905–908, IEEE, Boston, Massachusetts (1991). Google Scholar

15. A. Kienleet al., “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37(4), 779–791 (1998).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.37.000779 Google Scholar

16. E. Sultanet al., “Modeling and tissue parameter extraction challenges for free space broadband fNIR brain imaging systems,” Proc. SPIE 7902, 790223 (2011).PSISDG0277-786X http://dx.doi.org/10.1117/12.875618 Google Scholar

17. J. J. DuderstadtW. R. Martin, Transport Theory, John Wiley & Sons, New York (1979). Google Scholar

18. J. B. Fishkinet al., “Gigahertz photon density waves in a turbid medium: theory and experiments,” Phys. Rev. E 53(3), 2307–2319 (1996). http://dx.doi.org/10.1103/PhysRevE.53.2307 Google Scholar

19. D. ContiniF. MartelliG. Zaccanti, ‘‘Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,’’ Appl. Opt. 36(19), 4587–4599 (1997).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.36.004587 Google Scholar

20. S. T. FlockB. C. WilsonM. S. Patterson, “Monte Carlo modeling of light propagation in highly scattering tissues. II. Comparison with measurements in phantoms,” IEEE Trans. Biomed. Eng. 36(12), 1169–1173 (1989).IEBEAX0018-9294 http://dx.doi.org/10.1109/10.42107 Google Scholar

21. J. Heinoet al., “Anisotropic effects in highly scattering media,” Phys. Rev. E 68(3), 031908 (2003).PLEEE81063-651X http://dx.doi.org/10.1103/PhysRevE.68.031908 Google Scholar

22. C. ZhuQ. Liu, “Validity of the semi-infinite tumor model in diffuse reflectance spectroscopy for epithelial cancer diagnosis: a Monte Carlo study,” Opt. Express 19(18), 17799–17812 (2011).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.19.017799 Google Scholar

23. A. Sassaroli, “Fast perturbation Monte Carlo method for photon migration in heterogeneous turbid media,” Opt. Lett. 36(11), 2095–2097 (2011).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.36.002095 Google Scholar

24. K.-W. Guanet al., “A two-layer model of laser interaction with skin: a photothermal effect analysis,” Opt. Laser Technol. 43(3), 425–429 (2011).OLTCAS0030-3992 http://dx.doi.org/10.1016/j.optlastec.2009.12.007 Google Scholar

25. M. Schweigeret al., “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).MPHYA60094-2405 http://dx.doi.org/10.1118/1.597634 Google Scholar

26. K. Mansetaet al., “Development challenges of brain functional monitoring using untethered broadband frequency modulated fNIR system,” in IEEE Topical Meeting on Microwave Photonics (MWP), pp. 354–357, IEEE, Montreal, Quebec (2010). Google Scholar

27. K. Mansetaet al., “Untethered helmet mounted functional near infrared (fNIR) biomedical imaging?,” in IEEE MTT-S Int. Microwave Symposium Digest (MTT), pp. 1–4, IEEE, Baltimore, Maryland (2011). Google Scholar

28. J. ZhouJ. Bai, “The most favorable path method for the propagation of light in scattering media: a finite element solution,” in Proc. 23rd Ann. Int. Conf. of the IEEE Eng. Med. Biol. Soc., Vol. 4, pp. 3203–3206, IEEE, Beijing, China (2001). Google Scholar

29. P. Di NinniF. MartelliG. Zaccanti, “The use of India ink in tissue-simulating phantoms,” Opt. Express 18(26), 26854–26865 (2010).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.18.026854 Google Scholar

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Ebraheem Sultan, Laleh Najafizadeh, Amir Gandjbakhche, Kambiz Pourrezaei, Afshin S. Daryoush, "Accurate optical parameter extraction procedure for broadband near-infrared spectroscopy of brain matter," Journal of Biomedical Optics 18(1), 017008 (15 January 2013). http://dx.doi.org/10.1117/1.JBO.18.1.017008
JOURNAL ARTICLE
8 PAGES


SHARE
KEYWORDS
Photons

Finite element methods

Brain

Inverse optics

Absorption

Modulation

Optics manufacturing

Back to Top