Optical and thermal simulations for the design of optodes for minimally invasive optogenetics stimulation or photomodulation of deep and large cortical areas in non-human primate brain

The MC simulations were performed using the GATE software (release 8.0) [24]. The Geant4 Application for Emission Tomography (GATE) open-source simulation platform, based on the Geant4 toolkit, has been developed since 2001 [25–27] by the OpenGATE collaboration (www.opengatecollaboration.org) and is currently widely used by the research community involved in molecular imaging, radiotherapy. A module dedicated to optical transport in tissues is available within GATE since release 6.2 [25]. Absorption, Mie scattering, refraction and reflection at boundaries are simulated. The simulation of photon transportation in biological tissues relies on a commonly used and validated implementation [26]. A set of parameters including refractive index (n), anisotropy (g), the absorption length (L_a), and scattering length (L_s) are defined for each type of tissues. In GATE, the mean absorption length (L_a), and mean scattering length (L_a) are defined as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{L}_{{\rm a}}}=\frac{1}{{{\mu}_{{\rm a}}}},\quad {{L}_{{\rm s}}}=\frac{1-g}{{{\mu}_{{{{\rm s}}^{\prime}}}}}.\nonumber \end{align}$$

For the evaluation of the thermal effects we used a hybrid approach based on the bioheat model developed by Pennes [27] and validated previously [28]. Our previous in vivo and post mortem comparisons of the temperature increase of brain tissues undergoing optical stimulation have shown only a slight difference presumably due to blood flow in the living brain (see supplemental figure (stacks.iop.org/JNE/15/065004/mmedia) in [29]). As this difference of the order of a fraction of Celsius degree even for a sustained optical stimulation, we chose to use the simplified Pennes model that accounts for the heat diffusion but not for the convection due to heat transfer in blood vessels.

For NHP cortical stimulation, the targeted cell layers are typically located between 2 and 4 mm beneath the brain surface, which corresponds to the average thickness of NHP cortex. These deep layers should receive sufficient fluence to activate the opsins. The commonly admitted threshold fluence is between 1 and 5 mW mm⁻² depending on the type of opsin used and the transfection efficiency [1]. The fluence in mW mm⁻² was chosen as the parameter of reference as it is directly related to the opsin activation. It incorporates the influence of the light source power but also its coupling to a fiber, and the dimension of the fiber section. Table 1 summarizes the geometrical and optical parameters and their range considered in the MC simulations.

We first assessed the influence of the numerical aperture (NA values of 0.1, 0.37, 0.58, 0.8 and 0.9) for fibers with a fixed diameter of 500 µm. The NA describes the angular emission of fiber. It is the product of the refractive index of the medium where light exits the fiber with the sinus of the half maximum angle of light emission. We then evaluated the impact on the fluence distribution of the fiber diameter for a fixed NA of 0.58. Fibers with diameters of 0.1, 0.2, 0.5, 0.75 and 1 mm have been successively simulated. Figure 1(a) shows the set-up for the single fiber simulations carried out to study the influence of the fiber NA and diameter. The optical fiber axis was orthogonal to the tissue and the brain was tangent to the cortical tissue. In a third set of simulations, we studied the use of several fibers to increase the volume of activated tissue. We studied the effect of the distance between fibers on the fluence distribution homogeneity and the superimposition of fluence from each fiber. We have simulated distance between the fibers of 2, 3, 4 and 5 mm. The distance between the fibers is defined as the distance between their optical axes. For these simulations, fibers with a diameter of 1 mm and an NA of 0.58 were chosen according to the results from the previous simulations. In order to maintain the dimensions of the optode mechanically feasible and to limit its bulkiness we did not carry out simulation for distances larger than 5 mm between the fibers optical axes (figure 1(b)).

Zoom In Zoom Out Reset image size

Since brain tissue strongly absorbs light in the blue/green region of the spectrum, optical stimulation of large and deep volume of brain tissues require red light that is less absorbed and scattered that other visible wavelengths. Regarding optogenetics in NHP, opsins sensitive to red light have been developed for deep tissue optical stimulation. The use of one of these red-shifted opsin, called Jaws, was recently validated in NHP cortex [30]. Accordingly, we simulated photons with a 674 nm wavelength. The optical properties of NHP cortex have been measured previously at this wavelength and used as reference optical properties in the context of optogenetics in NHP [7]. The optical properties of tissues are the aoptical absorption coefficient µ_a, the optical scattering coefficient µ_s, and the optical anisotropy coefficient g. The reduced optical scattering coefficient is a practical coefficient often used in optical transport simulation as it allows to pool the optical scattering coefficient and the optical anisotropy coefficient in a single parameter. It is defined as $\mu _{{\rm s}}^{\prime}$  =  µ_s (1  −  g). The optical [31] and thermal [32] properties for NHP cortex used in the simulations are shown in table 2. We have chosen to use the optical parameters measured on cortical human tissues [31] as they have been obtained in vivo on non-processed tissues. Other data have been obtained for human brain properties [33, 34] but they have been obtained in vitro on processed (cryogenized or saline immersed tissues). Saline immersion leads to partial blood removal and thus decrease the absorption while cryogenization has an strong impact especially on the reduced scattering coefficient [35, 36]. The optical properties use in the present studies have been measured at 674 nm are therefore relevant for the wavelengths used to activate red shifted opsins.

Table 2. Cortical tissues optical and thermal properties used in GATE simplified simulations.

µs	g	$\mu _{{\rm s}}^{\prime}$  =  µs(1  −  g)	µa	n	k	C	ρ	K  =  k/ρC
Optical scattering coefficient	Optical anisotropy coefficient	Reduced optical scattering coefficient	Optical absorption coefficient	Refractive index	Thermal conductivity	Heat Capacity	Density	Thermal diffusivity
mm−1	—	mm−1	mm−1	—	W · m−1 °C−1	J · kg−1 °C−1	kg · m−3	mm2 · s−1
12.44	0.91	1.12	0.017	1.4	0.51	3630	1046	0.134

To further refine the simulations to design NHP optodes, we have used a realistic geometry based on anatomic MRI images of a primate brain segmentated into grey matter, white matter, and cerebrospinal fluid. The simulation of the motor cortex with an optode composed of four fibers was carried out. The optical and thermal properties used in the simulation are summarized in table 3.

Table 3. Optical and thermal properties used in the segmented primate brain geometry.

	µs	g	$\mu _{{\rm s}}^{\prime}$  =  µs (1  −  g)	µa	n	k	C	ρ	K  =  k/ρC
	mm−1	—	mm−1	mm−1	—	W · m−1 °C−1	J · kg−1 °C−1	kg · m−3	mm2 · s−1
White matter	9.0	0.85	1.35	0.0165	1.4	0.51	3630	1046	0.134
Gray matter	9.4	0.9	0.94	0.08	1.4	0.51	3630	1046	0.134
CSF	2.5	0.9	0.25	0.004	1.4	0.51	3630	1046	0.134

The GATE simulations outputs were the coordinates of photons' absorption within a 5  ×  5  ×  2 cm³ area beneath the fiber tip. For each simulation, 10⁸ optical photons were launched from each fiber and propagated through the tissues. For each set of optical and geometrical parameters, we derived the relative fluence isocontours considering the radial distance to the fiber for absorbed photons and normalizing by the total number of launched photons. Fluence profiles along the direction of fiber axis (z profile) and perpendicular to the fiber axis (y profile at 0.5 mm depth) were both calculated with a spatial resolution of 10 µm. The simulation results were normalized to the real number photon emitted by the fiber.

Fluence in each voxel was derived from the number of photon counts passing through the voxels C_photons (photons per mm³) using the following equation (1).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F={{C}_{photons}}\times \frac{c}{n}\times {{E}_{photon}}\nonumber \end{align} \tag{ 1 }$$

where F is the fluence in mW mm⁻², C_photons is the number of photons per mm³, c is the speed of light in vacuum in mm.s⁻¹, n is the media optical index, i.e. 1.4 for biological tissues, and E_photon is the energy of a photon at 674 nm.

We used a hybrid MC/analytical approach to simulate the dynamic increase of temperature during an optical stimulation. Optical absorption was dynamically recorded in a voxelized geometry and used to derive the dynamic temperature increases in the tissues considering the diffusion modeled on the basis of the classical bioheat model [27]. Details of the implementation of the heat generation and diffusion from optical absorption can be found in [28].

The relative fluence isocontours obtained for a single 500 µm diameter fiber with increasing NA are shown in figure 2(A). Figure 2(B) shows that lower NA lead to a modest increase of light transmission in the axial direction within the first millimeter deep into the cortex. Similarly, the lateral distribution changes with the NA are modest and located close to the light entry point. Higher NA values produce a larger lateral spread of the photons. Yet, due to the strong optical diffusion of brain tissues, the fluence quickly becomes homogeneous as depth increases. Beyond 2 mm from the entry point, the fluence distribution and isocontours are identical for all NA (see insert in figure 2(B)). For all NA, about 1% of the light reaches a 2 mm depth in the phantom, 0.5% and 0.2% for 3 and 5 mm depths respectively. The fluence immediately at the fiber tip is larger for low NA fibers, potentially resulting in stronger local thermal effects. A NA of 0.5 appears to be a good compromise to obtain an optical stimulation of a large volume of cortical tissue at 1 mm depth from the surface. It allows minimizing the fluence for the tissues around the light entry point, while maintaining appropriate light penetration. Consequently, in the following sections, the NA was set to 0.58 corresponding to common available multimode fibers.

Zoom In Zoom Out Reset image size

The relative fluence isocontours for fibers with increasing core diameters are presented in figure 3(A). The corresponding axial and lateral fluence profiles are shown in figures 3(B) and (C). The diffusive nature of cortex leads to fluence homogenization beyond the first millimeter of tissues surrounding the fiber tip. The fluence near the surface of the cortex is significantly higher for fibers with lower diameters. Fibers with larger core result in much lower local fluences due to the distribution of the photons over a larger surface. Considering the lateral fluence profiles, the fluence peaks are much sharper for lower diameter fibers. Consequently, larger diameter fibers should be chosen to minimize heat effects at the brain surface. Still, one should consider the availability of large core fibers and the bulkiness and mechanical flexibility of each light guide (including fiber and ferrule), especially for optodes with multiple fibers arrangement, necessary for the stimulation of large brain tissue volumes. Clearly, several fibers are requested to illuminate a 200–500 µl volume of cortical tissue with sufficient fluence to actuate light controlled channels.

Zoom In Zoom Out Reset image size

Two fibers with a NA of 0.58 and a core diameter of 1 mm were simulated in order to evaluate the influence of the distance between fibers on the fluence distribution. The fluence isocontours shown in figure 4 correspond to the coronal plane through the axis of both fibers and the transverse plane at 3 mm below the cortical surface. The fluence decreases rapidly with depth. Interestingly, the use of fibers allows obtaining a significant number of photons down to several millimeters inside the tissue. For all the simulated configurations with two fibers, about 1% of the total number of photons are absorbed locally at 2 mm depth which is similar to what we observed for a single 1 mm diameter, 0.58 NA fiber (figure 3(A), extreme right). Since both simulations are normalized by the total number of photons, this means that twice more photons are absorbed in the two fibers configuration. In addition, the homogeneity of the fluence distribution depends on the distance between the fibers. Thus, for a 5 mm distance between fibers, the individual fluence patterns of each fiber can still be observed at a depth of 3 mm. The fluence becomes homogeneous for an inter-fiber distance below 4 mm. The lateral extension of the illuminated volume also increases with the inter-fiber distance. To illuminate a large brain volume deep down to several millimeters, a distance of 4 mm between the fibers seems adequate for 1 mm diameter fibers with NA 0.58 because it induces a large and homogeneous illumination coverage.

Zoom In Zoom Out Reset image size

Based on the simulation results of the previous sections, four fibers of 1 mm diameter and 0.58 NA with an inter-fiber distance of 4 mm, were used to simulate an optimized optode. A practical implementation is presented in supplementary data 1. A fluence of 150 mW mm⁻² at the fiber tips was simulated. Figure 5(A) shows the relative fluence maps at different depths starting from the cortical surface down to 3 mm depth. The illumination pattern for each individual fiber can be seen down to 2 mm depth. For deeper planes, the fluence is homogeneous due to the diffusion of light in tissues. Figure 5(B) shows the relative fluence maps of the central coronal plane. Deeper than 3 mm below the surface, the fluence is homogeneous over a lateral extent of at least 8 mm. The fluence decreases with the depth, yet it remains above the opsin stimulation threshold of 5 mW mm⁻² down to 5.2 mm depth.

Zoom In Zoom Out Reset image size

Figure 6 shows the temperature increase in the axial planes at depths corresponding to the fluence maps of figure 5(A). As expected, the temperature increase is maximum at the tip of each fiber, though it is limited by the use of large diameter and large NA fibers. Temperature increase of 0.65, 0.35 and 0.1 °C are obtained respectively at depths of 0, 1 and 2 mm below each fiber while the temperature increase below the center of the optode is limited to 0,1 °C (profile on figure 6(B)).

Zoom In Zoom Out Reset image size

In order to illustrate the ability of the developed simulation tools to provide insights for translational application of optical stimulation, we have used a more realistic input geometry based on segmented MRI images of a non-human-primate brain. The four fiber optode is positioned over the cortex targeting motor cortex tissues assuming to contain neurons transfected with a red shifted inhibitory opsin. Figure 7 shows coronal slices of the brain superimposed with the fluence and temperature increase distributions.

Zoom In Zoom Out Reset image size

The fluence profile on figure 7(B) is close to the profile observed in the homogeneous geometry though the 5 mW mm⁻² threshold is valid only down to 4.65 mm almost 1 mm above the same limit for the homogeneous simulations (figure 6(C)). This threshold corresponds to a volume of activated tissues of 285 mm³ composed 117 mm³ of GM, 151 mm³ of WM and 17 mm³ of CSF (values of V₁ are provided in supplemental table 1). The fluence profile is smooth contrarily to the temperature increase profile that reflect the inhomogeneities of photon absorption through the different tissues types. The absorption is very low in CSF, then, as the photons enter the brain they encounter layers of gray and white matters.

GATE was used in the present study to obtain fluence and resulting temperature increase maps. Other public MC code are available for the simulation of photon transport Monte Carlo multi layered (MCML), MCextreme or tMCimg [26, 42, 43], but do not allow to simulate thermal effects simultaneously. MCML has been used to simulate the transport and detection of photons in the context of brain neurosurgery for the safe and accurate positioning of electrodes [44] and is limited to semi-infinite layer geometries. tMCimg and MCextreme allow to simulate more realistic geometries based on segmented MRI or CT images but require data to be formatted into meshed files, which are not straightforward to handle for neurobiologists or neurosurgeons. Recently, OptogenSIM, an open-source platform based on mcxyz MC code, has been developed for optogenetic application. However, this platform is limited so far to optogenetics in mice [23]. GATE offers several key features: new materials with optical and thermal properties can be easily added and analytical as well as voxelized phantoms and sources can be used. We believe that the current limitation of MC simulation results mainly from the limited accuracy of the tissues optical and thermal properties. In vivo human brain optical properties measurements are notably scarce in the literature [31, 33, 34, 45]. Here, we used the brain optical properties taken from the only study, to our knowledge, that was carried out in vivo on human brain tissues [31]. Very promising direct local measurement approaches have been recently developed and could enable to improve the accuracy of the simulations [46]. A previous study used MC simulations and a similar bioheat model to quantify the irradiance dosimetry and heat effects during optogenetics in NHP [47]. In this study the stimulation at 473 and 561 nm using an optrode implanted in the brain tissues was considered. Two different fiber diameters were simulated (10 µm and 200 µm) with very short and low power optical stimulation (1 mW, 1 ms). Although carried out at lower wavelength and in a different geometry (single fiber implanted in brain tissues versus four fibers at the brain surface) the results regarding the influence of the fiber diameter are similar to our results. The larger diameter fiber led to lower local heat effects while maintaining the same volume of stimulated tissues. In the configuration simulated (5 mW of 561 nm light injected into a 200 µm diameter fiber) the volume of activated tissues is estimated to be of about 3 mm³ thus supporting the need of red shifted opsins and higher power and several fibers to stimulate large volumes of brain tissue. Supplementary table 1 provides the volume of stimulated tissues for the two fibers and four fibers configuration simulated showing that volumes of tissues of at least 300 mm³ of brain tissues can be stimulated using the four fiber optode.

The data shown on figure 7 show the ability to perform simulations in realistic geometries based on segmented MRI data. The simulation provides with 3D maps of photon absorption and subsequent thermal effect including heat diffusion. Yet, the optical properties available in the literature for primate brain tissues are scarce and scattered over one order of magnitude [31, 33, 34, 47], and the variability of thermal coefficient for specific tissues is not known. Consequently, MC simulations should be considered as general guidance for NHP optodes design and stimulation protocols rather than as absolute quantitative indicators. More experimental work, including the in vivo measurement of tissue dependent opto-thermal properties is required to facilitate the translation of photomodulation techniques to clinical applications.

The impact of NA and core diameter of fiber were evaluated. The results indicate that larger fiber diameters and larger NA allow minimizing the fluence at the cortical surface and subsequent heating. These results are in agreement with previous calculations carried out for optical stimulations in the rodent brain [12, 23, 48]. Plastic optical fibers of 1 mm diameter with a NA around 0.5 are commercially available. They also have the advantage of being mechanically robust and easy to manipulate. Yet, the high-power sources used to provide the required fluence into deep cortical layers will produce a large amount of heat at the laser-fiber interface that would damage plastic fibers. Silica core fibers are a better choice though they are more fragile and usually have smaller diameters and NA between 0.1 and 0.3. To our knowledge, at least one multimode silica core fiber with a diameter of 1 mm and a NA of 0.5 is commercially available off the shelf (FP100ERT Thorlabs) and several fiber manufacturers offer custom fiber design capabilities. Its maximum bend radius of 40 mm without loss is well adapted to the overall geometry of brain tissue stimulation. In addition, the maximum damage threshold of the fiber materials is of about 250 kW cm⁻², which is several orders of magnitude above what is required for the stimulation of brain tissues. Finally, such fibers can be easily sterilized and are commonly used in preclinical rodent studies.

The use of a single fiber is not adequate to optically modulate a large volume of brain tissue. It would require a too high fluence to ensure that the opsin threshold is reached millimeters away from the fiber tip. The total volume of optically modulated tissue is limited by the diffusion and absorption of light. The normalized fluence isocontours on figure 3(A) show that, for a fiber with a 1 mm diameter and 0.58 NA, the fluence at the fiber tip has to be 1000 mW mm⁻² to obtain a fluence of 5 mW mm⁻² down to 3 mm deep in the brain. Large diameter fibers (above 5 mm) could allow reaching the fluence threshold in larger volumes while minimizing the fluence at the brain surface. However, this would require coupling a laser diode with several hundreds of mW into a non-standard and fragile large diameter glass fiber. A much cost-effective approach is to consider an optode composed of several fibers. In this geometry, photons from each fiber can add up locally to ensure that the fluence is above the opsin activation threshold. Basically, this allows dividing the fluence required at the fiber tip by the number of fiber used. Our results indicate that a distance of 3–4 mm between the fibers allows obtaining a broad and homogeneous fluence down to 3–4 mm depths in cortical tissues. A volume of about 300 mm³ is then appropriately stimulated above the 5 mW mm⁻² threshold if a fluence of 150 mW mm⁻² is injected in each fiber. We considered four fibers, each one coupled to a laser diode. This allows controlling the power output and optical stimulation pattern (duration, frequency, duty cycle) for each individual fiber. This opens the way for complex simultaneous or sequential optical stimulation patterns. In order to validate the feasibility of rapidly designing such multifiber optode we have designed a four fiber optode using off the shelf optical fibers, a surgical grid used in routine neurosurgery in humans and custom designed relay fibers (supplemental figure 1).

Recent studies have provided simulations and experimental data on the thermal effects following optical stimulation of brain tissues using optogenetics. Most of these data have considered solely the excitation in the blue-green range [23, 48] in the rodent brain. In a previous work, we studied experimentally in vivo in rodents the minimally invasive photostimulation using a fiber (0.4 mm diameter and a 0.48 NA) with 5 ms pulses of 638 nm light with a power varying from 100 to 600 mW mm⁻² and a frequency varying from 20 to 60 Hz [49]. We observed a maximum temperature rise of about 0.1 °C to 2.5 °C depending on the experimental parameters. Histological analysis did not show light induced damage at the cellular level when a 200 mW mm⁻² continuous power was applied. We also did not observe light induced but opsin-independent neural activity contrarily to others [49]. Considering safety issues for translational research, it is highly desirable to minimize temperature increase which could have more subtle effects such as local vasodilation of blood vessels, or opsin independent increase of neuronal firing rates leading to behavioral alterations [50, 51]. Our simulation data show limited temperature increase at the brain surface with a fluence up to 150 mW mm⁻² at the output of each fiber. These temperature increases are within the range of normal temperature increase due to environmental changes or physiological activation [52, 53]. Experimental data showed a linear relationship between the total energy deposition and temperature increase [49]. On this basis, MC simulations can easily be scaled up to reflect the stimulation protocol, thus providing an efficient method to optimize irradiance levels and stimulations parameters such as duration and pulses frequencies for the optical stimulation of large and deep cortical volumes in the primate cortex. Simulations based on segmented MRI images of the brain show the influence of the local inhomogeneities of the optical properties on the heat effects.