ABSTRACT

Objectives: Enhanced infrared neural stimulation (EINS) using nanoparticles is a new research hotspot. In this paper, the numerical modeling of the interaction between a light source and brain tissue during EINS is studied. Materials and Methods: This model is built with the finite element method (FEM) to mimic the propagation and absorption of light in brain tissue with EINS. Only the thermal change is considered in this model since the photothermal effect is the main mechanism of EINS. The temperature response of brain irradiation is governed by the extensively used Pennes’ bio-heat equation in a multilayer model. Results: The temperature distribution in the brain under laser irradiation is determined. And the relationships between the brain tissue temperature and the three factors (the laser pulse time, the laser energy and the enhanced absorption coefficient of the tissue caused by the nanoparticles) are analyzed. Conclusions: The results indicate that the brain tissue is easier to warm up with the enhancement of nanoparticles and parameters of the laser can alter the temperature increase of the brain tissue. These findings offer a theoretical basis for future animal experiments.

Keywords:

Enhanced Infrared Neural Stimulation, Pennes Bioheat Equation, Finite Element Method, Temperature Distribution, Nanoparticles

1. Introduction

Neural functional disorder is a serious disease affecting millions of patients worldwide. All the motor and sensory functions of the body are the process of the muscle or nerve cells changing from a resting potential to an action potential (AP) after receiving a variety of stimuli. Thus, a nerve functional disorder can be repaired by applying an artifact neural interface with an external stimulus such as with electrical, optical, chemical and mechanical methods. Among these methods, electrical stimulation, which has already entered the clinical stage, is the most popular method [1]. However, electrical techniques are challenged due to low data and energy transmission efficiency, poor spatial resolution, and short service time. In addition, people often suffer from muscle fatigue, pain and other negative effects from electrical stimulation because coarse fibers with lower excitation thresholds are usually excited earlier than fine fibers in an electrical stimulus, which is opposite to normal physiology. In this case, high spatial resolution to realize selective stimulation is needed for treating neurological disorders.

Infrared neural stimulation (INS) with a noncontact interface, high spatial selectivity and no stimulation artifacts has advantages in achieving selective stimulation [2]. Some wavelengths of lasers are found succeeded in excite the neuron, such as 1450 and 1540 nm [3], 1850 and 1870 nm [4] [5] [6], 1940 nm [7] and 2120 nm [6]. Some scientists believed that optical sources with wavelengths between 1450 - 2200 nm were the best choice for nerve stimulation because the tissue absorption coefficient of lasers with those wavelengths is higher [6]. However, the higher the absorption coefficient, the shorter the propagation length of light in the biotissue is. Then, only the skin layer of the tissue can be stimulated since the deep layer is hard to reach. Since the photothermal effect is the main mechanism of INS, any wavelength, not only those confined to a 1450 - 2200 nm laser, that can cause sufficient thermal effects may succeed in INS. Therefore, in this paper, an 808 nm laser is tried for INS since this wavelength light has the ability to reach deeper layers. However, the tissue absorption coefficient of the 808 nm laser is too low to excite the nerve. Last several years, scientists found that nanometer particles have very high absorptivity for light shining at their resonance frequency [8] [9]. Therefore, with the help of nanometer particles, laser at wavelength of 808 nm might activate the deep nerves.

Heat injury is a difficult issue of INS that should not be neglected. When the laser energy exceeds two times the required minimum energy, the tissue produces thermal damage [9]. Heat injury is more likely to happen when the laser is used to inhibit the neurons. Therefore, in INS, the most important issue is to control the temperature increase and distribution over the living tissue. The tissue temperature during INS mostly relies upon the light and the thermal performance of tissue. In other words, the thermal response of different tissues will be different even if they are irradiated with the same light. In INS, prediction of the thermal response of the tissue, especially the human brain, is very important. In this paper, a simulation model for INS on human brain tissue is built to accurately calculate the thermal changes of the biological tissue caused by infrared laser irradiation to maintain safe stimulation and to test the possibility of INS by an 808 nm laser combined with nanoparticles.

