Optical Properties and Fluence Distribution in Rabbit Head Tissues at Selected Laser Wavelengths

Abstract

The accurate estimation of skin and skull optical properties over a wide wavelength range of laser radiation has great importance in optogenetics and other related applications. In the present work, using the Kubelka–Munk model, finite-element solution of the diffusion equation, inverse adding-doubling (IAD), and Monte-Carlo simulation, we estimated the refractive index, absorption and scattering coefficients, penetration depth, and the optical fluence distribution in rabbit head tissues ex vivo, after dividing the heads into three types of tissues with an average thickness of skin of 1.1 mm, skull of 1 mm, and brain of 3 mm. The total diffuse reflectance and transmittance were measured using a single integrating sphere optical setup for laser radiation of 532, 660, 785, and 980 nm. The calculated optical properties were then applied to the diffusion equation to compute the optical fluence rate distribution at the boundary of the samples using the finite element method. Monte-Carlo simulation was implemented for estimating the optical fluence distribution through a model containing the three tissue layers. The scattering coefficient decreased at longer wavelengths, leading to an increase in optical fluence inside the tissue samples, indicating a higher penetration depth, especially at 980 nm. In general, the obtained results show good agreement with relevant literature.

1. Introduction

The use of lasers in biomedical diagnostics and therapy requires spectroscopic measurements of the beams that propagate through the examined tissue or organ. Optical diagnosis methods have a high impact because they are safe, minimally invasive, and non-destructive. In neurology and brain examination, light travels through the scalp, skull, and brain. Despite the important effects of the upper tissue layers on measuring brain hemodynamics, the effect of the scalp and skull on the sensitivity and intensity of laser radiation has not been well characterized. Innovation headways in photonics have made massive advancements toward the development of creative strategies and systems for clinical practical optical imaging, laser surgery procedures, biomedical diagnostics, and phototherapy. The progress of optical biomedical techniques and methods focuses on the investigation of the optical properties of tissues, which characterize the viability of tissue optical testing and light activity on the tissue and, when known (estimated), allow for predicting the exact photon spread directions and the fluence rate distribution inside irradiated tissues [1]. Tissue removal, coagulation, laser cutting, and many other applications depend on the spectroscopic properties of tissues for their accuracy. Therefore, information regarding the optical properties of tissues is essential for the interpretation and evaluation of the diagnostic data and for the prediction of the absorbed light and energy distribution for therapeutic and surgical use. Various publications connected with the assurance of brain tissue optical properties and demonstrating the fluence rate distribution are accessible in the literature [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. The examination of these publications shows that the optical properties of human and different animal brain tissues have been explored in the visible and NIR spectral ranges. However, these analyses exhibit an assortment of cutting edge biomedical advances that require exact information on the optical properties of the tissues [10,11,12,13,14,15,16].

The optical properties of tissues that assume a significant role in tissue characterization are the absorption coefficients (μ_a), scattering coefficients (μ_s), and the reduced scattering coefficients (μ_s′), including the refractive index (RI) and the anisotropy factor (g). These parameters are high-wavelength-dependent and give practical data about the tissue, for example, the hemoglobin content, tissue oxygenation, and water fraction. However, assessing the optical properties of any tissue requires estimations regarding the diffuse reflectance (R_d), diffuse transmittance (T_d), and collimated transmittance (T_c). Experimentally, there are several methods to obtain these estimates, either by integrating spheres or using techniques in light with distant detector arrays or with other devices. Optical differentiation is a useful tool in biomedical diagnosis, for the most part, as a result of its safety. As per histopathology, the values of the tissue optical properties will contrast for different tissues and consequently could be utilized for differentiation of the norm and pathology. The optical fluence rate distribution inside the tissue limits depends on the optical boundaries. Therefore, presenting the parameters of such distributions can give an optical method for biomedical diagnosis [17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]. This paper presents the results of the measurements and estimation of tissue optical properties, where an experimental setup was implemented to measure the diffuse optical reflectance and transmittance of the ex vivo samples of the head skin, skull, and brain at 532, 660, 785, and 980 nm laser wavelengths. Using the measured values, the optical parameters of the samples were determined utilizing an inverse adding-doubling (IAD) method and the Kubelka–Munk model. The assessed values of the optical parameters were utilized in the diffusion equation to simulate the fluence rate at the tissue surface using the finite element method. The results were confirmed with the Monte-Carlo simulation, which included modeling of the in-depth fluence rate distribution for the three-layer head tissue.

2. Materials and Methods

2.1. Sample Preparation

Ten adult rabbits were acquired from a local breeder who sold rabbits for consumption. Then, the animals were sacrificed one at a time on different days. All of the animal studies were approved by the Animal Care and Experimentation Committee of the university. Skin (scalp), skull, and brain samples were manually removed from the heads using appropriate scalpels for ex vivo investigation. The average thickness of the skin was 1.1 mm, skull was 1 mm, and brain was 3 mm. For the RI measurements, several samples of rabbit brain tissue were prepared with a 3 mm thickness. For the spectral measurements, 10 samples were prepared with a 0.5 mm thickness. The skin samples were shaved for accurate measurements (see Figure 1b). Transmittance and reflectance and measurements were done in the fresh native state and the rest of the samples were placed in the refrigerator at −20 °C. The studied samples were positioned between two 1 mm thick microscope glass slides for the spectroscopic measurements.

2.2. Experimental Arrangements

An optical setup based on an integrating sphere (McPherson, KS, USA) coupled with a commercial digital fiber spectrometer (STDFSM, Touptek Photonics Co., Ltd., Hangzhou, China) was utilized for measuring the tissue transmittance and reflectance. The sphere had two output ports with diameters of 25 and 15 mm, while the diameter of the radiation inlet port was 10 mm. The measurements were obtained using four CW laser sources (PGL-DF, Changchun New Industries Optoelectronics Tech. Co., Ltd., Changchun, China) at wavelengths of 532 nm (~100 mW power, 2 mm beam diameter, near TEM₀₀), 660 nm (~130 mW power, 3.5 mm beam diameter, near TEM₀₀), 785 nm (~126 mW power, 4 mm beam diameter, near TEM₀₀), and 980 nm (~108 mW power, 3 mm beam diameter, multimode) [22]. An analysis of the results was completed using the spectrometer software and Matlab R2018a program.

2.3. Tissue Diffuse Light Measurements

Tissue diffuse transmittance and reflectance were measured using a single integrating sphere and CCD (TCD1304AP) connected to STDFSM digital fiber spectrometer through a data collecting optical fiber used in the experiment. Two different geometrical arrangements were implemented for measuring the total transmittance T_t and diffuse reflectance R_d, as illustrated in Figure 2. Because there was no hole to allow the collimated beam to escape the integrating sphere, the collimated transmittance T_c was measured in a separate optical arrangement (outside the sphere) using a laser source, pinholes, and detector, as illustrated in Figure 2e. Accordingly, the diffuse transmittance T_d is calculated as T_d = T_t − T_c. Finally, the collected data were then applied in the Kubelka–Munk model for estimation of the absorption and scattering coefficients.

2.4. Refractive Index Measurements

The refractive index of the tissue samples was measured using the multi-wavelength Abbe refractometer (Atago, Tokyo, Japan) (Figure 3). The RI was measured for samples of rabbit head tissues ex vivo after dividing the head into three types of tissues with an average thickness of the skin of 1.1 mm, skull of 1 mm, and brain of 3 mm. The multi-wavelength refractometer Abbe permits measuring the RI in the wavelength range of 450 to 1550 nm with an exactness of 0.0002. The functioning principle of the refractometer technique depends on deciding the critical angle of the total reflectance, where the occurrence light waves are totally reflected at a 90-degree angle to the ordinary position. The total inside reflectance method was applied for estimation of the RI of biological tissue, which is described by high light scattering and absorption. We carried out the measurements to calculate the optical parameter RI of the rabbit head tissue at wavelengths of 532, 660, 785, and 980 nm.

2.5. Estimation of Optical Properties

Kubelka–Munk’s (KM) model permits two inside radiation fluxes through the tissue, one flux in the direction of the incident beam J₁ and the other in the backscattered direction J₂; the essential flux spreads in a similar direction as the incident radiation and the other flux engenders conversely, as displayed in Figure 4. Two coefficients ($A_{\mathrm{KM}}$ and $S_{\mathrm{KM}}$) are proposed to show the absorption and scattering of the diffuse radiation, individually [23,24,25,26,27,28,29]:

$$\frac{d{J}_1}{dz}=-S_{\mathrm{KM}}{J}_1-A_{\mathrm{KM}}{J}_1+S_{\mathrm{KM}}{J}_2,$$

(1)

$$\frac{d{J}_2}{dz}=-S_{\mathrm{KM}}{J}_2-A_{\mathrm{KM}}{J}_2+S_{\mathrm{KM}}{J}_1,$$

(2)

where z alludes to the primary direction of the incident radiation, as indicated by Kottler [23]. The Kubelka–Munk coefficients are connected so as to measure the diffuse transmittance $T_{\mathrm{d}}$ and diffuse reflectance $R_{\mathrm{d}}$ as [23,24]:

$$R_{\mathrm{d}}=\frac{\sinh (S_{\mathrm{KM}}yd)}{x\cosh \left({\mathrm{S}}_{\mathrm{KM}}yd\right)+y\sinh (S_{\mathrm{KM}}yd)},$$

(3)

$$T_{\mathrm{d}}=\frac{y}{x\cosh \left(S_{\mathrm{KM}}yd\right)+y\sinh (S_{\mathrm{KM}}yd)},$$

(4)

where $d$ is the optical thickness of the piece to be considered, and parameters $x$ and $y$ are expressed as follows:

$$x=\frac{1+R_{\mathrm{d}}^2-T_{\mathrm{d}}^2}{2R_{\mathrm{d}}},$$

(5)

$$y=\sqrt{x^2-1}.$$

(6)

The Kubelka–Munk model is a unique instance of the supposed multi-flux theory, where the transport equation is changed into a matrix differential equation by considering the radiance at numerous discrete angles [25]. The convenience of applying the Kubelka–Munk method arises from the fact that the absorption coefficient ${\mathsf{\mu }}_{\mathrm{a}}$ and reduced scattering coefficient ${\mathsf{\mu }}_{\mathrm{s}}^{\prime }$ can be directly calculated from the measured diffuse transmittance $T_{\mathrm{d}}$ and diffuse reflectance $R_{\mathrm{d}}$, i.e., as follows:

$$S_{\mathrm{KM}}=\frac{1}{yd}\ln \left[\frac{1-R_{\mathrm{d}}\left(x-y\right)}{T_{\mathrm{d}}}\right],$$

(7)

$$A_{\mathrm{KM}}=\left(x-1\right)S_{\mathrm{KM}}.$$

(8)

Then, the relation between $S_{\mathrm{KM}}$ and $A_{\mathrm{KM}}$ and the scattering ${\mathsf{\mu }}_{\mathrm{s}}$ and absorption ${\mathsf{\mu }}_{\mathrm{a}}$ coefficients of the sample can be expressed as follows:

$$A_{\mathrm{KM}}=2{\mathsf{\mu }}_{\mathrm{a}},S_{\mathrm{KM}}=\frac{3}{4}{\mathsf{\mu }}_{\mathrm{s}}\left(1-\mathrm{g}\right)-\frac{1}{4}{\mathsf{\mu }}_{\mathrm{a}},$$

(9)

where g represents the scattering anisotropy factor and ${\mathsf{\mu }}_{\mathrm{s}}^{\prime }={\mathsf{\mu }}_{\mathrm{s}}\left(1-\mathrm{g}\right)$.

In addition, the optical penetration depth $\mathsf{\delta }$ can be estimated using the following equation, which is valid in the diffusion approximation: [28]

$$\mathsf{\delta }=\frac{1}{\sqrt{3{\mathsf{\mu }}_{\mathrm{a}}\left({\mathsf{\mu }}_{\mathrm{a}}+{\mathsf{\mu }}_{\mathrm{s}}^{\prime }\right)}}.$$

(10)

Regarding the determination of the refractive index (RI) of the studied tissues, the dispersion is represented as follows:

$$n^{*}\left(\mathsf{\lambda }\right)=n\left(\mathsf{\lambda }\right)+ik\left(\mathsf{\lambda }\right),$$

(11)

where $n\left(\mathsf{\lambda }\right)$ is the real part of the complex RI $n^{*}\left(\mathsf{\lambda }\right)$ and $k\left(\mathsf{\lambda }\right)$ is its imaginary part which is related to the tissue absorption coefficient by the Kramers–Kronig relations, as follows [29]:

$$k\left(\mathsf{\lambda }\right)=\frac{\mathsf{\lambda }{\mathsf{\mu }}_{\mathrm{a}}\left(\mathsf{\lambda }\right)}{4\mathsf{\pi }}.$$

(12)

Additionally, the Fresnel equation connects the specular reflectance at a normal incidence $R\left(\mathsf{\lambda }\right)$ and the tissue RI as follows [30,31,32,33,34,35,36,37,38,39]:

$${\mathrm{R}}_{\mathrm{F}}\left(\mathsf{\lambda }\right)=\frac{{\left(\mathrm{n}\left(\mathsf{\lambda }\right)-1\right)}^2+\mathrm{k}{\left(\mathsf{\lambda }\right)}^2}{{\left(\mathrm{n}\left(\mathsf{\lambda }\right)+1\right)}^2+\mathrm{k}{\left(\mathsf{\lambda }\right)}^2}.$$

(13)

In general, the RI of the tissues samples is measured at discrete wavelengths using multi-wavelength Abbe refractometers [40] or by utilizing the total internal reflectance method with different lasers [41]. When the RI values are estimated, the tissue dispersion for that wavelength range is calculated by apply the experimental RI data with equations such as the Cauchy (Equation (14)), the Conrady (Equation (15)), or the Cornu (Equation (16)) equations [41,42,43]:

$$n_{\mathrm{tissue}}\left(\mathsf{\lambda }\right)=\mathrm{A}+\frac{\mathrm{B}}{{\mathsf{\lambda }}^2}+\frac{\mathrm{C}}{{\mathsf{\lambda }}^4},$$

(14)

$$n_{\mathrm{tissue}}\left(\mathsf{\lambda }\right)=\mathrm{A}+\frac{\mathrm{B}}{\mathsf{\lambda }}+\frac{C}{{\mathsf{\lambda }}^{3.5}},$$

(15)

$$n_{\mathrm{tissue}}\left(\mathsf{\lambda }\right)=\mathrm{A}+\frac{\mathrm{B}}{\left(\mathsf{\lambda }-\mathrm{C}\right)},$$

(16)

where, A, B, and C are the Cauchy, the Conrady, and the Cornu parameters, respectively. Using the Kramers–Kronig relations, which were developed for non-scattering materials, a broader dispersion can be calculated from μ_a(λ) if the calculated tissue dispersion is not available for the entire wavelength range of interest. Equation (12) is the first one to be used to obtain the imaginary part of tissue dispersion (κ(λ)) [41,44]. When known, κ(λ), the dispersion that corresponds to the real part of the refractive index (RI), can be calculated from [41,44]:

$$n_{\mathrm{tissue}}\left(\mathsf{\lambda }\right)=1+\frac{2}{\mathsf{\pi }}{{\displaystyle \int }}_0^{\infty }\frac{{\mathsf{\lambda }}_1}{\mathsf{\Lambda }}\frac{{\mathsf{\lambda }}_1}{{\mathsf{\Lambda }}^2-{\mathsf{\lambda }}_1^2}K\left(\mathsf{\Lambda }\right)\mathrm{d}\mathsf{\Lambda }$$

(17)

where Λ is the integrating variable over the wavelength space and λ₁ is a decent wavelength that can be adapted for better matching of the determined scattering to the one obtained from the discrete experimental information. When the expansive band tissue scattering is determined through Equations (12) and (17), we can choose discrete values from it to use in the inverse adding-doubling (IAD) simulations.

Running (IAD) simulations will generate ${\mathsf{\mu }}_s^{\prime }$ values are then given with the equation [28,42]:

$${\mathsf{\mu }}_s^{\prime }\left(\mathsf{\lambda }\right)=a\left[f_{\mathrm{Ray}}{\left(\frac{\mathsf{\lambda }}{500\mathrm{nm}}\right)}^{-4}+\left(1-f_{\mathrm{Ray}}\right){\left(\frac{\mathsf{\lambda }}{500\mathrm{nm}}\right)}^{-b_{Mie}}\right],$$

(18)

where a is the scaling factor that addresses the value of ${\mathsf{\mu }}_s^{\prime }$ at 500 nm, f_Ray is the Rayleigh scattering fraction, and b_Mie represents the mean size of the Mie scatterers. These parameters can be achieved when the discrete ${\mathsf{\mu }}_s^{\prime }$ values that were produced through IAD are given, using Equation (18).

Such a spectrum can be joined with the μ_s spectrum that is determined using the measurements of collimated transmittance $T_{\mathrm{c}}\left(\mathsf{\lambda }\right)$ = T_t − T_d and thickness d [28,42]:

$${\mathsf{\mu }}_{\mathrm{s}}\left(\mathsf{\lambda }\right)=\frac{\ln \left[T_{\mathrm{c}}\left(\mathsf{\lambda }\right)\right]}{d}-{\mathsf{\mu }}_{\mathrm{a}}\left(\mathsf{\lambda }\right).$$

(19)

To obtain the wavelength reliance of the tissue dispersing anisotropy factor (g(λ)), both scattering coefficients are used: [45]

$$g\left(\mathsf{\lambda }\right)=1-\frac{{\mathsf{\mu }}_s^{\prime }\left(\mathsf{\lambda }\right)}{{\mathsf{\mu }}_s\left(\mathsf{\lambda }\right)}$$

(20)

In general, g(λ) gives a rising exponential behaviour expanding wavelength in the UV–NIR range [1], which can be described mathematically using Equation (21) [45] or Equation (22) [46].

$$g\left(\mathsf{\lambda }\right)=\mathrm{a}+\mathrm{b}\left[1-\exp \left(\frac{\mathsf{\lambda }-\mathrm{c}}{\mathrm{d}}\right)\right],$$

(21)

$$g\left(\mathsf{\lambda }\right)=\mathrm{a}\exp \left(\mathrm{b}\mathsf{\lambda }\right)+\mathrm{c}\exp \left(\mathrm{d}\mathsf{\lambda }\right),$$

(22)

where a, b, c, and d are the parameters that are acquired during the gives of g data. This estimation system depends just on IAD simulations to obtain discrete (${\mathsf{\mu }}_{\mathrm{s}}^{\prime }$) values [47].

2.6. Monte-Carlo Simulation for Modeling of Light Diffusion

Monte-Carlo (MC) is a very common computational modeling algorithm used to simulate light spread in one- or multi-layer biological slabs [31,32,33]. As an advanced method, the MC execution needs input data regarding the analyzed piece. Consequently, each layer is defined by its RI, scattering anisotropy factor, scattering coefficient, absorption coefficient, and thickness. The simulation begins by submitting the fitting number of photons and the tracing process starts as per the predefined step size and limit conditions. The MC simulation process is summarized in the following flowchart (Figure 5).

To use the simulation, five parameters with respect to the sample must be known. The simulation output gives the worth of the diffuse reflectance that ought to be obtained from the experimental estimations. In our experimental, the achieved thickness and optical parameters of the samples were introduced to MC Multi-Layered (MCML) [33], expecting matched boundary conditions in order to achieve the values of the diffuse reflectance in each case, hence validating our results. The simulation was applied using the Matlab R2018a program.

Where all the measurements of the parameters that we obtained from our experiment were entered in the Matlab (2018) program, and after applying the Monte Carlo program code, the photons were distributed—approximately 100,000 photons in the tissues—which were then measured after interacting with them, according to the data and the angles of deviation in the photon path when scattering occurred. This statistic depends on calculating the spread of a large number of photons. As a result, this method requires a long computational time.

2.7. Simulating Fluence Distribution

The distribution of the optical fluence of the head tissue layers was determined by utilizing the finite element solution of the diffusion equation. The fluence rate is determined from the following equation [26]:

$$\frac{\partial \mathsf{\Phi }\left(\stackrel{\rightharpoonup }{r},t\right)}{c\partial t}+{\mathsf{\mu }}_{\mathrm{a}}\mathsf{\Phi }\left(\stackrel{\rightharpoonup }{r},t\right)-\nabla \left[D\nabla \mathsf{\Phi }\left(\stackrel{\rightharpoonup }{r},t\right)\right]=S\left(\stackrel{\rightharpoonup }{r},t\right),$$

(23)

where $D=\frac{1}{3\left({\mathsf{\mu }}_{\mathrm{a}}+{\mathsf{\mu }}_{\mathrm{s}}^{\prime }\right)}$ is the tissue diffusion coefficient, $S\left(\stackrel{\rightharpoonup }{r},t\right)$ represents the source term, and $\mathsf{\Phi }\left(\stackrel{\rightharpoonup }{r},t\right)$ is the fluence rate.

To solve diffusion equations numerically, which is a forward problem, the finite element method (FEM) can be used [26]. In this case, the fluence distribution $\mathsf{\Phi }\left(\stackrel{\rightharpoonup }{r},t\right)$ is gained in the domain as a function F of optical properties μ(r), as follows:

$$\mathsf{\Phi }=F\left(\mathsf{\mu }\right),$$

(24)

where $\mathsf{\mu }=\left[{\mathsf{\mu }}_{\mathrm{a}},D\right]$.

In diffuse optical imaging, the image reconstruction process requires the solution of the inverse problem in which the distribution of $\mathsf{\Phi }\left(\stackrel{\rightharpoonup }{r},t\right)$ at the boundary is given, while optical properties of the domain are unknown [27]; this can be represented by the following:

$$\mathsf{\mu }=F^{-1}\left(\mathsf{\Phi }\right)$$

(25)

Solving the inverse problem requires minimization of the error function χ², which can be calculated as follows:

$${\mathsf{\chi }}^2=\frac{1}{2}{\displaystyle \sum }_{J=1}^M{\left({\mathsf{\Phi }}_{\mathrm{meas}}-\mathrm{F}\left(\mathsf{\mu }\right)\right)}^2,$$

(26)

where ${\mathsf{\Phi }}_{\mathrm{meas}}$ is the fluence measurements at the boundaries and $\mathrm{F}\left(\mathsf{\mu }\right)$ is the calculated fluence (${\mathsf{\Phi }}_{\mathrm{calc}}$) using the forward model. In the present work, the finite element method (FEM) using the diffusion equation was implemented to get up fluence rate distribution images at the boundary of the sample surfaces using Matlab R2018b and COMSOL Multiphysics 5.4 program software.

In COMSOL Multiphysics, the diffusion of Equation (23) in the steady-state can be presented using the Helmholtz equation:

$$\nabla \left(-c\nabla u\right)+au=f$$

(27)

Identifying the parameters from Equation (23) yields the following:

$$\begin{array}{c}\mathrm{u}=\mathsf{\Phi },\mathrm{a}={\mathsf{\mu }}_{\mathrm{a}},\\ \mathrm{c}=\mathrm{D},\mathrm{f}=\mathrm{S}.\end{array}$$

(28)

The executed model was a rectangle of 3 cm × 2 cm in width. A point source addressing the laser source was placed at a location of (1.5, 1) to simulate the laser source position in the integrating sphere configuration, as illustrated in Figure 6.

3. Results and Discussion

3.1. Determination of Scattering and Absorption Coefficients

The tissue sample total transmittance T_t, collimated transmittance T_c, and diffuse reflectance R_d were measured using an integrating sphere and two-pinhole arrangement, as illustrated in Figure 2. For the RI measurements, several samples of rabbit brain tissue were prepared with a 3 mm thickness. For the spectral measurements, 10 samples were prepared with a 0.5 mm thickness. The results of the measurements for the freshly native rabbit brain at wavelengths of 532, 660, 785, and 980 nm are presented in Figure 7. Accordingly, reduced scattering and absorption coefficients were extracted for the freshly native rabbit brain, skin, and skull at wavelengths of 532, 660, 785, and 980 nm. The obtained values were compared with previously published data, as summarized in

Table 1. The obtained optical parameters were compared with the relevant data from the literature.

Literature	Tissue	λ, nm	Optical Parameters
			μa, cm−1	${\mathsf{\mu }}_{\mathit{s}}^{\prime }$, cm−1	g	n
Soleimanzad et al. [34]	Mice skull (ex vivo)	532	14 ± 2	25 ± 1.5	-	-
		660	5 ± 0.5	25 ± 2.5	-	-
		705	4.7 ± 0.7	22.9 ± 1.2	-	-
Firbank et al. [35]	Pig skull (ex vivo)	660	0.3 ± 0.025	23.7 ± 1.2	~0.93	-
		780	0.25 ± 0.025	18.5 ± 1.2	~0.93	-
Beek et al. [36]	Rabbit skin (in vitro)	630	0.94 ± 0.13	40 ± 2.2	0.812 ± 0.017	-
		790	0.7 ± 0.07	18.4 ± 0.5	0.940 ± 0.003	-
Pitzschke et al. [37]	Rabbit brain (in vivo)	635	1.1 ± 0.2	12.8 ± 0.1	0.92	1.37
		671	0.8 ± 0.2	10.7 ± 0.9	0.92	1.37
		808	0.5 ± 0.1	6.2 ± 0.5	0.92	1.37
Gonçalves et al. [29]	Rabbit brain (cerebral cortex, ex vivo)	401.4	~9.5	~20	~0.72	1.3850 ± 0.0053
		534.6	~3.9	~14	~0.79	1.3736 ± 0.0058
		626.6	~2.1	~11	~0.82	1.3680 ± 0.0045
		782.1	~0.9	~8	~0.84	1.3611 ± 0.0050
		820.8	~0.8	~7	~0.85	1.3597 ± 0.0066
		850	~0.8	~6	~0.85	1.3589 ± 0.0076
P. Sawosz et al. [39]	Human scalp	750	0.125	13.0	-	1.4
		850	0.180	13.0	-	1.4
	Human skull	750	0.082	12.0	-	1.4
		850	0.118	12.0	-	1.4
	Human CSF	750	0.025	10.0	-	1.4
		850	0.041	10.0	-	1.4
	Human brain	750	0.156	12.5	-	1.4
		850	0.200	12.5	-	1.4
This study	Rabbit skin (ex vivo)	532	4.4 ± 0.6	15.0 ± 0.2	0.80	1.3769 ± 0.002
		660	2.8 ± 0.4	14.1 ± 0.2	0.82	1.3691 ± 0.002
		785	1.5 ± 0.8	11.4 ± 0.3	0.85	1.3631 ± 0.003
		980	0.8 ± 0.6	6.0 ± 0.2	0.90	1.3551 ± 0.002
	Rabbit skull (ex vivo)	532	4.7 ± 1.5	16.0 ± 0.5	0.80	1.3784 ± 0.005
		660	3.0 ± 0.4	14.8 ± 0.1	0.81	1.3695 ± 0.001
		785	1.7 ± 0.1	11.5 ± 0.4	0.90	1.3637 ± 0.004
		980	0.9 ± 0.4	7.2 ± 0.2	0.90	1.3596 ± 0.002
	Rabbit brain (ex vivo)	532	4.0 ± 0.3	14.3 ± 0.9	0.80	1.3761 ± 0.002
		660	1.2 ± 0.1	10.5 ± 0.3	0.83	1.3676 ± 0.003
		785	1.0 ± 0.7	8.7 ± 0.2	0.85	1.3626 ± 0.002
		980	0.8 ± 0.6	5.9 ± 0.2	0.90	1.3541 ± 0.009

.

The data for the (RI) refractive index of all rabbit head tissue types measured independently are presented in Table 1, and Figure 8 shows the (RI) refractive index of the rabbit brain tissue with the standard deviation.

The data for the (RI) refractive index compared with the data from the literature are presented in Figure 9.

Figure 10 illustrates the variation in the calculated optical properties and penetration depth at the studied wavelengths.

In our implementation, we measured the anisotropy factor (g) of all of the tissue types. These data are presented in Figure 11.

The proportion of the contrast factor in Figure 11 differs slightly in the scalp compared with in the skull, as well as the brain for the rabbit head, as well as with the different wavelengths used in our experiment. Therefore, little variation was noted between them, while in Figure 12 which compared our data with data from the literature, there was a difference between this and the variance factor measured for [36] in contrast with [29], which was closer to our measurements.

In the present study, ex vivo optical parameters of the rabbit skin, skull, and brain were determined and compared with data from the literature (Figure 13). The parameters were obtained for four laser wavelengths ranging from visible to NIR optical regions. The proposed results showed similarly reduced scattering coefficients of the skull at 532 nm to that published by Soleimanzad et al. [34] for mice skull, and Firbank et al. [35] for pig skull at 660 and 780 nm, while there were some differences in the absorption coefficient values. Regarding the skin samples, the optical parameters obtained at 660 were a little higher than at 780 nm, which followed almost the same behavior compared with those obtained by Beek et al. [36] at 630 and 790 nm. Moreover, the obtained optical parameters of the rabbit brain were comparable to those obtained by Pitzschke et al. [37] at near wavelengths with some variations within the acceptable range. In principle, the sample preparation method, sample thickness, experimental technique, and the employed mathematical model could be reasons for the relatively spread values for μ_a and ${\mathsf{\mu }}_{\mathrm{s}}^{\prime }$ in the literature [38].

3.2. Modeling the Fluence Rate Distribution

Using the MCML source code, the fluence rate distribution inside each separate tissue sample was simulated at each studied wavelength, as demonstrated in Figure 14. The fluence rate distributions through a block containing the three tissue layers (skin, skull, and brain) were obtained for all four wavelengths.

As shown in Figure 14, the fluence rate distribution inside the tissue increased at longer wavelengths for each tissue type. The maximum fluence rate was obtained at 980 nm because of the low scattering and better transmittance. The simulation results showed a clear increase in the fluence distribution inside the three-layered model at 980 nm, followed by 785 nm, 660 nm, and 532 nm.

The fluence rate at the outer surface of the sample was based on the finite element method (FEM) solution of the diffusion equation via COMSOL multiphysics software. The obtained results are presented in Figure 15. The fluence rate at the tissue surface became less diffused for longer wavelengths as a result of the low scattering and better penetration.

The obtained optical parameter values of the rabbit head tissue samples were introduced, through Equation (23), to the diffusion of light propagation in biological tissue to investigate the change in the fluence rate distribution at the tissue surface as a result of the change in optical parameters. Figure 14 and Figure 15 show the difference in the optical fluence rate images in a model containing three layers (skin, skull, and brain) of rabbit head tissue samples at 532, 660, 785, and 980 nm wavelengths.

The difference in fluence rate values at each wavelength resulting from the change in optical parameter values introduced to the diffusion equation were used to obtain the fluence rate distribution at the sample’s surface. Varying the sample condition affected the water and blood contents, and hence affected the values of the optical parameters. Moreover, it can be observed from Figure 15 that the fluence rate distribution at the surface of the samples changed with the tissue type due to the change in their optical properties.

4. Conclusions

Recent technological advances in photonics have prompted extensive advancements in the improvement of imaginative strategies and systems for clinical and functional optical imaging, phototherapy, and laser surgery. Intensive studies and the development of biomedical optical methods and techniques have stimulated great interest in the quantitative assessment of the optical properties of tissues. The optical properties of the head tissue depend on the laser wavelength, and the delivery of radiation to the target in the depth of the tissue depends on the fluence rate on the surface and the thickness of the head tissue layers. The Kubelka–Munk model and Monte-Carlo modeling were used in this work to evaluate the optical properties of the skin, skull, and brain tissues of the rabbit head from measurements of diffuse reflectance and the transmittance of tissue samples using a single integrating sphere. The obtained data for the optical properties were close to the previously published data and significantly supplement them.

As the attenuation of radiation at the considered wavelengths is quite strong, it is necessary to use methods of optical clearing of the upper layers of the head tissue. The experimental and theoretical methods developed in this study should be used for analyzing the distribution of the laser fluence rate in the case of optical clearing, which is used to reduce the radiation density in the upper layers of the head in order to avoid strong heating when exposed to intense laser beams.

In conclusion, a method based on measuring the optical properties of biological tissues using integrative sphere and optical fluorescence distributions on tissue surfaces has been applied as a new diagnostic tool. In this method, a combination of the Kubelka–Munk model and Monte-Carlo model was used to calculate the optical parameters of different samples of biological tissues from the measured values of diffuse reflectance and transmittance using the diffusion equation. The measurements were carried out using an experimental setup based on the refractometer and integration sphere in rabbit head tissues at selected laser wavelengths. The currently proposed method is suitable for in vitro measurements; however, it can be upgraded via the use of optical fiber probes to be suitable for in vivo study. The results obtained here are promising because the fluency rate distribution images have discriminatory features between different tissue types. Therefore, this method can be used in the diagnosis and differentiation of biological tissues.