INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 Published online 20 November 2013 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/cnm.2606
3D Numerical modeling and its experimental verifications for an inhomogeneous head phantom using broadband fNIR system E. Sultan 1 , K. Pourrezaei 2 , A. Ghandjbakhche 3 and A.S. Daryoush 1, * ,† 1 Department
of Electrical and Computer Engineering, Drexel University, Philadelphia, PA 19104, U.S.A. of Biomedical Engineering, Science, and Health Systems, Drexel University, Philadelphia, PA 19104, U.S.A. 3 Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, Bethesda, MD 20892, U.S.A.
2 School
SUMMARY Modeling behavior of broadband (30–1000 MHz) frequency modulated near infrared photons through a multilayer phantom is of interest to optical bio-imaging research. Photon dynamics in phantom are predicted using three-dimension (3D) finite element numerical simulation and are related to the measured insertion loss and phase for a given human head geometry in this paper based on three layers of phantom each with distinct optical parameter properties. Simulation and experimental results are achieved for single, two, and three layers solid phantoms using COMSOL (COMSOL AB, Tegnérgatan 23, SE-111 40, Stockholm, Sweden) (for FEM) simulation and custom-designed broadband free space optical transmitter (Tx) and receiver (Rx) modules that are developed for photon migration at wavelengths of 680, 795, and 850 nm. Standard error is used to compute error between two-dimension and p 3D FE modeling along with experimental results by fitting experimental data to the functional form of a f reque ncy C b. Error results are shown at narrowband and broadband frequency modulation. Confidence in numerical modeling of the photonic behavior using 3D FEM for human head has been established here by comparing the reflection mode’s experimental results with the predictions made by COMSOL for known commercial solid brain phantoms. Copyright © 2013 John Wiley & Sons, Ltd. Received 3 May 2013; Accepted 23 September 2013 KEY WORDS:
fNIR; TBI; FEM; DE; COMSOL; optical transmitter Tx; optical receiver Rx; triwavelength; CW; TD; FD; PDW; VCSEL; IL; IP
1. INTRODUCTION The future of non-invasive medical imaging depends on many aspects such as accurate diagnosis, cost, mobility , and transportability. Imaging using near infrared (NIR) photons reveals accurate assessments of functionality, which can be achieved by accurate extraction of optical parameters [1, 2]. NIR imaging uses the electromagnetic spectrum of the NIR (600nm–2500 nm) as a base of spectroscopic modality [3]. The general form of this spectroscopy depends on photons that travels through biological media and interact with different particles of human cells and organs. These photons will either be absorbed or scattered while traveling in the most favorable path through the media [4]. Accurate diagnosis will depend on knowing the amount of photons entering the media and the amount leaving the media, and amount of phase shift that occurs due to multi-scattering [5, 6]. Photon behavior in biological high-scattering media is of great interest to researches in the field of bio-optics imaging, where the relationship between NIR photons and biological media reveals information of absorption or scattering that is related to hemodynamic process [7]. Modeling the *Correspondence to: A.S. Daryoush, Department of Electrical and Computer Engineering, Drexel University, Philadelphia, PA 19104, U.S.A. † E-mail: [email protected] Copyright © 2013 John Wiley & Sons, Ltd.
354
E. SULTAN ET AL.
Figure 1. Untethered frequency modulated photon migration in brain using helmet mounted structure; a) block diagram; b) differential detection; and c) multi optical Tx and Rx mounted on human head.
hemodynamic process through either absorption or scattering of photons in the biological media using either analytical or numerical methods helps in accurate medical diagnosis. In previous publication [1, 8], we have introduced the standard diffusion equation (DE) as an analytical modeling tool to model the frequency modulated photons in homogenous highly scattering media and extend the solution through FEM for two-dimension (2D) finite element modeling. One of the main purposes of this paper is to accurately model the modulated photon behavior when traveling in multilayer media and inhomogeneous complexity. Confidence in the FEM modeling of the homogenous media has been established [8] by comparing the analytical DE to the 2D FEM model simulation along with experimental result in homogenous phantoms. More complex biological media such as the human head is considered important inhomogeneous multilayers that contain one of the most important organs in the human body, the brain. The human brain is considered the control element of all physiological activities, and it happens to be a very fragile and sensitive organ when compared to other organs, and its functionality depends on knowing the amount of oxygenated and deoxygenated blood flows circulating in order to maintain normal neuro-activities [7, 9]. This paper will examine the modulated photon traveling in multilayer inhomogeneous brain phantom that resembles the human head. A conceptual representation of human head and propagation of diffused density waves are depicted in Figure 1(c), where optical transmitter (Tx) and receivers (Rx) are strategically located on head. This study is part of early detection of traumatic brain injury using a helmet mounted functional NIR (fNIR) device, as shown in Figure 1(a). Both wireless connectivity and differential detection of insertion loss (IL) and phase through brain are conducted for a helmet mounted broadband radio frequency (RF) electronics. A three-dimension (3D) FEM modeling is key in a high spatial resolution modeling the modulated photon migration through the inhomogeneous human head. Therefore, this paper will begin with demonstrating confidence in accuracy of the FEM-based modeling by identifying practical mesh sizes as the 2D FEM modeling is being extended to 3D. In particular, accuracy of optical parameter extraction is quantified at the shortest and the most challenging wavelength of 680 nm. Knowing a practical mesh size and the achieved FEM modeling accuracy, then we present accuracy of the two and three layers of phantom stacked to resemble human head by making comparison between the experimental measurement and numerical modeling results at three different wavelengths of 680, 795, and 850 nm. Moreover, the accuracy of 3D modeling to measurements results are compared over various broadband frequency bandwidths. These results will set understanding of how the broadband measurements could be applied for optical parameter extraction of an inhomogeneous brain medium without resorting to analytically challenging methods [10]. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
355
Figure 2. Frequency modulated optical system used on biological phantom.
2. NUMERICAL MODELING AND EXPERIMENTAL METHODS 2.1. Frequency domain fNIR system Functional imaging based on optical photon absorption and scattering related to hemodynamic changes in the human body is developed for different system topologies such as continuous wave, time domain (TD), and frequency domain (FD) [6]. Continuous wave is based on unmodulated photons, while TD systems depend on propagation loss and delay of tunable short NIR pulses. FD is based on frequency swept diffused photon density waves, DPDW, as both optical absorption and scattering parameters of biological tissue are extracted by monitoring amplitude and phase change of modulated photons [11]. Many techniques have been demonstrated in TD to obtain quantitative optical properties such as absorption (a / and scattering (s / coefficients [12]. While the TD spectroscopy instrument has shown good sensitivity, the FD method can be achieved more economically and is more suitable for real time monitoring in clinical settings [13]. The FD method, in which the light source intensity is modulated at either single or broadband frequency, has been applied as an applicable method to study photon migration in a multilayer biological media and therefore monitor and extract the optical properties. As light travels outward from the photon source, optical transmitter (Tx), in a highly scattering medium, the flux energy density decreases more than exponentially with increasing distance from the source. In FD instrument, the intensity of the light source is sinusoidally modulated as shown in Figure 2. The photon density wave can be said to propagate through the biological phantom outward from the source. In a frequency modulated system, the source point is modulated and can be described by the equation S.r, t / D SAC 1 C me j!t (1) where SAC is the amplitude of unmodulated light, and m is optical modulation index for RF modulation frequency of !. 2.2. Finite element method-based numerical modeling The DE is widely used for analytical analysis to explain photon migration dynamics in homogenous medium, where the anisotropic factor, g, is introduced for radiance in a weakly anisotropic medium. Even though g was studied in different homogenous biological media, accuracy of optical parameter extraction and the complexity of the inverse problem solution for inhomogeneous media Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
356
E. SULTAN ET AL.
are still under investigation [14]. A variety of approximation, such as Monte Carlo simulation of the photon behavior, has been extended to two layer diffusion, but it would require high amount of computational power [15]. The analytical model of photon transport in tissue using FEM has been first introduced by Arridge et al. in 1993 [16] to demonstrate the ability of FEM to model the photon density inside an object with photon flux along the geometric boundaries. The FEM modeling manifestation is based on the assumption that photon transport is the diffusion approximation to a radiative transfer equation. The advantage of the FEM modeling is its high accuracy for a practical convergence speed than any other numerical modeling. Other researchers have applied FEM to construct methods of imaging the photon distribution and extract the optical parameters of bio-media. Commercial multiphysics-based FEM software, such as COMSOL, is used in this paper as a robust numerical tool to numerically predict the photon behavior in a highly scattering media, as we have demonstrated in 2D modeling [8]. The Helmholtz equation is used for FEM numerical analysis [16, 17] in COMSOL as r.Dr;.r, t // C a ;.r, t / D S.r, t /
(2)
where ;.r, t / is the photon flux, D D 1 0 , and in a frequency modulated NIR system S.r, t / D 3s SO 1 C me j!t where m is amplitude modulation index at time harmonic modulation frequency of !. The main objective of this paper is to validate the use of a commercially available tool (COMSOL) in 3D to help the process of solving the inverse problem and to extract optical parameters of the homogenous and inhomogeneous media, such as absorption and reduced scattering coefficients for broadband frequency modulated photons (from 30 to 1000 MHz). Accuracy of FEM simulation results depends on resolution of mesh structures in comparison to wavelength of light and modulating RF frequencies and any associated points in space, where DPDW solution is to be computed using partial difference equation of Helmholtz equation. Higher resolution meshing structure is required to accurately simulate the photon dynamics, which naturally results in longer computation times or even cause computation systems to crash; therefore, a compromise in accuracy and mesh resolution has to be made in order to find the optimal simulated environment. It is our experience that having high resolution around optical source and receiver would predict an accurate result for the specific amplitude and phase change while meeting a typical computation time of 4 s for single 2D simulation and 26 s for single 3D simulation using an Intel Core i3 processor PC with a reasonable error of less than 5%—a mesh element sizes of 80 and 1000 nm are considered in regions around optical Tx and Rx and everywhere else, respectively. Different mesh topologies have been reviewed at the wavelength of 680 nm for high meshing of 80 nm selection around the boundaries between layers and the result shown to be the same as if selecting the same meshing around only the Tx and Rx. The difference between these two topologies is the computational time, where applying 80 nm around boundaries takes 94.6 s for single simulation. These meshing element sizes are also recommended when inhomogeneous phantoms resembling head tissue is considered [18]. In COMSOL, the boundary condition is set either through Dirichlet or Neumann boundary conditions, where they are set to express the fluence rate, u, at desired boundaries. Dirichlet and Neumann boundary conditions are chosen appropriately for air-dielectric (or dielectric–dielectric) and radiation condition in region surrounding phantom. The mathematical representation is expressed as n ..cru/1 ..cru/2 / C q u D g huDr
(3) (4)
where h, r, q, and g are the boundary coefficients for the phantom modeling. For the diffused photon flux, h D 1 and r D 0. For the matching boundary of air dielectric and dielectric and dielectric interfaces q D 0 and g D 0. It is important to mention that extrapolated boundary conditions take the phantom-air boundary into account as a radiation boundary by setting the fluence rate to zero using Eq. 3 without any impact on error achieved. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
357
Table I. Manufacture provided optical parameters for the phantoms at 680 nm. a cm1
0 s cm1
0.35 0.15 0.12
15 13 8
Cortex Scalp Skull
Figure 3. Optical system set up showing optical transmitter and receiver measuring insertion loss and insertion phase for three stacked layers of phantoms.
2.3. Experimental measurement procedures Accurate modeling requires experimental verification of change in amplitude and phase of the incident photon wave that is related to absorption and scattering of the biological media is presented in the FD as IL and insertion phase (IP), where IL is any change in amplitude in units of dB, and IP is any change in phase in units of degree. This measurement of IL and IP requires certain hardware and testing equipment that would perform high-frequency modulation [19, 20]. The overall system would consist of Automatic Network Analyzer (Anritsu MS4623B) as an RF source and sensitive RF receiver and broadband optical Tx and sensitive optical Rx. The ANA acts as a frequency modulator and measurements tools, while the optical Tx and Rx modules act as modulated photons source and detector. To a phantom resembling human head, three layers of phantoms were designed with different optical properties and physical geometry as shown in Table I [20, 21], and all three layers can be stacked on top of each other with different combination as shown in Figure 3.
3. EXPERIMENTAL AND SIMULATION RESULTS 3.1. 3D modeling of biological phantom A transition to the 3D is required in order to model the multilayer media with high-spatial resolution. The method used in this section is based on the comparison between results obtained in previous publication [8], and with new 3D numerical model simulation results. This section will focus on accuracy comparison of homogenous phantoms that are modeled in 2D and 3D model at the shortest NIR wavelength of 680 nm in the broadband frequency modulation of the photons. This modeling comparison for an optimum meshing structure will provide confidence in the 3D FEM modeling at a reasonable computation time using our current computational resources for 3D Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
358
E. SULTAN ET AL.
Figure 4. Insertion loss COMSOL simulation for a) two-dimension (2D) scalp and skull, b) three-dimension (3D) scalp and skull, c) 2D reflection mode cortex, d) 3D reflection mode cortex, e) 2D semi-transmission mode, and f) 3D semi-transmission mode.
Figure 5. Insertion phase COMSOL simulation for a) two-dimension (2D) scalp and skull, b) three-dimension (3D) scalp and skull, c) 2D reflection mode cortex, d) 3D reflection mode cortex, e) 2D semi-transmission mode, and f) 3D semi-transmission mode. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
359
frequency modulated photon traveling through multilayer inhomogeneous media. Three different modes of photon traveling are used to make the transition between 2D and 3D FEM modeling; these three modes of transmission, semi-transmission, and reflection can been used separately to extract the optical parameters of the phantoms [1, 8] using COMSOL. A broadband simulation of photon flux is computed for both cases of 2D and 3D, and comparison is made to the experimental results to quantify optical parameter extraction accuracies. Photon flux predictions for transmission mode is applied for scalp and skull phantom, while semi-reflection and reflection mode is applied to cortex phantom as shown in Figures 4 and 5. In the case of transmission and reflection mode, the transmitter is located at the center of the phantom, while for the semi-transmission mode, the optical Tx is located 1 cm away from the edge of the phantom. Location of the optical Rx can be chosen as desired on various positions of phantom, and for accurate optical parameter extraction without being hampered by optical Tx and Rx interface errors, two locations are considered for the differential detection procedures established earlier [1, 8]. The concept of the homogenous modeling depends on identifying an extrapolated boundary region layer and choosing an isotropic point of source limited to the media in sense of space and frequency. The source location is set to be very close to the extrapolated region with certain depth of one transport mean free path ltr beneath the surface of the phantom [18]. Comparison of the experimental results with the FEM numerical result is shown in Figure 6 for broadband frequency modulation. These result in Figure 6 shows that IL between low frequency to high frequency for transmission mode in scalp is about 6 dB and in skull is about 5 dB, while for reflection mode of cortex is about 35 dB and for semi-transmission mode of cortex is about 16 dB. IP between low frequency and
Figure 6. Broadband frequency modulation measurements and simulation at 680 nm for a) insertion loss (IL) scalp and skull, b) insertion phase (IP) scalp and skull, c) IL reflection mode cortex, d) IP reflection mode cortex, e) IL semi-transmission mode, and f) IP semi-transmission mode. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
360
E. SULTAN ET AL.
Table II. Error percentage between FEM using COMSOL and analytical diffusion equation using MATLAB (MathWorks, 3 Apple Hill Drive, Natick, MA 01760-2098, United States). Scalp transmission Error %
30–1000 MHz 30–300 MHz 100–500 MHz 500–1000 MHz 30–500 MHz 100–1000 MHz
Skull transmission Error %
Cortex semi-trans Error %
Cortex reflection Error %
IL
IP
IL
IP
IL
IP
IL
IP
2.32 2.13 3.18 3.56 2.83 2.94
1.92 1.86 2.86 3.04 2.78 2.17
2.55 2.49 3.22 3.16 2.82 2.83
2.12 2.05 2.92 3.44 2.65 2.42
2.62 2.57 3.15 3.57 3.24 2.98
2.87 2.79 3.54 3.68 3.18 3.05
2.8 2.67 3.47 3.76 3.28 3.11
2.51 2.38 3.28 3.42 3.15 2.78
IL, insertion loss; IP, insertion phase.
Table III. Error percentage between two-dimension FEM and three-dimension FEM using COMSOL . Scalp transmission Error %
30–1000 MHz 30–300 MHz 100–500 MHz 500–1000 MHz 30–500 MHz 100–1000 MHz
Skull transmission Error %
Cortex semi-trans Error %
Cortex reflection Error %
IL
IP
IL
IP
IL
IP
IL
IP
1.32 1.58 1.65 1.77 1.64 1.51
1.12 1.24 1.37 1.42 1.32 1.24
1.55 1.63 1.61 1.72 1.75 1.56
1.18 1.31 1.41 1.48 1.35 1.23
1.76 1.85 1.84 1.95 1.87 1.81
1.21 1.44 1.35 1.52 1.42 1.3
1.62 1.75 1.81 1.86 1.76 1.71
1.19 1.35 1.34 1.38 1.34 1.22
IL, insertion loss; IP, insertion phase.
Figure 7. Insertion loss COMSOL simulation for two layer of scalp/cortex or skull/cortex and three layers of scalp/skull/cortex.
high frequency for transmission mode in scalp is about 5ı and in skull is about 3ı , while for reflection mode of cortex is about 14ı and for semi-transmission mode of cortex is about 12ı . Accuracy through computing error between these result is shown in Tables II and III and is discussed in the last section of this paper. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
361
Figure 8. Insertion phase COMSOL simulation for two layer of scalp/cortex or skull/cortex and three layers of scalp/skull/cortex.
Figure 9. Broadband frequency modulation measurements and simulation for two and three layers of phantom of a) 680 nm insertion loss (IL) , b) 680 nm insertion phase (IP), c) 795 nm IL, d) 795 nm IP, e) 850 nm IL, and f) 850 nm IP. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
362
E. SULTAN ET AL.
3.2. 3D multilayers frequency modulated photon travel Because confidence has been established in 3D numerical modeling of the photon evolutions in terms of amplitude and phase losses as it travels through homogenous media resembling scalp, skull, and cortex, then the rest of this paper addresses on modeling the modulated photon evolution as it travels in multilayer media resembling human head. First, the FEM modeling results of two layers of scalp/cortex and skull/cortex are compared against experimental results and then extend for the first time to the challenging problem of multilayer (e.g., a three layer phantom) of scalp/skull/cortex. All of the photon flux analyses are based on three wavelength of 680, 795, and 850 nm as shown in Figures 7 and 8, while comparison between experimental and numerical FEM modeling result for broadband modulation is depicted in Figure 9. IL between two layer of skull/cortex and cortex is about 3–4 dB, two layers of scalp/cortex and cortex is about 5–6 dB, and between three layers and cortex is 7–8 dB. IP between two layer of skull/cortex and cortex is about 1–2ı , two layers of scalp/cortex and cortex is about 3–4ı , and between three layers and cortex is 5–6ı .
4. DISCUSSION AND CONCLUSIONS The confidence in 3D FE modeling depends first on the accuracy between the analytical DE and the numerical simulated result performed by COMSOL. The transition between 2D and 3D layers is performed by identifying the error between these dimensions. The error between 2D and 3D modeling analysis was based on calculating the standard error between the two broadband plots of IL and IP using narrowband and broadband frequencies as shown in Tables II and III. The error between the 2D and 3D FE numerical analysis is no more than 2%, which means the difference between the twodimensional models is insignificant. The next analysis deals with the two and three layers modeling and experimental results. Measurements of broadband frequency modulation along with the numerical simulation are consistent as shown in Figure 9, and it renders that two layers of scalp/cortex and skull/cortex have poor sensitivities as values of IL and IP at low frequencies are very close to one another, but at high frequencies, it shows gradual deviation emerging, which is expected because the shallow layers have different optical properties. The higher the optical properties of the top layer phantom are, the higher IL and IP are expected at higher frequencies. This behavior is shown when at three wavelengths, the IL and IP are higher at high frequencies for the scalp/cortex than the skull/cortex. The important part of this modeling and experimental result is the understanding of the modulated photon behavior in both two and three layers media in order to set the stage for inverse problem solving for medical applications. It is observed that at two layers of scalp/cortex and skull/cortex, the slope of the curve of IL and IP at low frequencies is very close, while at higher frequencies, it changes, and slope change could be used as a point of identification for high optical activity. The last analysis is done to compute the error between the experimental raw data and 3D COM SOL FEM simulation at different band of frequencies. This analysis is performed by curve fitting p the experimental data to a f reque ncy C b and then computing the standard error between curve fitted and 3D COMSOL simulation result as shown in Table IV. The transition between analytical DE and 2D finite element shows an error less than 4%, and error between 2D FEM to 3D FEM applied to homogenous media shows to be insignificant with less than 2%. The last error computed between experimental curve fitted data of two and three multilayers and 3D FEM shows to be less than 5%. This paper addressed two important aspects in the process of 3D modeling of FD functional imaging applied to human head. One is the transition between 2D and 3D FEM for frequency modulation photon traveling in homogenous media and two modeling the modulated photon behavior in inhomogeneous human head media using 3D FEM numerical method. These two important analyses are used to build up confidence in 3D numerical modeling and set the stage for future modeling the modulated photon behavior in inhomogeneous layers addressing different medical cases. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
Copyright © 2013 John Wiley & Sons, Ltd.
IL, insertion loss; IP, insertion phase.
850 nm
795 nm
2.42 3.89 4.51 4.32 2.39 3.72 4.48 4.27 2.5 3.83 4.63 4.38
680 nm
30–1000 MHz 30–300 MHz 300–500 MHz 500–1000 MHz 30–1000 MHz 30–300 MHz 300–500 MHz 500–1000 MHz 30–1000 MHz 30–300 MHz 300–500 MHz 500–1000 MHz
IL
Wavelength 2.25 2.87 3.15 3.08 2.18 2.78 3.12 3.10 2.15 2.75 3.1 3.07
IP
Cortex phantom Error %
2.52 4.21 4.75 4.52 2.43 4.18 4.59 4.44 2.52 4.27 4.67 4.53
IL 2.22 2.75 3.13 3.04 2.25 2.68 3.17 3.13 2.31 2.74 3.22 3.19
IP
Two layer scalp/cortex Error %
2.31 3.75 4.43 4.26 2.26 3.64 4.31 4.15 2.3 3.68 4.34 4.19
IL
2.19 2.81 3.12 3.05 2.16 2.72 3.08 3.04 2.19 2.75 3.13 3.07
IP
Two layer skull/cortex Error %
Table IV. Error percentage between 3D FEM and curve fitted experimental data.
2.39 4.0 4.53 4.41 2.32 3.94 4.26 4.32 2.28 3.9 4.23 4.28
IL
2.2 2.79 3.17 3.12 2.18 2.68 3.15 3.07 2.25 2.75 3.2 3.14
IP
Three layer scalp/skull/cortex Error %
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
363
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
364
E. SULTAN ET AL.
REFERENCES 1. Sultan E, Manseta K, Khwaja A, Najafizadeh L, Gandjbakhche A, Pourrezaei K, Daryoush AS. Modeling and tissue parameter extraction challenges for free space broadband fNIR brain imaging systems. Proceedings of SPIE, 7902, Imaging, Manipulation, and Analysis of Biomolecules, Cells, and Tissues IX, 790223, 2011, DOI: 10.1117/12.875618. 2. Zhao J, Ding HS, Hou XL, Zhou CL, Chance B. In vivo determination of the optical properties of infant brain using frequency-domain near-infrared spectroscopy. Journal of Biomedical Optics 2005; 10(2). 3. Ferrari M, Quaresima V. A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application. NeuroImage 2012; 63(2). 4. Contini D, Martelli F, Zaccanti G. Photon migration through a turbid slab described by a model based on diffusion approximation I. Theory. Applied Optics 1997; 36:4587–4599. 5. Farrell TJ, Patterson MS, Wilson B. A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo. Medical Physics 1992; 19(4):879–88. 6. Mansouri C, Kashou NH. New window on optical brain imaging; medical development, simulations and applications. Selected Topics on Optical Fiber Technology 2012. Chapter 11, pp. 271–288. ISBN: 978-953-51-0091-1. 7. Izzetoglu K, Bunce S, Onaral B, Pourrezaei K, Chance B. Functional optical brain imaging using near-infrared during cognitive tasks. International Journal of Human-Computer Interaction 2004; 17(2). 8. Sultan E, Najafizadeh L, Gandjbakhche AH, Pourrezaei K, Daryoush A. (2013). Accurate optical parameter extraction procedure for broadband near-infrared spectroscopy of brain matter. Journal of Biomedical Optics 2013; 18(1):017008–017008. 9. Hoshi Y. Functional near-infrared optical imaging: utility and limitations in human brain mapping. Psychophysiology 2003; 40:511–520. 10. Ripoll J, Ntziachristos V, Culver JP, Pattanayak DN, Yodh AG, Nieto-Vesperinas M. Recovery of optical parameters in multiple-layered diffusive media: theory and experiments. Journal of the Optical Society of America A-Optics Image Science And Vision 2001; 18(4):821–830. 11. Pham TH, Coquoz O, Fishikin JB, Anderson E, Tromberg BJ. Broad bandwidth frequency domain instrument for quantitative tissue optical spectroscopy. Review of Scientific Instruments 2000; 71:2500–2513. 12. Grosenick D, Moesta KT, Wabnitz H, Mucke J, Stroszczynski C, Macdonald R, Schlag PM, Rinneberg H. Timedomain optical mammography: initial clinical results on detection and characterization of breast tumors. Applied Optics 2003; 42:3170–3186. 13. Yang Y, Liu H, Li X, Chance B. Low-cost frequency-domain photon migration instrument for tissue spectroscopy, oximetry, and imaging. Optical Engineering 0001; 36(5):1562–1569. 14. Bevilacqua F, Piguet D, Marquet P, Gross JD, Tromberg BJ, Depeursinge C. In vivo local determination of tissue optical properties: applications to human brain. Applied Optics 1999; 38(22):4939–50. 15. Zhu C, Liu Q. Validity of the semi-infinite tumor model in diffuse reflectance spectroscopy for epithelial cancer diagnosis: a Monte Carlo study. Optics Express 2011; 19:17799–17812. 16. Arridge SR, Schweiger M, Hiraoka M, Delpy DT. A finite element approach for modeling photon transport in tissue. Medical Physics 1993; 20:299. 17. Schweiger M, Arridge SR, Hiraoka M, Delpy DT. The finite element method for the propagation of light in scattering media: boundary and source conditions. Medical Physics 1995; 22:1779–1792. 18. Zhou J, Bai J. The most favorable path method for the propagation of light in scattering media: a finite element solution, Engineering in Medicine and Biology Society, 2001. Proceedings of the 23rd Annual International Conference of the IEEE, vol.4, Istanbul, Turkey, 2001; 3203–3206. vol. 4. 19. Manseta K, Sultan E, Khwaja AM, Pourrezaei K, Joshi A, Najafizadeh L, Gandjbakhche A, Daryoush AS. Development challenges of brain functional monitoring using untethered broadband frequency modulated fNIR system. IEEE Microwave Photonics Proceeding, Montreal, Canada, MWP 2010; 295–299. 20. Manseta K, Khwaja AM, Sultan E, Daruwalla P, Pourrezaei K, Najafizadeh L, Gandjbakhche A, Daryoush AS. Untethered helmet mounted functional near infrared (fNIR) biomedical imaging? Microwave Symposium Digest (MTT), 2011 IEEE MTT-S International 2011:1, 4, 5–10. 21. Hoshi Y, Shimada M, Sato C, Iguchi Y. Reevaluation of near-infrared light propagation in the adult human head: implications for functional near-infrared spectroscopy. Journal of Biomedical Optics 0001; 10(6):064032–064032–10.
Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D Numerical modeling and its experimental verifications for an inhomogeneous head phantom using broadband fNIR system E. Sultan 1 , K. Pourrezaei 2 , A. Ghandjbakhche 3 and A.S. Daryoush 1, * ,† 1 Department
of Electrical and Computer Engineering, Drexel University, Philadelphia, PA 19104, U.S.A. of Biomedical Engineering, Science, and Health Systems, Drexel University, Philadelphia, PA 19104, U.S.A. 3 Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, Bethesda, MD 20892, U.S.A.
2 School
SUMMARY Modeling behavior of broadband (30–1000 MHz) frequency modulated near infrared photons through a multilayer phantom is of interest to optical bio-imaging research. Photon dynamics in phantom are predicted using three-dimension (3D) finite element numerical simulation and are related to the measured insertion loss and phase for a given human head geometry in this paper based on three layers of phantom each with distinct optical parameter properties. Simulation and experimental results are achieved for single, two, and three layers solid phantoms using COMSOL (COMSOL AB, Tegnérgatan 23, SE-111 40, Stockholm, Sweden) (for FEM) simulation and custom-designed broadband free space optical transmitter (Tx) and receiver (Rx) modules that are developed for photon migration at wavelengths of 680, 795, and 850 nm. Standard error is used to compute error between two-dimension and p 3D FE modeling along with experimental results by fitting experimental data to the functional form of a f reque ncy C b. Error results are shown at narrowband and broadband frequency modulation. Confidence in numerical modeling of the photonic behavior using 3D FEM for human head has been established here by comparing the reflection mode’s experimental results with the predictions made by COMSOL for known commercial solid brain phantoms. Copyright © 2013 John Wiley & Sons, Ltd. Received 3 May 2013; Accepted 23 September 2013 KEY WORDS:
fNIR; TBI; FEM; DE; COMSOL; optical transmitter Tx; optical receiver Rx; triwavelength; CW; TD; FD; PDW; VCSEL; IL; IP
1. INTRODUCTION The future of non-invasive medical imaging depends on many aspects such as accurate diagnosis, cost, mobility , and transportability. Imaging using near infrared (NIR) photons reveals accurate assessments of functionality, which can be achieved by accurate extraction of optical parameters [1, 2]. NIR imaging uses the electromagnetic spectrum of the NIR (600nm–2500 nm) as a base of spectroscopic modality [3]. The general form of this spectroscopy depends on photons that travels through biological media and interact with different particles of human cells and organs. These photons will either be absorbed or scattered while traveling in the most favorable path through the media [4]. Accurate diagnosis will depend on knowing the amount of photons entering the media and the amount leaving the media, and amount of phase shift that occurs due to multi-scattering [5, 6]. Photon behavior in biological high-scattering media is of great interest to researches in the field of bio-optics imaging, where the relationship between NIR photons and biological media reveals information of absorption or scattering that is related to hemodynamic process [7]. Modeling the *Correspondence to: A.S. Daryoush, Department of Electrical and Computer Engineering, Drexel University, Philadelphia, PA 19104, U.S.A. † E-mail: [email protected] Copyright © 2013 John Wiley & Sons, Ltd.
354
E. SULTAN ET AL.
Figure 1. Untethered frequency modulated photon migration in brain using helmet mounted structure; a) block diagram; b) differential detection; and c) multi optical Tx and Rx mounted on human head.
hemodynamic process through either absorption or scattering of photons in the biological media using either analytical or numerical methods helps in accurate medical diagnosis. In previous publication [1, 8], we have introduced the standard diffusion equation (DE) as an analytical modeling tool to model the frequency modulated photons in homogenous highly scattering media and extend the solution through FEM for two-dimension (2D) finite element modeling. One of the main purposes of this paper is to accurately model the modulated photon behavior when traveling in multilayer media and inhomogeneous complexity. Confidence in the FEM modeling of the homogenous media has been established [8] by comparing the analytical DE to the 2D FEM model simulation along with experimental result in homogenous phantoms. More complex biological media such as the human head is considered important inhomogeneous multilayers that contain one of the most important organs in the human body, the brain. The human brain is considered the control element of all physiological activities, and it happens to be a very fragile and sensitive organ when compared to other organs, and its functionality depends on knowing the amount of oxygenated and deoxygenated blood flows circulating in order to maintain normal neuro-activities [7, 9]. This paper will examine the modulated photon traveling in multilayer inhomogeneous brain phantom that resembles the human head. A conceptual representation of human head and propagation of diffused density waves are depicted in Figure 1(c), where optical transmitter (Tx) and receivers (Rx) are strategically located on head. This study is part of early detection of traumatic brain injury using a helmet mounted functional NIR (fNIR) device, as shown in Figure 1(a). Both wireless connectivity and differential detection of insertion loss (IL) and phase through brain are conducted for a helmet mounted broadband radio frequency (RF) electronics. A three-dimension (3D) FEM modeling is key in a high spatial resolution modeling the modulated photon migration through the inhomogeneous human head. Therefore, this paper will begin with demonstrating confidence in accuracy of the FEM-based modeling by identifying practical mesh sizes as the 2D FEM modeling is being extended to 3D. In particular, accuracy of optical parameter extraction is quantified at the shortest and the most challenging wavelength of 680 nm. Knowing a practical mesh size and the achieved FEM modeling accuracy, then we present accuracy of the two and three layers of phantom stacked to resemble human head by making comparison between the experimental measurement and numerical modeling results at three different wavelengths of 680, 795, and 850 nm. Moreover, the accuracy of 3D modeling to measurements results are compared over various broadband frequency bandwidths. These results will set understanding of how the broadband measurements could be applied for optical parameter extraction of an inhomogeneous brain medium without resorting to analytically challenging methods [10]. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
355
Figure 2. Frequency modulated optical system used on biological phantom.
2. NUMERICAL MODELING AND EXPERIMENTAL METHODS 2.1. Frequency domain fNIR system Functional imaging based on optical photon absorption and scattering related to hemodynamic changes in the human body is developed for different system topologies such as continuous wave, time domain (TD), and frequency domain (FD) [6]. Continuous wave is based on unmodulated photons, while TD systems depend on propagation loss and delay of tunable short NIR pulses. FD is based on frequency swept diffused photon density waves, DPDW, as both optical absorption and scattering parameters of biological tissue are extracted by monitoring amplitude and phase change of modulated photons [11]. Many techniques have been demonstrated in TD to obtain quantitative optical properties such as absorption (a / and scattering (s / coefficients [12]. While the TD spectroscopy instrument has shown good sensitivity, the FD method can be achieved more economically and is more suitable for real time monitoring in clinical settings [13]. The FD method, in which the light source intensity is modulated at either single or broadband frequency, has been applied as an applicable method to study photon migration in a multilayer biological media and therefore monitor and extract the optical properties. As light travels outward from the photon source, optical transmitter (Tx), in a highly scattering medium, the flux energy density decreases more than exponentially with increasing distance from the source. In FD instrument, the intensity of the light source is sinusoidally modulated as shown in Figure 2. The photon density wave can be said to propagate through the biological phantom outward from the source. In a frequency modulated system, the source point is modulated and can be described by the equation S.r, t / D SAC 1 C me j!t (1) where SAC is the amplitude of unmodulated light, and m is optical modulation index for RF modulation frequency of !. 2.2. Finite element method-based numerical modeling The DE is widely used for analytical analysis to explain photon migration dynamics in homogenous medium, where the anisotropic factor, g, is introduced for radiance in a weakly anisotropic medium. Even though g was studied in different homogenous biological media, accuracy of optical parameter extraction and the complexity of the inverse problem solution for inhomogeneous media Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
356
E. SULTAN ET AL.
are still under investigation [14]. A variety of approximation, such as Monte Carlo simulation of the photon behavior, has been extended to two layer diffusion, but it would require high amount of computational power [15]. The analytical model of photon transport in tissue using FEM has been first introduced by Arridge et al. in 1993 [16] to demonstrate the ability of FEM to model the photon density inside an object with photon flux along the geometric boundaries. The FEM modeling manifestation is based on the assumption that photon transport is the diffusion approximation to a radiative transfer equation. The advantage of the FEM modeling is its high accuracy for a practical convergence speed than any other numerical modeling. Other researchers have applied FEM to construct methods of imaging the photon distribution and extract the optical parameters of bio-media. Commercial multiphysics-based FEM software, such as COMSOL, is used in this paper as a robust numerical tool to numerically predict the photon behavior in a highly scattering media, as we have demonstrated in 2D modeling [8]. The Helmholtz equation is used for FEM numerical analysis [16, 17] in COMSOL as r.Dr;.r, t // C a ;.r, t / D S.r, t /
(2)
where ;.r, t / is the photon flux, D D 1 0 , and in a frequency modulated NIR system S.r, t / D 3s SO 1 C me j!t where m is amplitude modulation index at time harmonic modulation frequency of !. The main objective of this paper is to validate the use of a commercially available tool (COMSOL) in 3D to help the process of solving the inverse problem and to extract optical parameters of the homogenous and inhomogeneous media, such as absorption and reduced scattering coefficients for broadband frequency modulated photons (from 30 to 1000 MHz). Accuracy of FEM simulation results depends on resolution of mesh structures in comparison to wavelength of light and modulating RF frequencies and any associated points in space, where DPDW solution is to be computed using partial difference equation of Helmholtz equation. Higher resolution meshing structure is required to accurately simulate the photon dynamics, which naturally results in longer computation times or even cause computation systems to crash; therefore, a compromise in accuracy and mesh resolution has to be made in order to find the optimal simulated environment. It is our experience that having high resolution around optical source and receiver would predict an accurate result for the specific amplitude and phase change while meeting a typical computation time of 4 s for single 2D simulation and 26 s for single 3D simulation using an Intel Core i3 processor PC with a reasonable error of less than 5%—a mesh element sizes of 80 and 1000 nm are considered in regions around optical Tx and Rx and everywhere else, respectively. Different mesh topologies have been reviewed at the wavelength of 680 nm for high meshing of 80 nm selection around the boundaries between layers and the result shown to be the same as if selecting the same meshing around only the Tx and Rx. The difference between these two topologies is the computational time, where applying 80 nm around boundaries takes 94.6 s for single simulation. These meshing element sizes are also recommended when inhomogeneous phantoms resembling head tissue is considered [18]. In COMSOL, the boundary condition is set either through Dirichlet or Neumann boundary conditions, where they are set to express the fluence rate, u, at desired boundaries. Dirichlet and Neumann boundary conditions are chosen appropriately for air-dielectric (or dielectric–dielectric) and radiation condition in region surrounding phantom. The mathematical representation is expressed as n ..cru/1 ..cru/2 / C q u D g huDr
(3) (4)
where h, r, q, and g are the boundary coefficients for the phantom modeling. For the diffused photon flux, h D 1 and r D 0. For the matching boundary of air dielectric and dielectric and dielectric interfaces q D 0 and g D 0. It is important to mention that extrapolated boundary conditions take the phantom-air boundary into account as a radiation boundary by setting the fluence rate to zero using Eq. 3 without any impact on error achieved. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
357
Table I. Manufacture provided optical parameters for the phantoms at 680 nm. a cm1
0 s cm1
0.35 0.15 0.12
15 13 8
Cortex Scalp Skull
Figure 3. Optical system set up showing optical transmitter and receiver measuring insertion loss and insertion phase for three stacked layers of phantoms.
2.3. Experimental measurement procedures Accurate modeling requires experimental verification of change in amplitude and phase of the incident photon wave that is related to absorption and scattering of the biological media is presented in the FD as IL and insertion phase (IP), where IL is any change in amplitude in units of dB, and IP is any change in phase in units of degree. This measurement of IL and IP requires certain hardware and testing equipment that would perform high-frequency modulation [19, 20]. The overall system would consist of Automatic Network Analyzer (Anritsu MS4623B) as an RF source and sensitive RF receiver and broadband optical Tx and sensitive optical Rx. The ANA acts as a frequency modulator and measurements tools, while the optical Tx and Rx modules act as modulated photons source and detector. To a phantom resembling human head, three layers of phantoms were designed with different optical properties and physical geometry as shown in Table I [20, 21], and all three layers can be stacked on top of each other with different combination as shown in Figure 3.
3. EXPERIMENTAL AND SIMULATION RESULTS 3.1. 3D modeling of biological phantom A transition to the 3D is required in order to model the multilayer media with high-spatial resolution. The method used in this section is based on the comparison between results obtained in previous publication [8], and with new 3D numerical model simulation results. This section will focus on accuracy comparison of homogenous phantoms that are modeled in 2D and 3D model at the shortest NIR wavelength of 680 nm in the broadband frequency modulation of the photons. This modeling comparison for an optimum meshing structure will provide confidence in the 3D FEM modeling at a reasonable computation time using our current computational resources for 3D Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
358
E. SULTAN ET AL.
Figure 4. Insertion loss COMSOL simulation for a) two-dimension (2D) scalp and skull, b) three-dimension (3D) scalp and skull, c) 2D reflection mode cortex, d) 3D reflection mode cortex, e) 2D semi-transmission mode, and f) 3D semi-transmission mode.
Figure 5. Insertion phase COMSOL simulation for a) two-dimension (2D) scalp and skull, b) three-dimension (3D) scalp and skull, c) 2D reflection mode cortex, d) 3D reflection mode cortex, e) 2D semi-transmission mode, and f) 3D semi-transmission mode. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
359
frequency modulated photon traveling through multilayer inhomogeneous media. Three different modes of photon traveling are used to make the transition between 2D and 3D FEM modeling; these three modes of transmission, semi-transmission, and reflection can been used separately to extract the optical parameters of the phantoms [1, 8] using COMSOL. A broadband simulation of photon flux is computed for both cases of 2D and 3D, and comparison is made to the experimental results to quantify optical parameter extraction accuracies. Photon flux predictions for transmission mode is applied for scalp and skull phantom, while semi-reflection and reflection mode is applied to cortex phantom as shown in Figures 4 and 5. In the case of transmission and reflection mode, the transmitter is located at the center of the phantom, while for the semi-transmission mode, the optical Tx is located 1 cm away from the edge of the phantom. Location of the optical Rx can be chosen as desired on various positions of phantom, and for accurate optical parameter extraction without being hampered by optical Tx and Rx interface errors, two locations are considered for the differential detection procedures established earlier [1, 8]. The concept of the homogenous modeling depends on identifying an extrapolated boundary region layer and choosing an isotropic point of source limited to the media in sense of space and frequency. The source location is set to be very close to the extrapolated region with certain depth of one transport mean free path ltr beneath the surface of the phantom [18]. Comparison of the experimental results with the FEM numerical result is shown in Figure 6 for broadband frequency modulation. These result in Figure 6 shows that IL between low frequency to high frequency for transmission mode in scalp is about 6 dB and in skull is about 5 dB, while for reflection mode of cortex is about 35 dB and for semi-transmission mode of cortex is about 16 dB. IP between low frequency and
Figure 6. Broadband frequency modulation measurements and simulation at 680 nm for a) insertion loss (IL) scalp and skull, b) insertion phase (IP) scalp and skull, c) IL reflection mode cortex, d) IP reflection mode cortex, e) IL semi-transmission mode, and f) IP semi-transmission mode. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
360
E. SULTAN ET AL.
Table II. Error percentage between FEM using COMSOL and analytical diffusion equation using MATLAB (MathWorks, 3 Apple Hill Drive, Natick, MA 01760-2098, United States). Scalp transmission Error %
30–1000 MHz 30–300 MHz 100–500 MHz 500–1000 MHz 30–500 MHz 100–1000 MHz
Skull transmission Error %
Cortex semi-trans Error %
Cortex reflection Error %
IL
IP
IL
IP
IL
IP
IL
IP
2.32 2.13 3.18 3.56 2.83 2.94
1.92 1.86 2.86 3.04 2.78 2.17
2.55 2.49 3.22 3.16 2.82 2.83
2.12 2.05 2.92 3.44 2.65 2.42
2.62 2.57 3.15 3.57 3.24 2.98
2.87 2.79 3.54 3.68 3.18 3.05
2.8 2.67 3.47 3.76 3.28 3.11
2.51 2.38 3.28 3.42 3.15 2.78
IL, insertion loss; IP, insertion phase.
Table III. Error percentage between two-dimension FEM and three-dimension FEM using COMSOL . Scalp transmission Error %
30–1000 MHz 30–300 MHz 100–500 MHz 500–1000 MHz 30–500 MHz 100–1000 MHz
Skull transmission Error %
Cortex semi-trans Error %
Cortex reflection Error %
IL
IP
IL
IP
IL
IP
IL
IP
1.32 1.58 1.65 1.77 1.64 1.51
1.12 1.24 1.37 1.42 1.32 1.24
1.55 1.63 1.61 1.72 1.75 1.56
1.18 1.31 1.41 1.48 1.35 1.23
1.76 1.85 1.84 1.95 1.87 1.81
1.21 1.44 1.35 1.52 1.42 1.3
1.62 1.75 1.81 1.86 1.76 1.71
1.19 1.35 1.34 1.38 1.34 1.22
IL, insertion loss; IP, insertion phase.
Figure 7. Insertion loss COMSOL simulation for two layer of scalp/cortex or skull/cortex and three layers of scalp/skull/cortex.
high frequency for transmission mode in scalp is about 5ı and in skull is about 3ı , while for reflection mode of cortex is about 14ı and for semi-transmission mode of cortex is about 12ı . Accuracy through computing error between these result is shown in Tables II and III and is discussed in the last section of this paper. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
361
Figure 8. Insertion phase COMSOL simulation for two layer of scalp/cortex or skull/cortex and three layers of scalp/skull/cortex.
Figure 9. Broadband frequency modulation measurements and simulation for two and three layers of phantom of a) 680 nm insertion loss (IL) , b) 680 nm insertion phase (IP), c) 795 nm IL, d) 795 nm IP, e) 850 nm IL, and f) 850 nm IP. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
362
E. SULTAN ET AL.
3.2. 3D multilayers frequency modulated photon travel Because confidence has been established in 3D numerical modeling of the photon evolutions in terms of amplitude and phase losses as it travels through homogenous media resembling scalp, skull, and cortex, then the rest of this paper addresses on modeling the modulated photon evolution as it travels in multilayer media resembling human head. First, the FEM modeling results of two layers of scalp/cortex and skull/cortex are compared against experimental results and then extend for the first time to the challenging problem of multilayer (e.g., a three layer phantom) of scalp/skull/cortex. All of the photon flux analyses are based on three wavelength of 680, 795, and 850 nm as shown in Figures 7 and 8, while comparison between experimental and numerical FEM modeling result for broadband modulation is depicted in Figure 9. IL between two layer of skull/cortex and cortex is about 3–4 dB, two layers of scalp/cortex and cortex is about 5–6 dB, and between three layers and cortex is 7–8 dB. IP between two layer of skull/cortex and cortex is about 1–2ı , two layers of scalp/cortex and cortex is about 3–4ı , and between three layers and cortex is 5–6ı .
4. DISCUSSION AND CONCLUSIONS The confidence in 3D FE modeling depends first on the accuracy between the analytical DE and the numerical simulated result performed by COMSOL. The transition between 2D and 3D layers is performed by identifying the error between these dimensions. The error between 2D and 3D modeling analysis was based on calculating the standard error between the two broadband plots of IL and IP using narrowband and broadband frequencies as shown in Tables II and III. The error between the 2D and 3D FE numerical analysis is no more than 2%, which means the difference between the twodimensional models is insignificant. The next analysis deals with the two and three layers modeling and experimental results. Measurements of broadband frequency modulation along with the numerical simulation are consistent as shown in Figure 9, and it renders that two layers of scalp/cortex and skull/cortex have poor sensitivities as values of IL and IP at low frequencies are very close to one another, but at high frequencies, it shows gradual deviation emerging, which is expected because the shallow layers have different optical properties. The higher the optical properties of the top layer phantom are, the higher IL and IP are expected at higher frequencies. This behavior is shown when at three wavelengths, the IL and IP are higher at high frequencies for the scalp/cortex than the skull/cortex. The important part of this modeling and experimental result is the understanding of the modulated photon behavior in both two and three layers media in order to set the stage for inverse problem solving for medical applications. It is observed that at two layers of scalp/cortex and skull/cortex, the slope of the curve of IL and IP at low frequencies is very close, while at higher frequencies, it changes, and slope change could be used as a point of identification for high optical activity. The last analysis is done to compute the error between the experimental raw data and 3D COM SOL FEM simulation at different band of frequencies. This analysis is performed by curve fitting p the experimental data to a f reque ncy C b and then computing the standard error between curve fitted and 3D COMSOL simulation result as shown in Table IV. The transition between analytical DE and 2D finite element shows an error less than 4%, and error between 2D FEM to 3D FEM applied to homogenous media shows to be insignificant with less than 2%. The last error computed between experimental curve fitted data of two and three multilayers and 3D FEM shows to be less than 5%. This paper addressed two important aspects in the process of 3D modeling of FD functional imaging applied to human head. One is the transition between 2D and 3D FEM for frequency modulation photon traveling in homogenous media and two modeling the modulated photon behavior in inhomogeneous human head media using 3D FEM numerical method. These two important analyses are used to build up confidence in 3D numerical modeling and set the stage for future modeling the modulated photon behavior in inhomogeneous layers addressing different medical cases. Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
Copyright © 2013 John Wiley & Sons, Ltd.
IL, insertion loss; IP, insertion phase.
850 nm
795 nm
2.42 3.89 4.51 4.32 2.39 3.72 4.48 4.27 2.5 3.83 4.63 4.38
680 nm
30–1000 MHz 30–300 MHz 300–500 MHz 500–1000 MHz 30–1000 MHz 30–300 MHz 300–500 MHz 500–1000 MHz 30–1000 MHz 30–300 MHz 300–500 MHz 500–1000 MHz
IL
Wavelength 2.25 2.87 3.15 3.08 2.18 2.78 3.12 3.10 2.15 2.75 3.1 3.07
IP
Cortex phantom Error %
2.52 4.21 4.75 4.52 2.43 4.18 4.59 4.44 2.52 4.27 4.67 4.53
IL 2.22 2.75 3.13 3.04 2.25 2.68 3.17 3.13 2.31 2.74 3.22 3.19
IP
Two layer scalp/cortex Error %
2.31 3.75 4.43 4.26 2.26 3.64 4.31 4.15 2.3 3.68 4.34 4.19
IL
2.19 2.81 3.12 3.05 2.16 2.72 3.08 3.04 2.19 2.75 3.13 3.07
IP
Two layer skull/cortex Error %
Table IV. Error percentage between 3D FEM and curve fitted experimental data.
2.39 4.0 4.53 4.41 2.32 3.94 4.26 4.32 2.28 3.9 4.23 4.28
IL
2.2 2.79 3.17 3.12 2.18 2.68 3.15 3.07 2.25 2.75 3.2 3.14
IP
Three layer scalp/skull/cortex Error %
3D FE NUMERICAL MODELING AND ITS EXPERIMENTAL VERIFICATION
363
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
364
E. SULTAN ET AL.
REFERENCES 1. Sultan E, Manseta K, Khwaja A, Najafizadeh L, Gandjbakhche A, Pourrezaei K, Daryoush AS. Modeling and tissue parameter extraction challenges for free space broadband fNIR brain imaging systems. Proceedings of SPIE, 7902, Imaging, Manipulation, and Analysis of Biomolecules, Cells, and Tissues IX, 790223, 2011, DOI: 10.1117/12.875618. 2. Zhao J, Ding HS, Hou XL, Zhou CL, Chance B. In vivo determination of the optical properties of infant brain using frequency-domain near-infrared spectroscopy. Journal of Biomedical Optics 2005; 10(2). 3. Ferrari M, Quaresima V. A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application. NeuroImage 2012; 63(2). 4. Contini D, Martelli F, Zaccanti G. Photon migration through a turbid slab described by a model based on diffusion approximation I. Theory. Applied Optics 1997; 36:4587–4599. 5. Farrell TJ, Patterson MS, Wilson B. A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo. Medical Physics 1992; 19(4):879–88. 6. Mansouri C, Kashou NH. New window on optical brain imaging; medical development, simulations and applications. Selected Topics on Optical Fiber Technology 2012. Chapter 11, pp. 271–288. ISBN: 978-953-51-0091-1. 7. Izzetoglu K, Bunce S, Onaral B, Pourrezaei K, Chance B. Functional optical brain imaging using near-infrared during cognitive tasks. International Journal of Human-Computer Interaction 2004; 17(2). 8. Sultan E, Najafizadeh L, Gandjbakhche AH, Pourrezaei K, Daryoush A. (2013). Accurate optical parameter extraction procedure for broadband near-infrared spectroscopy of brain matter. Journal of Biomedical Optics 2013; 18(1):017008–017008. 9. Hoshi Y. Functional near-infrared optical imaging: utility and limitations in human brain mapping. Psychophysiology 2003; 40:511–520. 10. Ripoll J, Ntziachristos V, Culver JP, Pattanayak DN, Yodh AG, Nieto-Vesperinas M. Recovery of optical parameters in multiple-layered diffusive media: theory and experiments. Journal of the Optical Society of America A-Optics Image Science And Vision 2001; 18(4):821–830. 11. Pham TH, Coquoz O, Fishikin JB, Anderson E, Tromberg BJ. Broad bandwidth frequency domain instrument for quantitative tissue optical spectroscopy. Review of Scientific Instruments 2000; 71:2500–2513. 12. Grosenick D, Moesta KT, Wabnitz H, Mucke J, Stroszczynski C, Macdonald R, Schlag PM, Rinneberg H. Timedomain optical mammography: initial clinical results on detection and characterization of breast tumors. Applied Optics 2003; 42:3170–3186. 13. Yang Y, Liu H, Li X, Chance B. Low-cost frequency-domain photon migration instrument for tissue spectroscopy, oximetry, and imaging. Optical Engineering 0001; 36(5):1562–1569. 14. Bevilacqua F, Piguet D, Marquet P, Gross JD, Tromberg BJ, Depeursinge C. In vivo local determination of tissue optical properties: applications to human brain. Applied Optics 1999; 38(22):4939–50. 15. Zhu C, Liu Q. Validity of the semi-infinite tumor model in diffuse reflectance spectroscopy for epithelial cancer diagnosis: a Monte Carlo study. Optics Express 2011; 19:17799–17812. 16. Arridge SR, Schweiger M, Hiraoka M, Delpy DT. A finite element approach for modeling photon transport in tissue. Medical Physics 1993; 20:299. 17. Schweiger M, Arridge SR, Hiraoka M, Delpy DT. The finite element method for the propagation of light in scattering media: boundary and source conditions. Medical Physics 1995; 22:1779–1792. 18. Zhou J, Bai J. The most favorable path method for the propagation of light in scattering media: a finite element solution, Engineering in Medicine and Biology Society, 2001. Proceedings of the 23rd Annual International Conference of the IEEE, vol.4, Istanbul, Turkey, 2001; 3203–3206. vol. 4. 19. Manseta K, Sultan E, Khwaja AM, Pourrezaei K, Joshi A, Najafizadeh L, Gandjbakhche A, Daryoush AS. Development challenges of brain functional monitoring using untethered broadband frequency modulated fNIR system. IEEE Microwave Photonics Proceeding, Montreal, Canada, MWP 2010; 295–299. 20. Manseta K, Khwaja AM, Sultan E, Daruwalla P, Pourrezaei K, Najafizadeh L, Gandjbakhche A, Daryoush AS. Untethered helmet mounted functional near infrared (fNIR) biomedical imaging? Microwave Symposium Digest (MTT), 2011 IEEE MTT-S International 2011:1, 4, 5–10. 21. Hoshi Y, Shimada M, Sato C, Iguchi Y. Reevaluation of near-infrared light propagation in the adult human head: implications for functional near-infrared spectroscopy. Journal of Biomedical Optics 0001; 10(6):064032–064032–10.
Copyright © 2013 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2014; 30:353–364 DOI: 10.1002/cnm