Many articles have tried to model the temperature changes in biotic tissue under laser irradiation [2] [10] [11]. There are generally two parts of content for the simulation, the light transmission in the tissues and the temperature changes of tissue after absorption of light energy. There are several mathematical models for learning light transmission in the tissues [12] [13] [14], and Monte Carlo method is widely used. However, the Monte Carlo model is only appropriate for optical propagation. For the thermal distribution in biological tissue, which is the other step to simulate the physical interaction of light with biological tissue, the Pennes’ bio-heat equation [15] is the most widely known approach. The Pennes’ equation is quite good for microvessels as it considers the effects of the blood perfusion and metabolism on the biological tissue [10]. In addition, the Pennes’ bio-heat equation is simple and accurate, so it is used to calculate the thermal distribution in brain tissue in this model.

There are two popular numerical methods to implement photothermal model, a finite element method (FEM) and a simple explicit one-step finite difference time domain (FDTD) method [10] [11]. FDTD method can reach convergence in a suitable small time and space steps, but the amount of calculation increases rapidly when used to solve the model with a very small step. Moreover, if used to solve the continuity condition, the space mesh boundaries might be ineffective because of the virtual boundaries. However, the FEM shows high accuracy when dealing with more intricate geometry [10] [11]. Therefore, in this paper, the FEM is selected to solve the temperature distributions using Pennes’ bio-heat transfer equation in Enhanced infrared neural stimulation (EINS), and a Gaussian function is used instead of the Monte Carlo model to calculate the heat source in the Pennes’ bio-heat equation. And in this work, the relationships between the brain tissue temperature and the three factors (the laser pulse time, the laser energy and the enhanced absorption coefficient of the tissue caused by the nanoparticles) are analyzed.

2. Modeling Method

This model simulates the temperature distribution of the human brain tissue under enhanced near-infrared laser irradiation. In this work, a laser with wavelength of 808 nm is chosen as an example and the model is built using the FEM simulation software COMSOL Multiphysics 3.5a.

2.1. Bioheat Equation

In biological tissue, heat transfer is a complex process including energy conduction and radiation, air convection and heat evaporation and the metabolism of biological tissue. In addition, the blood flow also affects the heat transfer. Therefore, the Pennes’ bio-heat equation is used to calculate the temperature of the tissue irradiated by laser exposure:

${\delta }_{\text{ts}}\rho C\frac{\delta T}{\delta t}+\nabla \cdot \left(-k\nabla T\right)={\rho }_{\text{b}}{C}_{\text{b}}{\omega }_{\text{b}}\left(T_{\text{b}}-T\right)+Q_{\text{met}}+Q_{\text{laser}}$ (1)

where:

δ_ts = time ratio (the default value is 1)

ρ = density of tissue (kg/m³)

C = specific heat of tissue (J/(kg∙K))

T = temperature (˚C)

K = thermal conductivity of tissue (W/(m∙K))

ρ_b = density of blood (kg/m)

C_b = specific heat of blood (J/(kg∙K))

ω_b = volumetric perfusion rate (kg/(s∙m))

T_b = temperature of arterial blood (˚C)

Q_met = metabolic heat source (W/m³)

Q_laser = laser heat source (W/m³)

2.2. Laser Propagation Calculation

In Equation (1), Q_laser is derived from the absorption of the laser energy by the tissue. Only a portion of the laser propagating in the tissue is absorbed; the remaining portion will be scattered, reflected or refracted. Here, we suppose the light source is Gaussian laser, so Q_laser of Equation (2) is defined such as [16] :

$S\left(r,z\right)={\mu }_{\alpha }\left(1-R\right){\varphi }_0\exp \left[-0.5\frac{r^2}{{\omega }_0^2}\exp \left(-{\mu }_{\text{s}}z\right)-{\mu }_{\text{t}}z\right]$ (2)

S(r, z) is the Q_laser of Equation (1) and where:

μ_α = absorption coefficient (1/m)

R = specular reflectance

μ_s = scattering coefficient (1/m)

μ_t = μ_s + μ_α, Attenuation coefficient of tissue (1/m)

