Modeling time resolved light propagation inside a realistic human head model.

BACKGROUND: Near infrared spectroscopy imaging is one of the new techniques used for investigating structural and functionality of different body tissues. This is done by injecting light into the medium and measuring the photon intensity at the surface of the tissue.
METHODS: In this paper the different medical applications, various imaging and simulation techniques of NIRS imaging is described. Each method is introduced and discussed. Then, the optimized model is prepared for numerical simulations. In this paper, the finite element method is used for solving the diffusion equation numerically.
RESULTS: Diffusion equation was solved for realistic human head model using finite element approach for a point light source and time resolved case. The photon intensity distribution in different head layers has been shown and the intensity orientation via the CSF layer has been illustrated.
CONCLUSION: Simulating the photon transformation inside the tissue is essential for investigating the NIRS imaging technique. The finite element approach is a fast and accurate method for simulating this fact. The time resolved approach of this technique could illustrate the photon migration and intensity orientation in the tissue for time dependent light sources in tissues.

Introduction

One of the greatest advances in recent history of medical science is the rapid progress in medical imaging techniques development. These techniques make the disease diagnosis procedure more accurate, facile, and quick. A novel method in medical imaging applicable to investigate different body parts with high temporal and partly high spatial resolution is Near InfraRed Spectroscopy (NIRS). This method uses the near infrared light (with wavelength of 650-900nm) to penetrate the biological tissue and presents measurements to study the tissue structure and/or functional. In this technique, the light sources are placed on the surface of tissue, injecting light into the medium and an array of detector fibers also placed on the same tissue surface measuring the transilluminated light. These measurements can be used to investigate the tissue structure and functionality. Nowadays this imaging method is widely used in breast and prostate tomography, and brain imaging, particularly suitable for young infants.

There are several effective properties that this imaging method profits: low operating cost, tolerance of movement, being noninvasive, capability of investigating real-time in vivo functionality of tissue, ability of constructing three dimensional images of biological tissues and determining the amount of water and oxy and de-oxy hemoglobin in the medium. These are rarely come together in other medical imaging techniques. All these properties make the near infrared imaging a powerful, evolutionary and widely applicable medical imaging technique.

In the next section, medical applications of this imaging technique are presented. The NIRS imaging techniques are described in third section.  The fourth section is dedicated to techniques that are used for simulating the light propagation in tissue. In the fifth section the Finite Element Method (FEM) is introduced and presented for simulating the Diffusion equation numerically in a medium. The simulation results are shown in sixth section. The discussion and conclusion is presented in the final section.

NIRS applications

Cognitive developmental neuroscience and psychology

Because of different molar extinction coefficient for oxy and de-oxy hemoglobin in near infrared wavelengths, each media has its own attenuation (figure 1). Using the difference of measured attenuations for different light frequencies, the amount of oxy and de-oxy hemoglobin’s concentration could be determined inside the tissue. The functionality of brain could be inquired by capturing these measurements. Lots of efforts have been done in using near infrared imaging to investigate the functional and development of infant’s brain [1] such as language development study [2, 3], speech notwithstanding and language-specific reactions [4, 5], response to upright and inverse faces [6], recognizing faces seen from different angles [7], eye gaze, eye, and mouth movements [8, 9] and modality-specific effects of joint auditory and visual stimulation [10-13]. The use of near infrared imaging would lead to more acquaintance of human brain functionality and development from the birth to childhood.

Figure 1

Different molar extinction coefficient for oxy and de-oxy hemoglobin in near infrared wavelengths lead to different attenuation for each media [1].

Breast tomography

The breast cancer is the most frequent malignancy in women worldwide [14-17]. The diagnostic mammography has been shown to have a sensitivity of 75% [18] in cancer detection. The MRI can solve this problem and boost the sensitivity up to 95%, but its specificity is significantly lower than mammography [19] and also it is more expensive. Thus there would be demands on a complementary imaging technique to improve the diagnostic accuracy in breast cancer detection.

The tumor tissues need more vessel and blood to grow and survive than healthy ones. On the other hand, the near infrared imaging is capable of monitoring the total value of hemoglobin inside the medium. This imaging technique could identify tumors inside the breast.  Furthermore, the breast surface is fully accessible to place the source and detector fibers, it has small size and there are no bony structures inside it. The influence of using this imaging technique as a complementary to mammography on diagnostic accuracy has been studied in [18].

Prostate tomography

The prostate cancer is one of the most diagnosed cancers in men [20]. This kind of cancer is screened using several kinds of tests including measurement of serum prostate-specific antigen (PSA) [21], digital rectal examination (DRE) and usually a combination of both techniques are used to detect the prostate cancer [22]. It is difficult to differentiate between benign and malignant case using just the measurement value of PSA [23, 24]. On the other hand, the palpation during a DRE, is subjective, insensitive, and more than half of all prostate cancers detected are not palpable [22]. The difference between blood and parenchymal tissue is very difficult to obtain using conventional imaging techniques but the near infrared attenuation through tissue demonstrates a significant contrast gradient between these two [25-30]. In multi-spectral case of near infrared imaging even the malignancy is quantified directly [25-30]. The studies had clarified the vascular density gradient difference between the malignant and benign tissues in prostate [31]. Furthermore, the different water concentration in cancerous and benign tissue has been investigated [32]. There are lots of efforts done based on non-invasive near infrared imaging for prostate tumor detection [32-34]. There are reports for surface measurements of implanted prostate tumor in [35, 36]. The studies demonstrate the potential of using near infrared to detect and characterize prostate cancer [37].

