Computational Modeling of the Photon Transport, Tissue Heating, and Cytochrome C Oxidase Absorption during Transcranial Near-Infrared Stimulation

Transcranial near-infrared stimulation (tNIRS) has been proposed as a tool to modulate cortical excitability. However, the underlying mechanisms are not clear where the heating effects on the brain tissue needs investigation due to increased near-infrared (NIR) absorption by water and fat. Moreover, the risk of localized heating of tissues (including the skin) during optical stimulation of the brain tissue is a concern. The challenge in estimating localized tissue heating is due to the light interaction with the tissues’ constituents, which is dependent on the combination ratio of the scattering and absorption properties of the constituent. Here, apart from tissue heating that can modulate the cortical excitability (“photothermal effects”); the other mechanism reported in the literature is the stimulation of the mitochondria in the cells which are active in the adenosine triphosphate (ATP) synthesis. In the mitochondrial respiratory chain, Complex IV, also known as the cytochrome c oxidase (CCO), is the unit four with three copper atoms. The absorption peaks of CCO are in the visible (420–450 nm and 600–700 nm) and the near-infrared (760–980 nm) spectral regions, which have been shown to be promising for low-level light therapy (LLLT), also known as “photobiomodulation”. While much higher CCO absorption peaks in the visible spectrum can be used for the photobiomodulation of the skin, 810 nm has been proposed for the non-invasive brain stimulation (using tNIRS) due to the optical window in the NIR spectral region. In this article, we applied a computational approach to delineate the “photothermal effects” from the “photobiomodulation”, i.e., to estimate the amount of light absorbed individually by each chromophore in the brain tissue (with constant scattering) and the related tissue heating. Photon migration simulations were performed for motor cortex tNIRS based on a prior work that used a 500 mW cm−2 light source placed on the scalp. We simulated photon migration at 630 nm and 700 nm (red spectral region) and 810 nm (near-infrared spectral region). We found a temperature increase in the scalp below 0.25 °C and a minimal temperature increase in the gray matter less than 0.04 °C at 810 nm. Similar heating was found for 630 nm and 700 nm used for LLLT, so photothermal effects are postulated to be unlikely in the brain tissue.


Introduction
Near-infrared (NIR) light has been reported to be able to penetrate the extra-cranial layers such as scalp, skull, and cerebrospinal fluid, and reach the superficial layers of the cerebral cortex due to the optical window. It has been hypothesized that interaction of NIR light with cytochrome c oxidase (CCO) can potentiate the CCO in the mitochondria, a component of the electron transport chain and key complex in energy production [1]. CCO is the primary chromophore in the mitochondria (near-infrared laser) [24]. Moreover, the use of transcranial NIR laser (810 nm) in low-level light therapy showed improvement in patients suffering from anxiety and depression [25]. Since increased oxygen consumption occurs during increased neural activity [26], which leads to increased CCO activity, so an assessment and modulation of CCO activity can open a pathway to monitor and modulate neuronal activity [27,28]. Here, redox state-dependent changes in the NIR spectrum is an essential tool for near-infrared spectroscopy of the oxidation state of CCO [2].
In this paper, investigation of the interaction of light with the chromophores that are responsive to photons in the red-NIR spectral region has been performed. Since there is an increased absorption of longer wavelengths by water in the tissue, we postulate that NIR light interaction with neural tissue may have effects of photobiomodulation as well as photothermal neurostimulation which needs consideration for rational dosing of tNIRS due to the biphasic dose response. Although it has been reported that tNIRS is a modulator of cortical excitability in a healthy human brain [1] which forms the basis of this paper, the exact mechanisms of the neuromodulation has been elusive. In this paper, we apply computational modeling to dissociate photothermal effects from photobiomodulation during tNIRS with 810 nm while comparing that with the red spectrum (630 nm and 700 nm) used for LLLT [16] to investigate the mechanisms underlying neuromodulation [1]. Here, the primary aim is to better understand the extent of optically induced tissue heating (primarily due to water and fat absorption) during tNIRS based on the experimental results by Chaieb et al. [1].