φ₀ = the flux density of the incident light ( ${\phi }_0=2P/\text{π}{\omega }_0^2$ , W/m²)

P = laser power (W)

ω₀ = Gaussian beam radius (m)

When the biological tissue is divided into multiple layers, the optical properties of each layer are different. The vertical coordinates of the end boundary of each layer are defined as Z₁, Z₂, Z_3, so the heat source of the second layer is:

$S_2\left(r,z\right)={\mu }_{\alpha 2}{\varphi }_0\exp \left\{-0.5\frac{r^2}{{\omega }_0^2}\exp \left[-{\mu }_{\text{s}1}{z}_1-{\mu }_{\text{s}2}\left(z-z_1\right)\right]-{\mu }_{\text{t}1}{z}_1-{\mu }_{\text{t}2}\left(z-z_1\right)\right\}$ (3)

The heat source of the third layer is:

$\begin{array}{l}{S}_3\left(r,z\right)={\mu }_{\alpha 3}{\varphi }_0\exp \{-0.5\frac{r^2}{{\omega }_0^2}\exp \left[-{\mu }_{\text{s}1}{z}_1-{\mu }_{\text{s}2}\left(z_2-z_1\right)-{\mu }_{\text{s}3}\left(z-z_2\right)\right]\\ \begin{array}{c}\\ \end{array}-{\mu }_{\text{t}1}{z}_1-{\mu }_{\text{t}2}\left(z_2-z_1\right)-{\mu }_{\text{t}3}\left(z-z_2\right)\}\end{array}$ (4)

Since the heat source varies with time, the heat source term of each layer should be multiplied by a time factor:

$f\left(t\right)=\exp \left[-4\frac{{\left(t-\tau \right)}^2}{{\tau }^2}\right]$ (5)

where t is the computation time and τ is the laser irradiation time.

2.3. FEM Model: Geometry, Material Properties and Boundary Setting

The geometry of FEM software COMSOL Multiphysics 3.5a is shown in Figure 1. The model was built with the bioheat equation in the heat transfer module in COMSOL. The range of laser irradiation is extremely limited, so the brain tissue under exposure is equivalent to a homogeneous cylinder in order to simplify the model. To further increase the computional speed, the cylinder is established using a 2D axisymmetric coordinate, in which r axis represents the radius of the cylinder (brain tissue) and z axis represents the central axis of rotation. The laser light transmitted by optical fiber can propagate through the skin and skull to directly illuminate the cerebral cortex with the dura torn up in practice. This geometry model consists of three layers: white matter, gray matter and cerebrospinal fluid (CSF). Figure 1 indicates the size of each layer. The thickness of the human cerebral cortex is generally 2.5 - 3 mm, therefore, the gray matter and white matter is set to 2.5 mm respectively. The optical fiber is set 0.5 mm away from the gray matter, which is covered by CSF and the laser spot center is the coordinate origin.

And the subdomain setting of the geometry model is calculated by Bioheat Equation, Equation (1). The thermal properties of the brain layers for Equation (1) are summarized in Table 1. In the simulation examples, the blood water content is 80% so the density of blood is 1060 kg/m [3] [24]. And the blood heat capacity is 3770 J/(kg∙K) [25] while the arterial blood temperature is assumed to be 37˚C and remains constant. The brain is exposed to a continuous-wave laser with wavelength λ = 808 nm and Gaussian beam radius ω₀ = 2e−4 m. For the 808 nm laser, the optical parameters of the brain layers are summarized in Table 2. In Equation (2), the specular reflectance $R={\left(n-1\right)}^2/{\left(n+1\right)}^2$ , where n is the refractive index of the tissue and $n=1.53-0.2\ast w$ [16], where w is the water content in the organization. In this model, CSF is the first layer that the laser irradiates directly, and the water content of CSF is defined as 99%, so n of CSF will be 1.53 − 0.2*0.99 = 1.33, then $R={\left(1.33-1\right)}^2/{\left(1.33+1\right)}^2=0.02$ .