NIRS imaging techniques

Continuous light near infrared imaging

This imaging technique is the most common one that is widely used in brain investigations [1]. In this method the monochromatic near infrared light is induced into the medium and using the measurements of light attenuation, the amount of the change in oxy and de-oxy hemoglobin could be determined inside the tissue. In fact when light is propagating through the medium, it could be absorbed, scattered or transmitted without any interference with the media. The influence of each event on the light is totally based on the medium light properties at that particular wavelength. The concentration of certain chromophors in the biological tissues can be determined by calculating the absorption and scattering coefficients of tissue using the light energy measurements.

Frequency resolved near infrared imaging

Another approach to near infrared imaging is the frequency resolved one. In this method the near infrared light is modulated at a specific frequency and then it is injected into the medium. At the light detectors, the attenuation and the change in modulation phase of light caused by the scattering events are measured. These measurements can lead to calculate the absolute amount of oxy and de-oxy hemoglobin, total hemoglobin, and water within the tissue. The frequency resolved near infrared imaging was first introduced and tested by [38]. This kind of imaging is usually used in tomography applications. In [39] the medium reconstruction by solving the inverse problem in finite element scheme is described.

Time resolved near infrared imaging

In time resolved near infrared imaging the single light pulses are used instead of continuous wave. When a very short pulse of light is injected into the tissue, it travels through multiple scattering events and is detected by detectors. By measuring the difference between light injection and detection time and multiplying it by light speed, the mean path length of light can be determined inside the tissue. These systems were introduced by [40] in 1988 and since then lots of systems have been developed [41-43] using this technique.

NIRS simulation techniques

According to the applications of NIRS, simulation of the imaging procedure would be very essential in designing and developing such systems. On the other hand, studying near infrared imaging technique is based on scrutiny of light propagation modeling in tissue. There are two main approaches in modeling this phenomenon: Monte-Carlo and partial differential equation based approach.

The Monte-Carlo approach

This method has been first introduced by Wilson and  Adam [44] in the field of light-tissue interaction. With improvements made in this method in recent years [45-49], it is widely used to simulate the light transillumination within tissues. These models are particularly based on the probabilistic interpretation of photon existence and movement in the medium. The results of Monte-Carlo simulations acquired after launching large number of photons into the medium and tracking each photon’s path until it is observed by or exits the tissue. This technique is capable of checking each photon’s transformation path individually. This is important when there is a problem defined based on tracking light path from a particular light source to some specific light detector.

This technique results in a very high accuracy but it suffers from significant low computational efficiency. For example for a typical domain size in the human head over an hour or more of computational time is required to simulate the light propagation [50]. In recent years, using massively powerful graphic accelerator processors and entering the parallel processing techniques in Monte-Carlo modeling, the simulation time is decreased significantly up to 300 times faster [51]. However, these hardwares are expensive and conflict with low-cost property of this imaging technique when imaging devices need to accomplish photon migration simulations during imaging procedure.

Partial Differential Equation (PDE) based approach

In these techniques the light propagation in tissue is modeled using partial differential equations. In this approach, the light is considered as following two main cases:

1.   Light as electromagnetic wave:

In this case the light is considered as electromagnetic energy. Therefore the Maxwell’s equations accurately describe its propagation behavior in a medium. However, with assumption of tissue to be a complex medium containing multiple scattering events, the computational efficiency of methods solving Maxwell’s equations is decreased massively [52].Thus this model is not widely applied to cases contain light migration in tissue.

2.   Light as particles:

In this case the light is considered to be made of massless neutral particles called photons. The equation describing the particle like nature of light is called the radiative transfer equation (RTE). This equation is precisely describing photon migration in any well-defined medium. In the following paragraphs, some definitions are given before introducing RTE.

Specific intensity u: The photon distribution can be expressed in terms of number of photons dn crossing a surface of unit area, per unit time, per unit frequency dv , per unit stradian dΩ. This quantity is simply the ratio between radiant intensity and the energy of one photon, and is given as follows:

Absorption µa: The probability of a photon being absorbed while traveling a distance ds is µ. The coefficient µa (m -1 ) is defined as absorption and  is the true loss of a photon from count.

Scattering µ's: The definition of scattering coefficient µ's is like the definition of absorption. It means that the probability of a photon being scattered while traveling a distance ds is µ's (r, v, Ω, t). A scattering event changes the direction and/or frequency of the photon during the propagation in medium. So it is convenient to define a differential scattering coefficient µs:

which is the probability of scattering from direction Ω' to direction Ω and frequency v' to v while traveling distance ds. The scattering coefficients µ's and µs are related as follows:

