1.

Introduction

There has been tremendous interest in recent years in implementing new techniques for optical stimulation of neurons. Optical stimulation is often considered the superior technology, compared with electric stimulation, since the inherent parallelism of optics allows researchers to excite or inhibit the activity of cells in large-scale networks of the brain. Among all developed techniques, optogenetics has become the most popular and widely used method for the study of the dynamics of the brain and circuitry of neurological and psychiatric disorders. In optogenetic brain stimulation, specific cell-types of interest, on the surface or inside the brain tissue, are genetically targeted to express certain light-gated ion channels or ion pumps.^1–8 Once these proteins are produced in a cell, the activity of the cell can be increased or suppressed simply by exposing the cell to specific wavelengths, mostly in the visible range of the spectrum. Obviously, delivering an adequate amount of light to stimulate the target area within the brain is an important step in the design of all optogenetic experiments. An insufficient amount of light cannot generate effective stimulation, whereas excessive light intensity can potentially cause damage to the tissue or penetrate deeper and stimulate other brain objects. In optogenetic stimulation, the amplitude of the induced photocurrent in a neuron depends on multiple factors, including biological variables such as the kinetics of the expressed opsin⁹ and optical properties of the tissue, as well as system design variables such as the intensity and wavelength of light. Two separate approaches have been explored for light delivery in optogenetics. For cortical stimulation and superficial areas, spatial light modulators, such as MEMS devices¹⁰ or liquid crystals,^11,12 are employed to generate complex and even three-dimensional (3-D) patterns of light distributions inside tissue. However, to target deep brain objects, optical fibers are implanted to guide and deliver laser pulses to the region of interest.^3,4,13–15

Predicting the light intensity delivered by an optical fiber to a target area in the brain is necessary for proper in vivo optogenetic stimulation and has been briefly studied before. A simplified Kubelka-Munk model,¹⁶ which neglects the absorption coefficient, was proposed by Aravanis et al.¹⁵ and later by Stark et al.¹⁷ To develop this model, they conducted a sequence of experiments and measured the intensity of the transmitted light through mouse brain slices of different thicknesses, then used this experimental data to formulate a set of equations which estimate the axial variations of light intensity inside the brain tissue. More recently, Al-Juboori et al.¹⁸ used the fiber punch-through method to measure the effective attenuation coefficient of different areas of the mouse brain and employed the Beer-Lambert law to measure the effective attenuation coefficient in tissue slices of 600 μm thickness. Then, based on the extracted optical properties, the diffusion equation in one-dimension was adapted to determine the required input light intensity needed to deliver the desired optical power to target a deep brain object. Although this approach offers a systematic pathway to measure the tissue effective attenuation coefficient, it is still unable to clarify other optical properties that are essential for precise estimation of light distribution in any turbid medium. Moreover, these models do not incorporate the impact of the tissue heterogeneity on light distribution in the medium.

Obviously, an initial requirement for the design of all biophotonic experiments (including optogenetic stimulation) or related theraputic protocols is obtaining reasonable estimates of the optical properties of tissue.¹⁹ In most published literature, tissue optical properties are parameterized by three variables: the absorption coefficient, ${\mu }_{\mathrm{a}}$, the scattering coefficient, ${\mu }_{\mathrm{s}}$, and the anisotropy factor, $g$.¹⁶ The mathematical models that are developed to emulate light transport in turbid media use these variables to estimate the distribution and penetration depth of light within the medium.¹⁹ These phenomenological properties of biological tissue can be measured in vivo or in vitro. For in vivo measurements, multiple time-domain or frequency-domain methods have been developed with instrumentation such as confocal microscopy and optical coherence tomography used to collect the raw data.²⁰ Then, the values of tissue optical properties are extracted from the raw data by solving a set of inverse problems.^21–23 Data acquisition in vitro is performed mostly by using double-integrating-sphere setups,²⁴ which are used to measure the total diffuse reflection and transmission and the intensity of the transmitted ballistic beam. By making these three measurements and using mathematical algorithms, e.g., inverse adding doubling (IAD) method,²⁵ the optical properties of the tissue under test are estimated.

In this paper, we present an approach to measure the three variables that define optical properties for the rat brain tissue. To extract the values of these parameters, 500-μm brain slices were prepared and scanned by two customized optical setups and the diffuse transmittance, diffuse reflectance, and ballistic transmittance for each slice were measured with high spatial sampling rate of 3000 points per square centimeter. These measurements were made for three different cuts (transverse, sagittal, and coronal) and three different wavelengths, 405, 532, and 635 nm, to cover the visible range of the spectrum. Then, the IAD method²⁵ was employed to reconstruct the optical properties of the tissue using the collected experimental data. By repeating this procedure for all slices, a 3-D database of the tissue optical properties was developed. Next, this database was connected to our 3-D Monte Carlo toolbox^26,27 to simulate light-tissue interactions and find a reasonable estimation of the distribution of light inside the tissue. The developed software incorporates the effect of tissue heterogeneities, physical parameters of the source, including optical fiber core diameter, numerical aperture (NA), or beam type, such as uniform or Gaussian,²⁸ and the wavelength of the source to offer a more accurate model for light distribution in the brain. The process of extracting tissue optical properties is shown in Fig. 1.

Fig. 1

The algorithm adapted for extracting tissue optical properties. The process starts with sample preparation in which brain slices of 500 μm thickness are produced and scanned by the customized optical setups shown in Fig. 2. Next, the inverse adding doubling (IAD) reconstruction algorithm is applied to the collected data to extract the value of parameter, which determine the optical properties of the tissue.

2.

Material and Methods

2.1.

Sample Preparation

In these experiments, 27 female Sprague-Dawley rats weighing from 250 to 300 g were used. All animal procedures were approved by the University of Wisconsin-Milwaukee Institutional Animal Care and Use Committee and were conducted in accordance with National Institute of Health standards on the humane treatment of laboratory animals. Animals were anesthetized by isoflurane (1.5 to 2% oxygen) before decapitation by a lab guillotine. The intact brain tissue was extracted from each animal and fixed on an agarose block for slicing. Next, transverse, sagittal, and dorsal brain slices of 500 μm thicknesses were cut using a vibratome (Lancer vibratome series 1000, St. Louis, Missouri) while the brain was submerged in phosphate buffered saline solution (Sigma-Aldrich, St. Louis, Missouri) and cooled by dry ice. Then, fresh slices were loaded into the optical setup displayed in Fig. 2 and scanned one slice at a time. In each round, one slice was cut and scanned during the period the vibratome was cutting the next slice. To minimize the timing of the experiment during which the brain is kept in the freezing solution to preserve the tissue and its optical properties, it is important to synchronize the slicing and the scanning procedures. In our experience, following this protocol, the optical properties did not have significant change during the scanning of each slice and it took $\sim 2\mathrm{h}$ to complete the whole brain scan.

Fig. 2

Schematic of the experimental setups used to measure diffuse reflected and transmitted light (setup 1) and transmitted ballistic light (setup 2).

2.2.

Measurement Procedures

2.2.1.

Measuring diffuse reflectance and transmittance

Based on the IAD model, three different measurements (diffuse reflectance, diffuse transmittance, and ballistic transmittance) are required to extract the complete set of tissue optical properties, including ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}$, and $g$.²⁵ Diffuse reflectance and diffuse transmittance were measured simultaneously using optical setup 1 and the ballistic transmittance measurements were performed by optical setup 2, which are both illustrated in Fig. 2. Diffuse reflectance and diffuse transmittance measurements were acquired using the double-integrating-sphere technique. For this experiment, after preparing a brain slice, the slice was placed between two integrating spheres (one two-port, 6-in. CVI Melles Griot BPS integrating sphere with 1.5-in. sample port diameter, and one four-port, 6-in. Spectraflect Integrating Sphere, Newport Corporation, Irvine, California, with 1-in. sample port diameter) and was scanned using a focused laser beam (with three different wavelengths: blue laser Thorlabs 405-nm laser diode, green laser Ultralaser 532-nm diode-pumped solid-state, and red laser Thorlabs 635-nm laser diode) with a 200-μm spot size at the focal point on the sample. The diffuse reflected light intensity, $R(r_{\mathrm{s}}^{\text{direct}},r_{\mathrm{s}})$, was measured using a photodetector (Thorlabs, Newton, New Jersey, PDA36A, Si switchable gain detector) that is mounted on the top sphere, and diffuse transmitted light intensity, $T(t_{\mathrm{s}}^{\text{direct}},r_{\mathrm{s}})$, was measured by a second photodetector (Thorlabs, APD110A2/M, avalanche photodetector) that was installed on the bottom sphere. The variables reflectance ($M_R$) and transmittance ($M_T$), defined by the amount of light reflected by and transmitted through the sample normalized to the intensity of the incoming light, respectively, were calculated using the following set of equations:²⁵