Q_laser of Equation (1) for the CSF is defined as S(r, z) which is calculated by Equation (2). Q_laser for Gray is defined as S₂(r, z) which is calculated by Equation (3). And Q_laser for White is defined as S₃(r, z) which is calculated by Equation (4).

In this COMSOL model, four boundary conditions are set: 1) the axial symmetry is set to thermal insulation; 2) the CSF surface is set as a natural convection boundary considering the convection with the surrounding air, $q_{\text{conv}}=h\left(T_{\text{inf}}-T\right)$ [10], where q_conv is the heat flux due to convection and T is the surface CSF temperature; the heat transfer coefficient, h, is 10 W/(m²∙K) and the external temperature, T_inf, is 25˚C; 3) the interior boundaries are set to continuity; and 4) other boundaries are set as Dirichlet constant temperature, where the initial temperature of both gray matter and white matter was fixed at 37.2˚C and the CSF temperature was fixed at 37.1˚C [17].

2.4. Nanoparticle

Gold nanorods are selected as an example in this paper. As shown in Figure 2(a), the size ratio of the nanorod is 4.1 with 79 nm length and 19 nm width on average. Figure 2(b) provides the absorption spectra of nanorods with a concentration of 3e−5 g/l. As we can see form the diagram, the strongest absorption is at 810 nm, which is very closed to 808 nm laser that we used. The absorbance of the nanorods at the 808 nm is approximately 2, so the absorption coefficient of the nanorod at this concentration is about 460 1/m (Appendix A). This absorption coefficient can be easily adjusted by changing the concentration of nanorod.

3. Results and Discussions

Simulation results on the temperature distribution of the brain model at different laser pulse duration and varying laser power are displayed in Figure 3, and Figures 4(a)-(d). In Figure 3, the laser energy is 10 W with 2 ms pulse duration. Figure 3(a) is the temperature distribution of the multilayer model at Time = 8 ms

(laser pulse opens at Time = 0 ms). Although the CSF is in front of the gray matter, the temperature of the CSF is quite lower than that of the gray matter. Furthermore, almost all the energy is gathered in the gray matter, and little energy from the laser can propagate into the white matter. Figure 3(b) shows the temperature distribution along the central axis of the laser irradiation that is z axis (r = 0 m) at Time = 8 ms. We can see that the highest temperature point is between z = 5e−4 m and z = 1e−3 m which is in the superficial layer of the gray

matter (Figure 3(b)). The above phenomenon is because the gray matter has the strongest absorption of the 808 nm laser compared to CSF and white matter. Figure 3(c) shows the temperature versus time for the highest temperature point in Figure 3(b) (z = 6e−4 m, r = 0 m). As shown in Figure 3(c), the temperature does not stop rising when the laser pulse stops due to the thermal hysteresis. In addition, the temperature rises rapidly to the highest point and then drops very slowly after reaching the peak value.

For laser irradiation at different pulse durations and powers, Figures 4(a)-(d) provides the temperature distribution of the brain model while the temperature rises to the peak value. All of the results from the above simulations are summarized in Figure 5(a) and Figure 5(b), which clearly show the varied effects of pulse duration and laser energy on the maximum temperature values. It is clear that with the increase of the pulse duration, the maximum temperature of the brain also increases. In addition, when laser power is increased by keeping the exposure duration constant, the maximum temperature of the brain is increasing linearly, which is consistent with Vitkaya’s [10] simulation result.

According to Fribance’s simulation result, local temperature needs to be rapidly increased by 6.6˚C - 11.2˚C in 0.1 - 2.6 ms to activate the axon [26]. Wells et al. also demonstrated that the axon temperature needs to be increased by 6˚C - 10˚C in 0.005 - 5 ms to activate action potentials. However, Figure 4(b) shows that the temperature increase is only (41.177 − 37.2) = 3.977˚C within 2 ms even when the laser power reaches 30 w, which is very high for an 808 nm laser. Therefore, it is very unlikely to be able to activate the nerve using the 808 nm laser.