Head Model Selection
To develop the computational model of light interaction with the chromophores in the human head by non-invasive approach, a digital brain phantom based on high-resolution brain atlas [29] was used in the study. From the Colin27 head atlas, the different layers of the brain were segmented to form layered tissues of the head model [30]. Volume mesh was created from each layer after surface smoothing using CGAL surface mesh toolbox [31] with each surface having its own mesh criteria and density. After the multi-layered surface mesh was generated, the volume mesh was generated using the Delaunay tetrahedralization algorithm [32]. Figure 1 shows the multi-layered head mesh generated after Delaunay tetrahedralization.
The details of the mesh components of the complete four-layered Colin27 head model is given in Table 1.

Geometry and Domain Assignment
The tetrahedral mesh was imported to COMSOL Multiphysics software, a Finite Element Method(FEM) simulation software, as a Computer Aided Design (CAD) model and each layer was designated as a domain and corresponding optical properties were assigned to each domain. The domains are as follows: •

Simulation of Radiative Transfer Equation Using Diffusion Approximation
Photon transport in scattering media, such as biological tissues, is generally modeled using the radiative transfer equation (RTE) [33] due to its more accurate solution for highly scattering medium as in the case of brain tissues [34] and higher computational efficiency for complex medium [35,36]. The tetrahedral mesh generated was converted into a CAD file and imported to COMSOL. After importing the mesh to COMSOL, the computation of photon propagation was solved through the diffusion approximation of the RTE. The second order partial differential equation (Equation (1)) describes the time behavior of photon fluence rate distribution in a low-absorption high-scattering medium.
Here, µ a , µ s , and µ t are the absorption, reduced scattering, and total attenuation coefficients, respectively, L(r, S, t), the radiance at position r with direction of propagation S , v the velocity of light through the medium (v = c/n where c is the velocity of light in vacuum and n the refractive index of the medium), Q(r, S, t) the source term, and S, S , r the phase function for scattering. A standard approximation method for the RTE assumes that the radiance in tissue can be represented by an isotropic fluence rate, φ(r,t), plus a small directional flux, J(r,t), where: The final diffusion approximation of RTE, i.e., diffusion equation, is derived as: where D(r) is defined as: The µ s is the reduced scattering coefficient and is obtained from Equation (6): where g is the anisotropy factor. The diffusion equation is solved by using the COMSOL Multiphysics software using the Partial Differentiation Equation (PDE) toolbox (comparison with Monte Carlo simulation shown Figure S1). The entire head model had four domains: scalp and skull combined, Cerebrospinal Fluid (CSF) , gray matter, and white matter, as listed in Table 2. The physics was applied to each domain at steady state with the initial condition being zero. The source term was taken from the published literature [1] where the power density was 500 mW/cm 2 at the scalp surface as presented by Chaieb and colleagues [1]. The head model was assumed to be surrounded by air at room temperature (25 • C). We placed our sources at the air-tissue interface, which is at the scalp, following Chaieb and colleagues [1]. The boundary condition here is as follows: The optical properties, namely scattering coefficient at these wavelengths have been reported in various prior works [37,38]. The absorption coefficients of the tissues are calculated as the summation of the absorption coefficient due to the contribution of each component of interest in the corresponding tissue [39]. The optical properties for the whole tissues at the three wavelengths used in the study are given in Table 2. The reduced scattering coefficients are calculated based on the scattering coefficient and the anisotropy factor (Equation (6)). The anisotropy factor, g = 0.89 has been assumed for all the tissue layers. Although the literature has shown that diffuse reflection occurs at the skin surface, in this paper, the reflection effects have been excluded.

Optically Induced Thermal Effects Modeled Using Bioheat Transfer Mechanism
The thermal effect due to the absorbed incident light is modeled using the bioheat transfer mechanism [40]. The algorithm analyzes the temperature distribution and heating profile when the heat is applied to the tissue. The Penne s Bioheat equation (Equation (8)) is used to model this phenomenon for localized and distributed energy source.
Here, ρ (kgm −3 ) is the tissue density, c is specific heat of the tissue (kJ/kg/K), K is thermal conductivity, cb (3664 J/kg. • C) is blood specific heat, ω b is blood volumetric perfusion rate, T a is the arterial blood temperature (37 • C), ρ b (1050 kgm −3 ) is the blood density and Q met and Q r are the volumetric metabolic heat and the external spatial heating respectively. The heat source term is related to the local fluence rate and tissue absorption coefficient [41] as follows: The bioheat physics is applied with the different tissue components with their respective thermal and blood perfusion properties. The thermal and the blood perfusion parameters are taken from [42] to be used for the computation of the bioheat transfer, as shown in Tables 3 and 4. The boundary condition was applied at the skin surface. It was assumed that there is heat loss at the skin surface by convection to ambient [42]. For the whole scalp(skin), the convective heat flux value was assumed to be 4 W/m 2 . • C [43].