Eq. (1)

$$\begin{gathered} M_R\equiv r_{\mathrm{std}}·\frac{R(r_{\mathrm{s}}^{\text{direct}},r_{\mathrm{s}})-R(0,0)}{R(r_{\mathrm{std}},r_{\mathrm{std}})-R(0,0)} \\ {M}_T\equiv \frac{T(t_{\mathrm{s}}^{\text{direct}},r_{\mathrm{s}})-T_{\text{dark}}}{T(0,0)-T_{\text{dark}}}. \end{gathered}$$

In these equations, parameters $R(0,0)$, $T(0,0)$, and $T_{\text{dark}}$ are calibration measurements, and $r_{\mathrm{std}}$ is the wall reflectance coefficient of the integrating sphere. More details of the calibration and measurement process are available in Ref. 25.

By using a galvo scanner (Te-Lighting, China, PT-30K laser scanner galvo system kits) for each brain slice, $M_R$ and $M_T$ variables were measured with a spatial resolution of 3000 points per square centimeter as discussed later. In our measurements, we utilized a CCD camera to obtain the geometry of the tissue and used this information to avoid scanning the empty areas around or within the slice to speed up the scanning process. Therefore, it was necessary to register the galvo scanner on the camera output. The registration was carried out by formulating the relationship between the voltages applied to the scanning mirrors of the galvo system and the coordinates of the laser beam on the tissue, which was recorded by the camera. For this purpose, four different voltages were applied to the scanning mirrors while the camera was monitoring the position of the beam, and the system used this information to map the position of the scanning mirrors on the pixels of the CCD camera. The registration needs to be done only once during the whole experiment. Then, for each slice, the geometrical shape was obtained by capturing a picture of the brain slice by the CCD camera, and this image was used as a mask to separate the tissue from the background. The developed binary mask defines the region of interest (ROI) which includes just the tissue area. This complex ROI also excludes the edges of the tissue where measurements are not accurate and also the footprints of ventricles within each slice. The process of scanning and data acquisition for each slice took $\sim 1\min$., and it was controlled through a custom software developed using LabVIEW (LabVIEW 8.2, National Instruments Corp., Austin, Texas). The procedure for system calibration and data evaluation are discussed in Appendix A.

2.2.2.

Measuring ballistic transmittance

In the ballistic transmittance test, we measured the percentage of the beam that passes through a sample with no scattering or absorption. Traditionally, for this test, a photodetector is used to measure the intensity of the beam that passes directly through the sample in combination with an aperture placed in front of the detector that rejects the scatterer rays. For heterogeneous samples, such as brain slices, this test can be performed either by moving the sample or possibly the aperture and the photodetector, which makes the process quite time consuming. To avoid this problem, we used a telecentric lens system instead of the aperture, which only passes and projects the collimated portion of the transmitted beam onto a CCD rather than the photodetector.²⁹ In these experiments, each brain slice was placed on a sample holder between the galvo mirrors and the camera and once again the slice was scanned with the same spatial resolution of 3000 points per square centimeter (see Fig. 2). Similar to the previous measurements, only the tissue area was scanned. By processing the captured frame, the intensity of the ballistic beam was computed by integrating the values of all the pixels within a circular area that has the same diameter as the reference beam-waist. Then, the transmitted ballistic beam ratio was calculated using the following equation:

Eq. (2)

$$M_U=\frac{I_B-I_N}{I_R-I_N},$$

