Simulating Light Diffusion in Human Brain Tissues Using Monte-Carlo Simulation and Diffusion Equation

Medical diagnosis with optical techniques is favorable due to its safe and painless features. Every tissue type can be distinguished by its optical absorption and scattering properties that are related to many physiological changes and considered to be very important signs for tissue heath. Characterizing light propagation in the human brain tissues is a vital issue in many diagnostic and therapeutic applications. In this work, light propagation in different brain tissues in normal and coagulated state was investigated. A Monte-Carlo simulation model was implemented to obtain spatially resolved steady state diffuse reflectance profiles of the examined tissues. Furthermore, the diffusion equation was solved to create images presenting the optical fluence rate distribution at the tissue surface using the finite element method. The proposed diffuse reflectance curves and fluence rate images show different features regarding tissue type and condition that promises to be effective in medical diagnosis.


Introduction
Light propagation in any biological tissue is described by its optical parameters. The primary optical parameters of a tissue are the absorption coefficient, scattering coefficient and the anisotropy, however, these parameters are highly dependent on the wavelength of the illuminating light [1]. The examined tissue is excited by an appropriate light then, some detector measurements are recorded. Furthermore, for optical parameters estimation, a suitable mathematical model is utilized to relate the experimental measurements with the unknown optical parameters [2]. The possible experimental measurements are the diffuse reflectance, the diffuse transmittance and the collimated transmittance. The recorded light is called diffuse due to the turbid nature of the biological tissues [3].
Different mathematical models have been applied to represent the propagation of light in biological tissues, these models classified into forward and inverse models. In forward models, a known optical parameters tissue is examined to obtain the detector recordings. While, the experimental measurements is introduced to the inverse models in order to predict the unknown optical parameters [4]. Monte-Carlo simulation [5], diffusion equation [6] and adding doubling [3,7] method are examples of the forward models. While, Kubelka-Munk [8], inverse Monte-Carlo [9] and invers adding doubling [10] are inverse models. Monte-Carlo model is a very common numerical model used to simulate the light propagation in turbid multilayer media [11] and during photodynamic therapy [12]. This model assumes an infinitely narrow pencil beam that is perpendicularly incident on a multi-layered scattering medium, each layer has different absorption coefficient, scattering coefficient and anisotropy [13]. Monte-Carlo simulation has been implemented to detect abnormalities in tissues at different locations [14] and also to model light propagation in Coral tissue [15].
The diffusion approximation to the radative transport equation is usually applied to illustrate the light propagation in biological tissues [16], the diffusion approximation simplify the required mathematics assuming some boundary conditions to abridge the difficulties in calculations resulting from the in-homogeneities of biological tissues [17]. The results obtained by the diffusion equation usually compared with Monte-Carlo simulation for verifications [18].
Investigating light behaviour and propagation in human brain is considered an important subject in many diagnosis and therapeutic processes regarding the brain activities and physiological changes such as Oxygenation and Haemoglobin concentrations [19], this behaviour depends on the optical absorption and scattering parameters of examined tissue. Optical parameters of variety of mammalian brain tissues have been estimated using both ex vivo and in vivo measurement techniques [20]. Furthermore, characterizing the optical properties of brain tissue with high numerical aperture optical coherence tomography has been introduced in [21], the obtained optical properties provided unique information to differentially identify laminar structures in primary visual cortex and various nuclei in the midbrain.
In this work, Monte-Carlo Model of Light propagation in biological tissues 'MCML' was implemented to investigate some known optical parameters human brain tissues including cerebellum, thalamus, gray-matter and whitematter in normal and coagulated condition. The output of the simulation obtain spatially resolved steady state diffuse reflectance profiles of the examined tissues. Furthermore, the optical parameters were then introduced to the diffusion equation to obtain fluence rate distribution images at the surface of the sample using finite element method.

Methods
The Monte-Carlo Model for Light Transport in biological tissue has been implemented to simulate the spatially resolved steady state diffuse reflectance of different brain tissue types, the diffuse reflectance profile is considered a figure print for each tissue type. Moreover, diffusion approximation of the radative transport equation has been employed to create images presenting the fluence rate distribution at the tissue surface according to the different optical parameters values of each examined sample.

