Eight reconstructed morphologies of layer 5 pyramidal neurons of rat somatosensory cortex were downloaded from an online database¹ (Ascoli, 2006; Hay et al., 2013; Cohen et al., 2020); a model cell is shown in Figure 1A. Soma diameter was set to 20 μm in all model PCs. Axonal trees were split into hillock, AIS, unmyelinated axon, nodes of Ranvier, and myelinated axon (see Table 1). The hillock (L = 0–2 μm) was attached to the soma followed by the AIS (L = 35–48 μm) and an unmyelinated axon section up to the first axonal bifurcation (100–150 μm from the soma). Following axonal branches were split into nodes of Ranvier (0.5–1 μm) and internodes (L = 100 × section diameter; Rushton, 1951). Nodes of Ranvier were placed at the beginning of each branch as well as at the end of myelin sections except at terminal branches. The length of each node of Ranvier was dependent on the length of the adjacent myelin section; this was done to prevent self-spiking that could result from a too large nodal area.

34 reconstructed morphologies of mouse alpha RGCs were taken from a previous study (Werginz et al., 2020); Figure 1B shows a representative model cell. Each model neuron was divided into dendrites, soma (D = 14–24 μm), hillock (L = 10–47 μm), AIS (L = 12–33 μm), and distal (unmyelinated) axon (L ∼ 1000 μm).

For each cell, the soma was either modeled as a single spherical compartment or as a sphere approximated by 41 cone-shaped compartments. For the multi-compartment soma approach (Figure 1C), the main axis of the soma was always pointing toward the stimulating electrode to study the gradient of the applied electric field across the soma (Fellner et al., 2019). Responses to intracellular and extracellular stimulation in both cell types were performed in NEURON 7.8 (Carnevale and Hines, 2006); Python 3.8² was used to control the simulations. Compartment length was between 1 and 2 μm in the axon and below 10 μm in dendritic sections. Extracellular stimulation was modeled via NEURON’s “extracellular” mechanism. A monophasic cathodic pulse with a duration of 0.1 ms was used in the majority of all experiments. Charge-balanced biphasic pulses with different cathodic/anodic ratios were used in a subset of simulations.

For PCs, the biophysical properties (i.e., ion channel kinetics and densities, Table 1) of the soma and dendrites were specified based on experimental data (see Table 1 and cell 5 data of Almog and Korngreen, 2014) and the axon kinetics were adapted from Mainen and Sejnowski (1996). Leak conductivity, intracellular resistivity, and specific membrane capacitance were set to 0.39 pS/μm², 120 Ω.cm, and 0.6 μF/cm², respectively, except for the myelinated part of the axon. The nodes of Ranvier had an increased leak conductivity of 200 pS/μm², whereas the myelinated internode had a lowered membrane capacitance of 0.04 μF/cm². Model temperature was set to 34°C. Model baseline function, including action potential initiation and (back)propagation as well as spike shape, was tested by intracellular current injections into the PC soma and compared to Almog and Korngreen (2014).

Biophysical properties of RGCs were slightly modified from previous studies (Fohlmeister et al., 2010; Werginz et al., 2020; Table 2) without changing baseline function, i.e., AP initiation in the AIS, forward and backward propagation of spikes into the distal axon and somatodendritic compartment, respectively. Leak conductivity, intracellular resistivity, and specific membrane capacitance were set to 2.5 pS/μm², 143 Ω.cm, and 1 μF/cm², respectively. Model temperature was set to 33°C.

Stimulating electrodes were placed in close vicinity of the soma within distances of 15, 30, 45, 60, 100, and 200 μm to the soma center (Figures 1A,B, green circles). The extracellular potential V_e was calculated for a point source in a homogeneous infinite medium. In addition, disk electrode stimulation was applied to RGCs in a semi-infinite medium. The extracellular resistivity ρ_e was smaller for cortical tissue (300 Ω.cm; Rattay and Wenger, 2010) than for retinal tissue (1000 Ω.cm; Werginz and Rattay, 2016).

For point source stimulation the extracellular potential was calculated as (Rattay, 1999)

with I_Stim being the stimulus current applied to the electrode and r representing the Euclidian distance from the compartment center to the electrode. In the multi-compartment soma configuration, r was computed as the distance to the compartment surface (Figure 1D).

The electric potential generated by a disk electrode was calculated as (Newman, 1966; Wiley and Webster, 1982; Rattay, 1988)

with a, r, and z being the electrode radius, the radial, and the axial distance to the electrode, respectively.

Comparison between two groups were performed by the Kruskal-Wallis test. Significance levels were set as follows: p < 0.05 ^∗, p < 0.01 ^∗∗, p < 0.001 ^∗∗∗. Boxplots use standard notation (1st Quartile, Median, 3rd Quartile). All statistical analyses were performed in Python 3.8.

To examine the contribution of the soma to LT and UT, model neurons were either equipped with a soma consisting of a single sphere or of 41 truncated cone compartments (Figures 1C, 2A left). The right part of Figure 2A shows the activating function, a predictor for excitation along the neural membrane (Rattay, 1999), for a cathodic electrode current of 1 μA. For the multi-compartment soma, the activating function is positive (red, indicating depolarization) in compartments close to the electrode, whereas compartments further away from the electrode have a negative activating function (blue) indicating hyperpolarization. In contrast, for the single-compartment soma, the activating function is small (white), resulting from no reflection of the extracellular potential gradient across the soma.

LTs and UTs in PCs (n = 8) and RGCs (n = 34) were computed at point source distances of 15, 30, 45, 60, 100, and 200 μm to the soma center (Figures 2B,C). As expected LTs and UTs increased with electrode distance; LTs were in similar ranges in RGCs and PCs for small distances (∼1 μA for 15 μm distance). However, in PCs LTs were higher but still within an order of magnitude compared to RGCs for larger distances. In contrast, UTs were almost 10-fold larger in PCs vs. RGCs.

No significant difference between both soma configurations, neither in LTs nor in UTs, was observed in PCs for all distances indicating little contribution of the trans-somatic electric field to UTs (Figure 2B, orange vs. gray). In contrast, for RGCs, LTs in the multi-compartment configuration were significantly lower than in the single-compartment soma model for electrode distances ≤60 μm. Also, UTs were significantly lower in the multi-compartment soma configuration for electrode distances of 15 and 30 μm (Figure 2C, orange vs. gray).

Threshold ratios (UT/LT) in PCs were calculated for the single- and multi-compartment soma configuration (Figure 2D). The ratios were smallest for the shortest investigated electrode distance, increased for distances up to 60 μm and decreased again for electrode-to-soma distances of 100 and 200 μm. No significant difference was found in ratios between both soma configurations at any examined electrode distance. Figure 2E displays threshold ratios for RGCs for single and multi-compartment soma configurations, respectively. Threshold ratios increased monotonically with electrode distance and a significant difference between both soma configurations could be observed in RGCs for the closest electrode distance suggesting a significant contribution of the soma in determining the stimulus window.

We were also interested in whether single morphological features influence LTs, UTs as well as threshold ratios for a given distance. We computed linear correlations between single anatomical parameters (e.g., AIS length, soma diameter, etc.); r² values for each correlation are shown in Supplementary Tables 1,2. We did not observe a clear trend for any feature tested, with an exception of PC dendritic area vs. LT that resulted in r² values >0.4 for all distances (Supplementary Table 1). Additionally, AIS length was somewhat correlated to thresholds and threshold ratios for small electrode distances, similar to results obtained by previous studies that showed a correlation between AIS length and LTs (Jeng et al., 2011; Werginz et al., 2020). These results indicate that LTs, UTs as well as threshold ratios are not dependent on single anatomical properties but to a combination of anatomy, biophysics as well as axonal geometry relative to the electrode location.

In summary, we found up to 300-times higher UTs in PCs compared to their LTs whereas threshold ratios were below 50 for all electrode distances in RGCs. A significant soma contribution was observed in both LTs and UTs for RGCs, particularly for small electrode distances, in contrast to a negligible soma contribution in PCs regardless of electrode distance.

In our initial experiments, spiking activity was detected in the distal axon in RGCs as it is the sole output pathway in these cells. PCs, however, have complex axons with numerous collaterals, and we decided to detect APs at the first node of Ranvier approximately 100–150 μm distant to the soma where the axon starts to bifurcate. In a second set of simulations, we questioned whether partial spiking, i.e., APs in one part of the axon but no spiking in other portions of the cell, could occur in PCs. Surprisingly, a complete block did not occur in any of the eight investigated PCs, but depending on the electrode distance to the soma, at least small portions of the cell generated an AP, even at amplitudes above previously determined UTs. Figure 3 demonstrates this partial block phenomenon for one model cell (PC7) when the electrode was at a distance of 15 μm. At LT, the axon and soma generated an AP which actively backpropagated into the dendritic tree [Figure 3A, note the broad dendritic calcium spike (Almog and Korngreen, 2014)]. When the pulse amplitude was set to 50% of the UT the axonally initiated AP did not propagate across the soma and backpropagation into the dendrites failed (Figure 3B). At stimulus levels equal to and higher than UT, some parts along the axon still generated APs while most portions of the cell were in a blocking condition (Figures 3C,D).

To quantify the partial spiking in axon collaterals, we computed the percentage of spiking nodes of Ranvier at stimulus amplitudes above LT (Figure 3E). For the smallest electrode distance of 15 μm, similar to the results from Figure 2D (see Figure 3E, triangles), the percentage of spiking nodes dropped rapidly when stimulus amplitude exceeded ∼50 × LT. For larger distances the spiking percentage decreased but eventually plateaued at approximately 25% indicating partial spiking even at stimulus amplitudes >500 × LT.

Taken together, no complete block was observed in PCs for all examined distances and AP blockage was mostly observed in regions with strong hyperpolarization in response to high-amplitude stimulation. The partial block is mostly specific to extracellular stimulation of the cell and originates from an inhomogeneous reflection of electric field on the transmembrane voltage of axon collaterals and inhomogeneous activating function.

In the next set of experiments, we were interested in the location where APs are initiated during stimulation close to the soma. The site of spike initiation (SSI) was defined as the compartment whose membrane voltage first crossed 0 mV. The SSI was examined in the multi-compartment soma configuration. When applying a monophasic cathodic pulse close to the soma (15 μm distance), following the activating function (cf. Figure 2A), a distinct polarization pattern of the somatic compartments could be observed. For both cell types, somatic compartments close to the stimulating electrode were strongly depolarized during the pulse whereas compartments on the far side of the soma were hyperpolarized (Figures 4A,C, red and blue arrows). Interestingly, the strong depolarization at the soma did not result in an AP directly but compartments along the proximal axon initiated the spike (Figures 4A,C, light blue). This region is called the axon initial segment and, due to its high density of sodium channels, has been shown to be the SSI during a variety of stimulation conditions in both PCs and RGCs (Rattay and Wenger, 2010; Werginz et al., 2014). The AIS-induced AP subsequently propagated bidirectionally into the distal axon as well as backward into the soma and dendrites (Figures 4A,C, thick black line).

The experiment was repeated for all model neurons at electrode distances of 15, 30, 45, 60, 100, and 200 μm. The SSI was confined to the AIS in all cases, with small variations for different electrode distances. This is shown by plotting the distance between the SSI and the soma vs. the sum of the distance between the AIS and the soma and the length of the AIS (Figures 4B,D). The solid gray lines show the best-fit logistic functions with r² values of 0.74 and 0.91 for RGCs and PCs, respectively.

Our results show that APs are always initiated within the AIS for micro-stimulation in the vicinity of the soma, regardless of electrode-to-soma distance. For short AIS lengths, the AP initiated close to the distal end of the AIS, whereas in cells with longer AIS which are located far from the soma the SSI was shifted to the center of the AIS.

So far, our simulations were restricted to stimulation with a point source and monophasic pulses from defined electrode locations. However, parameters such as electrode size and pulse duration cannot be chosen freely in neural implants due to hardware, software, and energy constraints. The placement of stimulating electrodes is also limited by constraints of the surrounding biological tissue. In the retina, for example, implants can be placed on the epiretinal surface to stimulate RGCs directly. In such implants an array of disk electrodes interfaces with the targeted RGCs. Therefore, the impact of the size of disk electrodes on LTs, UTs and the UT/LT ratios was examined. The electrode-to-soma distance was set to 15 μm as in our previous simulations as the strongest effect of the trans-somatic electric field on LTs and UTs was observed at smallest distances. Figure 5A shows UT/LT ratios for electrode diameters ranging from 10 to 200 μm and compares them to the point source electrode results. For electrode diameters in the range of 10–50 μm threshold ratios increased, however, for the largest electrode tested we found decreasing threshold ratios that were in a similar range as for point source stimulation.

To prevent potentially tissue damage due to charge accumulation, most studies use charge-balanced biphasic pulses. Therefore, we studied threshold ratios when a biphasic cathodic-first pulse was applied. The cathodic phase duration was held constant at 0.1 ms, whereas the following anodic pulse duration was varied from 0.1 up to 1 ms with its amplitude adjusted for charge-neutral stimulation. Threshold ratios are shown for phase duration ratios (D_cat/D_ano) ranging from 0.1 to 1 and compared with the monophasic pulse used in our previous simulations (Figure 5B). Our results indicate up to threefold higher threshold ratios for symmetric and close to symmetric pulse shapes (pulse ratio > 0.6), whereas this effect decreased for longer charge balancing pulses. For the longest anodic pulse of 1 ms the biphasic stimulus became pseudo-monophasic and its UT/LT ratio was close to the monophasic case.

Targeted stimulation of PCs with neural implants is challenging as the precise electrode location cannot be determined during surgery. The inserted micro-electrodes are placed within the targeted cortical layer(s) and thus can be located randomly in space relative to PC somata. In order to quantify the effect of the electrode-to-soma positioning, LTs, UTs as well as threshold ratios were calculated in the eight PC model neurons for six electrode locations around the soma. Electrode distance was set to 15 μm to the center of the multi-compartment soma (Figure 6A). LTs and UTs for all investigated positions are shown for each cell in Figure 6B. Interestingly, for each of the eight model neurons tested, one electrode location resulted in significantly lower LTs and UTs (outliers in Figure 6B). More detailed analysis revealed that these outliers were linked to the electrode distance to the hillock (Figure 6C). Similar LTs and UTs were observed for all positions, except for the closest electrode location to the axon hillock, which resulted in the lowest LTs and UTs (red ellipse). Figure 6D shows threshold ratios of the six investigated electrode locations for each cell. Median ratios ranged from 40 to 80 similar to threshold ratios for an electrode distance of 15 μm (cf. Figure 2D).

In a final set of experiments, we explored how the UT affects the activation pattern of a population of RGCs. Figure 7A shows the transition from LT to UT for 100 RGCs distributed on the epiretinal surface when stimulated by a 50 μm disk electrode. The pulse was biphasic with a pulse ratio of 0.1, i.e., a 0.1 ms cathodic stimulus followed by a 1 ms balancing pulse. Small amplitudes led to focal activation of RGCs which were located slightly offset from the stimulating electrode (Figure 7A, 1 and 2 μA). This can be explained by the low-threshold region at the distal AIS in RGCs (Fried et al., 2009; Jeng et al., 2011; Werginz et al., 2020). Increasing the stimulus amplitude led to an enlarged region of activation (Figure 7A, 4 and 8 μA). At high amplitudes the UT resulted in a different activation pattern with RGCs that had lowest thresholds being prevented from firing APs first (Figure 7A, 1 vs. 16 μA). Highest amplitudes generated a ring of activated RGCs around the stimulating electrode creating a non-responding region around the stimulating electrode (Figure 7, 32 and 64 μA). We were also interested in the total number of activated RGCs for different stimulus amplitudes. Figure 7B shows a monotonic increase of activated cells when the stimulus amplitude was increased from 1 to 64 μA. Interestingly, the number of activated cells still increased when RGCs close to the electrode were already blocked (gray shading). However, despite a larger number of activated cells at higher amplitudes, the different activation pattern (circular versus ring-shaped, Figure 7B, top) is not expected to generate well-defined visual percepts.