Recently, researchers have used different materials nanoparticle, such as gold, silica and silver to enhance the absorption of light because of their localized surface plasmons (LSPs) capability [8] [9]. Essentially, nanoparticles are the main suction light components after injected into the tissue since their absorption coefficient is much higher than the surrounding tissue. Therefore, the absorption of light for the injection parts can be similar to the absorption of light by the nanoparticles. Nanoparticles cause temperature rise after suction light, then heat transfer to the surrounding nerves, so as to make the temperature rise of nerve. According to paragraph “2.4 Nanoparticle”, the measured μ_α of the nanorod at concentration of 3e−5 g/l is about 460 1/m, which is 18.4 times of μ_α of gray. And this absorption coefficient can be easily adjusted by changing the concentration of nanorod. Then we analyzed the effect of the tissue absorption coefficient (which are 5 times and 10 times of μ_α of gray) on the temperature distribution, results are provided in Figure 4(e) and Figure 4(f). The maximum temperature values of different absorption coefficients are summarized in Figure 5(c). It can be clearly inferred that as the absorption coefficient is increased by 4 times or 9 times by keeping the pulse duration and laser energy constant which is easily implemented by adding nanoparticles, maximum temperature of the brain is increased, and it can easily rise above 6˚C, and the Neurons will be excited by near infrared laser combined with nanoparticles. Therefore, near-infrared laser will be able to activate the deeper nerves with the help of nanoparticles.

4. Conclusions

INS receives more and more attention due to its advantages including high spatial selectivity and no need for contact. The local transient temperature gradient induced by the absorption of the laser is the main mechanism of INS. Therefore, the study and prediction of the temperature distribution in biological tissues in INS is frequently required.

In this work, the FEM numerical algorithm is proposed to solve the Pennes’ bio-heat equation to get the temperature distribution of the brain irradiated by an 808 nm laser. The results obtained from this study indicate that the parameters of the laser that clearly affect the temperature distribution of brain are pulse duration, laser power and the absorption coefficient. The tissue absorption of the near infrared laser increases by adding nanoparticles around the tissue. By changing the absorption coefficient in this model, we found that the tissue temperature could rise above 6˚C with exposure to the near infrared laser, proving that it is possible to excite the nerve by near-infrared laser enhanced by nanoparticles.

The results of the brain photothermal model can provide a suitable parameter reference for laser treatment. It is possible to improve the model accuracy by finding a more accurate absorption and scattering coefficient of brain for the light, as well as more accurate thermal parameters for the tissue.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (grant number: 31500796).

Conflicts of Interest

The authors have stated explicitly that there are no conflicts of interest in connection with this article.

Cite this paper

Zhou, R., Chen, H.L. and Mou, Z.X. (2019) FEM Model of the Temperature Distribution in the Brain during Enhanced Infrared Neural Stimulation Using Nanoparticles. Journal of Applied Mathematics and Physics, 7, 381-393. https://doi.org/10.4236/jamp.2019.72029

References

Appendix A

According to Beer-Lambert Law,

$AU=\lg \left(1/T\right)=KCl$ (A1)

where,

T = the transmittance (the ratio of transmitted light intensity to incident light intensity),

K = the extinction coefficient of the solution (L/g/cm),

C = the concentration of the solution (g/L),

l = the optical path of the solution (cm).

In this paper, μ_α = 1/l_α

μ_α = the absorption coefficient (1/cm),

l_α = the absorption length (cm),

l_α is also the optical path when T = 1/e,

where $A{U}_{\alpha }=\lg \left(1/T\right)=\lg \text{e}=0.43429448$ . According to Equation (A1), AU_α = KCl_α, then

${\mu }_{\alpha }=1/l_{\alpha }=KC/A{U}_{\alpha }$ (A2)

From the measurement results we can see that AU is about 2 at λ = 808 nm, where C = 3e−5 g/l, and l = 1 cm. So according to Equation (A1),

$AU=KCl=2$ (A3)

where,

K and C is equal to those in Equation (A2), and l is the equal the thickness of the cuvette (l = 1 cm). Then solve Equation (A2) and Equation (A3), μ_α = 4.6 (1/cm) = 460 (1/m).