/mo> R ) ϕ 0 exp [ 0.5 r 2 ω 0 2 exp ( μ s z ) μ t z ] (2)

S(r, z) is the Qlaser 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)

φ0 = the flux density of the incident light ( φ 0 = 2 P / π ω 0 2 , W/m2)

P = laser power (W)

ω0 = 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 Z1, Z2, Z3, so the heat source of the second layer is:

S 2 ( r , z ) = μ α 2 ϕ 0 exp { 0.5 r 2 ω 0 2 exp [ μ s 1 z 1 μ s 2 ( z z 1 ) ] μ t 1 z 1 μ t 2 ( z z 1 ) } (3)

The heat source of the third layer is:

S 3 ( r , z ) = μ α 3 ϕ 0 exp { 0.5 r 2 ω 0 2 exp [ μ s 1 z 1 μ s 2 ( z 2 z 1 ) μ s 3 ( z z 2 ) ] μ t 1 z 1 μ t 2 ( z 2 z 1 ) μ t 3 ( z z 2 ) } (4)

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

f ( t ) = exp [ 4 ( t τ ) 2 τ 2 ] (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 ω0 = 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 = ( n 1 ) 2 / ( n + 1 ) 2 , where n is the refractive index of the tissue and n = 1.53 0.2 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 = ( 1.33 1 ) 2 / ( 1.33 + 1 ) 2 = 0.02 .

Figure 1. The geometry model of the human brain in COMSOL.

Table 1. Thermal parameters of the material in the model [17] - [23].

Table 2. Optical parameters of the material in the model [13].

Qlaser of Equation (1) for the CSF is defined as S(r, z) which is calculated by Equation (2). Qlaser for Gray is defined as S2(r, z) which is calculated by Equation (3). And Qlaser for White is defined as S3(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 conv = h ( T inf T ) [10], where qconv is the heat flux due to convection and T is the surface CSF temperature; the heat transfer coefficient, h, is 10 W/(m2∙K) and the external temperature, Tinf, 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

(a) (b)

Figure 2. Characterization of gold nanorods. (a) TEM image of gold nanorods with a geometric dimensions of length = 79 nm and width = 19 nm; (b) A graphs of absorption spectrum of gold nanorods showing a maximum absorption peak at λ = 810 nm.

(a)(b) (c)

Figure 3. (a) Temperature distribution for P = 10 W, τ = 2 ms; (b) Temperature distribution along the central axis of the laser irradiation (T = 8 ms); (c) Temperature versus time for the highest temperature point in (b) (z = 6e−4 m, r = 0 m).

(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

Figure 4. Temperature distribution when the temperature rises to the peak value for different parameters. (a) P = 20 W; τ = 2 ms; (b) P = 30 W, τ = 2 ms; (c) P = 10 W, τ = 20 ms; (d) P = 10 W, τ = 50 ms; (e) P = 10 W, τ = 2 ms, μα of gray matter = 125 1/m; (f) P = 10 W, τ = 2 ms, μα of gray matter = 250 1/m. μα of gray matter in Figures 4(a)-(d) is 25 1/m.

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.

(a)(b) (c)

Figure 5. Maximum temperature values for varying (a) laser power at exposure duration τ = 2 ms, μα of gray matter = 25 1/m; (b) Exposure durations at laser power P = 10 W, μα of gray matter = 25 1/m; (c) Absorption coefficients at pulse duration τ = 2 ms, P = 10 W.

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.


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

Appendix A

According to Beer-Lambert Law,

A U = lg ( 1 / T ) = K C l (A1)


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 α = lg ( 1 / T ) = lg e = 0.43429448 . According to Equation (A1), AUα = KClα, then

μ α = 1 / l α = K C / A U α (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),

A U = K C l = 2 (A3)


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

Cite this paper
Zhou, R. , Chen, H. and Mou, Z. (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. doi: 10.4236/jamp.2019.72029.
[1]   Ackermann, D.M., Foleds, E.L., Bhadra, N. and Kilgore, K.L. (2010) Conduction Block of Peripheral Nerve Using High Frequency Alternating Currents Delivered through an Intrafascicular Electrode. Muscle & Nerve, 41, 117-119.

[2]   Thompson, A.C., Wade, S.A., Brown, W.G. and Stoddart, P.R. (2012) Modeling of Light Absorption in Tissue during Infrared Neural Stimulation. Journal of Biomedical Optics, 17, Article ID: 075002.

[3]   McCaughey, R.G., Chlebicki, C. and Wong, B.J.F. (2010) Novel Wavelengths for Laser Nerve Stimulation. Lasers in Surgery and Medicine, 42, 69-75.

[4]   Izzo, A.D., Walsh, J.T., Jansen, E.D., Bendett, M., Webb, J., Palph, F. and Richter, C.P. (2007) Optical Parameter Variability in Laser Nerve Stimulation: A Study of Pulse Duration, Repetition Rate, and Wavelength. IEEE Transactions on Biomedical Engineering, 54, 1108-1114.

[5]   Ritcher, C.-P., Bayon, R., Izzo, A.D., Otting, M., Suh, E., Goyal, S., Hotaling, J. and Walsh, J.T. (2008) Optical Stimulation of Auditory Neurons: Effects of Acute and Chronic Deafening. Hearing Research, 242, 42-51.

[6]   Wells, J., Kao, C., Konrad, P., Milner, T., Kim, J., Mahadevan-Jansen, A. and Jansen, E.D. (2007) Biophysical Mechanisms of Transient Optical Stimulation of Peripheral Nerve. Biophysical Journal, 93, 2567-2580.

[7]   Izzo, A.D., Walsh, J.T., Ralph, H., Webb, J., Bendett, M., Wells, J. and Ritcher, C.P. (2008) Laser Stimulation of Auditory Neurons: Effect of Shorter Pulse Duration and Penetration Depth. Biophysical Journal, 94, 3159-3166.

[8]   Eom, K., Kim, J., Choi, J.M., Kang, T., Chang, J.W., Byun, K.M., Jun, S.B. and Kim, S.J. (2014) Enhanced Infrared Neural Stimulation Using Localized Surface Plasmon Resonance of Gold Nanorods. Small, 10, 3853-3857.

[9]   Wang, Y. and Guo, L. (2016) Nanomaterial-Enabled Neural Stimulation. Frontiers in Neuroscience, 10, 1-7.

[10]   Vitkaya, A.Y. and Seker, S.S. (2016) Study of Temperature Distribution in Light-Tissue Interaction Using the FEM. Turkish Journal of Electrical Engineering and Computer Sciences, 24, 807-819.

[11]   Fanjul-Velez, F., Romanov, O.G. and Arce-Diego, J.L. (2009) Efficient 3D Numerical Approach for Temperature Prediction in Laser Irradiated Biological Tissues. Computers in Biology and Medicine, 39, 810-817.

[12]   Yoon, G., Welch, A.J., Motamedi, M. and Gemert, M.V. (1987) Development and Application of Three-Dimensional Light Distribution Model for Laser Irradiated Tissue. IEEE Journal of Quantum Electronics, 23, 1721-1733.

[13]   Vo-Dinh, T. (2003) Biomedical Photonics Handbook. CRC Press, Boca Raton, FL, USA.

[14]   Wang, L., Jacques, S.L. and Zheng, L. (1995) MCML-Monte Carlo Modeling of Light Transport in Multi-Layered Tissues. Computer Methods and Programs in Biomedicine, 47, 131-146.

[15]   Pennes, H.H. (1948) Analysis of Tissue and Arterial Blood Temperatures in the Resting Human Forearm. Journal of Applied Physiology, 1, 93-122.

[16]   Li, X. (2004) Numerical Analysis and Experimental Research on Laser Induced Thermal Effect in Bio-Tissues. Tianjin University Press, Tianjin.

[17]   Lu, G.H., Maekawa, T. and Kimura, H. (2006) An Integrated Model of Thermodynamic-Hemodynamic-Pharmacokinetic System and Its Application on Decoupling Control of Intracranial Temperature and Pressure in Brain Hypothermia Treatment. Journal of Theoretical Biology, 242, 16-31.

[18]   Xu, F., Seffen, K.A. and Lu, T.J. (2008) Non-Fourier Analysis of Skin Biothermomechanics. International Journal of Heat and Mass Transfer, 51, 2237-2259.

[19]   Dexter, F. and Hindman, B.J. (1994) Computer Simulation of Brain Cooling during Cardiopulmonary by Pass. The Annals of Thoracic Surgery, 57, 1171-1179.

[20]   Xu, X., Tikuisis, P. and Giesbrecht, G. (1999) A Mathematical Model for Human Brian Cooling during Cold Water near Drowning. Journal of Applied Physiology, 86, 265-272.

[21]   Sukstanskii, A.L. and Yablonskiy, D.A. (2004) An Analytical Model of Temperature Regulation in Human Head. Journal of Thermal Biology, 29, 583-587.

[22]   Olsen, R.W., Hayes, L.J., Wisser, E.H., Nikaidoh, H. and Eberhart, R.C. (1985) Influence of Hypothermia and Circulatory Arrest on Cerebral Temperature Distribution. ASME’s Journal of Biomechanical Engineering, 107, 354-360.

[23]   Nelon, D.A. and Nunneley, S.A. (1998) Brain Temperature and Limits on Transcranial Cooling in Humans Quantitative Modeling Results. European Journal of Applied Physiology, 78, 353-359.

[24]   Duck, F.A. (1990) Physical Properties of Tissue: A Comprehensive Reference Book. Academic Press, Cambridge.

[25]   Torvi, D.A. and Dale, J.D. (1994) A Finite Element Model of Skin Subjected to a Flash Fire. Journal of Biomechanical Engineering, 116, 250-255.

[26]   Fribance, S., Wang, J., Roppolo, J.R., Groat, W.C. and Tai, C. (2016) Axonal Model for Temperature Stimulation. Journal of Computational Neuroscience, 41, 185-192.