OptogenSIM: a 3D Monte Carlo simulation platform for light delivery design in optogenetics

: Optimizing light delivery for optogenetics is critical in order to accurately stimulate the neurons of interest while reducing nonspecific effects such as tissue heating or photodamage. Light distribution is typically predicted using the assumption of tissue homogeneity, which oversimplifies light transport in heterogeneous brain. Here, we present an open-source 3D simulation platform, OptogenSIM, which eliminates this assumption. This platform integrates a voxel-based 3D Monte Carlo model, generic optical property models of brain tissues, and a well-defined 3D mouse brain tissue atlas. The application of this platform in brain data models demonstrates that brain heterogeneity has moderate to significant impact depending on application conditions. Estimated light density contours can show the region of any specified power density in the 3D brain space and thus can help optimize the light delivery settings, such as the optical fiber position, fiber diameter, fiber numerical aperture, light wavelength and power. OptogenSIM is freely available and can be easily adapted to incorporate additional brain atlases.


Introduction
New advances in optics and molecular genetics have spurred the rapid advances in the field of optogenetics to modulate neural activities in specific cell-types of interest using light [1][2][3][4][5]. Optogenetics has become a powerful tool for neuronal control and stimulation, and has recently demonstrated great promise to investigate the neuro-circuitry of a number of neurological diseases, such as Parkinson's diseases, autism, Schizophrenia, anxiety, depression, etc [6][7][8]. As the applications for optogenetics expand, strategies that can ensure adequate light delivery precision for modulating selected neurons while minimizing undesired side effects are critically needed [4].
In order to best design the optimal light delivery strategy, it is important to evaluate or estimate the effect of each of the parameters of the light delivery and of the brain tissue to be interrogated. Light distribution in brain can be affected by optical properties such as scattering coefficient, absorption coefficient, and scattering anisotropy factor, light source, and fiber geometry and positioning [2,4]. Monte Carlo methods [9,10] and analytical models including those related to diffusion theory or radiation transport equation [2,11] have already found their applications in predicting the light power density in optogenetic applications. However, these methods do not take into account the non-uniform 3D distribution of optical properties of the brain or the brain tissue heterogeneity as a whole. Examples of assumptions made in previous studies include: the whole brain was considered as all gray matter (GM) [2,9], different brain structures were treated separately and the light distribution was estimated within each of the structures separately using simplified one dimensional diffusion theory [11], only two dimensional Monte Carlo was applied [10]. In addition, there is no intuitive way to illustrate the light power distribution in three dimensional views of the brain.
Hence, to resolve the challenges of brain tissue heterogeneity in terms of spatially varying optical properties [12][13][14][15][16][17], complex light delivery strategies [4,18,19], and the complicated spatial anatomy of the mammalian brain [20,21], a complete 3D Monte Carlo modeling approach is needed. While a 3D Monte Carlo modeling platform will clearly improve the capabilities of optogenetics by offering powerful simulations of tissue light interactions, two practical issues must first be explored and addressed: 1) reliable optical properties of brain tissue may not be available for all wavelengths or tissue structures of interest; and 2) there has been no publicly available 3D rodent brain atlas designed for the purpose of conducting Monte Carlo simulation. For the first issue, optical properties of brain tissue vary across different regions and are wavelength dependent, and the reported measured optical properties are subject to variations in the measurement conditions [11][12][13][14]17,[22][23][24][25]. Any effective modeling approach has to take this variability into account. For the second issue, the rodent brain has been widely used in optogenetics and other researches in neuroscience, and both in vivo and ex vivo 3D atlases are available. However, the atlases are constructed based on neuronal anatomy, not optical properties, and are not able to be directly applied to the Monte Carlo simulation. For example, the brain was segmented into 671 unique brain structures in Allen mouse brain atlas [26], 62 structures in an average MRI mouse brain atlas [27], 20 regions in 3D MRM(Magnetic Resonance Microscopy) mouse brain atlases [28,29], and 29 brain regions in a 3D DTI(diffusion tensor imaging) adult rat brain atlas [30]. The structures or regions in the atlases are unable to be classified based on the optical properties since optical properties may have significant change in different regions of the same neuroanatomy structure. Generally, a neuronal structure is a mixture of gray matter (GM) and white matter (WM) that are significantly different in optical properties [14,[23][24][25] and thus these matter types need to be classified first in order to create an effective Monte Carlo applicable database.
To resolve the challenges described above, an open-source 3D simulation platform, OptogenSIM was developed that integrates a 3D Monte Carlo model, generic models of brain tissue optical properties, and a well-defined brain tissue atlas of structure and tissue type. Specifically, four major goals were accomplished in our platform. First, we created a 3D Monte Carlo method that accommodates brain tissue heterogeneity and different light delivery approaches in optogenetic applications. 3D Monte Carlo code has been developed by several groups including voxel-based tMCimg [31] and mesh-based 3D Monte Carlo simulation tools [32,33]. The "gold standard" Monte Carlo program MCML [34] is for simulation of 3D photon propagation in planar multi-layered tissues. Our current work builds on our previous 3D Monte Carlo work named "mcxyz", which is freely available (http://omlc.org/software/mc/mcxyz/, which presents a description of the method) and uses the same light transport model as in MCML [34] with only minor exceptions. This Monte Carlo system reads an input text file in a standard format, which can be created by a variety of editors routinely used to create input files such as MATLAB (The Mathworks, Inc., Natick, MA). The outputs are standard binary files of the fluence rate within the heterogeneous tissue. We previously used mcxyz for light power estimation [35,36]. In OptogenSIM, the fiber beam modeling has been enriched and incorporated, as shown in the appendix. Second, generic brain tissue models were created in the wavelength range of 400nm-1000nm for brain gray matter and brain white matter based on reported values in literature as well as our own experimental measurements [35,36]. These previous measurements yielded the optical properties of GM and WM for a generic brain model. Third, a mouse brain tissue type atlas was defined for the purpose of conducting 3D simulation by segmenting a public freely available mouse brain anatomy atlas [28] into GM, WM and cerebrospinal fluid (CSF) regions using a statistical parametric mapping tool [37,38]. The light delivered to arbitrary brain structures or brain atlas locations can be simulated, which takes into account brain tissue heterogeneity while avoiding a dramatic increase in the tissue model complexity. Fourth, a graphical user interface (GUI) based on MATLAB (The Mathworks, Inc., Natick, MA) was developed for this simulation platform, in which the user can intuitively adjust the simulation settings and check the light power density at any desired 3D position.
With tested applications of this 3D Monte Carlo platform, a significant to moderate impact of brain heterogeneity is demonstrated and the iso-fluence contour of any specified light density can be displayed in any region of the 3D brain. The 3D Monte Carlo platform presented here provides a useful tool for the design of light delivery strategies in optogenetics applications. This open source platform is named OptogenSIM and is freely available from www.loci.wisc.edu/software/optogensim.

Construction of 3D Monte Carlo simulation platform
The simulation platform has a MATLAB-based GUI and integrates a 3D Monte Carlo model, a generic brain tissue optical properties library, and a well-defined 3D mouse brain tissue type atlas.

3D Monte Carlo model
The 3D Monte Carlo method used in this study is based on previous work [34,39] and is a Cartesian-grid voxel-based method, which allows assignment of optical properties, including , to each individual voxel. The photon transport differs from the MCML [34] as follows: (1) Each voxel can have a different optical property, so a heterogeneous tissue can be modeled; (2) each photon step size (s) is calculated as s = -ln(random number)/μ s , and after each step photon weight (W) is deposited as W(1-exp(-μ a s)). When a photon step will cross multiple voxels, the photon takes partial steps across each encountered voxel and drops a portion of its photon weight in each voxel; (3) in addition to the standard collimated orthogonal point light source used in MCML, other light delivery mechanisms commonly used in optogenetics can also be implemented, such as optical fiber based stimulation in which physical parameters of the fiber such as diameter and numerical aperture (NA) can be defined as well; (4) a three-dimensional Cartesian coordinate system of voxels is used to record fluence rate rather than the cylindrical coordinates of MCML. A fluence rate is recorded for each individual voxel located in the 3D tissue atlas. The 3D code was validated against previous code in MCML using a 2-layer skin tissue structure that worked with both code bases (data is shown in the OptogenSIM website mentioned above). A more detailed description about the implementation of this model can be seen in the appendix.

3D brain atlas
The representative ex vivo adult male C57BL/6J mouse brain atlas used here was selected from a comprehensive digital atlas database of magnetic resonance images [28]. This 3D atlas was manually delineated into 20 brain anatomical structures and was used as a reference to segment other brain images. To integrate this atlas into our simulation platform, this atlas was segmented into GM, WM, and CSF regions by using a statistical parametric mapping based tool named SPMMOUSE [37,38]. An affine registration was first employed to reorient both the original volume atlas and the annotation atlas in order to align them with the prior template. Then the reoriented atlas was segmented and smoothed in its native space. The defined individual single mouse brain tissue atlas has 256 × 512 × 256 x × y × z voxels, and each voxel has a dimension of 47µm × 47µm × 47µm.

Generic brain tissue models
The generic tissue model is based on the tissue constitution and can provide a good estimate of optical properties for various tissue types at any desired wavelength [17]. The reduced scattering coefficient (μ s ') is modeled as where the scaling factor a equals the value of ' s μ at 500 nm wavelength. The b factor is the scattering power, typically ~0.
where S is the hemoglobin (HGb) oxygen saturation of the mixed arterio-venous vasculature, B is the average blood volume fraction, and W is water content. Each tissue type is specified by characteristic values of B, S, W, a and b. A choice of wavelength (λ) is made for a particular Monte Carlo simulation. The absorption spectrum for each absorber can be found at the Oregon Medical Laser Center website (http://omlc.org/). For the generic reduced scatting properties of GM and WM, the optical parameters were inferred from previous measurements [35] by a least-square fitting approach with g set to be 0.9. Two typical sets of the average reduced scattering coefficient μ s ' were extracted from a coronal slice and dorsal slice, respectively, at the wavelengths of 405nm, 532nm and 635nm by registering the optical properties map into the bright field slice image, and then manually segmenting the gray matter and white matter based on the image intensity as measured with MATLAB and Fiji [40].
To properly mimic the blood of in vivo GM and WM tissues, the mean blood content (B) and oxygen saturation (S), from Table 3 in the review [17] were used for the absorption model. The water content (W) was assumed to be 0.65.
The values for the average generic models are summarized in Table 1. For the CSF, its absorber was simplified as only water. Its reduced scattering coefficient at 500nm was set to 0.24 mm −1 with g set to 0.9. The region outside the brain is non-scattering ~non-absorbing medium with wavelength-independent optical properties of µ a = 0.01mm −1 , µ s = 1mm −1 , anisotropy of scattering, g = 1.0. The g value of 1.0 avoids any scatter. The very low µ a value allows negligible photon weight deposition, which is sufficient to later allow calculation of the fluence rate escaping the tissue, φ = (absorbed density)/μ a . The model in its current version does not include skull or skin, which will be incorporated in a future version. This version is appropriate for modeling the φ in the tissue around an implanted optical fiber. Figure 1 shows the optical properties of gray matter and white matter from the generic models and the ex vivo measurements from both our previous study [35] and others [14,41].  Table 1). Symbols are experimental data. (A) Reduced scattering coefficient shows a significant difference between gray matter(GM) and white matter(WM); (B) Absorption coefficient of model assumes an in vivo blood content (cyan curve), which differs from literature values based on ex vivo samples (symbols).

Investigate light delivery strategy
In this paper, the light delivery approaches were evaluated by creating the light density distribution in terms of fluence rate map, contour, and attenuation relative to the position of the optical fiber. First, the impact of tissue heterogeneity was simulated by launching photons at locations where photons will most likely transport across more than one tissue types and the results were compared with a homogeneous tissue model (gray matter). Second, the light distribution was compared to that yielded by launching photons at an adjacent position by stimulating areas close to or within the neocortex. Then, an application of the simulation platform was demonstrated to show its potential uses in practice by using different stimulating parameters.

Investigate impact of brain heterogeneity
Two demonstrations are presented to show that the brain heterogeneity can affect the light distribution. The first demonstration in Fig. 2 shows that the heterogeneous brain model has a different light distribution pattern than the homogeneous brain model if the fiber is implanted near the boundary of gray matter and white matter. The second demonstration in Fig. 3 shows the fluence change within a heterogeneous brain by implanting the light delivery fiber at two adjacent positions.  Table 1.  Figure 4 shows the application of OptogenSIM to study the effect of different stimulation parameters. The parameters can be set from the user graphic interface, including wavelength, fiber position, fiber diameter, fiber NA, beam profile type, etc. Four built-in iso-fluence contours or a single contour with a specified value can be chosen to be overlaid on the original brain atlas. The fluence rate at any 3D position of the brain can be checked by moving the mouse icon. The output window can show three views (zy, zx, and yx planes) with iso-fluence contours, and these views can be either the original atlas, or the segmented tissue type atlas (GM, WM, CSF). The display of iso-φ contours allows the user to see the region of brain that receives a fluence rate equal to or greater than that particular φ. The user can adjust the power, which rescales the Monte Carlo simulation, and redisplays the iso-φ contours until only the desired target region of brain is contained within the iso-φ contour. In this way, the user can determine the proper amount of power to deliver through the optical fiber. Other output files include a .csv file (comma-separated values) as well as a Matlab .mat file that the user can use for further post-processing of the fluence distribution.

Discussion
The significant difference in the optical scattering properties of gray matter and white matter is shown in Fig. 1. This key difference was the initial motivation for developing optogenSIM to estimate light distribution in a heterogeneous brain model in order to assist the optimal light delivery design in optogenetics applications. The generic models presented here can be used to estimate the optical properties of brain gray matter and brain white matter at any wavelength between 400 and 1000 nm. The models use scattering properties consistent with the literature and absorption properties consistent with estimates of in vivo tissue blood contents. The reported optical properties of brain tissue have a large variation, which can be due to measurement artifacts, variation in tissue constitution, variable tissue preparation, experimental technique, and inaccuracy of calculation models [12,17]. We assume the same B and S values for gray matter and white matter, as it is difficult to infer the in vivo blood content from ex vivo experiments. Certainly, there is a need for more experimental work on in vivo measurements of brain tissue optical properties, and the results can be readily incorporated in the OptogenSIM platform.
The simulations in Fig. 2 show an up to ~2 fold difference in the fluence rate of the heterogeneous model compared to the homogeneous model. The maximum difference occurs distant from the fiber tip where the fluence rate φ has dropped below 5% of the maximum φ near the tip. If the light source power is set such that the target region within a specific iso-φ is close to the fiber tip (i.e., φ/φmax > 5%), the effect of heterogeneity may become less important. However, Fig. 3 indicates that a slight movement of the fiber position can cause a significant change in the axial profile of φ. The absolute value of φ near the fiber tip can be significantly affected by the local heterogeneity. For example, the maximum fluence is 4246/cm 2 for P1, while only 3642/cm 2 for P2 and fluence rate drops faster versus distance from the fiber tip when the fiber is implanted at P1 than at P2. The maximum fluence rate difference is due to the average mean free path which leads to different energy loss before the fluence reaches its maximum value [42]. The position of the implanted optical fiber and the light wavelength used for the stimulation could significantly affect how much light reaches the target. In addition, the fluence rate attenuates very rapidly, for example by 95% within 500-800µm from the fiber tip. The fluence change due to the heterogeneity is dependent on the distance from the photon launching position. Combining the fact that the difference between white matter and gray matter is also wavelength dependent, the effect of tissue heterogeneity is also wavelength dependent.
Other atlases can be incorporated into OptogenSIM. We have tested the simulation on an average mouse brain atlas [38] and an average rat brain atlas [43] but due to copyright issues for direct distribution of the atlases with our tool, these atlases have not been included in the current OptogenSIM platform. For other 3D volumes or 3D atlases not designed for Monte Carlo simulations, segmenting the 3D volume into the tissue types that are included in the tissue library of the platform is needed. There are a number of options to segment the customized 3D volume. Most of the segmentation methods are intensity based statistical classification methods [44][45][46]. One of the most popular atlases, Allen mouse brain "reference atlas" was segmented into GM, WM, and CSF for other purposes [47]. There are also software packages available to help the 3D segmentation. For example the prior atlases discussed earlier can be used as the template for SPM software [37,48] to do the 3D segmentation.
Even though tissue heterogeneity is taken into account here by using the atlases with the segmented GM, WM, and CSF regions, finer structure segmentation would help improve the accuracy of the light density estimation. If the tissue optical properties at regions of interest are known, it would be beneficial to replace the tissue properties of this region by the known properties based on the associated tissue segmentation. For example, generic models of reduced scattering coefficient for cortex (frontal lobe), cortex (temporal lobe), astrocytoma of optic nerve, normal optic nerve, cerebellar white matter, medulloblastoma can be found in [17], and if these brain structures can be located in the 3D volume, the simulation would yield a more accurate prediction. The measured map of the optical properties from our previous study [35] was used to conduct the simulation because of its high spatial resolution. However, so far, the voxel resolution of the measured volume is limited and it is hard to apply this map to other popular 3D atlases without further processing. We are planning to make full use of the measurements by developing a robust 3D reconstruction, registration, and segmentation approach in the future. Figure 4 shows that light source wavelength, optical fiber diameter, numerical aperture (NA), and light source beam profile all play a role in light transport. Light at longer wavelength can penetrate deeper (Fig. 4(B)), larger fiber diameter significantly reduces the maximum fluence ( Fig. 4(C)), smaller NA seems to help light transport in the axial direction but its effect is moderate (Fig. 4(D)), and beam profile type could dramatically affect the maximum fluence ( Fig. 4(E)). However, the effect of each parameter should be evaluated by taking into account other parameters and what really matters is the combined effect of all the parameters.
This platform can potentially be adapted to characterize the light stimulation of high spatial precision methods including multi-site deep brain stimulation devices, such as 3Dmicro LED [9,50,] multi-array silicon probes [50], and 3D glass optrode arrays [51].The anisotropic light diffusion may occur in some anisotropic brain tissue structures [52] and our future versions may address this issue to yield more accurate fluence calculation.

Conclusion
A 3D Monte Carlo simulation platform has been developed that integrates a novel 3D Monte Carlo model, generic brain tissue models, and a 3D brain tissue atlas to assist in designing proper light delivery strategy. The 3D simulation platform can handle tissue heterogeneity and yields the light density distribution throughout the brain. Future work will (1) construct more high-resolution mammalian brain tissue atlases and generic tissue models based on direct measurements and/or publically available neuron anatomy atlases, (2) model light delivery strategies in advanced light stimulation such as multi-site multi-wavelength stimulation, time-sensitive stimulation, (3) quantitatively analyze light delivery under more specific optogenetic stimulation conditions such as mismatched tissue boundary, (4) address the issue of anisotropic tissue structure, and (5) accelerate the computation speed for largescale simulation. This 3D simulation platform is freely available to the public on the LOCI website for non-commercial use.

Appendix: fundamentals of 3D Monte Carlo implementation for Optogenetics application (gomxzyzOGS.c)
gomcxyzOGS.c was adapted from mcxyz.c for optogenetics application, specifically for OptogenSIM software by adding fiber based beam modelling as well as other custom simulation parameters.
mcxyz.c implements a voxel-based 3D Monte Carlo based on gold standard multi-layered MCML code. During the photon transport in the 3D voxels, mcxyz uses an algorithm to determine whether a photon is still within the same voxel after taking a step. If it is still within the same voxel, then it drops some of its weight. If not within the same voxel, the photon has crossed the voxel boundary, and an algorithm finds the face of the voxel where the photon is crossing, and calculates the distance that the photon travels immediately after it enters into the new voxel. The photon will drop partial weight based on the absorption along the above travel distance and then update the remaining step. Three boundary conditions are designed after the photon hits the volume boundary including no boundary (infinite medium assumption), escaping at boundaries, and only escaping at the top surface (i.e. z boundaries).
For the fiber-based light source, both de-focused uniformly distributed beam and two types of de-focused Gaussian beams were implemented to estimate the photo trajectory after it exits from the fiber. The half angle of the fiber is associated with NA using the ( ) The refractive index is assumed to be identical for all voxels, i.e., mismatched voxel boundaries have not yet considered in the current implementation. However, we are planning to implement this in the future by adapting the strategy of MCML for handling air-tissue mismatched boundary. In addition, a future version of the program will allow the angle of orientation of the optical fiber delivering light to be changed.