Abstract
Single-photon optogenetics enables precise, cell-type–specific modulation of neuronal circuits, making it a crucial tool in neuroscience. Its miniaturization in the form of fully implantable wide-field stimulator arrays enables long-term interrogation of cortical circuits and bears promise for brain–machine interfaces for sensory and motor function restoration. However, achieving selective activation of functional cortical representations poses a challenge, as studies show that targeted optogenetic stimulation results in activity spread beyond one functional domain. While recurrent network mechanisms contribute to activity spread, here we demonstrate with detailed simulations of isolated pyramidal neurons from cats of unknown sex that already neuron morphology causes a complex spread of optogenetic activity at the scale of one cortical column. Since the shape of a neuron impacts its optogenetic response, we find that a single stimulator at the cortical surface recruits a complex spatial distribution of neurons that can be inhomogeneous and vary with stimulation intensity and neuronal morphology across layers. We explore strategies to enhance stimulation precision, finding that optimizing stimulator optics may offer more significant improvements than the preferentially somatic expression of the opsin through genetic targeting. Our results indicate that, with the right optical setup, single-photon optogenetics can precisely activate isolated neurons at the scale of functional cortical domains spanning several hundred micrometers.
Significance Statement
Sensory features, such as the position or orientation of a visual stimulus, are mapped onto the surface of the cortex as functional domains. Their selective activation, which may enable eliciting complex percepts, is intensively pursued for basic science and clinical applications. However, delivering light into one functional domain in an optogenetically transfected cortex results in complex neuronal activity, spreading beyond the targeted domain. Our computational study reveals that neuron morphology contributes to this diffuse response in a cortical layer- and intensity-dependent manner. We show that enhancing stimulator optics is more effective than soma targeting of the opsin in increasing stimulation precision. Our simulations provide insights for designing optogenetic stimulation protocols and hardware to achieve selective activation of functional domains.
Introduction
Optogenetic stimulation offers neuroscientists a unique tool for manipulating neuronal activity with precision, providing cell type specificity and high spatial and temporal resolution (Deisseroth, 2015; Emiliani et al., 2022). Until recently, long-term optogenetic wide-field manipulations in higher mammals were limited but now become feasible through implantable stimulator arrays using single-photon–based stimulation (Rajalingham et al., 2021; Hee Lee et al., 2022). This advancement opens a path toward a sophisticated interrogation of sensory or motor code in meaningful volumes of the cortex (Huang et al., 2014; Nassi et al., 2015; Chernov et al., 2018; Rajalingham et al., 2021; Wang et al., 2022), as well as the development of neuroprosthetic systems for vision and hearing loss remediation (Chernov et al., 2018; Kleinlogel et al., 2020; Rajalingham et al., 2021; Sahel et al., 2021). However, these applications require spatially precise stimulation to selectively engage functional cortical representations at a spatial scale of cortical columns (Roe et al., 2020; Seidemann, 2023) requiring submillimeter precision (Ohki et al., 2005), which in experiments proved challenging due to the spatial spread of stimulation-evoked responses (Huang et al., 2014; Chernov et al., 2018; Wang et al., 2022). The observed spatial spreading and complex distribution of responses to stimulation may be attributed to the propagation and transformation of the stimulation-evoked responses in the network (Li et al., 2019; Luboeinski and Tchumatchenko, 2020; Mahrach et al., 2020; Antolik et al., 2021). Alternatively, studies on neuron excitability suggest that optogenetic activation can recruit neurons via distal parts of their morphology (Foutz et al., 2012; Schoeters et al., 2023), raising the question whether direct activation of distant neurons via their spatially extending morphology additionally drives spread-out and spatially complex responses to optogenetic stimulation?
Here, we approached this question using computational methods to reveal how optogenetic responses of a layer-2/3 and a layer-5 pyramidal neuron depend on the lateral displacement of an optical fiber placed on the cortical surface. We validated the model's basic optogenetic response dynamics against data from primate retinal ganglion cells (RGCs) recorded across a broad range of stimulation intensities, inducing threshold, sustained, and depolarization block responses.
Our simulations demonstrate that since distant parts of neuron morphology contribute to activation, optogenetic stimulation directly drives neurons distributed in a complex manner in the cortical network beneath the stimulator. We further find that nonlinearities in the optogenetic mechanism and neuronal response function cause a stimulation-intensity–dependent redistribution of the activated neurons. Furthermore, we find that due to neuronal morphology, different spatial patterns of cells will be directly optogenetically activated in layer-2/3 than in layer-5. These findings reveal that neuronal morphology causes a diverse stimulation intensity- and layer-dependent spatial distribution of optogenetic responses that can explain the experimentally observed variation of stimulation outcomes.
The simulations predict that surface optogenetic stimulation significantly activates isolated layer-2/3 neurons within 100 µm and layer-5 neurons within 200 µm lateral distance to the stimulator, which suggests that surface stimulation may enable selectively recruiting functional neuronal representations encoded at a submillimeter scale of cortical columns. Yet, higher precision may be necessary, since evoked activity could spread beyond the directly activated cortical volume through synaptic transmission in the network. We therefore investigated whether genetics-based or optical improvements to the optogenetic setup could further enhance precision. We simulated soma-targeted channelrhodopsin expression derived from fluorescence measurements in cortical brain slices, finding that spatial precision was only improved for the layer-5 neuron and high channelrhodopsin expression. On the other hand, increasing the collimation of the optical setup improved spatial precision irrespective of the neuron type and maximized spatial stimulation precision for the tightest light beam explored regardless of whether channelrhodopsin targeting to cell bodies was involved. These results imply that optical improvements of the stimulator are more effective in increasing the precision of surface-based stimulation than genetics-based targeting of the opsin to the cell body.
Materials and Methods
Animals
All experiments were conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. The protocol was approved by the local animal ethics committees and was conducted in accordance with Directive 2010/63/EU of the European Parliament. Cynomolgus macaques (Macaca fascicularis) of foreign origin and unknown sex were used for the study.
Experimental design
Experimental data from cell-attached recordings of primate retinal ganglion cells have been partially published previously with a detailed description of the methodology (Chaffiol et al., 2017). Electrophysiological data are available from the corresponding author upon reasonable request.
Adeno-associated virus production and eye injection
Adeno-associated virus (AAV) backbone plasmid contained a human codon–optimized calcium-translocating channelrhodopsin (CatCh) sequence. Gamma-synuclein promoter was cloned into the AAV backbone plasmid. The constructs further contained woodchuck hepatitis virus post-transcriptional regulatory element and bovine growth hormone poly (A).
Animals were anesthetized with ketamine/xylazine (10:1 mg/kg). The vitreous was injected with 100 µl of viral vector holding 1012 viral particles. After injection, the cornea was treated with an ophthalmic steroid and antibiotic ointment.
Electrophysiology
AAV-treated macaque retinas were recorded in oxygenized (95% O2, 5% CO2) Ames’ medium (Sigma-Aldrich), to which a selective group III metabotropic glutamate receptor agonist, ʟ-(+)-2-amino-4-phosphonobutyric acid (ʟ-AP4, 50 mM, Tocris Bioscience) was added.
Cell-attached recordings were conducted with an Axon MultiClamp 700B amplifier, using borosilicate glass capillaries (BF100-50-10, Sutter Instrument), pulled to 7–9 MΩ. Retinal ganglion cells were recorded in current-clamp configuration (current zero) with electrodes filled with Ames’ solution. Retinal ganglion cells were dark-adapted for 1 h prior to the recording.
Full-field photostimulation was performed with a polychrome V monochromator (Olympus), and output light intensities and wavelengths were calibrated.
Simulations and computational model
Simulation software
We simulated neuronal dynamics with the NEURON simulator, version 8.0.0 (Carnevale and Hines, 2006), through its Python interface (Hines et al., 2009) using Python 3.6. We used the Snakemake workflow manager (Mölder et al., 2021) to organize the simulation workflow, parallelize, and manage runs with various parameter sets.
Neuron models
We simulated neurons with active dendrites using existing multicompartment models based on reconstructions of a layer-5 and layer-2/3 pyramidal cell (Mainen and Sejnowski, 1996) and electrophysiological dynamics previously identified based on experimental recordings for the layer-5 cell (Hu et al., 2009). Somatodendritic electrophysiological dynamics included voltage-gated sodium (na) and potassium channels (kv), slow noninactivating (km) and calcium-dependent potassium channels (kca), and high-voltage activated calcium channels (ca). Axonal electrophysiology included low- and high-threshold sodium channels (na12/na16) as well as voltage-gated potassium channels (kv). Channel conductance throughout the neuron's morphology and reversal potentials are listed in Table 1. Channels are implemented in the corresponding NMODL files (channel abbreviation + suffix “.mod”). To model the conductance through channelrhodopsin-2 channels, we used an NMODL implementation of the channel by Foutz et al. (2012); see below, Channelrhodopsin model, for details. As a default, we assumed a uniform density of 130 channels/µm2 throughout the neuron morphology (Nagel et al., 1995; Foutz et al., 2012). This yielded 10,354,945 ChR2 channels for layer-5 cells and 2,588,736 for layer-2/3 cells (a fourth of the channels in the layer-5 cells). To investigate the effects of preferential expression of ChR2 at the soma (“soma targeting”), we set the somatic density to the otherwise uniformly applied density and decreased ChR density along the neuronal processes according to distributions derived from experimental data (see below, Derivation of spatial channelrhodopsin expression distributions). In addition to our default expression density of 130 channels/µm2, we explored a 10-fold stronger expression (1,300 channels/µm2) along with the simulations comparing somatic targeting and no targeting conditions (Figs. 7, 8).
Ion channel distribution and properties
Channelrhodopsin model
The channelrhodopsin dynamics is based on a four-state Markov model describing the changes in the conformation of the molecule (Nagel et al., 2003; Hegemann et al., 2005; Nikolic et al., 2009) and has been previously implemented in NMODL (NEURON; Foutz et al., 2012). Model parameters are summarized in another study (Foutz et al., 2012) and in the original model publication (Nikolic et al., 2009). The model accounts for two open, i.e., ion-conducting, and two closed states. Under light illumination, channelrhodopsin can transfer from the closed-1 to the open-1 state. From there, it can transfer back to closed-1 or transfer to another open state, open-2, in which its conductance is lower than that in open-1. From the open-2 state, the molecule can transfer back to open-1 or to closed-2. The closed-2 state can thermally decay into closed-1 or under light illumination photoexcite back into open-2.
Optical fiber light model
We simulated the propagation of light emitted by an optical fiber inside the cortical tissue accounting for absorption and scattering of photons [absorption coefficient, 0.125 mm−1; scattering coefficient, 7.37 mm−1 (Gradinaru et al., 2009; Foutz et al., 2012); index of refraction, n = 1.36 (Vo-Dinh, 2014)] according to the Kubelka–Munk general theory of light propagation (Kubelka and Munk, 1931; Aravanis et al., 2007; Foutz et al., 2012; Vo-Dinh, 2014). The model enables the simulation of various fiber models by adapting the parameters diameter and numerical aperture. The diameter defines the size of the output surface, and the numerical aperture defines the divergence of the emitted light beam. In our simulations, we assumed that the optical fiber is placed on top of the cortical surface with its output surface directly connected to the cortical medium. We did not account for light losses involved in the coupling of the light into the cortical tissue. We used a fiber with a 200 µm diameter and numerical aperture of 0.22 as default and explored variations of the diameter of 25–400 µm and of the numerical aperture of 0.1–0.5 specifically (Fig. 7).
Derivation of spatial channelrhodopsin expression distributions
We derived the spatial distribution of ChR inside the neuronal morphology from experimental data obtained by Shemesh et al. (2017). In this study, mouse cortical neurons were transfected with CoChR (Chloromonas oogama channelrhodopsin) with and without using somatic targeting motifs to restrict the channelrhodopsin expression to the soma. They attached green fluorescent protein (GFP) to CoChR to derive the relative surface density of CoChR through fluorescence measurements, which were taken in 10 µm steps along neurites and normalized to the fluorescence level at the soma.
In the condition without somatic targeting, the fluorescence measurements along the neurites had a slightly decreasing trend with distance from the soma. Considering their uncertainty, they were compatible with a uniform expression throughout the morphology, which we chose as distribution to model the “no targeting” condition (Fig. 6A, left). A uniform distribution in the “no-targeting” case is also supported by similar fluorescence measurements for a ChR2-YFP compound in cultured neurons from rats (Grossman et al., 2010).
Somatic targeting of CoChR resulted in an exponential drop of the relative fluorescence with distance to the soma. However, it was not clear whether the exponential decay plateaued or further decayed beyond their most distant data point at 100 µm distance from the soma. To account for these two alternatives, we fitted two distribution models to their soma targeting data (Fig. 6A, right): Distribution model 1 represented a relative somatic targeting of ChR. The density decayed exponentially with the distance from soma until it plateaued at ∼50 µm distance to soma at a level of ∼20% of the somatic density. Since this model represented only an overexpression of ChR at the soma compared with distant dendrites, we referred to it as “soft soma targeting.” Distribution model 2 accounted for a complete decay of the ChR density. It described a similar exponential decay as distribution model 1 but had an additional linear decay which yielded zero ChR density at ∼200 µm distance to the soma. We referred to this case as “strict soma targeting.”
Analysis of simulated optogenetic responses
To obtain spatial profiles of optogenetic responses exhibiting their dependence on the stimulator placement on top of the cortical surface, we simulated optogenetic stimulation with a constant intensity of 200 ms duration. We recorded the membrane voltage at the neuron's soma, interpolated the data at 0.1 ms, and thresholded it at 0 mV to count spikes. We converted the spike count to the firing rate by dividing it by the stimulation duration. We repeated this procedure at stimulator locations starting at the central location above the neuron's soma and displacing the stimulator in lateral steps of 25 µm up to 975 µm. We additionally varied the stimulator location along an angular direction in 16 steps along the full circle.
We decided to compare various conditions requiring the same target response level. We determined the response level as the maximum of the optogenetic response across the stimulator's radial distance coordinate and averaged over its angular coordinate values. Since the response level is dependent on the stimulation intensity, finding a specifically targeted response level requires a search algorithm. To find this target response level, we first mapped the optogenetic response for the layer-2/3 and layer-5 neurons on a logarithmic and coarse scale to determine the stimulation intensity thresholds for response and depolarization block. Using these thresholds as limits, we used a bisection algorithm to find the stimulation intensity resulting in the desired target response level. The algorithm used as a first test value the center between the limits on a logarithmic scale and repeated this procedure while updating the upper/lower limit with the previously tested value depending on whether the resulting response level exceeded or fell below the targeted response level. During this procedure, the target response level was determined from the stimulator-distance–dependent optogenetic response profile only up to a radial distance of 200 µm to accelerate the search. We terminated the optimization procedure when the simulated response level was within 10% tolerance of the targeted one or terminated and discarded the data point if the target response level was not found within 31 steps.
Finally, we measured the spatial precision of the stimulation by determining the radial stimulator distance at which the optogenetic response averaged over the stimulator's angular coordinate reached half of the response level. We used linear interpolation if the stimulator distance yielding half of the response level was not covered by sampling.
Code accessibility
All original simulation code has been deposited at https://github.com/CSNG-MFF/osmorph/ and is publicly available as of the date of publication. Data and any additional information required to reanalyze the data reported in this paper are available from the corresponding authors upon reasonable request.
Results
We simulated optogenetic neuronal responses of pyramidal neurons from cat primary visual cortex layers 2/3 and 5, expressing channelrhodopsin-2 (ChR2). We disregarded any synaptic inputs coming from the remaining cortical network, hence focusing on studying the direct influence of optogenetic input on individual cortical neurons. Neuronal electrophysiology was simulated with the NEURON simulator (Carnevale and Hines, 2006) based on previously developed multicompartment neuron models (Mainen and Sejnowski, 1996; Hu et al., 2009; Foutz et al., 2012; Fig. 1A,B). Our model further uses an optical fiber light model to calculate light intensities along the neuronal membrane accounting for attenuation through scattering and absorption (Aravanis et al., 2007; Fig. 1A) as well as a model of channelrhodopsin light activation (Nikolic et al., 2009) to predict optogenetically induced membrane conductance (Fig. 1C). We approximated the distribution of channelrhodopsin throughout the cell morphology as uniform in agreement with experimental data (Grossman et al., 2010; Shemesh et al., 2017). Key variables of the simulation workflow are illustrated in Figure 1D. Details of the model are described in Materials and Methods.
Optogenetic response model. A, Simulated light intensity distribution of an optical fiber (Aravanis et al., 2007; Foutz et al., 2012) and reconstructions of the layer-5 and layer-2/3 pyramidal neuron morphology (Mainen and Sejnowski, 1996; Foutz et al., 2012) in proportion to the cortical volume with layers indicated on the left. B, Hodgkin–Huxley-type neuron model with membrane mechanisms used in this study (Hu et al., 2009; Foutz et al., 2012). C, Four-state Markov model used to simulate the activation of channelrhodopsin-2 (ChR2; Nikolic et al., 2007; Foutz et al., 2012). D, Behavior of key variables during an example stimulation of the layer-5 neuron: the light propagates from the source through the cortex, reaching the cell membrane, where it results in conductance of ions through ChR2 channels, which causes a neuronal response recorded at the neuron's soma. Plots of light power and ChR2 conductance show sums across the cell membrane. Stimulation intensity was 9.5 mW/mm2 measured at the stimulator output surface. E, Agreement of intensity dependence of model optogenetic responses with biological data. The first and third columns (from the left) show cell-attached recordings of two primate retinal ganglion cells (RGCs) undergoing optogenetic stimulation at stepwise increasing stimulation intensity starting from 4 × 1014 photons/cm2/s (corresponding to a relative stimulation intensity of 1). Model optogenetic responses of the layer-5 pyramidal neuron (second and fourth column from the left) matched biological response dynamics for a ChR-2 expression of 130 channels/µm2 for RGC Cell 1 and 195 channels/µm2 for RGC Cell 2. To account for differences between the experimental and simulated computational setups, absolute intensities were increased by a factor of 14.75 for the computationally modeled responses.
Model optogenetic responses match cell-attached recordings of optogenetic responses in primate RGCs
We verified that our model dynamics of optogenetic responses are biologically plausible by comparing the optogenetic response of the simulated layer-5 neurons to cell-attached recordings of transduced primate retinal ganglion cells (RGCs; Fig. 1E). Optogenetic responses of RGCs were recorded for absolute stimulation intensities ranging between 4 × 1014 photons/cm2/s to 3 × 1017 photons/cm2/s measured at the neuron's soma. To account for variations in opsin expression, efficiency, and light attenuation, we compared stimulation intensity levels on a normalized scale. Biological responses were recorded at the multiplications of 1, 50, 200, 250, and 750 of the lowest stimulation intensity. This intensity range induced optogenetic responses covering threshold response at lowest, sustained response at 50-fold, and depolarization block at 200–750-fold normalized stimulation intensity. This behavior was observed in approximately half of the cells recorded from one retina. The remaining cells had varied responses but did not exhibit depolarization block within the recorded range of stimulation intensities, which may have been caused by lower opsin expression, and therefore insufficient stimulation intensity to reach depolarization block. Model responses simulated at an expression density of 130 channels/µm2 and the same (normalized) stimulation intensity steps captured response onset, spike count, and amplitude modulation of the biological responses of cells exhibiting depolarization block. For example, RGC 1 responded with a single spike at ∼30 ms latency at the lowest and a sustained response of 35 spikes at 50-fold stimulation intensity (Fig. 1E, left). At 250- and 750-fold normalized intensity, the cell exhibited few spikes of decaying amplitude within 10–20 ms after stimulation onset. Model responses qualitatively matched the latency of the single spike response (∼20 ms) at the lowest intensity, the spike count at the sustained response (33 spikes) appearing at 50-fold intensity, and a complete decay of spike amplitude within 10–20 ms in the depolarization block response at 250- and 750-fold intensity (Fig. 1E, left). Model responses also agreed with optogenetic responses recorded in another RGC (RGC 2) and further captured the modulation of spike amplitudes occurring during sustained response at 50-fold normalized stimulation intensity (Fig. 1E, right). In summary, our model exhibits intensity-dependent changes in response dynamics, including depolarization block, which we validated based on matching observations in a major number of neurons recorded for a range of stimulation intensities from a single macaque retina.
Morphology causes heterogenous spatial distribution of optogenetic activation beneath the light source
Keeping stimulation intensity constant, we placed the stimulator at various positions along the cortical surface (Fig. 2A) and simulated optogenetic responses of the layer-5 neuron to a 200 ms duration constant light pulse (Fig. 2B, top). We counted spikes in somatic membrane voltage through thresholding at 0 mV and quantified the firing rate by dividing it by stimulation duration. By mapping the optogenetic responses with respect to the position of the stimulator relative to the cell body projected onto the cortical surface (Fig. 2B, bottom), we found that optogenetic responses varied in a heterogenous way with the stimulator position. At a stimulation intensity of 16 mW/mm2 at the stimulator surface, the neuron's response was weaker for a central stimulator location above the soma compared with a 200 µm lateral displacement. While exhibiting a decaying trend toward more lateral positions from the center, isolated stimulator locations at >300 µm distance from the neuron's soma caused responses of similar strength compared with central positions. These results highlight that spatially extending neuronal morphology of layer-5 neurons could contribute to ∼500 µm spread of optogenetic activity along the cortical surface, considering the average full-width at half-maximum of the optogenetic response. Further, individual and spatially isolated responses may introduce noise into spatial stimulation patterns at a larger scale of up to 1 mm.
Heterogeneous distribution of stimulator locations induces strong optogenetic responses. A, Illustration of the overlap of light intensity throughout cortical space with neuron morphology for three optical fiber locations increasing in radial offset from center. B, Optogenetic responses to a 200 ms duration stimulation at the example optical fiber locations shown in (A). The membrane voltage is measured at soma and thresholded at 0 mV to count elicited spikes. Spike counts are encoded with color and mapped spatially at the respective stimulator position to reveal the impact of neuron morphology (blue) on the dependence of optogenetic response on stimulator location. Repetition of the optogenetic response simulation at various optical fiber locations reveals the full spatial response profile. Stimulation intensity was 16 mW/mm2 at the stimulator surface. Stimulator was an optical fiber with a diameter of 200 µm and numerical aperture of 0.22.
Increasing stimulation intensity moves the strongest ChR conductance deeper along the neuron's morphology
Next, we investigated the impact of varying stimulation intensity on the optogenetic responses of the layer-5 neuron. We varied the stimulator position along the cortical surface while stimulating at low, intermediate, and high stimulation intensities (Fig. 3A). Recording the ChR conductance across all compartments, we found that the origin of strongest ChR conductance shifted with increasing stimulation intensity from the apical tuft to the apical shaft (Fig. 3B). Contributions of single basal dendritic, somatic, or axonal compartments to total ChR conductance were negligible. The reasons for the shift were twofold. First, ChR channels in the apical tuft saturated with increasing light intensity at stimulator output earlier than in the apical shaft as apical tuft compartments were exposed to comparably higher light intensity being located closer to the stimulator. Second, a higher surface area in the apical shaft (Fig. 3C) accommodated a higher number of channels, which could in total sum up to a larger ChR conductance per compartment although being exposed to less light (Fig. 3D). The shift of strongest ChR conductance (per compartment) from apical tuft at low to the apical shaft at higher stimulation intensity resulted in a change of the spatial distribution of the neuron's spiking response (Fig. 3E). At low stimulation intensity (1.6 mW/mm2), spiking responses appeared only if the stimulator was located where the apical tuft reached close to the cortical surface, i.e., in isolated and off-center stimulator locations (Fig. 3E, left). At higher light intensities (16 mW/mm2 and more) for which the apical shaft contributed the strongest ChR conductance per compartment (Fig. 3B middle), spiking responses were highest for central stimulator positions, for which the apical shaft received the highest light exposure (Fig. 3E, middle). These results suggest that for spatial precision, stimulation of layer-5 pyramidal neurons with a wide branching apical tuft is unfavorable at low stimulation intensities, as neurons located laterally to the light source might receive stronger optogenetic excitation compared with neurons located right below. Increasing the stimulation intensity could mitigate this effect, as it recruits more ChR conductance in the apical shaft, which resulted in an optogenetic response spatially centered around the neuron's soma in our simulations.
Optogenetic responses redistribute with increasing stimulation intensity. A, Apical tuft (purple), apical shaft (orange), and remaining neuron morphology (blue) of the layer-5 pyramidal neuron and stimulator light profile inside the cortex viewed from the side (top) and neuron morphology viewed from the top (bottom). B, Dependence of maximum channelrhodopsin conductance in a single compartment on stimulator location for low (1.6 mW/mm2), intermediate (32 mW/mm2), and high stimulation intensity (160 mW/mm2). Color encodes compartment type where maximum conductance is reached. Basal dendritic, somatic, and axonal compartments had negligible contributions and therefore were not displayed for visualization purposes. C, Area per compartment tends to be larger in the apical shaft compared with the apical tuft. D, Total conductance through ChR increases more strongly in the apical shaft than in the apical tuft and therefore dominates the response at higher stimulation intensities. E, Optogenetic spiking response depending on stimulator location at low, intermediate, and high stimulation intensity. F, Strongest ChR conductance observed across compartments for varying stimulator locations (central location marked by a red circle, otherwise not encoded) and stimulation intensity (color-coded). The depolarization block is induced only for the highest stimulation intensity condition and is strongest in the central stimulator location.
Depolarization block preferentially occurs at central stimulator locations causing a ring-shaped response
At high stimulation intensity, depolarization block additionally reshaped the spatial characteristics of the optogenetic response (Fig. 3B,E,F). Central stimulator locations, which recruited the strongest ChR conductance per compartment, induced depolarization block with increasing stimulation intensity before other locations. As recruited maximal ChR conductance decreased with the stimulator offset from the center, depolarization block was induced inside a ring around the neuron's soma, which increased in radius with stimulation intensity. These observations indicate that neurons below the optogenetic stimulator can be suppressed due to depolarization block while neurons located laterally to the stimulator still respond with strong sustained responses at high stimulation intensities.
Stimulation drives neurons in layer 5 and layer 2/3 differently
Since our simulations indicated that morphology causes complex optogenetic responses, we were interested in whether this implies distinct layer-specific optogenetic activation patterns, as a neuron's morphological type is linked to its layer origin. We therefore simulated the stimulator location–dependent optogenetic responses of a layer-2/3 pyramidal neuron to compare them to the layer-5 neuron's responses (Fig. 4A). In contrast to the layer-5 neuron, optogenetic response maxima were evoked for centered stimulator locations also at low and not only at intermediate stimulation intensities. The depolarization block was, as in the layer-5 neuron, first induced at centered stimulator locations when stimulation intensity was increased. However, the thresholds for response and depolarization block were at different stimulation intensities compared with the layer-5 neuron. Consequently, the spatial response characteristics of the layer-2/3 and layer-5 neurons differed when stimulated at the same stimulation intensity (Fig. 4A). We found inverse response dependencies on radial stimulator offset at low stimulation intensity (3.2 mW/mm2) as response maxima were evoked by off-center stimulator locations in the layer-5 while evoked by central stimulator locations in the layer-2/3 neuron. At increased stimulation intensity (9.6 mW/mm2), response maxima were evoked at central stimulator positions for both neuron types. A further increase (32 mW/mm2) induced depolarization block at central stimulator positions only in the layer-2/3 neuron, while the layer-5 neuron retained a sustained response. At high stimulation intensity (160 mW/mm2), depolarization block occurred in both neurons for central stimulator locations, however, affecting the layer-2/3 neuron within a wider area of stimulator locations. These observations indicate that differences in morphology due to layer origin of a neuron give rise to a spatially different distribution of optogenetic excitation in populations of (isolated) layer-5 and layer-2/3 neurons, which we visualized accounting for different neuronal densities of these neurons in cat visual cortex (Binzegger et al., 2004) in Figure 4B. In summary, our results demonstrate that stimulation intensity not only controls response strength but also modulates spatial response distribution (Fig. 4C).
Spatial distribution of optogenetic excitation is dependent on the layer origin of the neuron. A, Optogenetic response dependence on stimulator location for layer-5 (left) and layer-2/3 pyramidal neuron (right) at varying stimulation intensity (rows). B, Activated layer-5 (left) and layer-2/3 pyramidal cells inside 1 mm2 of the cortex around a stimulator (optical fiber) located in the center (black circle) and accounting for neuron-type–dependent density (Binzegger et al., 2004). C, Qualitative illustrations of the different spatial distributions of optogenetic excitation in layer-5 (brown) and layer-2/3 (pink) neurons beneath the stimulator and depending on stimulation intensity (rows).
Quantification of neuron-type–specific spatial stimulation precision
Since optogenetic responses varied in their spatial extent with varying stimulation intensity, and therefore response strength, we quantified stimulation precision in layer-5 versus layer-2/3 neurons dependent on their activation level (Fig. 5). We varied stimulation intensity starting from below the response threshold all the way to a level that exceeded the threshold for depolarization block. For each simulated stimulation intensity step, we first averaged the optogenetically induced responses across the angular stimulator coordinate and then found the maximum response in the resulting distance-dependent profile, which we termed “peak response” and used as an indicator for the neuron's response strength. Next, we calculated the stimulator distance at which the neuron's response dropped to half of the peak response, which we defined as “response space constant” and utilized as a measure of stimulation precision (Fig. 5A).
Comparison of spatial stimulation precision of layer-5 and layer-2/3 neurons at equal target response levels. A, Examples of optogenetic response profiles and response space constant at a target response level of 50 Hz (top). Determination of “peak response” (maximum) and “response space constant” (radial fiber offset at half-maximum) from the stimulator-offset–dependent optogenetic response. B, Peak response of layer-2/3 (red) and layer-5 neuron (blue) plotted over stimulation intensity at the stimulator output (top). Response space constant for both neurons plotted over stimulation intensity (bottom). C, Response space constant plotted against targeted peak response, same color-coding as in (B). The response space constant is not directly dependent on peak response; instead, both variables are dependent on the underlying variable stimulation intensity. Diverging of the response space constant at a certain peak response level is occurring due to depolarization block at high underlying stimulation intensities.
The peak response increased monotonically with stimulation intensity in both neuron types (Fig. 5B, top). However, layer-5 and layer-2/3 neurons had different response thresholds, gains, and saturation points, which meant that at lower stimulation intensity, the layer-2/3 cell responded with a higher firing rate than the layer-5 cell and vice versa at higher intensities. Crucially, the response space constant increased with stimulation intensity in both neuron types (Fig. 5B, bottom). Hence, higher stimulation intensity resulted in a stronger response but also decreased spatial precision. Finally, we plotted the response space constants over the peak response, to provide a comparison of stimulation precision at equal response levels for the two examined neuronal types (Fig. 5C). We found that throughout the response range below depolarization block, layer-2/3 neurons can only be engaged above half-maximum within ∼100 μm distance between stimulator and soma while layer-5 neurons can be engaged within ∼200 μm to this strength.
Dendritic excitability drives responses at lateral stimulator positions
Dendritic excitability may differ between neurons raising the question of how its variance affects optogenetic responses. Lowering the excitability of the layer-5 model neuron by reducing dendritic sodium conductance suppressed the responses at lateral stimulator locations, but the spatial heterogeneity of the responses prevailed (Fig. 6A, “low Na cond.”). Likewise, increasing dendritic sodium conductance strengthened optogenetic responses at lateral stimulator locations (Fig. 6A, “high Na cond.”). Manipulating dendritic excitability via increasing or decreasing leak conductance leads to qualitatively similar but weaker effects (Fig. 6A, “low leak cond.,” “high leak cond.”). Higher leaks resulted in marginally suppressed, and lower leaks resulted in slightly increased lateral responses. Overall, spatial optogenetic responses behaved qualitatively similar across varying dendritic excitability conditions except when dendritic sodium conductance was decreased (Fig. 6B–D): Instead of a monotonous increase of the response radius with stimulation intensity as in the other conditions, the response radius followed a U-shaped curve with a dip at intermediate stimulation intensities, implying a gain in precision for this stimulation intensity range (Fig. 6C,D, red dash-dotted curve). This precision gain is caused by an increase in response strength at central stimulator locations between low and intermediate light intensities (Fig. 6A, “low Na cond.,” 3 vs 27 mW/mm2), which did not occur in the other conditions.
Impact of dendritic sodium and leak conductance on optogenetic response. A, Optogenetic responses of layer-5 pyramidal neuron depending on stimulator location for varying dendritic excitability (rows) and different stimulation intensities (columns). Rows represent (1) default, (2) 50% decreased sodium (Na) conductance, (3) 50% increased sodium (Na) conductance, (4) 50% decreased leak conductance, and (5) 50% increased leak conductance. B, Peak response and, C, response space constant of layer-5 neuron response plotted over stimulation intensity at the stimulator output depending on excitability conditions from (A). D, Response space constant plotted against targeted peak response, same color-coding as in (B) and (C).
Restricting ChR expression to the soma improves spatial stimulation precision only in the layer-5 neuron
An intuitive measure to prevent the activation of neurons from distant stimulator positions due to the illumination of their dendritic arbors is restricting the expression of optogenetic constructs to the cell body. Shemesh et al. (2017) have experimentally demonstrated that an adapted version of the CoChR protein, which selectively traffics to the cell body of neurons, improves the spatial stimulation precision of holographic optogenetic stimulation in cell cultures and cortical slices. To understand if such soma targeting of the opsin is likely to improve the precision when the stimulation is applied with a divergent light source from the cortical surface, we replicated their measured soma-targeted opsin expression distribution in addition to the default uniform distribution representing no targeting in our simulations (Fig. 7A,B; see Materials and Methods). We simulated two alternative somatic targeting conditions as their data were not clear on the presence of residual ChR transfection in the distal dendrites. First, we fitted an exponential decay plus y-offset plateauing at ∼20% of the density expressed at the soma (“soft soma targeting”), and second, an exponential and linear decay plus y-offset which decays to zero density at ∼200 µm distance from the soma (“strict soma targeting”). In addition, we explored the impact of different absolute ChR expression levels in addition to different relative ChR distributions. The default ChR expression density used in our simulations was 130 channels/µm2, which matched an empirical estimate of bacteriorhodopsin density in Xenopus oocytes (Nagel et al., 1995). We additionally tested a 10-fold higher expression density (1,300 channels/µm2).
Preferential somatic ChR expression improves spatial precision only in the layer-5 pyramidal neuron. A, Normalized channelrhodopsin expression distributions (solid lines) fitted to GFP fluorescence measurements (discrete data, mean, and STD; Shemesh et al., 2017). B, Expression distributions from (A) applied to the layer-5 pyramidal neuron. C, Comparison of response space constants observed at 50 Hz target peak response for the layer-5 and layer-2/3 pyramidal neuron type at channelrhodopsin expression levels of 130 channels/µm2 and 1,300 channels/µm2. D, Detailed results of the effect of somatic targeting for the layer-5 neuron and an expression level of 1,300 channels/µm2 for which strict somatic targeting achieved the strongest improvements in spatial precision. The 3D plot (top left subpanel) illustrates the dependence of peak response and response space constant on stimulation intensity. The bottom projection shows the dependence of the peak response on stimulation intensity, which is also plotted in 2D in the top right subpanel. The stars in the top right subpanel mark the data points at the 50 Hz target peak response, from which the bar plot in (C) was generated. Spatial response profiles corresponding to these data points are plotted below (bottom right subpanel). The bottom left subpanel shows the behavior of the response space constant with peak response, which is also plotted as a projection in the 3D plot (top left subpanel).
We compared spatial precision between conditions by comparing the response space constant at the same target peak response level. Somatic targeting improved spatial precision only for the layer-5 neuron and required strict somatic targeting of ChR, as shown for an example target response level of 50 Hz in Figure 7C. Higher ChR expression (1,300 channels/µm2) further enhanced the precision. As notable precision improvement was restricted to the layer-5 neuron with strict somatic targeting of ChR at high expression density (1,300 channels/µm2), we included detailed simulation results only for this condition in Figure 7D. Across conditions (neurons, targeting, ChR expression level), somatic targeting required an increase of stimulation intensity, which amounted to two orders of magnitude in the layer-5 neuron and one order of magnitude in the layer-2/3 neuron (Fig. 7D). Symmetry of the stimulator-location–dependent response profiles was improved with strict somatic targeting in both neurons (Fig. 7D). In conclusion, our results suggest that somatic targeting of ChR can improve spatial precision only in the layer-5 neuron, but symmetry in both neurons, at the cost of a strong increase of light irradiance. Its effect further requires zero residual ChR expression in distal dendrites.
To understand whether the reported effects of ChR targeting generalized across optical fiber parameterizations, we repeated the simulations for the “no targeting” and “strict targeting” conditions and evaluated the response space constant at a peak response of 50 Hz. Indeed, the substantial improvement of spatial precision through strict targeting of ChR in the layer-5 neuron and the lack of improvement in the layer-2/3 neuron were present for most simulated optical fiber parameterizations [fiber diameters d = (25, 400) µm, numerical apertures NA = (0.1, 0.5); Fig. 8].
Narrowing optical fiber diameter and numerical aperture improves stimulation precision. A, Light intensity levels throughout cortical space with layer-5 and layer-2/3 pyramidal neuron superimposed for several example optical fiber parameterizations. B, Dependence of the response space constant at a target peak response of 50 Hz on optical fiber diameter and numerical aperture for the layer-5 and layer-2/3 neurons and using no (red) and strict somatic targeting of ChR2 (blue).
Narrowing the stimulator light profile improves precision independently of neuron type
Narrowing the stimulator's light profile decreased the spatial extent of the optogenetic response in both neuron types (Fig. 8). The relative improvement of their response space constants between the widest (d = 400 µm, NA = 0.5) and narrowest (d = 25 µm, NA = 0.1) optical fiber conditions was approximately fourfold for the layer-2/3 neuron and three- or fourfold for the layer-5 neuron, depending on whether somatic targeting was involved or not. With the optical fiber narrowed to the tightest light beam profile explored (d = 25 µm, NA = 0.1), stimulation precision reached its maximum regardless of whether somatic targeting of channelrhodopsin was involved or not.
Discussion
Cortical implants using optogenetic stimulation to engage functionally meaningful neuronal populations are posed to play a key role in unraveling how cortical circuits process information (Huang et al., 2014; Nassi et al., 2015; Roe et al., 2020; Wang et al., 2022; Seidemann, 2023), develop functional architecture (Mulholland et al., 2024a,b), and in clinical applications as neuroprosthetic devices for restoration of sensory function, such as vision (Kleinlogel et al., 2020; Antolik et al., 2021; Sahel et al., 2021). For the success of many applications, the stimulation must be precise at a scale of hundreds of micrometers, at which many functional neuronal representations, for example, orientation preference in the primate visual cortex (Ohki et al., 2005), are encoded. Our simulations show that optogenetic stimulation directly drives neurons across several hundred micrometers, depending on the shape, size, and cortical depth of neuronal morphology. We find direct optogenetic impact on layer-5 neurons within 250 µm and on layer-2/3 neurons within 100 µm distance from the stimulator. These predictions align with experiments in the primary visual cortex of tree shrew, which found significant direct optogenetic activation within 150 µm and none beyond 300 µm when network inputs were abolished via excitatory synaptic blockers (Huang et al., 2014). Crucially, our simulations predict that the impact of neural morphology on direct stimulation precision depends on stimulation intensity in a nonlinear fashion, as increased intensity results in both elevated or decreased response depending on secondary factors such as the absolute intensity level and neuron type, which has practical implications for optogenetic experiments. Finally, while subcellular targeting of channelrhodopsin improves precision under restricted conditions, optical changes to the stimulator enhance precision more generally in our simulations.
Our findings can help interpret past and future experiments and guide the design of future optogenetics-based implants. For example, we find that high light intensity can cause suppression by depolarization block centrally below the stimulator while surrounding areas remain activated (Fig. 3). This unwanted activation pattern can lead to misinterpretation of experiments and may be difficult to detect. Predicting the intensity triggering depolarization block is challenging due to variable transfection levels between cells and experiments. Hence, when increasing stimulation intensity to recruit more neurons, it is important to monitor neurons in the center to detect suppression by depolarization block. Alternatively, pulsed instead of continuous stimulation could be used to prevent depolarization block (Herman et al., 2014), which we observed in our experiments and simulations within 10 ms after stimulation onset, hence requiring pulse frequencies larger than 100 Hz to avoid its onset. Interestingly, depolarization block could also be used intentionally to suppress neurons beneath the stimulator while activating those around it.
The second implication concerns the spatial distribution of optogenetic excitation, which varies with stimulation intensity in pyramidal neurons of layers 2/3 and 5 (Fig. 4A). Modulating intensity not only affects response strength but also allows optimization of the distribution of excitation (Fig. 4B,C). However, achieving central activation of both neuron populations beneath the stimulator requires restricting the intensity to a narrow range, limiting flexibility in controlling response levels.
Third, we found that increasing stimulation intensity not only raises the response level but also extends the distance at which a neuron is stimulated. Similarly, optogenetic excitation of inhibitory neurons in mice showed a concomitant increase in the inhibitory neurons’ response strength and spatial response extent with higher intensity (Li et al., 2019). Since the synaptic activity was not blocked and the observed response extent (>1 mm) exceeded our findings for isolated neurons, it suggests that network mechanisms may amplify this effect. Therefore, while increasing stimulation intensity amplifies neuronal responses, it also broadens their spatial reach. Vice versa, lowering intensity could lead to more spatially specific responses. This poses important constraints for designing stimulation devices and protocols, which must balance activation strength with precision, for which our simulations could provide guidance.
Fourth, we observed that somatic targeting of ChR could more than double the stimulation precision in layer-5 but had minimal effect in layer-2/3 neurons. This improvement was only seen with high ChR expression (1,300 channels/µm2) and no residual expression in the apical tuft, highlighting that precision enhancement requires a potent, well-targeted optogenetic construct. However, somatic targeting required up to two orders of magnitude higher light intensities to reach the same response (Fig. 7B,D), raising concerns about phototoxicity in superficial layers. Similarly, narrowing the stimulator's light emission profile to improve precision also raises the risk of phototoxicity as it increases light intensity close to the output (Fig. 8). Using more light-sensitive opsins may help alleviate this risk by enabling stimulation at lower intensities (Sridharan et al., 2022).
Our study predicts direct optogenetic activation at distances greater than 200 µm from the neuron's soma. As shown in Figure 6, these contributions depend on dendritic excitability and can be suppressed or enhanced. While more excitable dendrites do not change the overall pattern of increased response levels and reduced spatial precision with higher stimulation intensity, we find that neurons with less excitable dendrites may offer better precision at intermediate compared with low intensities. Additionally, since lateral responses are suppressed in neurons with less excitable dendrites, this has a similar effect in improving precision to somatic targeting of the opsin.
Various computational studies investigated optogenetic stimulation of single neurons (Grossman et al., 2011, 2013; Foutz et al., 2012; Jarvis et al., 2018; Schoeters et al., 2023), but only a few (Foutz et al., 2012; Schoeters et al., 2023) examined how the stimulator's position in relation to the neuron affects the optogenetic response. These studies focused on determining the stimulation intensity required to induce spiking, showing that dendrites contribute to activation and that altering ChR distribution across the cell can lower the stimulation threshold. Our results complement these studies by examining spatial optogenetic responses beyond threshold, demonstrating that response characteristics change with varied stimulation intensity, which is important for the design of spatially coordinated optogenetic interventions.
There are several limitations to this study that warrant further research. First, this study focuses on excitatory neurons, providing only limited insight into the response characteristics of inhibitory neurons. Our simulations revealed that proximity to the optogenetic stimulator increases a neuron's sensitivity to stimulation. Given that pyramidal cells which comprise most excitatory neurons typically have an apical dendrite extending toward the cortical surface, these excitatory neurons may, on average, be more sensitive to surface-delivered stimuli than inhibitory neurons located in deeper cortical layers, whose dendrites do not extend in the same way. However, less predictable factors, such as the illuminated ChR-transfected membrane area, intensity threshold for effective ChR conductance, and specific ChR transfection levels in inhibitory neurons, can significantly influence their response. Detailed computational studies are needed to explore these variables.
Second, since our work only examines single-neuron responses to optogenetic stimulation, it is limited to predictions on the direct optogenetic input to a neuronal network. Network-level effects, like those seen in previous studies showing broadening of activation from random connectivity (Luboeinski and Tchumatchenko, 2020) or functional sharpening in the orientation domain of the primary visual cortex (Antolik et al., 2021), must be explored using detailed spiking network simulations. This study paves the way toward such studies, which can integrate our single-neuron findings with network modeling to account for the impact of neural morphology and optogenetic stimulation at the network level.
Third, our model predicts depolarization block at high stimulation intensities, consistent with experimental data, but only as observed in approximately half of the ganglion cells from a single Macaque retina. The absence of depolarization block in the remaining cells may be attributed to lower opsin expression levels, which could not be experimentally verified. Computationally, we confirmed that lower opsin expression indeed requires higher light intensities to induce depolarization block. However, more experimental evidence is needed to answer whether depolarization block appears more universally at high stimulation intensities across various neuron types.
Fourth, to understand how large neuronal populations respond to stimulation, it is essential to characterize neurons with diverse morphologies and biophysical properties beyond the neurons examined here. For example, discovering different response characteristics between excitatory and inhibitory neurons could enable varied stimulation strategies, leveraging their distinct roles in network activation. Furthermore, future work could apply our analysis to computational models of a broader variety of the same neuron type, helping to quantify the effects of inter-neuron variations and potentially detect errors in our model parameterizations, which, while well-established, rely on a limited set of experimentally challenging measurements.
Finally, we omitted synaptic background activity, which could interact with the charge introduced through ChR photoactivation, potentially shifting the activation threshold or influencing response characteristics if the synaptic input is spatially biased. Accounting for the spatiotemporal properties of synaptic inputs accurately would require full network simulations, which is beyond the scope of this study. Furthermore, our study considers temporally steady light stimulation. Pulsed stimulation is often used in practice [see Grossman et al. (2011) for a comparison]; thus, studying how the parameters of such protocol impact the response characteristics is of great interest. Another limitation is that our neuron models did not include axonal arbors. Previous modeling by Foutz et al. (2012) suggests that large axon bundles, such as those in layer 2/3 of the primary visual cortex in cats (Stepanyants et al., 2008), could create additional clusters of stimulator positions that activate the neuron. However, such response clusters have not been observed in optogenetic experiments in the primary visual cortex of tree shrews (Huang et al., 2014). Future work incorporating neuronal reconstructions with axonal arbors could clarify their role in optogenetic activation.
Footnotes
We thank Thomas Foutz for his technical advice on the multicompartment neuron model and Lawrence Humphreys and Jose-David Celdran for their helpful discussions on the subcellular targeting of opsins. This work was supported by the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 861423 (Entrain Vision), the Charles University Primus Research Program 20/MED/006, and the ERDF-Project Brain Dynamics, No. CZ.02.01.01/00/22_008/0004643.
S.P. is a founder and consultant for GenSight Biologics and Pixium Vision.
- Correspondence should be addressed to Ján Antolík at jan.antolik{at}mff.cuni.cz.
