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 ( , 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:
The heat source of the third layer is:
Since the heat source varies with time, the heat source term of each layer should be multiplied by a time factor:
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  . And the blood heat capacity is 3770 J/(kg∙K)  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 , where n is the refractive index of the tissue and , 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 .
Figure 1. The geometry model of the human brain in COMSOL.
Table 1. Thermal parameters of the material in the model  - .
Table 2. Optical parameters of the material in the model .
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, , 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 .
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
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.
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  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 . 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.
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  . 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.
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).
According to Beer-Lambert Law,
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 . According to Equation (A1), AUα = KClα, then
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),
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).