where $I_N$ is the mean value of the background noise intensity when the laser beam is blocked, $I_B$ is the intensity of the transmitted ballistic beam, $I_R$ is the intensity of the reference beam, and $M_U\in [0,1]$ is the ballistic transmittance. The calibration process, including measuring the intensity of the reference beam ($I_R$) and the background noise ($I_N$), needs to be done only once during the whole experiment. To measure the intensity of the reference beam, a predefined voltage was applied to the laser and an image was captured by the CCD camera without any sample in the path. To prevent saturating the camera, a neutral density filter was positioned in front of the CCD sensor to reduce the laser intensity. As described earlier, the summation of all the pixel values inside the beam-waist was calculated to obtain the reference beam intensity ($I_R$).

In Fig. 3, a sample transverse cut of rat brain and the corresponding interpolated diffuse reflectance, diffuse transmittance, and ballistic transmittance measurements are shown for the green laser (532 nm) measurements. These data are used later to extract the optical properties of the tissue, including the absorption coefficient, scattering coefficient, and anisotropy factor.

Fig. 3

Results of the scanning process with the green laser wavelength of 532 nm. Data are presented in arbitrary units. (a) Image of a sample brain slice. (b) Interpolated reflectance measurements. (c) Interpolated transmittance measurements. (d) Interpolated ballistic transmittance.

2.2.3.

Mapping the experimental data

Since the measurements are made by two different setups, there is a possible displacement between the data collected by these two optical setups. Therefore, before any further process, this displacement should be compensated. We adapted the MATLAB® command “imregister” from the image processing toolbox (MATLAB® 7.14, The MathWorks Inc., Natick, Massachusetts) to register these two datasets. Once the dataset was prepared and displacement was compensated, these data were fed to the IAD code.

3.

Extracting Optical Properties

The IAD method is an iterative algorithm designed to extract the optical properties of a homogeneous sample from the measurements discussed in the previous section. The iterations are initiated with some random values assigned to the optical properties of the medium. Then, in each iteration, the algorithm solves the radiative transport equation ³⁰ to calculate the transmission and reflection of the sample and compares these results with the experimental data. Next, it readjusts the values of the optical properties to minimize the error between the calculations and experimental results and continues this process until convergence.

As mentioned before, in the IAD method, the assumption is that the sample is a slab with uniform optical properties. To satisfy the homogeneity assumption in IAD algorithm, the scattered light in the tissue should be confined to a small region in which the optical properties remain relatively uniform inside the illuminated area. Our Monte Carlo simulation for light propagation in a slab of turbid medium of thickness 500 μm, which has uniform optical properties close to typical values of biological tissue (${\mu }_{\mathrm{a}}=1{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}=31.6{\mathrm{mm}}^{-1}$, and $g=0.9$), reveals that the beam width of a laser, which is 200 μm in diameter before the slab, increases to $\sim 500\mu \text{m}$ after passing through this scattering medium, as shown in Fig. 4. Therefore, in our experiments, we collected data from discrete sample points and assumed that the optical properties are uniform within the illuminated area for each data point. Then, we used the interpolation to fill the gaps between discrete data points on the tissue.

Fig. 4

Lateral distribution of light transmitted through a slab with thickness of 500 μm and typical tissue optical properties when illuminated with a laser beam of 200 μm diameter.

Once the experimental data are prepared and displacement was compensated, these data were loaded to the IAD code. Since IAD is a time-consuming process, not all of the interpolated data points could be used for reconstruction. In our code, we used only a subset of the sampled data points to extract the optical properties via IAD. In the next step, the density of the extracted optical properties was increased by a two-dimensional bicubic interpolation algorithm.³¹ For the brain slice shown in Fig. 3(a), reconstructed values of absorption coefficient, reduced scattering coefficient, scattering coefficient, and anisotropy factor are shown in Fig. 5. Another example of such a reconstruction that illustrates the heterogeneity of the brain tissue regarding the optical properties is displayed in Fig. 6.

Fig. 5

Extracted optical properties of the rat brain slice produced by IAD algorithm, which is applied to the raw data shown in Fig. 3. (a) Reduced scattering coefficient. (b) Absorption coefficient. (c) Scattering coefficient. (d) Anisotropy factor.

Fig. 6

Extracted reduced scattering coefficient at 532 nm for a rat brain tissue: (a) transverse cut and (b) sagittal cut.