Investigation of Individual Chromophore Absorption in the Tissue
Absorption of red or near-infrared photons by cytochrome C oxidase (unit IV of the mitochondrial respiratory chain) has been established by prior works [44,45]. In this section, we investigated the absorption by CCO at 630 nm, 700 nm, and 810 nm wavelengths, which can cause its activation and may lead to photobiomodulation in the gray matter [44]. Besides CCO, we also investigated other major chromophores in the gray matter along with the investigation of water absorption. Chromophores present in the gray matter that were investigated in this section are as follows: Cytochrome c oxidase(reduced and oxidized state) • Lipid Figures 2 [46], 3 [47], 4 [48] and 5 [49] shows the absorption spectra of two states of hemoglobin, two states of cytochrome c oxidase and lipid respectively. In Figure 4, the main plot shows the absorbance due to 4.9 µM of CCO and the inset shows the absorbance due to five times the concentration of CCO [48].

Optical Parameters of Individual Chromophores
The three wavelengths, 630 nm, 700 nm, and 810 nm in the red and the NIR spectral regions, have been reported to be promising for photobiomodulation [44]. The two reported wavelengths were chosen from red (630 nm, 700 nm) and one in NIR (810 nm) spectral region.
The water has very low absorption in the red and near-infrared spectrum, although it increases with increasing wavelength. Calculation of tissue absorption specifically due to water was performed by obtaining the value of the absorption coefficient of pure water at the three wavelengths [46,50,51]. Since 75% water per unit volume (i.e., volume fraction) is present in brain tissues, hence, 0.75 mu a is the absorption coefficient of the tissue specifically due to water [52] (more details provided in the Supplementary Materials). The percentage of the dry weight of lipid in gray and white matter [53] (i.e., mass fraction), density of gray and white matter along with the specific absorbance of lipid [47] contribute to the absorption coefficient of the tissue specifically due to lipid. The absorption due to hemoglobin in the brain tissue is dependent on cerebral blood volume [54] and was calculated based on the volume fraction of the blood in the cerebral tissue, hemoglobin oxygen saturation of mixed arterio-venous vasculature, and the absorption coefficient of pure oxy and deoxyhemoglobin [55]. The molar concentration of oxidized and reduced CCO (in mM) were first obtained for the gray and white matter [56]. The absorption coefficient of 1 mM of CCO was obtained [2] based on which the tissue absorption coefficient due to oxidized and reduced CCO was calculated. The absorption coefficients of gray and white matter due to each specific chromophore and the chromophore concentration are shown in Table 5. The calculations of the absorption coefficients are provided in the Supplementary Materials.
The whole tissue absorption coefficient is given in Table 6.

Finite Element Analysis
For the simulation of the RTE (Equation (1)) coupled with the bioheat transfer equation (Equation (8)), the Finite Element Analysis used the PDE toolbox of COMSOL to solve the equations based on discretization. In this case, the head model with the four layers has been discretized into more than 917,075 tetrahedral elements forming a complete mesh. The computation of the partial differential equations is performed at each discrete unit, more precisely, at each node of the tetrahedral element of the mesh. The approximation of the solution on the entire three-dimensional head model is performed by interpolation of the data in the space between the nodes using quadratic Lagrangian shape function in-built in the COMSOL software. The source position Cz, according to the 10-20 EEG system, has been chosen as the stimulation site. A point source of light with power density 500 mW/cm 2 adapted from the literature [1] was assumed at the Cz position in the head model. The source was a point source of light and was placed at the scalp at the Cz position, thus being in direct contact with the skin (Figure 6). The CAD model of the adult Colin27 human head model was used with each domain (layers) assigned the optical and bioheat parameters given in Tables 2 and 3. For simulating the light interaction due to individual chromophores, we investigated the absorption in the brain tissues due to specific chromophores (based on its absorption coefficient in the brain tissue). Since scattering is a property attributed by the geometry of the medium (i.e., the brain tissue), we have assumed the same reduced scattering coefficient of the tissue during all the chromophore-specific simulation. We used both absorption and scattering properties of the skull and scalp and the CSF since we wanted to study the fluence rate at the brain tissues after the light traveled through the scalp, skull, and CSF. Thus, we initially simulated for the brain tissues' contribution to light absorption by taking the lumped or total (due to all chromophores) absorption coefficients of the gray and white matter. The results for this specific simulation are presented as 'Whole Tissue'. Please note that the attenuation coefficient is the sum of the absorption coefficient and the scattering coefficient. Here, scattered light fluence rate is expected to be a constant (for all chromophores in the tissue) while the fluence attenuation is due to all the chromophores (absorption coefficients lumped in the tissue). To determine the attenuation due to individual chromophores of interest, we kept the tissue scattering the same (as in 'Whole Tissue' simulation) while determining the fluence rate attenuation due to individual chromophores in the tissue. The results were all plotted along a line cut through Cz (Figure 7) crossing the layers, thereby, depicting the straight path of light into the tissue.