Emission f : Photons can be injected into the medium by external and/or internal sources. The number of photons emitted by volume dr at r in the direction Ω about the differential solid angle dΩ, at frequency v at the interval v to v+dv , between t and t+dt is:

The equation of transfer: It can be shown that for known medium properties, the intensity u can be calculated using the 4th equation.

After imposition of boundary and initial conditions, this equation is called the radiative transfer equation (RTE), which can be used for modeling the light propagation in tissues.

3.   Diffusion approximation of RTE:

The RTE is a complex equation and solving this equation analytically and/or numerically in medium with complex geometries, e.g. human brain, is difficult. However, with the following assumptions the RTE can be simplified: i) The medium is more scattering than absorbing, ii) The medium is isotropically scatters, and iii) The light propagation is far enough from sources and boundaries [53, 54]. Using these assumptions and expanding the RTE over spherical harmonics and using the first two terms of this expansion leads to a new equation named the diffusion equation that can be solved using finite element method (FEM); equation 5.

Where D(r,t) is the diffusion coefficient

The µa and µ's are absorption and scattering coefficients respectively, and c is the speed of light. The boundary condition for modeling the light propagation in tissues for diffusion equation is modified Robin boundary condition [39] which is given by

where u(r,t) is the energy fluence rate on the boundary, and A is the refractive index mismatch coefficient of the form

and

and

where n is the refractive index inside the boundary and n is the refractive index outside the boundary.

Finite Element (FE) solution of Diffusion equation

At 1993 S.R.arridge et al [55] solved the diffusion equation using FEM. They had also do the comparison of this approach with the Monte-Carlo method and showed that the error between these two methods is negligible. In order to solve diffusion equation using FEM, the first step is to discretize the domain. This can be performed using meshing procedures. At the next step, the weak form of the equation should be obtained. For this reason the equation (5) is multiplied by the weight function and integrated over each element. Then, the integration by parts over the term involving higher order derivatives using the gradient or divergence theorem results in equation 11 [56].

where the last term in the right hand side q contains the boundary terms for the element Ω. The parameter v is a function of only x and y. The integration by parts is applied only over spatial dimensions. At the next step, a finite element approximation for u is substituted considering that the time dependence can be separated from the spatial variation. This approximation is given by:

where uͤ denotes the value of u(x,y,t) at the location ( x,y) at time t. In finite element method, by substituting v=ψͤ in equation (11) and u by equation (12) in (11):

which can be rewritten in matrix form given by

where

The boundary condition (7) can be applied to the equation by its influence on stiffness matrix K [57]

where l is the length of the segment of the element along the boundary, λ=1/2A and A is the refractive index mismatch described in equation (8).

Transient Analysis: For transient analysis, the finite difference scheme for diffusion equation over t is used. Here by using the α family of approximation [56] we have

where u̇=𝜕u/𝜕t is the derivative with time. Using equation (19) and (14), the equation (13) can be rewritten using a set of algebraic equations at time t+1 [56]:

where

Equation (20) is solved for each time t after the matrix assembly and imposition of boundary condition. The initial value of u ({ u}) is used at the right hand side of the equation (22) for the first time step t = 0.

Methods and Results

Simulation results for diffusion equation using FEM

The model was applied to study the light propagation in the brain. A sagital brain slice was selected from simulated magnetic resonance image using Brainweb simulator. Meshing procedure was performed using COMSOL MULTIPHYSICS software package. In order to obtain accurate results, the path of light propagation was meshed in smaller mesh sizes in comparison with the entire model (figure 2).

Figure 2

Geometry of the applied model

Medium was divided into 12576 triangle meshes using 6338 nodes. The source was placed at a depth 1/𝜇ʹ from the boundary  and it was considered as a delta function in time f(x,y,t)=δ(x. The time steps used for this simulation was 0.1ns. The medium light properties is shown in table 1.

Table 1

Specified absorption and scattering coefficients for different tissue layers of brain.

Tissue type	Absorption coefficient(1/mm)	Scattering coefficient (1/mm)
Scalp and Skull	0.019	0.86
CSF	0.004	0.7
Brain	0.01	1.11

Problem was solved by FEM using discussed configuration and the results are shown in figure 3. The contour plot (figure 4) shows the propagation distribution at time 1.2ns where light is crossing the scalp and is entering the CSF and brain.

Figure 3

Simulation of light propagation in the brain using FEM

Figure 4

Propagation distribution at time 1.2ns where light is crossing the scalp and is entering the CSF and brain

Discussion

In this paper an overall introduction to NIRS imaging systems has been presented. The various applications, imaging techniques and simulation models have been discussed; and the FEM has been introduced in order to numerically solve the Diffusion equation. At the end the time resolved Diffusion equation has been solved in a head model created based on brain MRI. From results, it is obvious that the light can penetrate the scalp and skull and reach the CSF and brain. Furthermore, the propagation distribution of photons in different layers of head has been investigated. The results show that when light is entering the CSF, photon distribution is oriented along this layer of brain and it is because of the lower absorption coefficient of this layer of head.