The described calculation of ballistic transmittance is admittedly approximate, and a more sophisticated deconvolution of the point spread function in angle space of transmitted light versus variable sample thickness would more appropriately define the ${\mu }_{\mathrm{s}}$ and $g$ values. However, when $g\gg 0.5$, the spatial distribution of light is dominated by the lumped parameter ${\mu }_{\mathrm{s}}^{-1}={\mu }_{\mathrm{s}}(1-g)$, which is reliably specified by the integrating sphere measurements, even in the region of nondiffuse light that occurs near the probe tip. See Appendix B for an illustration of this similarity principle. So the properties of ${\mu }_{\mathrm{s}}$ and $g$ deduced here, even if slightly in error, will still yield reliable predictions of light distribution in the tissue, even near the light source.

4.

Results and Discussion

In order to investigate the impact of tissue heterogeneities and different light source configurations on the distribution of light inside tissue, a 3-D Monte Carlo toolbox was prepared and linked to the developed database of the rat brain tissue optical properties. Results of these simulations are presented in the following sections.

4.1.

Impact of Brain Tissue Heterogeneity on Light Distribution

For any inhomogeneous medium, such as the brain tissue, we logically expect to observe a different or possibly more complex light distribution pattern compared to the light distribution pattern in a uniform medium of the same volume. To investigate the effect of tissue heterogeneities on light distribution in the brain, we simulated the photon transport process inside the tissue using our 3-D Monte Carlo code and compared our results with the dual homogenous medium. In both simulations, all photons are launched from a uniform source with diameter of 150 μm, which is placed at the specified area marked in Fig. 7(a). For the dual homogenous model, the mean value of the optical properties of the whole tissue was assigned as the optical properties of a virtual uniform medium. On the other hand, for the heterogeneous test, we used the extracted database of the optical properties. Optical properties of the tissue in the region are displayed in panels Figs. 7(b) to 7(d). An example of the difference between the light distributions in these two structures is shown in Figs. 7(e) and 7(g), and the corresponding contour maps of the light distributions are displayed in Figs. 7(f) and 7(h). Furthermore, light distribution along an arbitrary axial line is displayed in Fig. 7(i). This result highlights the effect of tissue heterogeneity on light distribution, which, in this experiment, has caused some considerable asymmetry in the distribution pattern in a way that the intensity of light at a distance of $\sim 2.50\mathrm{mm}$ on one side is more than three times stronger than the light intensity at the same distance from the center of the source on the other side. The curve in Fig. 7(j) shows the difference between the two distributions that is caused by the heterogeneity in absorption, scattering, and anisotropy factor. Based on our observation, in short distances in the brain tissue, the difference between the light distribution in the homogeneous approximation and the actual heterogenous model is minimal. However, as the distance from the source increases, this difference grows. Therefore, the heterogeneity effect in the design of optogenetic experiments is more essential once highly sensitive opsins, such as step-function opsin,⁶ are used. The reason is that in these cases, even low-intensity light that reaches deeper regions of brain tissue can stimulate the cells.

Fig. 7

(a) An optical fiber of 150 μm diameter is placed on the marked area inside the tissue to launch a uniform beam. The three-dimensional Monte Carlo simulations are run by launching 10 million photons. Optical properties of the tissue in the region are displayed in (b) to (d) for a blue laser at 405 nm. (e) Two-dimensional representation of the light distribution in the $XZ$ plane ($F_{XZ}$) for the homogeneous brain tissue. (f) Contour map of the light distribution in the homogeneous brain tissue. (g) Two-dimensional representation for the distribution of light in the $XZ$ plane ($F_{XZ}$) for the inhomogeneous brain tissue. (h) Contour maps of the light distribution in the inhomogeneous brain tissue. (j) Lateral fluence rate of the light along the $X$ axis for the inhomogeneous and homogenous brain tissue, which shows considerable difference between the two distributions. Distribution of the light for homogenous tissue along the $Z$ axis shows almost an identical change in both directions far from the fiber position (solid curve) while the light distribution has become asymmetric as a result of the tissue heterogeneity (dashed curve). (i) The difference between distribution of light in the homogeneous and inhomogeneous tissue.

The impact of the tissue heterogeneity on the axial distribution of light is also significant. Figure 8 shows the simulation results in which the axial penetration depth at two different regions, marked by points A and B in the figure, are compared. In this case, at a depth of 2 mm from the fiber tip, light fluence rate for point B is almost five times larger than that for point A.

Fig. 8

(a) Optical fiber with 100 μm diameter is placed on the marked regions in the tissue. (b) Comparing the axial fluence rate along $z$ axis for points A and B. The difference between the attenuation coefficients in these two regions has caused significant difference between the axial penetration depth of light at these two positions for a blue laser at 405 nm.

4.2.

Impact of Physical Parameters of the Source on Light Distribution

Choosing the proper light delivery mechanism is important in optogenetic experiments, particularly when the application demands precise stimulation of some target areas. Optical fibers are widely used for optogenetic stimulation of deep brain objects in vivo.^3,13–15 Physical parameters of the optical fiber, such as core diameter and NA, have a certain impact on the beam profile and its distribution in the medium. Here, Monte Carlo simulation was used to investigate the light propagation along the tip of the optical fiber for different fiber parameters and different wavelengths. In Fig. 9(b), the axial fluence rate (light fluence rate along the $z$ direction) for a 405-nm wavelength is shown for different fiber diameters: 200, 150, 100, and 50 μm. The simulation results show that for the same optical power, the axial fluence rate decreases when increasing the fiber diameter. On the other hand, it seems that reducing the diameter of the fiber below 100 μm does not have significant impact on the axial fluence rate. Figure 9(c) demonstrates the axial fluence rate for a 100-μm-diameter optical fiber with different NAs, including 0, 0.1, 0.22, and 0.39. According to these results, light intensity drops faster along the axial direction for the optical fibers that have larger NA. Figure 9(d) represents the impact of the wavelength on light distribution. Here, for shorter wavelengths, the attenuation rate is higher as expected.¹⁶ For the blue light (dashed line), the attenuation is higher and the intensity drops much faster than the green and red lights. A comparison of the axial attenuation between the Gaussian and uniform light source is presented in Fig. 9(e). These data show that the Gaussian beam has a higher fluence rate.

Fig. 9

(a) An optical fiber is placed on the marked area and the Monte Carlo simulation software is used to investigate the effect of source parameters on light distribution inside the brain tissue for a blue laser at 405 nm. (b) Effect of the fiber diameter. (c) Effect of fiber numerical aperture. (d) Spectral response. (e) Effect of the beam profile on the axial distribution of light.

5.

Conclusion

In this paper, we present an approach to extract the major optical properties (absorption and scattering coefficients and anisotropy factor) of the rat brain tissue with spatial resolution of $\sim 3000$ sample points per square centimeter for three different cuts (transverse, sagittal, and coronal) and three different wavelengths: 405, 532, and 635 nm. A large database of the brain optical properties is produced by slicing the tissue and measuring these parameters in each slice followed by mapping the measurements to the atlas of rat brain.³² The 3-D simulation software is able to incorporate the effect of tissue heterogeneities, physical parameters of the source (e.g., optical fiber core diameter, NA, or beam type, such as uniform or Gaussian), and offers a more accurate model for prediction of light distribution in the brain. When designing an experiment to optogenetically stimulate an area inside the brain, some decisions should be made regarding the physical parameters of the fiber and the intensity of light required to generate effective stimulation in the tissue. The integration of the rat brain atlas and the 3-D Monte Carlo software can provide a useful toolbox and the simulations can help to choose these parameters more accurately, which certainly improves the quality of the results and the rate of success, particularly for in vivo optogenetic experiment. One aspect in tissue light interaction that is not discussed in the current manuscript is the effect of blood on light distribution and penetration. Previous experiments have shown that mammalian blood has relatively high absorption and scattering coefficients,³³ and incorporating the optical properties of the blood in the models potentially improves the precision of the predictions. One possible approach to include the effect of blood is to use advanced optical imaging methods, such as optical coherence tomography, to scan the tissue and generate angiograms, which provide insight to the details of blood distribution in the tissue under test. Next, this angiogram can be integrated into the developed brain atlas of optical properties, which is used for Monte Carlo simulations.

Although brain tissue and optogenetic neurostimulation was the main target in this study, the proposed approach for extracting optical properties and predicting the light distribution can be adapted as a general method for other tissue samples and different applications, such as laser surgery, laser therapy, optical tomography, etc.