Source
The cutline has been drawn through the source which has coordinates 92 mm, 104 mm, and 174 mm. The x-axis on the graphs show the values of z-coordinates of the points along the cutline. Thus, values of z-coordinates decrease as the light travels further from the source placed at the scalp surface(the grid of coordinates shown in Figure 6 where the z-axis is the depth). Variation of parameters have been shown in figures till the white matter at z-coordinate 92.45mm to better understand the penetration ability of light as it reaches the gray and white matter during tNIRS.

Photothermal Effects
The optical fluence rate due to absorption and scattering by each layer is obtained from the solution of the RTE Equation (7). Figure 8 shows the normalized (natural logarithm) fluence rate for the 'Whole Tissue' in the tissue layers for the three wavelengths used for the study along the straight line taken through Cz. It is seen that for wavelength 700 nm and 810 nm, fluence rate is comparatively higher at greater depths (less attenuation) when compared to that at 630 nm, i.e., a higher penetration depth near the NIR optical window. The fractional fluence rate from the scalp surface to gray matter is shown in Figure 9 where we see that minimal amount of light penetrates from the scalp through the skull and cerebrospinal fluid to the gray matter across a distance of more than 20 mm along the cutline (Figure 7) showing that minimal amount-around 0.2% NIR light is able to penetrate the skull.    It was seen that 810 nm comparatively shows a higher absorption of power at the gray matter, and thus we hypothesized that this wavelength a better choice for photothermal neuromodulation. We performed the bioheat simulation for all three wavelengths to verify our hypothesis.
The temperature along the line at the Cz location (10-20 EEG system) at different domains of the head model due to the 630 nm, 700 nm, and 810 nm optical stimulation was obtained from bioheat transfer solution, as shown in Figure 11.
The results showed a temperature rise, at the scalp surface as well as at the other layers, from the average body temperature of 37 • C. The increase of temperature at the scalp at Cz is less than 0.25 • C so well within the safety limit ( Figure 11), but the rise of temperature at the gray matter underlying Cz area was much lower less than 0.04 • C ( Figure 12).  The temperature plotted over the volume of gray and white matter and represented through color map further elucidates the temperature distribution over the entire volume of the two domains, as shown in Figure 13.
In Figure 13, it can be seen that at all the wavelengths, there is no considerable increase in temperature in the gray and white matter and temperature is very close to the average body temperature. In both cases, the photothermal effect leading to changes in neural excitability is not expected at such a small change in temperature. Hence, we investigated the other aspect of light interaction with the neural tissue, i.e., photobiomodulation.

Photobiomodulation
The light interaction with the chromophores in the gray and white matter (see Table 5) was performed to analyze how the three wavelengths (630 nm, 700 nm, and 810 nm) in different spectral regions are absorbed in the brain tissues that can lead to photobiomodulation. The fluence rate has been computed along the cutline shown in Figure 7.
The fluence rate distribution due to absorption by specific chromophores (at the gray and white matter) through the different layers of the head model is shown in Figures 14-19. Figures 14, 16 and 18. These show how the fluence rate varied in the layers along the surface normal through Cz. Our simulations support the prior findings that around 1-5% of light in the NIR spectral region reaches the gray matter [1]. We plotted the curve by taking the logarithm of the computed absolute fluence rate across different layers to better visualize the attenuation for each chromophore. An important finding when comparing between 630 nm, 700 nm, and 810 nm is that cytochrome c oxidase contributes significantly in the attenuation of light at the gray matter, thus showing lower flux. Moreover, 810 nm was found to have better depth penetration (Figure 9), and so this wavelength is seen to be more promising than 630 nm and 700 nm for non-invasive brain stimulation.

