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₁ and pyr₂ is mpg141208_B_idA and mpg141209_A_idA, respectively. The bistratified cell (int₁) has identification number 980513B and the basket cell (int₂) 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₁ vs. pyr₂ (int₁ vs. int₂). 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²⁺-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.

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 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 [μs′=μs(1-g)]. The refractive index (n) is 1.36. The simulations are run with 10⁷ 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.

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²) at a single section is determined as follows:

with g_ChR2 the specific conductance (S/cm²), G(V) the rectification function given by Equation 2, V the membrane potential (mV), I_rr the irradiance (mW/mm²), 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⁻¹ and 10^1.5 μS (for the uniform field nine additional points are included between 10⁻¹ and 10²).

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₀ is the pulse onset time, T_end = max(500 ms, t₀+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⁻¹, 10³] are used in the simulations with uniform light field).

The tests and metrics used to analyze the optogenetic response are described below.

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.

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₁ vs. int₂ (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₁ 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 non-accomodating 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₁ 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₁ the most excitable (p < 0.001). Aside from the axon in pyr₂, 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₁ vs. soma-int₂). 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₁) to the most (basal pyr₁) 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 ( R̄2TAC) value of 0.99978. The median rheobase (TAC₀) 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 R̄2tot = 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 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.

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 allsec-pyr₁ 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.036 mm], 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². 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 x-axis. This can be seen in the sudden changes at the highest I_fiber (>100mW/mm²). 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 single-sided, 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₁ combination has an increased (decreased) excitability under the π (0) fiber pitch. Also at pitch 0, the axon of pyr₂ has a lowered excitability with a SoFPAN significantly lower (p < 0.001) than soma-pyr₁. On the other hand, axon-int₂ has a higher SoFPAN than all-int₂ under the π/2 and π pitches. Basal-pyr₁ 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₂, however, has a significantly higher SoFPAN (p < 0.001) under pitch π than π/2. The median of the relative change in SoFPAN of the location-cell-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₁, pitch = 0) to most (basal pyr₁, 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₁ 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).

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². 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². 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 0-pitch with I_fiber = 1, 000 mW/mm² (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.

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₁. Finally, there is the channel distribution inside the cell itself. The ratio of depolarizing (e.g., Na⁺ and Ca²⁺) and hyperpolarizing (e.g., K⁺) channels defines the membrane threshold for AP initiation. This ratio is double in the axon of pyr₂ compared to pyr₁, motivating the low rank of axon-pyr₁. These observations are in agreement with the findings of Foutz et al. (2012).

Confinement to the basal dendrites of pyr₁ is the most excitable, while the highest I_th is required for the axon of int₁ (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₁ were imposed on pyr₂ (pyr₃) and vice versa (pyr₄). The mean relative difference of all pd and G_max combinations is shown in Table 5. The I_th of pyr₃ drops with ~30 and ~50% vs. pyr₁ and pyr₂, respectively, for all subcellular regions except for the axon (both cases) and for the basal dendrites in pyr₃ vs. pyr₁. Contrarily, the excitability of pyr₄ drops (higher I_th), except for the basal dendrites for pyr₄ vs pyr₂ or axon in pyr₄ vs. pyr₁. 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₃ vs. pyr₂. 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 thus 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₃ and pyr₄. Combined with rank 1 and 2 for pyr₁ and pyr₂, 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² 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.83 mV) = 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.77 10⁶ opsin complexes. Spread over the whole cell this is ±50 channels/μm² but confined only to the soma this rises toward >1,000 channels/μm². This value is higher than the estimated 130 channels/μm² based on bacteriorhodopsin expression (Nagel et al., 1995; Foutz et al., 2012) but lower than the indirectly estimated 4.4 10⁴ channels/μm² 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., 2017, 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 soma- and 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.

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.

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² (0.04 and 1.19 mm²) 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.

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, 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.83 mV) = 1 and G(V = −76.07 mV) = 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 red-shifted 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., 2017, 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². 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.