Monte-Carlo Modelling for Light Propagation in Biological Tissues
Monte Carlo model simulates the propagation of photon in turbid medium basing on the random walk concept. It is a numerical approach to the radative transport equation. In the model, a photon bundle is tracked over the examined tissue till it gets completely absorbed or exits the tissue boundary. Monte Carlo simulation can be made through eight main steps: (1) launching of the photons, (2) Photon's step size, (3) Moving the photon, (4) Internal Reflection, (5) Photon absorption, (6) Photon Scattering, (7) Boundary crossing of a photon, (8) Photon Termination [12].
After initializing the photon bundle, each photon moves a distance ∆ where it may be scattered, absorbed, propagated without interruption, internally reflected, or exit the tissue through reflection and/or transmission. If the photon left the tissue boundary, the amount of reflectance or transmittance is documented. In case of photon absorption, the absorption location is recorded. It is a repetitive process till tracing the pre-desired number of photons [13].
A typical Monte-Carlo model assumes an infinitely narrow photon beam, perpendicularly incident on a multi-layered turbid medium. Each layer is defined by its thickness t, refractive index n, absorption coefficient µa, scattering coefficient µs and anisotropy ց in layersmatrix array. The two basic sentences of the simulation are: Layersmatrix = [n µa µs ց t]; (1) S=MCML(filename,number_of_photons, Layersmatrix) (2) where, S is a MatLab structure that contains all simulation results the simulation input parameters as well.
The function that used in the simulation process is called MCML [22], the previously mentioned five parameters regarding the sample have to be known to use the simulation. In our simulation problem, the obtained optical parameters and thickness of the samples have been introduced to MCML assuming matched boundary conditions to obtain the diffuse reflectance profile at each case. All simulations were implemented under Matlab R2015a environment.

Diffusion Equation
The optical fluence rate distribution at the external surface of the tissue sample can be calculated from the diffusion equation [1].
where the constant ̀ is the diffusion coefficient and ̀ 1 " is the reduced scattering coefficient, " is the anisotropy factor.
, represents the source term and Ф , is the fluence rate.
The diffusion equation was solved using the Finite Element "FE" method to obtain fluence rate distribution images at the surface of the examined samples. A model of 2 cm width and 2 cm height square was created. Then, a point representing the laser source was located at (0.2, 1) to simulate the distant-detector experimental conditions [8]. The geometry of the model, and the produced finite element mesh is presented in Figure 1. Some known optical parameters tissues from different parts of the human brain were investigated [23]. Table 1.
Illustrates the optical parameters values used in the implemented Monte-Carlo simulation and in the diffusion equation as well, the all values were calculated at 800 nm laser irradiation.

Results and Discussion
From Monte-Carlo simulation output, the spatially resolved steady state diffuse reflectance profiles of cerebellum, thalamus gray-matter and white-matter in normal and coagulated were obtained using the optical parameters values presented in Table 1, these profiles are introduced in Figure 2. The curves present the diffuse reflectance as a function of the distance from the light source r (mm), therefore, it is called spatially resolved reflectance.
As conducted from Figure 2, the four tissue types showed observable higher values of the diffuse reflectance in coagulated rather than the normal state for cerebellum, thalamus and gray matter. However, the opposite occurs in normal and coagulated white matter. Though, both the normal and coagulated forms of the tissues show an exponential shape diffuse reflectance profile with different maximum values and very nearby minimum values.
The profile of the diffuse reflectance in biological tissues depends on the incident light wavelength. However, changing sample condition can cause many physiological changes regarding water and blood contents that typically affect absorption and scattering properties of the tissue. The same optical parameters values were then utilized in the diffusion equation to determine the fluence rate distribution at the surface of the examined tissues using finite element method. Figure 3 presents the optical fluence rate distribution images for each tissue type in normal and coagulated state. The fluence rate distribution at the surface of the examined tissues changes with the tissue type and condition due to the change in their optical properties as illustrated in Figure 3. It can be observed from the Figure that, the maximum value of log Ф in normal cerebellum was 5.29 and the minimum value was -28.5, while in coagulated cerebellum, the maximum value of log Ф changed to 5.9 and the minimum value was -76.5. In normal thalamus the maximum log Ф was 4.38 and the minimum value was -12.5, while in coagulated thalamus, log Ф maximum value was 4.55 and the minimum value was -17.2. Moreover, the maximum value of log Ф in normal gray matter was 3.76 while in the coagulated one was 5.08 and the minimum value was -5.24 for normal gray matter while the values to -17.2 in coagulated one. Finally, in normal white matter the maximum value 5.1 and the minimum value was -21.9 while, in coagulated white matter the maximum value was 5.15 and the minimum was -21.3.
The presented fluence rate images show observable differences between normal and coagulated cerebellum, thalamus and gray matter, however, in white matter the differences are not strong but still noticeable, and this can be acceptable as the optical parameters values in normal and coagulated white matter are very close as provided from

Conclusion
In conclusion, light propagation in different brain tissues in normal and coagulated form have been investigated using Monte-Carlo simulation and diffusion equation. With Monte-Carlo model, spatially resolved steady state diffuse reflectance profiles of the examined tissues has been obtained showing a maximum reflectance value in the coagulated state greater than in the normal state for cerebellum, thalamus and gray-matter, while, the opposite was occurred in case of while-mater. Images representing the optical fluence distribution at the surface of each examined sample have been obtained using the finite element solution of the diffusion equation. The resultant images provided some observable discriminations between normal and coagulated samples that make it suitable method for tissue characterization and differentiation. The proposed results are supposed to help in medical diagnosis and therapeutic procedures through many medical applications such as biostimulation and photodynamic therapy.