Systematic Model Errors
The simulated model shown here is a realistic head model based on the Colin27 head atlas. The mathematical model could show errors, mainly due to simplifications. The errors may be as follows: • In our model, the brain has been assumed as a highly scattering medium which is not true for CSF which is a low scattering medium where RTE can produce erroneous results [57] (Comparison between RTE and Monte Carlo has been shown in Supplementary Material).

•
Another error is related to computational limitation during the FEM modeling and discretization. There were limitations of accessible memory. Although enhancing resolution leads to better convergence of FEM results [58], but computational limitation restricted us to the standard COMSOL mesh refining process.

•
The reflection effects due to light interaction have been excluded from the simulations to focus on light interaction with tissues due to absorption and scattering.

•
The optical properties of the tissues in the head model vary significantly based on prior works. We selected a set of optical parameters from review literature [39,49,56]. We did not consider chromophores in the skin such as melanin, lipofuscin.

•
The simulation of the Bioheat Transfer assumed that the heat loss at the skin surface is due to convection and radiative heat loss was considered insignificant at that temperature.

Discussion and Conclusions
The computational pipeline aimed to investigate the temperature change induced by the light absorption at the three wavelengths, 630 nm, 700 nm, and 810 nm, as well as the absorption by the chromophores in the neural tissue. In this multiphysics modeling of light diffusion with bioheat transfer, we found that the temperature change in the scalp is well within 1 degree Celsius as reported by Chaieb and colleagues [1] for a light source of power density 500 mW/cm 2 at the scalp surface. As the light gets attenuated while propagating through the skull and cerebrospinal fluid (CSF) to reach the gray matter (0.2% reaches the gray matter; hence less than 1%), the low fluence rate leads to insignificant heating in the gray matter. Here, we assumed the initial body temperature at 37 • C, and the temperature increase at the Cz area was found to be around 0.033 • C at the gray matter-CSF interface. The sharp decrease in the fluence rate as the light propagated further into the gray matter is shown in Figure 12.
Prior works on brain temperature that was assessed on resting clinical patients showed an average brain temperature ranging around 36.9 ± 0.4 • C [59]. Brain activity has been shown to be associated with the rise in brain temperature. Studies have suggested that temperature changes of even less than 1 • C can result in functional alterations in the various areas of the nervous system [60], indicating the high thermal sensitivity of the brain. Thus, the temperature is an important active and dynamic variable that can modulate brain activity and needs to be monitored during stimulation. The brain, being at a typically higher temperature than the body, is cooled down by perfusing blood, which was considered in our modeling of thermal changes through bioheat transfer. Here, blood perfusion acts as a heat sink, thus cooling the brain down. Therefore, the temperature change in biological tissue is significantly dependent on the bioheat, which is further dependent on the heat source and the blood perfusion sink. The simulated results showed insignificant temperature change (0.033 • C) to cause photothermal neuromodulation. Hence, the chromophore simulations suggest a possible photobiomodulation effect of the NIR light interaction with the tissue. The results obtained from the simulation of the absorption by each chromophore, including lipid and water, elucidated the fact that besides water and lipid, light attenuation in the gray matter is due to the absorption of NIR light by the reduced and oxidized CCO. In fact, for 630 nm, 700 nm, and 810 nm wavelengths, we found that the two forms of CCO are the two major contributors to light attenuation besides water, lipid, and hemoglobin. Thus, we can conclude from the three categories of data that neuromodulation of the gray matter by photothermal effect is not significant with 500 mW cm −2 at the scalp surface at 630 nm and 700 nm (red spectral region) and 810 nm (near-infrared spectral region). However, the biochemical effects of CCO absorption need further investigation in conjunction with the heating effects since a small, steady state temperature change can affect the kinetics of photobiomodulation. Our simulation data comparing the fluence rate attenuation among 630 nm, 700 nm, and 810 nm also showed that 810 nm has higher penetration depth than the 630 nm and 700 nm, which supports the use of tNIRS for non-invasive brain stimulation.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3425/9/8/179/s1, Figure S1: Comparison between Monte Carlo and Diffusion Approximation. Acknowledgments: We thank the reviewers for their insightful comments and suggestions that improved the work presented in the paper. We also thank Steffy Rodriguez (undergraduate student of Biomedical Engineering at the University at Buffalo) for proofreading the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: