Quantitative analysis of the optogenetic excitability of CA1 neurons

Introduction Optogenetics has emerged as a promising technique for modulating neuronal activity and holds potential for the treatment of neurological disorders such as temporal lobe epilepsy (TLE). However, clinical translation still faces many challenges. This in-silico study aims to enhance the understanding of optogenetic excitability in CA1 cells and to identify strategies for improving stimulation protocols. Methods Employing state-of-the-art computational models coupled with Monte Carlo simulated light propagation, the optogenetic excitability of four CA1 cells, two pyramidal and two interneurons, expressing ChR2(H134R) is investigated. Results and discussion The results demonstrate that confining the opsin to specific neuronal membrane compartments significantly improves excitability. An improvement is also achieved by focusing the light beam on the most excitable cell region. Moreover, the perpendicular orientation of the optical fiber relative to the somato-dendritic axis yields superior results. Inter-cell variability is observed, highlighting the importance of considering neuron degeneracy when designing optogenetic tools. Opsin confinement to the basal dendrites of the pyramidal cells renders the neuron the most excitable. A global sensitivity analysis identified opsin location and expression level as having the greatest impact on simulation outcomes. The error reduction of simulation outcome due to coupling of neuron modeling with light propagation is shown. The results promote spatial confinement and increased opsin expression levels as important improvement strategies. On the other hand, uncertainties in these parameters limit precise determination of the irradiance thresholds. This study provides valuable insights on optogenetic excitability of CA1 cells useful for the development of improved optogenetic stimulation protocols for, for instance, TLE treatment.


. Introduction
Optogenetics is a technique that can be used to manipulate cellular activity with light (Boyden et al., 2005;Deisseroth, 2011).Specific cell types are sensitized to light via the introduction of gene constructs that encode for optogenetic actuators (Rost et al., 2017).Opsins are light-gated ion channels, pumps or receptors that when illuminated allow hyperor depolarizing currents through the cell membrane (Deisseroth, 2011).The expression of these opsins is controlled by promoters active in the target cells (Rost et al., 2017).Consequently, when applied to the brain, optogenetics provides the possibility to target specific neuronal populations which makes the technique highly promising to treat a variety of neurological disorders, such as epilepsy (Carrette et al., 2015;White et al., 2020).
Temporal lobe epilepsy (TLE) is one of the most common forms of epilepsy.30% of patients with TLE cannot be helped with anti-epileptic drugs (Cela and Sjöström, 2019).Based on the premise that temporal lobe seizures originate in the hippocampus, inhibition of the excitatory subpopulations is a possible strategy for TLE seizure suppression (Tønnesen and Kokaia, 2017;Walker and Kullmann, 2020).Using the anion selective opsin halorhodopsin (NpHR) to inhibit hippocampal pyramidal cells has resulted in a decrease of seizure activity (Tønnesen et al., 2009;Krook-Magnuson et al., 2013, 2015).Excitation of the inhibitory hippocampal interneurons using cation conducting opsins, such as channelrhodopsin-2 (ChR2) and variants, has also yielded promising results (Krook-Magnuson et al., 2013;Ladas et al., 2015;Assaf and Schiller, 2016;Lu et al., 2016).Furthermore, optogenetics can be employed in studies on the initiation and propagation of seizures, and to investigate the role different neuronal populations play in the seizure dynamics (Krook-Magnuson et al., 2015;Cela et al., 2019).
These studies are very promising but the development of clinical applications using optogenetics within the brain still faces some challenges (Cela and Sjöström, 2019;White et al., 2020).A first challenge is having a good understanding of the optogenetic effect, because aforementioned seizure-reduction strategies can have ictogenic effects as well (Assaf and Schiller, 2016;Wykes et al., 2016;Tønnesen and Kokaia, 2017).Another challenge is bringing sufficient light into the brain.The optic frequencies have a poor penetration in brain tissue necessitating in-vivo implantation of the illumination source.In order to minimize structural damage due to intrusion, the light source's dimensions are limited (Vandekerckhove et al., 2020).As a result only a small volume is illuminated with a single optical fiber, which could be insufficient for the human brain (Tønnesen and Kokaia, 2017).Furthermore, light absorption by the brain tissue causes heating.Therefore the light intensity should be limited in order to prevent brain damage (Stujenske et al., 2015;Shin and Kwon, 2016;Owen et al., 2019;Peixoto et al., 2020;Acharya et al., 2021).
Computational modeling can aid in developing a better understanding of the optogenetic response while avoiding invivo animal testing (Williams et al., 2013).In 2006, a first mathematical description of the channel dynamics of ChR2 was published (Nikolic et al., 2006).Since then, several studies have been published that use this and other, more accurate and efficient models in combination with neuron models to design in silico experiments (Nikolic et al., 2009;Grossman et al., 2011;Schneider et al., 2013;Williams et al., 2013;Gupta et al., 2019;Bansal et al., 2021;Schoeters et al., 2021).These experiments allow for easy parameter manipulation and exploration of the stimulation parameter space.Opsin expression can be spatially constrained, expression in different cell types can be tested, and interaction with the 3D light intensity profile can be evaluated.This optic field can be obtained via Monte Carlo based simulations (Wang et al., 1995) that simultaneously can be used to investigate the effect of tissue optical properties (Stujenske et al., 2015;Shin and Kwon, 2016;Acharya et al., 2021).
These optical field studies have shown that the light intensity spatial variation occurs on a neuronal scale.Studies combining optic field simulation with optogenetic response of neurons are scarce.Foutz et al. (2012) investigated this interaction with a multicompartment, neocortical layer V, pyramidal model.Arlow et al. (2013) performed a study on a myelinated axon MRG model.They identified a complex interplay of various simulation parameters on the activation threshold.On the other hand, Grossman et al. (2013) showed that opsin spatial confinement has an impact on action potential generation and propagation.This non-uniform opsin expression has been touched upon by Foutz et al. (2012) as well.In this study, the optogenetic excitability of cornu ammonis (CA1) cells is investigated with the aim of gaining insights that could help guide optogenetic experiments concerning the suppression and initiation of seizures.The novelty is the use of morphologically reconstructed and data-driven biophysical models of CA1 pyramidal neurons and interneurons (Migliore et al., 2018) that are extended to include ChR2(H134R) dynamics, and are subsequently subjected to a Monte Carlo simulated optic field.The effect of opsin expression level and spatial confinement on stimulation thresholds for multiple pulse durations is determined.Furthermore, the impact of various uncertain parameters, e.g., optical field properties, cell to optical fiber orientation and 3D structural cell morphology, is quantified.Based on these results, possible subcellular improvement strategies coupled with ideal optical fiber positioning are identified.

. Methods and materials
The CA1 cell models used to test the optogenetic excitability are described below.Next, the light intensity fields determined via the Monte Carlo (M.C.) method and the opsin model are elaborated on.Finally, the methods and metrics used to analyze the optogenetic response are explained.

. . Neuron models
Four models from Migliore et al. (2018) are used in this study: two models of CA1 pyramidal cells and two CA1 interneurons located in the stratum pyramidale, i.e., a pravalbumin positive basket cell and a bistratified cell.All models have a different three-dimensional structural morphology (see Figure 1A).The cell identification number of the pyramidal cells pyr 1 and pyr 2 is mpg141208_B_idA and mpg141209_A_idA, respectively.The ).This reference frame is applicable throughout the whole study, i.e., the soma is at z = mm, and the somato-dendritic axis is parallel to the z-axis.As depicted, three optrode pitch setups are tested.The color represents the di erent opsin expression locations.(B) The light intensity profile in gray matter with µ a = .mm − , µ s = .mm − , and g = .bistratified cell (int 1 ) has identification number 980513B and the basket cell (int 2 ) has number 060314AM2.Although the pyramidal cells have a different structural morphology, they are classified under the same morphological type (m-type, Markram et al., 2015;Romani et al., 2022).To avoid confusion, in this manuscript we will talk about difference in cell instead of morphology when comparing pyr 1 vs. pyr 2 (int 1 vs. int 2 ).Both pyramidal cells are continuous accommodating, while the electrical type or e-type of the interneurons is continuous non accommodating (Petilla convention; Ascoli et al., 2008;Romani et al., 2022).The models include 10 active conductances: a sodium current, three types of potassium (K DR , K M , K A ), three types of calcium (Ca T , Ca L , Ca N ), two Ca 2+ -dependent potassium currents (K Ca , K Cagk ), and a nonspecific current.Moreover, a simple calcium accumulation mechanism is included with a single exponential decay of 100 ms.The models were fitted to experimentally obtained voltage traces using a genetic algorithm.For more details of the models and the fitting procedure we refer to Migliore et al. (2018).The models themselves can be found on modelDB (https://modeldb.science/)under accession number 244688.The center of cell's soma is at the origin.The somato-dendritic axis of the neurons is aligned to the z-axis.This somato-dendritic axis is the first principal component determined via a principal component analysis on the 3D-positions of all compartments except the axon.In case of the pyramidal cells, the apical trunk is always directed in the positive z-direction.

. . Light field in gray matter
To activate the opsins, the neurons are subjected to light.The light intensity field produced by an optical fiber [100 µm radius, 0.39 numerical aperture (NA); after the optrode in the study of Default values and coefficient of variation (C.V.).Acharya et al., 2021] is determined via the Monte Carlo method.
The used method is based on the direct photon flux recording strategy of Shin and Kwon (2016).The Monte Carlo simulations are solved in cylindrical coordinates for a homogeneous medium as in Stujenske et al. (2015).Because the hippocampus is predominantly a gray matter structure, the in-vivo optical properties of gray matter are selected.The absorption coefficient (µ a ) at 470 nm is obtained by extrapolation of the data reported in Johansson (2010) via a third order polynomial fit on the points between 480 and 550 nm.The scattering coefficient (µ s ) and anisotropy factor (g) are obtained via interpolation from Yaroslavsky et al. (2002).This is often combined into the reduced scattering coefficient The refractive index (n) is 1.36.The simulations are run with 10 7 photons and a radial (dr) and axial discretization (dz) of 5 µm.The result of the Monte Carlo simulation with parameters given in Table 1 is shown in Figure 1B, with corresponding irradiance at the neuron level given in the inset.Throughout the study three pitch orientations with respect to the CA1 cells are tested as depicted in Figure 1A.

. . ChR (H R) opsin
The selected opsin is a genetically engineered Channelrhodopsin-2 (ChR2) variant: ChR2(H134R).It has an enhanced steady-state photocurrent with a slower closing mechanism than the wildtype ChR2.Its peak operation is at a wavelength of 470 nm (Berndt et al., 2011;Mattis et al., 2012).The opsin's kinetics are modeled with the double two-state opsin model, more specifically the RSRS-final reported in Schoeters et al. (2021), because of its improved computational efficiency compared to alternative four state Markov models.It consists of two independent two-state pairs: O (opening and closing of the channel) and R (change in conductance due to dark-light adaptation).The transmembrane current density (mA/cm 2 ) at a single section is determined as follows: (1) rr with g ChR2 the specific conductance (S/cm 2 ), G(V) the rectification function given by Equation 2, V the membrane potential (mV), I rr the irradiance (mW/mm 2 ), and E ChR2 = 0 mV the equilibrium potential.The photocurrents for an optical pulse with duration of 1 and 100 ms with increasing irradiance under voltage clamp recording of -70 mV is shown in Figure 1C.
The opsin is spatially confined to specific neuronal membrane compartments (in this manuscript further denoted as subcellular region).In case of the pyramidal cells, it is located in the soma, or distributed over the axon, basal dendrites (basal), apical dendrites (apic), all dendrites (dend = basal ∪ apic), or all sections (allsec = dend ∪ soma ∪ axon).In case of the interneurons, no distinction between apical and basal dendrites is made.The different regions are illustrated via the color code in Figure 1A.In each simulation, g ChR2 is uniformly spread over the subcellular region of interest.Its value is calculated from a preset total maximum conductance: with A k the total membrane surface area of a subcellular region k.The surfaces are summarized in Table 2. G max are eight points spaced evenly on a log scale between 10 −1 and 10 1.5 µS (for the uniform field nine additional points are included between 10 −1 and 10 2 ).
Finally, the total temporal averaged current at the excitation threshold (TAC; Williams and Entcheva, 2015) is calculated by: where pd is the pulse duration, t 0 is the pulse onset time, T end = max(500 ms, t 0 + pd + 100 ms) (to ensure to include channel closure) and j is every compartmental segment in region k.The pd varied over 5 points logarithmic evenly spaced between 0 and 100 ms (9 pds ǫ [10 −1 , 10 3 ] are used in the simulations with uniform light field).

. . Analyses
The tests and metrics used to analyze the optogenetic response are described below.
. /fncom. . . . .Surface of fiber positions for the activation of neurons As metric to assess the optogenetic excitability, the surface of fiber positions for the activation of neurons (SoFPAN) is determined.This metric is similar to the more well-known volume of tissue activation (VTA; Butson et al., 2007;Duffley et al., 2019).In VTA, the frame of reference is the position of the stimulation device, while here, the frame of reference is the stratum pyramidale with the center of cell's soma at origin.In other words, SoFPAN encompasses the positions at which the optical fiber can be located in order to activate the neuron of interest.This is chosen because of the layered structure of the hippocampus, constraining the cell bodies of the considered neurons to the stratum pyramidale.Because the optical fiber is limited to the plane shown in Figure 1A to reduce the number of simulations, only a 2D plane is explored resulting in a surface instead of a volume being determined.
The intensity threshold (I th ) to elicit an action potential (V>−10 mV) recorded in the soma for a given pulse duration and fiber position is determined first.This threshold value is the intensity at the fiber surface (I fiber ).It is obtained with a titration process using the bisection method for seven iterations (i.e., c i = (a i + b i )/2, with b 0 /a 0 = 10).One hundred and twenty-one fiber positions are evaluated with z taking on eleven linearly spaced values in [−400, 700] µm.For a pitch of π/2, x ǫ [−1,000, 4,000] µm with x = 500 µm.For the fiber pitch of 0 and π , x ǫ [0, 2,500] µm with x = 250 µm.For a given pulse duration, the SoFPAN is then determined as the union of the spatial points for which the threshold is lower than the fiber intensity (true: I th < I fiber ) multiplied by the discretization surface (i.e., 110×250 or 110×500 µm 2 depending on the fiber pitch).A lower bound is determined by counting only the discretization surfaces enclosed by four true-points.In case of the upper bound, the true-false field is dilated first with a 3 × 3 mask before the enclosed-by-four-truepoints regions are counted.For the pitches of 0 and π , the obtained SoFPAN is multiplied by two, such that the maximal SoFPAN for all pitches is equal to 5.5 mm 2 .The SoFPAN is calculated for nine logarithmic evenly spaced fiber intensities 0.1 and 1,000 mW/mm 2 .
If the cell's dimensions are assumed to be negligible with respect to the spatial variation of the intensity field, the irradiance is uniform over the whole neuron.The cell can thus be represented as a single point in space.The estimated SoFPAN under this assumption is denoted as SoFPAN uniform .

. . . Wilcoxon signed-rank test
A single-sided, paired Wilcoxon signed-rank test is used to test whether the results from two populations are significantly different.For the uniform field 20 classes (cell and opsin location: 4×{all, axon, soma, dend} + 2×{basal, apic}) are mutually compared.In case of the M.C. field there are 42 classes pitch, cell, and opsin location: 3× 4×{all, axon, soma} + 2×{basal} .A scoring system is used to identify the most excitable setup.Here, a 1, 0.1, and 0.01 is given if the SoFPAN is statistically greater with p-value below 0.001, 0.01, and 0.05, respectively.The maximum score is consequently 41.

. . . Regression
A two step fitting procedure is performed to analyze the effect of pulse duration and expression level on the intensity threshold as suggested by Williams and Entcheva (2015).First Lapicque's formulation is fit to the TAC.This gives the strength-duration relationship like with electrical direct current stimulation (Reilly, 1998).Second a linear regression is performed with independent variables the log 10 of TAC (µA) and G max (µS).

TAC =
TAC 0 1 − exp(−pd/τ TAC ) ( 6) Here, TAC 0 , τ TAC , a G , a pd , and c are the regression parameters.I th and TAC are the estimators of the threshold intensity and total averaged current, respectively.

. . . Optimal and worst fiber position
The optimal and worst fiber positions for a given pulse duration are defined as the z-position where the depth of activation (amount of positions along the x-axis with I th < I fiber ) is the highest or the lowest, respectively.As a tie-breaker, the z-position is selected where the average of TAC along x is, respectively, the lowest or the highest.

. . . Elementary e ects
The inputs to the simulations can be divided into two categories (Figure 1D).On the one hand, there are the parameters that are known (e.g., optrode radius, NA, and pitch) or controlled by the user (e.g., pd and I fiber ) in an optogenetic experiment.On the other hand, there are variables that cannot be controlled and contain some uncertainty, e.g., the tissue's optical properties (µ a , µ ′ s ), the neuron's structural morphology (cell: pyr/int 1 vs. pyr/int 2 ) and orientation (roll, yaw), and the opsin expression [level (G max ) and location].Because of the expected non-linear interactions and high computational time (±12 h for a 121 position sweep with five pulse durations), the elementary effects method is adopted for global sensitivity analysis (Saltelli et al., 2019).
Six influential factors are investigated, i.e., µ a , µ ′ s , G max , cell, opsin location, and roll.The measure used (µ * ) is the mean of the absolute value of r = 16 elementary effects.Also, the standard deviation of the elementary effects (σ ) is calculated to track interactions and non-linear effects.The r repetitions are sampled with the radial-based design using Sobol's numbers (Campolongo et al., 2011).µ a , µ ′ s , and G max are normally distributed with mean and coefficient of variation (C.V.) given in Table 1.The mean and C.V. of G max are 1 µS and 0.15.The roll is uniformly distributed between [-π , π [ radians.Cell and location are discretely uniformly distributed, with {pyr 1 , pyr 2 } ({int 1 , int 2 }) and {all, axon, soma, basal} ({all, axon, soma}) classes, respectively, for pyramidal (inter-) neurons.

. . Software
Simulations were done with NEURON 8.0.0 (Carnevale and Hines, 2006) and Python 3.9.12. on the HPC system with AMD Epyc 7552 processing units, provided by the Flemish Supercomputer Center.

. Results
The results of this study can be subdivided into three sections.First, there is the optogenetic response under a uniform field.This means that all cell sections receive the same irradiance.This facilitates the analysis of the importance of subcellular regions, pulse durations, and expression levels.Next, the Monte Carlo simulated light field is included.This allows us to indicate the effect of light propagation on optogenetic excitability and provides information on optimal or worst fiber positioning.Finally, an elementary effects study is performed to identify the most sensitive parameters.

. . The optogenetic excitability in a uniform light field
Figure 2A depicts the intensity threshold (I th ) and temporal averaged current (TAC) curves as function of pulse duration (pd) of several subcellular and cell combinations.Both decrease with increasing pd until a saturation value is reached.For the same subcellular region (i.e., allsec), the threshold curve is shifted down in case of pyr 1 vs. int 2 (blue solid vs. blue dashed-dotted line).Staying within the same cell, it can be seen that changing the subcellular region (all → basal or soma) can result in a downward shift, as well.The same effects can be seen on the TAC, but less pronounced.Increasing the expression level (G max ) from 1 to 10 µS (blue solid line with circle markers) results in an expected decrease of the threshold.On the other hand, the TAC remains the same.This is highlighted in Figure 2B where the TAC is constant for all pds for a G max > 1 µS in the allsec-pyr 1 setup.A linear relationship for I th can be observed for G max ≥ 0.518 µS.No threshold could be found for G max < 0.215 µS.This is because the photocurrent saturates at high intensities (see Figure 1C) and therefore cannot compensate for the low G max .At high pulse durations a decrease in TAC is already observed at higher G max (< 1 µS).For the corresponding I th , the photocurrent is biphasic with a peak and steady-state value.At these pulse durations (100 and 1,000 ms) the action potential (AP) occurs during the inactivation from peak to steady state.The occurrence of the peak is much more efficient in eliciting an AP.A large part of the total photocurrent occurs before the AP while this diminishes at higher G max with monophasic photocurrents.Additionally, during the AP, the current drops to zero due to channel shunting.Consequently, the net effect of this drop on the TAC is larger here than in case of no biphasic photocurrents (at low irradiances) or lower pulse durations, where the AP is elicited during the deactivation phase of the photocurrent.
The I th is determined for all pyramidal and continuous nonaccomodating interneuron models of Migliore et al. (2018), along with virtual clones optimized using HippoUnit (Saray et al., 2021).The result of these additional 38 pyramidal and 39 interneuron cell models is shown in Supplementary Figure 1.Overall, for a given subcellular opsin location, the variation is within one order of magnitude, except when the opsin is expressed in the axon.This observation holds for both pulse durations of 1 and 100 ms and maximum total conductances of 1 and 10 µS.The selected cells of this manuscript do not appear as outliers except for the axon of the bistratified cell.Therefore, it can be concluded that the selected cells are representative for the population.
Mutually paired Wilcoxon singed-rank tests are performed to identify the most excitable subcellular region.The p-values of a single-sided test are shown in Figure 2C of the threshold intensity and TAC on the left and right, respectively.Each compared population consists of 153 points (9 pds and 17 G max values).The basal-pyr 1 is the most excitable with a I th significantly (p < 0.001) lower than any other location-cell combination.The same is true for the TAC.The pyramidal cells are more excitable than the interneurons with pyr 1 the most excitable (p < 0.001).Aside from the axon in pyr 2 , the order of the excitability scores is the same in both cells (see Section 2.4.2).The apical dendrites is the least excitable subcellular region of the pyramidal cells (p < 0.001).On average a lower I th is required for the basket cell than the bistratified cell (p < 0.001).Overall a significantly lower threshold is accompanied with a significantly lower TAC.Some exceptions exist where the TAC is significantly greater (e.g., all-pyr 1 vs. soma-int 2 ).The median of the relative change in I th of the location-cell combinations as ranked in Figure 2C is shown in Supplementary Figure 2A.The average of the medians is −22.76%.Therefore, the I th drops every time on average with 22.76% when moving from the least (axon int 1 ) to the most (basal pyr 1 ) excitable location-cell combination.Optimal subcellular expression in a cell can result in I th drops of >75%.Compared to no subcellular specificity (opsin over the whole cell), I th reductions of >60% can be achieved in the pyramidal cells, increasing toward 83 and 92% for the basket and bistratified cell, respectively (cf.Table 3).The full analysis of the median relative differences of all different subcellular-cell combinations is shown in Supplementary Figure 3A.
The two step regression is performed on each location-cell combination separately.The combined results are shown in Figure 2D.The variability on TAC is captured by Lapicque's formulation resulting in a median adjusted coefficient of determination ( R2 TAC ) value of 0.99978.The median rheobase (TAC 0 ) is 3.17 nA and the median time constant (τ TAC ) is 42.95 ms.The variability of I th is well-explained by the two step regression model (median R2 tot = 0.96388).The median parameter values of G max and pd (a G and a pd ) are −1.47 and 0.53, respectively.Even though the latter is positive, I th decreases with pd due to the Lapicque's formulation.Based on these values, G max has thus a stronger impact on I th than the pd.
The input impedance of the cell's sections, measured after 100 ms after initialization at −70 mV, are shown in Figure 2E.The impedance is a proportionality factor on the input current resulting in a given voltage change.Therefore, the voltage change is Median relative change between best and worst subcellular location, and best and no subcellular specificity (all) of all pd and Gmax combinations.Excitability under a uniform light field (I th ) and with light propagation (SoFPAN).
higher for a higher input impedance given a constant input current.
The axons have on average the highest impedance but also the highest spread with many outliers.At the left, it can be seen that the impedance increases from axon hillock to the synapses.The soma has a low impedance compared to the median of the other subcellular regions. .

. The optogenetic response in the Monte Carlo light field
In this section, the cells are subjected to an intensity field produced by a 100 µm fiber.The 3D-profile is obtained via the Monte Carlo method in homogeneous gray matter tissue with the default optical parameters as summarized in Table 1.Unlike in the previous section, the irradiance at the different neuron sections will depend now on their relative position with respect to the fiber and its orientation.Consequently, I th now depends on the position in the reference frame as indicated in Figure 1A.An example is shown in Figure 3A, with the I th map on the left and corresponding TAC on the right, for an allsecpyr 1 with G max = 1.18 µS, pd = 10 ms, and a fiber pitch of π/2.A hotspot is observed around [x = 0.43 mm, z = −0.036mm], where the threshold is the lowest.Even when the neuron lies behind the fiber, i.e., x < 0 mm, it is still possible to excite the cell but at higher I th .Two islands can be observed for the TAC with a factor two difference between minimum and maximum.The farther away from the cell the lower the variation in TAC.
Based on these threshold maps, the surface of fiber positions for the activation of neurons (SoFPAN) can be estimated.This is shown in Figure 3B for the three investigated fiber pitches.The SoFPAN is a measure of the excitability of the subcellular optogenetic configuration in the light field.The uncertainty caused by the discretization on the SoFPAN is indicated by the shaded area (see Section 2.4.1).There is a larger uncertainty at lower intensities due to the rough simulation grid.The soma appears to be more excitable than the opsin in all sections (red vs. blue) with already a non-zero SoFPAN for a I fiber of 0.1 mW/mm 2 .This is true for all fiber pitches.The SoFPAN saturates due to the finite simulation domain.
An optimal and worst z-position for each fiber pitch can be determined, as well.The result for the two cases as above are shown in Figures 3C, D. It can be seen that the positions vary with increasing I fiber .At low intensities, the optimal position is near the subcellular region (cf.soma, red) or the most excitable region (cf.allsec, blue), which is the basal dendrites.At higher intensities, it is better to illuminate the whole simulation domain, explaining the optimal positions at −0.4 and 0.7 mm for the 0 and π pitches, respectively.Once the whole simulation field is excited, the optimal position is completely determined by the average TAC along the xaxis.This can be seen in the sudden changes at the highest I fiber (> 100 mW/mm 2 ).In green the density distribution of all 560 combinations (pd × G max × loc × cell) is displayed.The optimal position is more concentrated with a 0 or π/2 pitch.On the other hand, the worst position is more concentrated with the π/2 and π pitches.The optimal and worst fiber positions are most distant with the π/2 pitch position.
The excitability of the subcellular region is summarized in Figure 3E.The scores are based on the p-values of the singlesided, paired Wilcoxon signed-rank test, as described in Section 2.4.2.At each pitch individually, the subcellular excitability order remains, generally, the same as under the uniform field (cf.Supplementary Figure 4.Only the axons appear to be susceptible to the fiber pitch position.For instance, the axon-pyr 1 combination has an increased (decreased) excitability under the π (0) fiber pitch.Also at pitch 0, the axon of pyr 2 has a lowered excitability with a SoFPAN significantly lower (p < 0.001) than soma-pyr 1 .On the other hand, axon-int 2 has a higher SoFPAN than all-int 2 under the π/2 and π pitches.Basal-pyr 1 is under all pitches the most excitable combination, with a score of 41 at pitch π/2.Overall, pitch π/2 is significantly more excitable than π , which is in turn more excitable than 0 (p < 0.001).The axon-pyr 2 , however, has a significantly higher SoFPAN (p < 0.001) under pitch π than π/2.The median of the relative change in SoFPAN of the locationcell-pitch combinations according to this ranking, i.e., lowest to highest score, is shown in Supplementary Figure 2C.The average of medians is now 10.45%.Therefore, the SoFPAN increases every time on average with 10.45% percent when moving from least  (axon int 1 , pitch = 0) to most (basal pyr 1 , pitch = π/2) excitable location-cell-pitch combination.The same analysis is performed when restricted to a single pitch.The average of medians of the relative changes in SoFPAN increase toward 56.86, 25.89, 59.18% for a pitch = 0, π/2, and π , respectively.Matching the ideal fiber location to the subcellular expression of a single cell can result in doubling of the SoFPAN compared to the worst combination.Even a tenfold increase is observed in case of the bistratified cell.Compared to subcellular unspecificity (all), SoFPAN increases of 33-100% can be obtained by specifying subcellular expression (cf. Table 3).The full analysis of the median relative differences of all different subcellular-pitch combinations for each cell separately is shown in Supplementary Figure 3B.
The SoFPAN for a light field that is considered to be uniform over the whole cell (SoFPAN uniform , see bottom Section 2.4.1) is calculated, as well.The relative error compared to the SoFPAN under a Monte Carlo field, i.e., (SoFPAN uniform -SoFPAN M.C. )/SoFPAN M.C. , is shown in Figure 3F.The errors is calculated for all the G max , pd, pitch, and I fiber combinations.The error of the soma is negligible.The outliers are not shown.Of these, only 0.85% produce a relative error above 5%.These results validate the method as the soma should not depend on the M.C. field as its section is only one point in the 3D-space.The SoFPAN of basal dendrites gets overestimated, with 18.17% having a relative error above 20%.In case of the pyramidal cells and with the opsin distributed over the whole cell, the SoFPAN is predominantly underestimated and the interquartile range (IQR) of the relative error is 60% of the IQR for the basal dendrites.The estimation for the axon of the int 1 cell is the worst with a median relative error of −50%.At least in one test case of each subcellular region the error is either −100% or infinity (see Table 4).

. . Parameter uncertainties
There is uncertainty on various parameters used in this study.The optical parameters depend on multiple factors, e.g., tissue and wavelength.The absorption coefficient is extrapolated which introduces an uncertainty, as well.Moreover, the values of gray matter are used while different gradations exist.The effect of a change in optical parameters on the light field in homogeneous tissue, is illustrated in Figure 4A.On the left, the field as used in the section above is depicted.The effect of an increased absorption and decreased reduced scattering coefficient is shown in the middle.The result is a more conical field with higher degradation.A more round field is obtained when the reduced scattering coefficient is higher (cf.right).Also at the cell level, there are multiple sources of uncertainty.In experimental setting, the opsin expression G max will not be exactly known.Moreover, the subcellular location will probably not be discrete as used in these simulations.Finally, the morphology of the tested cells and its orientation (cf.roll) are fixed.To address the impact of these uncertainties on the output, a global sensitivity analysis is performed.The used approach is a screening method: the Elementary Effects test.
The influence of these six parameters, i.e., cell roll, µ ′ s , µ a , opsin subcellular location (loc), cell model (cell), and G max , on the SoFPAN for the three fiber pitches and two cell types (pyr and int) is investigated.For each fiber pitch and cell type, the elementary effects test (EET) is repeated for 5 pd and 9 I fiber combinations (these 45 combinations correspond with 100% in Figures 4B, D).The rank according to the µ * -measure is summarized in Figure 4B.For certain pd and I fiber no differentiation could be made based on the µ * -measure, explaining the bar height >100% at rank 0. The subcellular location has most frequently the highest impact on the excitability for all six test cases.This is followed by G max in second place (rank 4).The roll and µ a have on average the lowest impact.However, in case of the pyramidal cells and π fiber pitch, the roll has occupied the highest rank for some (pd, I fiber ) combinations.On the other hand, µ ′ s is more important when the fiber pitch is 0 for the pyramidal cells.The µ * and σ measures of three cases of the pyramidal cells with pd of 10 ms are shown in Figure 4C.The circles indicate the setup where the roll has the highest rank, i.e., pyramidal cell with π pitch and 1 mW/mm 2 .It can be seen that even though it has the highest rank, its measures are in the same range as the other two cases.The diamonds represent the EET of a π/2 pitch at I fiber = 1, 000 mW/mm 2 .Here, the reduced scattering coefficient has the highest impact.The non-linear and interaction effects are higher in this case, reflected by the higher σ .For the 0pitch with I fiber = 1, 000 mW/mm 2 (cross), the location is ranked highest with a clear difference in µ * .The effect of cell appears to be more linear than the other parameters indicated by its shift toward lower σ values.
The same analysis is performed on the optimal position.The rank is summarized in Figure 4D.For the interneurons, the subcellular location stays dominant.There is, however, a clear shift toward the optical parameters for the pyramidal cells.Here, the reduced scattering and absorption coefficient are most often ranked highest for the π and 0 pitch, respectively.Also, in case of the π/2 pitch, the absorption coefficient is more important for the optimal position than it is for the SoFPAN.The cell and roll appear to have the lowest effect.While, on the other hand, the cell is important in case of the interneurons.The optimal position for both cell types normalized over all simulations with allsec subcellular location of the EET study, is shown in Figure 4E.These exclude the preset subcellular selectivity.For the interneurons this is more smeared out and a focus toward the stratum pyramidale is observed.On the other hand, for the pyramidal cells there is a clear preference for a position such that the majority of the light reaches the stratum oriens region.At pitch π this is more smeared out due to the possibility to retract (z more positive) the fiber at higher intensities to illuminate a bigger region in the stratum oriens.This is limited for the 0 pitch explaining the high peak at −0.4 mm.

. Discussion
In this study we focused on the optogenetic excitability of CA1 cells.We attempted to not only gain more insight into the effect of the various stimulation and uncertain parameters but also to identify strategies for increased optogenetic efficiency.These insights are of interest for the development of better stimulation protocols that can be used as treatment for TLE.A broad view is adopted where both the excitability of pyramidal and interneurons is investigated.Even though excitation of inhibitory neurons is one of the two main investigated strategies as treatment of TLE, insights in the excitability of pyramidal neurons can be of interest as well.Like with electrical deep brain stimulation, the latter could be used as counter-irritation (Carrette et al., 2015) with various modes of action (Sprengers et al., 2020) that can be tested.Moreover, stimulation of both types could be beneficial for restoring the excitation-inhibition balance (Wykes et al., 2016).
Frontiers in Computational Neuroscience frontiersin.org

. . Excitability of spatially confined opsin expression
The results show that the optogenetic excitability of CA1 cells depends on various parameters.The irradiance threshold ranges over multiple order of magnitudes.As expected, an increase in expression level (G max ) or pulse duration (pd) results in a decrease of the intensity threshold (I th ).There is also a clear dependence on the subcellular region of opsin expression and variance among different cells.There is no single explanation for the relative excitability of the considered subcellular opsin locations, due to the complex interplay of many non-linear relationships.By comparing the membrane areas and impedances (cf.Table 2 and Figure 2E) of the subcellular regions some observations can be made.For a fixed G max , the specific channel conductance (g ChR2 ) is locally higher for regions with a lower total surface area.Thus, for the same I fiber , there will locally be a higher depolarizing current to elicit an action potential (AP).This combined with the fact that the AP is measured in the soma (therefore does not have to travel through the cell), explains why the soma-confinement is highly ranked in each cell.The observation concerning the locally raised channel conductance also holds when comparing basal dendrites with apical, all dendrites, and all sections.An argument for why confinement to the basal dendrites is more excitable than the soma could be found by comparing their input impedances.For a fixed depolarizing current, a higher impedance results in a larger membrane depolarization.Because the impedance of the basal dendrites is significantly higher than that of the soma (log-scale in Figure 2E) it will facilitate AP initiation.However, this contradicts the rank of axon-pyr 1 .Finally, there is the channel distribution inside the cell itself.The ratio of depolarizing (e.g., Na + and Ca 2+ ) and hyperpolarizing (e.g., K + ) channels defines the membrane threshold for AP initiation.This ratio is double in the axon of pyr 2 compared to pyr 1 , motivating the low rank of axon-pyr 1 .These observations are in agreement with the findings of Foutz et al. (2012).
Confinement to the basal dendrites of pyr 1 is the most excitable, while the highest I th is required for the axon of int 1 (cf.Figure 2C).Similar ranking is observed in the SoFPAN calculations (cf. Figure 3E).To identify the effect of endogenous channel distribution, the channel distributions of pyr 1 were imposed on pyr 2 (pyr 3 ) and vice versa (pyr 4 ).The mean relative difference of all pd and G max combinations is shown in Table 5.The I th of pyr 3 drops with ∼30 and ∼50% vs. pyr 1 and pyr 2 , respectively, for all subcellular regions except for the axon (both cases) and for the basal dendrites in pyr 3 vs.pyr 1 .Contrarily, the excitability of pyr 4 drops (higher I th ), except for the basal dendrites for pyr 4 vs pyr 2 or axon in pyr 4 vs. pyr 1 .Switching 3D structural morphology while maintaining endogenous channel distribution (top two) or vice versa (bottom two) impacts optogenetic excitability.It is clear that the interaction of structural morphology and channel distribution has an impact on optogenetic excitability.The axon subregion appears to be the most susceptible with tripled excitation threshold of pyr 3 vs.pyr 2 .Neuron degeneracy, i.e., the ability to perform the same functioning whilst being structurally different or having different ion channel distributions (Migliore et al., 2018), is something that should be taken into account in determining irradiance thresholds.This is also observed in the variability of I th in Supplementary Figure 1.Still, the basal dendrites region is also the most excitable in pyr 3 and pyr 4 .Combined with rank 1 and 2 for pyr 1 and pyr 2 , respectively, it can be concluded that this is the most effective subcellular target region for opsin expression in CA1 pyramidal cells.
This spatial dependence was also observed in the study of Grossman et al. (2013).With the specific conductance (g ChR2 ) as metric, they determined whole cell illumination to be most efficient, i.e., uniform opsin distribution over whole cell with a uniform light field compared to soma, axon initial segment or apical dendrite confinement.They determined that for a 20 ms pulse and irradiance of 1 mW/mm 2 the required g ChR2 when in all sections is only 6% of that when restricted to the soma.On the other hand, G max was 60% higher.These values are in agreement with our results where the ratio of the specific conductance of all sections to soma targeted expression is 2-4% under the same conditions in the pyramidal cells.However, in our study we argue that ranking should be based on G max , i.e., where the number of opsin complexes is fixed.This translates toward an equal comparison of total elicited photocurrent, while, on the other hand, for a fixed g ChR2 the total photocurrent is scaled by the surface area.As a result, the confinement to the soma is classified here to be more excitable.
After correction for the difference in rectification function [i.e., G(V = −68.83mV) = 1 in Grossman et al. (2011) vs. 0.07 in this study], the absolute values of G max were slightly lower but in the same order as reported in Grossman et al. (2013).For a G max of 1 µS (= 0.07 µS after correction) and with a single channel conductance of 40-100 fS this translates toward expression of 0.71-1.7710 6 opsin complexes.Spread over the whole cell this is ±50 channels/µm 2 but confined only to the soma this rises toward >1,000 channels/µm 2 .This value is higher than the estimated 130 channels/µm 2 based on bacteriorhodopsin expression (Nagel et al., 1995;Foutz et al., 2012) but lower than the indirectly estimated 4.4 10 4 channels/µm 2 by Arlow et al. (2013).With our tested G max values up to 100 µS especially when restricted to the soma, this could pose cellular toxicity problems, if these channel numbers would be achieved (Yizhar et al., 2011;Rost et al., 2017Rost et al., , 2022)).This can be avoided if single channel conductance is increased.
The opsin is in all conditions uniformly distributed but can be restricted to a spatial region.In-vivo this highly specific and discrete separation is not possible.Still, by merging the opsin with signaling and targeting constructs, localized enhancement can be obtained.For instance, the addition of the soma-targeting motif of the somaand proximal dendrite-localized voltage-gated potassium-channel Kv2.1 improves somato-dendritic expression (Rost et al., 2017;Mahn et al., 2018).The real distribution in those cases is not known.A normal distribution could be imposed but this would come with two more degrees of freedom.Therefore, the spatially restricted but uniform distribution is used.Our results encourage localized enhancement and advances in this research direction.A reduction of I th with more than 64% can be achieved via subcellular specificity.This is tempered toward, but still significant, increases in SoFPAN of 33-100%, when light propagation is included.Consequently, if made possible, spatial confinement of opsins to specific membrane compartments could significantly increase optogenetic efficiency.On the other hand, more than 76% reduction of I th between optimal and worst subcellular regions is possible.This is also reflected in SoFPAN, where an ideal subcellular-pitch combination can result in a 1.5-10 fold increase compared to the worst combination (cf.Table 3).A good knowledge of the optogenetic interaction at the subcellular level is therefore required in order to achieve the optimal configuration.

. . Optimal fiber position
Due to the finite size and discrete nature of the test grid (121 points), the SoFPAN, and optimal and worst positions could not be unambiguously determined.An upper and lower bound of the SoFPAN is defined (cf.Section 2.4.1),indicating a larger uncertainty at lower intensities.In case of the optimal fiber position, a tie-breaker based on the TAC is introduced.From a homeostatic point of view, it is ideal that the required result is achieved via the lowest perturbation of normal functioning.High/long transmembrane currents could lead to ion concentration imbalance.Especially when using ChR2(H134R), which has a high H + permeability, this can result in neuron acidification, which in turn can result in decreased neuron functioning or unexpected behavior (Schneider et al., 2013;Mahn et al., 2016;Vierock et al., 2017).Therefore, the positions that generate the lowest TAC are preferred.The results in Figure 3C show that the optimal position is at the depth of the region of subcellular expression or with a focus on the most excitable region (in case of opsin distribution over the whole or majority of the cell).Overall, the rank of excitability between uniform and M.C is unchanged field stimulation.Subcellular excitability appears to be dominant over spatial distribution.This spatial preference was also observed in Foutz et al. (2012).In their L5 pyramidal model, they found the apical tuft and soma to be most excitable.
For the π/2 pitch, the optimal and worst position are the most stationary for increasing intensities.For the other pitches (0 and π ), either the optimal or worst position is more smeared out while they are located more closely to each other for low intensities.Consequently, there is a higher risk for sub-optimal fiber positioning.Combined with the highest excitability according to the SoFPAN (cf. Figure 3E), it can be conclude that π/2 is the better fiber position.

. . Contribution of optical field simulation
This study combined simulations of light propagation and neuronal modeling.Light propagation is simulated using the Monte Carlo method for a uniform medium.The hippocampus is a predominantly gray matter structure.However, there is uncertainty on the exact values of the optical parameters (amplified by inter-and extra-polation).The effect of the uncertainty of these parameters is tested using the elementary effects method.The influence of µ ′ s on the excitability is ranked in the middle, while µ a is ranked lower.On the optimal position they were ranked much higher.The median and maximum µ * on SoFPAN are, respectively 0.14 and 1.71 mm 2 (0.04 and 1.19 mm 2 ) for µ ′ s (µ a ).These parameters have some influence, but are subordinate to the other uncertain parameters such as subcellular location, expression level and cell morphology.Moreover, the need to include the light intensity profile in the neuron simulation was addressed by calculating the SoFPAN from excitation thresholds under uniform illumination, as well.Deviations of more than 20% are observed in more than 25% of the tested setups (soma excluded).Confinement to the basal dendrites result in the lowest percentage (18.17%)while the highest is achieved in case of the axon subcellular expression (34.60%).Overall, investigating optogenetic excitability under uniform field conditions provides a good initial approximation, but, accuracy drops for larger and asymmetrical section distributions.The latter was also observed in the strong pitch dependence of axon excitability in the pyramidal cells.

. . Limitations and future work
This study focuses on single pulse excitation of CA1 cells.The occurrence of other spiking behavior is excluded.Likewise, when calculating the SoFPAN with high I fiber , a cell near the fiber end may exhibit bursting or depolarization block due to intense irradiance (Grossman et al., 2011(Grossman et al., , 2013)).Unlike with electrical stimulation, the photocurrent saturates for high light intensities (see Figure 1C).Therefore, extreme behavior is not expected when short pulses are applied.Additionally, the studied ChR2(H134R) opsin exhibits light-dark adaptation, i.e., the photocurrent is higher for a full dark adapted neuron but decreases toward a steady state value under prolonged illumination.The recovery time is in the order of seconds.During pulsed stimulation, the photocurrent of the first pulse will be higher than that of the subsequent pulses.Therefore, higher irradiances will be required to reliably elicit pulse-locked spiking (Foutz et al., 2012).Future work should test the capability of eliciting reliable spiking when opsin expression is confined to a specific subcellular region.Additionally, since this study focuses on isolated cells that are at rest prior to optogenetic stimulation, it is necessary to investigate if phase locking is possible in a network setting, quantify its impact on excitation thresholds, observe subcellular excitation's effects on cellular and network responses (Youssef et al., 2023), and evaluate its influence on synaptic plasticity.Therefore, the interaction in neuronal networks will be of interest in future work, particularly with a focus on hyperexcitable systems such as those in temporal lobe epilepsy.Clearly, the cell's optogenetic excitability depends on multiple factors.The results in Table 5 show both 3D structural morphology and endogenous channel distribution dependencies.Still, the individual impact and interaction effects are yet to be determined.Obtaining a better understanding will be of interest in future work.
The effective level of opsin expression in-vivo is uncertain.Furthermore, to account for potential improvements in plasma membrane expression (Rost et al., 2017;Mahn et al., 2018), G max is treated as a free parameter.However, due to the presence of inward rectification, it is unclear how the model parameter g ChR2 relates to the actual opsin expression level in-vivo.In our formulation of G(V), all proportionality factors are absorbed by g ChR2 (see Equations 1 and 2).While in the formulation of Grossman et al. (2011) and Williams et al. (2013), G(V = −68.83mV) = 1 and G(V = −76.07mV) = 1, respectively.These rectification functions cause a reduction of g ChR2 with a factor of 14.15 or 12.89 at those specific membrane potentials, compared to our formulation in Equation 2.
As aforementioned, there is still some uncertainty on the tissue optical properties.Different studies have reported values that can differ up to an order of magnitude (Yaroslavsky et al., 2002;Gebhart et al., 2006;Johansson, 2010;Genina et al., 2019).Furthermore, brain tissue is binary classified as either gray or white matter, while tissue gradation is more continuous.The impact of the uncertainty of these parameters on the optogenetic excitability is tested in this study.However, it is investigated only locally near the parameters' reported means (see Table 1), while the reduced scattering coefficient of white matter is reported to be 7-times higher than that of gray matter (Yaroslavsky et al., 2002).Additionally, tissue alterations due to foreign body reaction occur.A fibrous capsule is formed around the implanted fiber as reaction to blood-brain-barrier injury and gliosis caused by the presence of the implanted fiber itself (Vandekerckhove et al., 2020).In future work light propagation in a heterogeneous medium and its effect on excitability could be determined.Mesh and voxel based Monte Carlo algorithms exist that can accurately compute the light field distribution in complex tissues (Yan and Fang, 2020).Moreover, modern deep learning algorithms can be utilized to reduce the inherent stochastic noise of these Monte Carlo simulations (Ardakani et al., 2022).Exploring the propagation of light at different wavelengths, such as for the excitation of redshifted opsins, would also be of interest.The simulated fiber has a flat tip with fixed diameter and numerical aperture.Fiber tapering, flat tip patterning and alteration of the geometric properties can result in improved output coupling or broadened and multi-site illumination.However, this is out of the scope of this research (Kostovski et al., 2014;Vandekerckhove et al., 2020).
New opsins, either natural or genetically engineered, are discovered on a yearly basis (Yizhar et al., 2011;Klapoetke et al., 2014;Vierock et al., 2017Vierock et al., , 2022;;Oppermann et al., 2019).While this is generally beneficial, it can hinder the development of optogenetic tools.Dividing research among multiple opsins may limit the understanding of a single opsin's interactions with neurons and its capabilities.The question arises whether the gathered insights here are transferable to other opsins as well.Previous studies have demonstrated the influence of channel kinetics on factors like irradiance thresholds, spike reliability, and behavior (Grossman et al., 2011;Mattis et al., 2012;Klapoetke et al., 2014).The exact values of I th will thus differ for another opsin.However, these values are already uncertain due to multiple other uncertainties in other parameters like G max .Furthermore, these differences will affect all tested cases equally.Therefore, we expect that the observed trends and rankings regarding optogenetic excitability will be applicable across different opsins.In future work, a similar study could be performed focusing on optogenetic silencing with inhibitory opsins like GtACR2 and WiChr (Govorunova et al., 2015;Vierock et al., 2022).
Tissue illumination causes heating.To avoid permanent tissue damage, the local temperature increase cannot exceed 6 • C (Acharya et al., 2021).Moreover, behavioral changes are already possible at lower temperature changes (1 • C).Several neural parameters, e.g., capacitance, ion channel conductance, and transmitter release and uptake, have been shown to be temperature dependent (Kiyatkin, 2019;Fekete et al., 2020;Peixoto et al., 2020;Acharya et al., 2021).The reported SoFPAN values are for fiber intensities up to 1,000 mW/mm 2 .After extrapolation of the change in temperature results reported in Acharya et al. (2021) (Figure 3A), this intensity corresponds with an estimated temperature increase of 4.85 • C after 100 ms.Consequently, temperature-induced changes in optogenetic excitability should be included in future work.

. Conclusion
In this in-silico study, we examined the optogenetic excitability of four CA1 cells using ChR2(H134R).Our findings reveal that, for a fixed amount of opsin channels (G max ), confining the opsin to specific neuronal membrane compartments significantly enhances excitability.This confinement leads to threshold reductions exceeding 64% and up to 100% gains in the surface of fiber positions for the activation of neurons.Additionally, we determined that the perpendicular orientation of the fiber relative to the somatodendritic axis yields superior results.Furthermore, we observed substantial inter-cell variability, with differences in thresholds above 20%.The bistratified cell exhibited the least excitability, while pyramidal cell 1 demonstrated the highest excitability, especially when the opsin is confined to the basal dendrites.These findings highlight the importance of considering neuron degeneracy while developing optogenetic tools.By screening various uncertain parameters, we identified opsin location and G max having the greatest impact on simulation outcomes.Our study showed the advantages of computational modeling coupled with light propagation.An increased excitability is seen with optimal fiber positioning, i.e., perpendicular to the somatic-dendritic axis and focus on the most excitable cell region.Spatial confinement and enhancements of opsin expression levels are promoted strategies to improve optogenetic excitability.However, it should be noted that uncertainty in these parameters limits determining the exact irradiance thresholds.their affiliated organizations, or those of the publisher, the editors and the reviewers.Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

FIGURE
FIGUREThe simulation framework.(A) Graphical representation of the D reconstructed CA cell models used in this study: two pyramidal cells (pyr and pyr ), bistratified (int ) and basket cell (int ).The models are adopted fromMigliore et al. (   ).This reference frame is applicable throughout the whole study, i.e., the soma is at z = mm, and the somato-dendritic axis is parallel to the z-axis.As depicted, three optrode pitch setups are tested.The color represents the di erent opsin expression locations.(B) The light intensity profile in gray matter with µ a = .mm − , µ s = .mm − , and g = . .The irradiance at the cell, with soma µm bellow the fiber, is shown in the inset.(C) The opsin current density at a voltage clamp of − mV for increasing irradiance (light to dark: -→ mW/mm ) of a (left) and ms (right) optical pulse, indicated in gray.(Bottom) The peak (solid) and steady state (dashed) current densities as function of irradiance for a and ms optical pulse in green and blue, respectively.(D) The simulation flowchart with input and uncertain parameters.
. The irradiance at the cell, with soma µm bellow the fiber, is shown in the inset.(C) The opsin current density at a voltage clamp of − mV for increasing irradiance (light to dark: -→ mW/mm ) of a (left) and ms (right) optical pulse, indicated in gray.(Bottom) The peak (solid) and steady state (dashed) current densities as function of irradiance for a and ms optical pulse in green and blue, respectively.(D) The simulation flowchart with input and uncertain parameters.

FIGURE
FIGUREThe optogenetic response under a uniform light field.(A) The intensity threshold (I th ) and corresponding total temporal averaged current (TAC, top and bottom, respectively) as function of pulse duration.All lines shown are for a G max = µS except the circle markers (G max = µS).(B) The intensity threshold and corresponding total temporal averaged current (top and bottom, respectively) as function of G max for allsec-pyr .(C) The p-value of the mutually, paired Wilcoxon signed-rank tests.Single-sided test if class at the row is lower than class in the column for the I th (left) and TAC (right).The color code at the top and left indicate the cell model.The cell opsin location combinations are sorted on their excitability score (see Section . .). (D) The results of the two step regression.The goodness of fit is indicated with the adjusted R-squared measure: R TAC for the Lapicque fit to the TAC and R tot for the linear regression fit to the threshold.Box plots are shown for R -values and regression parameters, calculated for the di erent cell-types and opsin locations.(E) Summary of the input impedance at Hz, (left) the local impedance of each section in pyr , (right) the summarized impedances of the four tested cells.The gray bars are µm.

FIGURE
FIGUREThe optogenetic response in a homogeneous Monte Carlo simulated light field.(A) The threshold intensities (left) of the pyr cell, with opsin distributed over all sections, and corresponding total temporal averaged current (right).In the colormaps, the position of the optical fiber is varied with respect to the soma, which is fixed at (x, z) = ( , ). (B-D) The Surface of Fiber Positions for the Activation of Neurons, the optimal and worst z-position for neuron activation, respectively, as function of fiber intensity.The shaded area in (B) is enclosed by the upper and lower bound.The lines are the result of the pyr cell with opsin location shown in legend at the right top corner.G max = .µS and pd = ms in (A-D).In (C, D) the scatter plots overlay the density, averaged over all cell, location, G max and pd combinations.(E) Excitability score based on paired Wilcoxon signed-rank tests.Each population consists of ( pd × G max × I fiber values) combinations.The pitch of the fiber is illustrated by the orientation of the fiber icon on the left and valid for the whole row.(F) Relative error of SoFPAN, if the light intensity is considered uniform over the neuron: rel.error = (SoFPAN uniform -SoFPAN M.C. )/SoFPAN M.C. , the outliers are not shown.The errors are calculated for all the G max , pd, pitch, and I fiber combinations.

FIGURE
FIGUREResults of the elementary e ects study.(A) Intensity fields with di erent optical parameters; from left to right µ a = ., ., .mm − and µ ′ s = ., ., .mm − .(B, D) The influential parameters on SoFPAN and optimal fiber position, respectively, ranked for pds and fiber intensities (in %).The fiber pitch and cell type are shown on top and left, respectively.(C) The two measures of the elementary e ects of three pyramidal cell setups, i.e., pitch = π, π/ , and , and I fiber = , , , and , mW/mm indicated by circle, diamond and cross, respectively, for a pd = ms.The legend between (B, D) also applies on (C).(E) Normalized optimal z-position over all results with allsec subcellular location of the elementary e ects study.
TABLE Parameters used in the Monte Carlo (M.C.) simulations.
TABLE Total membrane surface area of subcellular regions in µm .
TABLE Summary relative error SoFPAN Monte Carlo vs. SoFPAN uniform (all cells pooled together).
TABLE Mean relative di erence in I th due to three-dimensional structural morphology (top) or by changing channel distribution (bottom; relative with respect to original).