1 Introduction

Optogenetics is a technique for exciting or silencing cells within living tissue via genetic photosensitization and remote optical activation (Boyden et al. 2005; Yizhar et al. 2011). This photostimulation technology allows the interrogation of neural circuits with high spatial and temporal resolution (Wang et al. 2007; Zhang et al. 2006, 2007a). In addition, it offers substantial prospects for the development of novel neuroprosthetic interfaces (Busskamp and Roska 2011; Degenaar et al. 2009). The photosensitization of cells is achieved by transfecting cells with an opsin, such as the light-sensitive algal protein channelrhopsin-2 (ChR2) (Nagel et al. 2003; Gunaydin et al. 2010), or with a light- gated chloride pump halorhodopsin NpHR (Zhang et al. 2007b) or proton pump Arch (Chow et al. 2010). Spatial selectivity can be achieved by genetically targeted expression of opsins, or by focusing the light beam onto targeted areas. For example, the activation light can be focused on a single soma using laser-coupled optical fibers (Han et al. 2009), or onto a large number of subcellular compartments using micro-LED arrays (Grossman et al. 2010) or amplitude (Shoham et al. 2005) and phase (Lutz et al. 2008) modulation, as well as two-photon excitation (Rickgauer and Tank 2009; Andrasfalvy et al. 2010). Despite its immense importance, it is still unclear how the spatial distribution of the stimulating light affects the spiking dynamics of the cells and the backpropagating APs (bAPs) in particular.

This study investigates, via modelling and computer simulations, the fundamental effects of different spatial illumination strategies on the neural spiking response. For that purpose we used rat layer V pyramidal neuron models. Although one of the best studied types of neurons, e.g. (Stuart et al. 1997; Schaefer et al. 2003; Letzkus et al. 2006) seemingly there is still not a unique model of this cell type with fixed parameter values which can accurately reproduce all experimental results. Here we use the model proposed by Hay et al. (2011), with parameters optimised for two cases: high and low dendritic tree excitability. However, when only the axon was illuminated we decided to use the model of Hu et al. (2009) which has a better description of the complete axon but shows an over excitability of the apical dendritic tree.

The pyramidal neurons were photosensitized by incorporating ChR2, but the modelling framework we establish here is applicable to other opsins. We simulated the photoconductivity of ChR2 with a new six-state model which is very similar to the branched four-state model we previously proposed (Nikolic et al. 2009). Another innovation in the present work is implementation of the six-state ChR2 model in NEURON (Hines and Carnevale 1997). The new model was needed for better representation of the ChR2 kinetics in the NEURON simulation environment. By implementing the model in NEURON, we enable its use for studying the dynamics and integration of optogenetic perturbations of three-dimensional morphologically realistic neurons and neuronal circuits.

The photosensetized cells were optically excited with four major types of illuminations: (1) somatic (2) whole cell (3) dendritic tree (common for in-vivo experiments) and (4) axon initial segment (AIS). We explored the impact of the illumination pattern on the waveform of the evoked APs and the frequency-irradiance response. Finally, we investigated the impact on the latency and magnitude of the bAPs which could have important implications in use of optogenetics in optical induction of synaptic plasticity, such as pairing of light stimulation with postsynaptic depolarization induced long-term potentiation (Zhang and Oertner 2007).

2 Methods

2.1 ChR2 model

The ChR2 current (\(I_{\rm ChR}\)) is given by a simple expression: \(I_{\rm ChR}=A\cdot g_{\rm ChR}\cdot (v-E_{\rm rev})~,\) where \(A\) is the effective area of the cell/compartment, \(g_{\rm ChR}\) is the conductance of ChR2 ion channels per unit area, and \(v\) is the membrane potential (typical in units of mV). The reversal potential for Chlamydomonas reinhardtii ChR2 is approximately zero for pH \(=\) 7.35 and the extracellular ion concentrations of sodium, potassium and calcium near their physiological values (the actual value could be in the range \(E_{\rm rev}\approx 0-8\,\)mV) (Lin et al. 2009; Bamberg et al. 2008). The conductance of ChR2 is best reproduced by separating its light and voltage dependences:

$$ g_{\rm {ChR}} (\phi,v,t)=\overline{g}_{\rm {ChR}}\cdot \psi(\phi,t)\cdot f(v)$$
(1)

where \(\phi \) is the photon flux in [photons\(\cdot \)s\(^{-1}\cdot \)cm\(^{-2}\)], \(v\) is the membrane potential, \(\psi \) a normalized light dependent function (\(\psi =1\) when all channels are in the maximal conductive state), \(f\) is the voltage dependent function and \(\overline {g}_{\rm ChR}\) is the maximum conductance (for \(f=1\)) in [nS/\(\mu \)m\(^2\)], which accounts for both the conductance of a single ChR2 molecules and the level of protein expression in the cell.

The voltage dependence of ChR2 was derived from photocurrent recordings under constant irradiance and changing clamped voltages (Bamberg et al. 2008; Feldbauer et al. 2009). The empirical equation was derived in our previous work (Grossman et al. 2011):

$$ f(v)=\frac{1-\exp(v/v_0 )}{v_1}$$
(2)

where \(v_0\) and \(v_1\) are fitting constants. The light dependence of the ChR2 conductance was derived from photocurrent recordings at fixed voltage (voltage clamp measurements, typically at \(-\)70 mV) and changing irradiances. We previously showed that the dynamics of the photocurrents can be accurately reproduced using a four-state kinetic model which consists of four functional states (Nikolic et al. 2006, 2009), see also (Hegemann et al. 2005; Nagel et al. 2003): dark adapted ground state (\(s_1\) in Fig. 1), high conductance open state (\(s_3\)), weak conductance open state (\(s_4\)) and close (non-ground) state (\(s_6\)). Here we introduce two additional, intermediate states to correctly account for the time needed for ChR2 to activate after trans-cis (\(s_2\)) and cis-trans (\(s_5\)) retinal transformations. In our previous model (Nikolic et al. 2009), we calculated this term analytically and introduced time-dependent transition rates between \(s_1\rightarrow s_3\) and \(s_6\rightarrow s_4\). However, in order to avoid having any rates explicitly time-dependent (some must implicitly depend on time through the light flux dependence) intermediate states \(s_2\) and \(s_5\) in Fig. 1 were introduced. In this way the model is better adapted for the NEURON simulation platform as well. The transition between the states is described here by a forward (\((a)\)) and a backward (\((b)\)) rates shown in Fig. 1.

Fig. 1
figure 1

The six-state model of the ChR2 photocycle

Full size image

The model is governed by a set of linear differential equations:

$$ \dot{s_1} = -a_1s_1 + b_1s_3 + a_6s_6$$
(3)
$$ \dot{s_2} = a_1s_1 - a_2s_2$$
(4)
$$ \dot{s_3} = a_2s_2 - (b_1+a_3)s_3 + b_2s_4$$
(5)
$$ \dot{s_4} = a_3s_3 - (b_2+a_4)s_4 + b_3s_5$$
(6)
$$ \dot{s_5} = -b_3s_5 + b_4s_6$$
(7)
$$ \sum s_i = 1$$
(8)

The normalized light dependent conductance of ChR2 is given by

$$ \psi(\phi)=(s_3+\gamma s_4)$$
(9)

where \(s_3\) and \(s_4\) are the population of the strong and weak open states, respectively, and \(\gamma \) is the ratio of the single-channel conductance in states \(s_4\) (\(G_{1,s4}\)) and \(s_3\) (\(G_{1,s3}\)): \(\gamma =G_{1,s4}/G_{1,s3}\) (the state \(s_4\) is usually much less conductive than \(s_3\) hence \(\gamma <1\)).

The set of equations (3)–(8) can be solved analytically at steady-state when \(\dot {s}_i=0\) (where \(i=1,\cdots , 6\)). The population of state \(s_i\) at steady state is given by:

$$ s_{i,\mathrm{steady-state}}=r_i/R, $$
(10)

where \(r_i\) is a compound populating rate (\(r_1=b_1 b_2 b_3\), \(r_2=a_1 b_2 b_3\), \(r_3=a_1 a_2 b_3\) and \(r_4=a_1 a_2 a_3\)) and \(R=\sum r_i\). The ratio between the population of the weak and the strong conducting states at steady state is \(s_4/s_3=a_3/b_2\). The resulting (normalized light dependent) conductance in steady state is therefore:

$$ \psi_{\mathrm{steady-state}}= \frac{r_3}{R}\left(1+\gamma\frac{a_3}{b_2}\right)$$
(11)

At high irradiance \(r_2\,,\,r_3 \gg r_1\,,\,r_4\) and \(\psi _{\mathrm {steady-state}}\) is governed solely by the rates between the open states:

$$\psi_{\mathrm{steady-state}}=\frac{r_2}{r_2+r_3} \left(1+\gamma\frac{a_2}{b_2}\right) = \frac{b_2+\gamma a_2}{a_2+b_2} ~. $$

The steady-state response of the model is therefore readily controlled by \(b_1\) and \(\gamma \) parameters.

At the transient peak \(\dot {s}_3=0\). At high irradiance, the transient peak conductance of ChR2 is approximately \(\psi _{\mathrm {peak}}\approx a_1/(a_1+ b_1+ a_2 )\). Since typically \(a_2<b_1\), the conductance peak is best controlled by the \(a_2\) and \(b_1\) parameters. The plateau-to-peak ratio, or the ChR2 adaptation ratio, \(\beta = \psi _{\mathrm {steady-state}}/\psi _{\mathrm {peak}}\) is therefore adjusted by the combination of the steady-state and peak parameters. These analytical results regarding the model equation (3)–(8), provide not only better insights into the behaviour of ChR2 but also help search in the parameter space for the best fit of experimental results.

2.2 ChR2 Model—Implementation

The empirical expressions for the light flux-dependent rates are given in the Table 1; they have the same form as in the previous model (Nikolic et al. 2009). We found it convenient to normalise the photon flux \(\phi \) to some value \(\phi _0\), which is below the threshold illumination, so that the transition rates which have weak logarithmic dependence on illumination have values equal to their spontaneous transition values in dark (i.e. when \(\phi \leq \phi _0\) we assume \(\phi =\phi _0\) and the log term is zero).

Table 1 Empirical expressions for the flux-dependent transition rates
Full size table

Parameter \(\overline {g}_{\mathrm {ChR}}\) can be deduced from a single photocurrent recording using saturating light and \(-\)70 mV clamped voltage (assuming the parameters in Eq. (2) are chosen so that \(f(-70)=1\)). A set of values for the model parameters is summarized in Table 2.

Table 2 Fitting parameters
Full size table

The ChR2 model was implemented in NEURON via a new mechanism called ChR2.mod. The ChR2 mechanism is implemented as a POINT_PROCESS module with an ELECTRODE_CURRENT. It is a shunt that is restricted to a small enough region so it can be described in terms of a net conductance and total current. In the localized ELECTRODE_CURRENT type point process positive currents depolarize the membrane while negative currents hyperpolarize it (Carnevale and Hines 2006). The six-state model scheme is described in a KINETIC block, in which a flow out of one state equals flow into another. The unknown of the model scheme, i.e. the relative population of the individual states, are declared in the STATE block, which cause the NMODL translator to convert it into a family of ODEs whose unknowns are the states (Hines and Carnevale 2000). The light dependent forward and reverse reaction rates are calculated in a separate PROCEDURE, which is called by the KINETIC block. The empirical model constant values are given in a PARAMETER block. The instantaneous flux of light is calculated in a self-events NET_RECEIVE block. The NET_RECEIVE block is essential to be used in order to define the different kinetics of the ChR2 and can only be implemented in a POINT_PROCESS allowing the current to change discontinuously. The kinetic model is integrated using the sparse method that separates the Jacobian evaluation from the calculation of the STATE derivatives, which is generally faster than computing the full Jacobian matrix (Carnevale and Hines 2006).

3 Results

3.1 ChR2 model validation

We sought validation of the six state model and its NEURON implementation against the same experimental dataset used to validate several previously described ChR2 models. Experimental procedures have been described in detail elsewhere (Nikolic et al. 2009; Grossman et al. 2011). We found the photokinetics of the six-state ChR2 model to be in good agreement with experimental data (Fig. 2). The peak and steady-state ChR2 currents and action potential latency were all well predicted over two orders of magnitude of irradiance (power per unit area). The recovery rate of the peak ChR2 current (due to slow transition between \(s_6\) and \(s_1\) states) was also well described by the six-state model for inter-pulse intervals between 1 and 10 s. It was also possible to accurately reproduce experimental photocurrent waveforms for both constant and pulsed illumination, see Fig. 2(e).

Fig. 2
figure 2

Experimental validation of the six-state ChR2 model implemented in NEURON simulation environment. Experimental results are shown with black lines (error bar represents standard deviation) and the simulation results are presented with red lines. (a) Photokinetics results of ChR2 conductance normalized to results at 90 mW/mm\(^2\), voltage clamp at \(-\)70 mV for: peak ChR2 current vs. light intensity (top, left panel), peak current recovery vs. inter-pulse interval (top-right) and plateau (steady-state) ChR2 current vs. light intensity (bottom-left). AP latency vs. irradiance is shown in bottom-right panel. (b) ChR2 photocurrent waveforms at constant and pulsed illuminations of various intensity. Scale bar: long pulse 50 ms, 0.25 nA; pulse train 200 ms, 0.5 nA

Full size image

3.2 ChR2 density and Threshold Irradiance

In order to compare the results and effects of illuminating different parts of a neural cell it is necessary to establish some meaningful metrics. For example, if constant concentration of ChR2 is assumed everywhere in the cell’s membrane, than by illuminating just the soma, much less photo-current is “injected” into the cell in comparison to the whole cell illumination, and then say the threshold illumination will be different and hence ChR2 adaptation will also differ making any comparison difficult. One possible way of comparing different spatial illumination patterns would be to compare the ChR2 densities (\(\overline {g}_{\rm ChR}\)) required for a fixed threshold illumination. Table 3 summarizes the results when the threshold was set to a single 20 ms light pulse of 1 mW/mm\(^2\) and then the parameter \(\overline {g}_{\rm ChR}\) was tuned to give the same threshold irradiance for spiking for each illumination case. Results for the case of targeting only AIS or other parts of the axon are given in a separate table in Discussion section in order to make it clear that a different model was used for that case. In this study we chose to calibrate the expression for the illumination threshold of 1 mW/mm\(^2\) which is equal to the half saturation irradiance of ChR2 (Grossman et al. 2011), thus allowing a good dynamic range of illuminations.

Table 3 ChR2 expression values for Layer V pyramidal neuron which would give the threshold irradiance of 1 mW/mm\(^2\) for a 20 ms pulse
Full size table

3.3 L5 neuron model optimised for both somaticand dendritic spiking (“High”)

This section uses the (Hay et al. 2011) layer V pyramidal neuron model that was optimized to replicate the experimental mean of both the somatic Na\(^+\) and the dendritic Ca\(^{2+}\) spiking characteristics. It has a high dendritic Na\(_{\rm t}\) density of 107 pS/\(\mu \)m\(^2\), but this should be biologically plausible (Stuart and Sakmann 1994; Kole et al. 2008). The model was modified to incorporate ChR2 based light sensitivity, as described above. The results are shown in Fig. 3. Figure 3(a)-Left shows the reconstructed morphology of the pyramidal cell used and the sites of recording in the soma and along the dendritic tree. Figure 3(a)-Right compares the neural response for a current injection pulse and a single short (5 ms) and saturating (40 mW/mm\(^2\)) light pulse on: the soma, the apical dendrites and the whole cell. Overall, they show a bursting response (multiple spikes long after the stimulus end), which is a manifestation of the interplay between the somatic and the dendritic spiking sources. The illumination that targeted the soma, first elicited a somatic spike that then propagated into the dendritic tree causing the generation of a local Ca\(^{2+}\) dendritic spike that then propagated back to the soma and triggered additional somatic AP. After the end of the spiking burst, both the soma and the dendrites repolarized quickly. In contrast, the illumination that targeted the apical dendrites evoked a stronger and more prolong depolarization of the dendrites, yielding larger Ca\(^{2+}\) spikes, which then propagated to the soma to evoke the first somatic action potential. This AP triggers another Ca\(^{2+}\) spike and this sequence repeats once more. After the end of the spiking burst, the dendritic tree remains highly deplolarized.

Fig. 3
figure 3

Model with both perisomatic and dendritic action potentials. (a) Left: Reconstructed morphology of a layer V neocortical pyramidal neuron. Recording sites are indicated by schematic electrodes at the soma (black), proximal apical dendrite (360 \(\mu \)m from the soma, red), middle apical dendrite (670 \(\mu \)m from the soma, green) and distal apical dendrite (1250 \(\mu \)m from the soma, blue). Right: Membrane potential at the recording sites evoked by an injection of 5 ms 2 nA current in the soma (top panel), or a single 5 ms 40 mW/mm\(^2\) light pulse on the soma, distal apical dendrite or whole cell (next three panels). (b) Left: Firing response recorded at the soma for a 40 mW/mm\(^2\) step illumination on the soma (red) or on the whole cell (black). Right: Corresponding instantaneous firing frequencies i.e. inverse inter-spiking times (red: soma, blue: apical dendrite, black: whole cell illumination). (c) Firing response for a 5 Hz train of 5 ms 40 mW/mm\(^2\) light pulses on the soma (red) or on the whole cell (black). Bar plot indicates the corresponding number of action potentials per light pulse. (d) F–E curve (where F is the initial spiking frequency and E is the light irradiance) for somatic (red) and whole cell (black) illuminations

Full size image

The illumination of the whole cell, in fact combined the characteristics seen in the somatic and in the dendritic illuminations. The synchronized excitation of the somatic and the dendritic sources yielded higher bursting frequency, and relatively fast repolarization of the neuron after the termination of the burst. Interestingly, the spiking dynamics that was evoked by illuminating the whole cell resembled the spiking output from a direct injection of current to the soma (5 ms pulse of 2 nA). The threshold voltage in which the neuron spiked during the direct current injection was slightly lower than the threshold voltage during the somatic and the whole cell illumination and was slightly higher than the threshold during the dendritic illumination.

Figure 3(b) compares the neuron’s response (recorded atthesoma)toastepofasaturating(40 mW/mm\(^2\)) light on the soma and on the whole cell. Overall, it shows that following a short period of bursting, the neuron switches to a steady non-bursting spiking pattern. The characteristics of the steady state spiking were affected by the location of the illumination. The illumination that targeted the soma yielded spiking with a constant interval between the spikes, while the illumination of the whole cell repeatedly evoked a pair of relatively close spikes. This feature is evident in the instantaneous frequency plot in Fig. 3(b)-Right.The frequency of the initial spikes bursting depends on the light irradiance, as can be seen in Fig. 3(e). The initial spiking frequency was calculated from the period between the first and second action potential. The analysis excludes depolarization spikelet transients below 0 mV. It shows that overall a somatic illumination resulted in slightly slower bursting at low irradiances and faster bursting at higher irradiances. Both types of illumination had similar half saturation irradiance \(E_{1/2} \approx 1.5\,\)mW/mm\(^2\), indicating a ChR2 limited spiking (Grossman et al. 2011). Figure 3(c) compares the neural response for a 5 Hz train of short (5 ms) and saturating (40 mW/mm\(^2\)) light pulses on the soma and on the whole cell. Overall, it shows that the first light pulse elicited a burst response that was then stabilized to a single spike output per light pulse input. In this case, there were no major differences between the illumination that targeted the soma and the one that targeted the whole cell.

3.4 L5 Neuron model with only somatic spiking (“Low”)

This section uses a modified Hay et al. layer V pyramidal neuron model which was optimised for accurately reproducing the perisomatic step current firing, but results in a strong attenuation of the bAP and a lack of Ca\(^{2+}\) spikes (Hay et al. 2011). In this case the apical dendrites have a lower fast inactivating Na\(^+\) (Na\(_t\)) conductivity density and higher fast inactivating (Kv3.1) and muscarinic (\(I_m\)) K\(^+\) densities, while maintaining somatic parameters in the same range. The results are shown in Fig. 4. Here, illuminations that targeted the apical dendritic tree yielded only subthreshold excitation, even with ChR2 density levels that were 3 orders of magnitude higher than the level for the 1 mW/mm\(^2\) somatic threshold (results not shown). Figure 4(a) compares the neural response for a single short (5 ms) and saturating (40 mW/mm\(^2\)) light pulse on the soma and the whole cell. In contrast to the previous model, here the neuron responded to the short light pulse with a single spike and no bursting. The illumination that targeted only the soma elicited an action potential that was followed by subthreshold excitation of the dendrites that then induced subthreshold after spike somatic depolarization. In comparison, the illumination that covered the whole cell yielded larger dendritic depolarization and hence more pronounced after spike depolarization at the soma. The spike waveform elicited by the illumination of the soma resembled the spike waveform that was created by a direct injection of current (5 ms pulse of 2 nA) to the soma, though with larger after spike depolarization. The spiking threshold was similar to the somatic current injection and was slightly lower than the whole cell illumination. Figure 4(b) compares the neuron’s response to a step of a saturating (40 mW/mm\(^2\)) light on the soma and on the whole cell. It shows that in this type of neuron, both types of illumination yielded a steady spiking, after the fifth action potentials, with a similar spiking frequency. The illumination that targeted the soma yielded a higher initial frequency associated with a short period of bursting. The dependency of the initial spiking frequency on the stimulating irradiance is plotted in Fig. 4(d). It shows that at low irradiances, illuminations of the soma and the whole cell elicited at the beginning a burst pattern. Again both types of illumination had similar half saturation irradiance (\(E_{1/2}\approx 1.5\,\)mW/mm\(^2\)), indicating a ChR2 limited spiking (Grossman et al. 2011). Interestingly, in this model the illumination of the whole cell at high irradiances (\(>10\) mW/mm\(^2\)) stopped the initial bursting. Figure 4(c) compares the neural response for a 5 Hz train of short (5 ms) and saturating (40 mW/mm\(^2\)) light pulses on the soma and on the whole cell. Overall, it shows that in this neural model, the pulse stimulation yields a low fidelity, with a slightly better spiking outcome when the whole cell is illuminated.

Fig. 4
figure 4

Results for L5 cell model constrained only by perisomatic step current firing and bAPs not guaranteed. (a) Left: As in Fig. 3(a), left. (a) Right: Comparison of the cell response to an injection of 5 ms 2 nA current in the soma, or a single 5 ms 40 mW/mm\(^2\) light pulse on the soma and whole cell. (b) Left: Firing response for a 40 mW/mm\(^2\) step illumination on the soma (red) or on the whole cell (black). Right: Corresponding instantaneous firing frequencies (red: soma, black: whole cell). (c) Firing response for a 5 Hz train of 5 ms 40 mW/mm\(^2\) light pulses on the soma (red) or on the whole cell (black). Bar plot indicates the corresponding number of action potentials per light pulse. (d) F–E curve for somatic (red) and whole cell (black) illuminations

Full size image

3.5 L5 Neuron: axonal excitation

This section uses the model proposed by Hu et al. (2009) which is based on the previous multicompartmental model of the full dendritic and somatic structures (Mainen and Sejnowski 1996) and the AIS and axon structures (Yu et al. 2008; Shu et al. 2006). Although the model has a very high (and probably infeasible) excitability of the apical dendrites (manifested for example in an increase in the amplitude of back propagating APs), it has a very good representation of the soma and axonal compartments. Since the level of the dendritic excitability does not have a strong effect on the response of the soma to a somatic step stimulation (Hay et al. 2011), it is reasonable to believe that it will have a minor effect as well on the spikes that are induced in the axon. The results of simulations are shown in Fig. 5. Figure 5(a) compares the neural response for a single short (5 ms, 5 mW/mm\(^2\)) light pulse on the soma and on the axon initial segment (AIS). It shows that APs generated by AIS illumination were elicited at lower threshold voltages and have a shorter latency, a larger after spike repolarisation and a smaller after spike depolarization in comparison to somatic illumination. A step illumination of the AIS fails to generate spikes trains, as can be seen in Fig. 5(b). In fact, in the case of AIS illumination, only a narrow range of irradiances (between 2 and 4 mW/mm\(^2\), Fig. 5(c)) was able to evoke more than one action potential under a constant illumination. Figure 5(d) compares the fidelity of pulse stimulation at different light irradiances. The fidelity was measured as the percentage of successful conversions from a light pulse to an AP. It shows that in both cases the fidelity of stimulation dropped at the same frequency in both AIS and somatic illuminations (though AIS has slightly lower fidelity at high frequency).

Fig. 5
figure 5

Illumination pattern affects bAPs. Left: Latencies (relative to time of somatic spike) and amplitudes of a backpropagating action potential along an apical dendrite (highlighted in red) using a 5 ms pulse of 10 mW/mm\(^2\). Right: Waveform of a bAPs train at three locations (490 \(\mu \)m, 913 \(\mu \)m and 1637 \(\mu \)m) along the apical dendrite, induced by a constant illumination of 50 mW/mm\(^2\). Threshold irradiance was 1 mW/mm\(^2\)

Full size image

3.6 Illumination impact on back-propagating APs

Figure 6(a) compares the latency (relative to the time of the somatic spike) and the amplitude of bAPs along an apical dendrite during somatic and whole cell illuminations (5 ms pulses, 40 mW/mm\(^2\)). Overall it shows that illuminations of the whole cell yielded shorter latencies and higher amplitudes of bAPs. The neuron model with a reduced excitability of the apical dendritic tree resulted in faster propagation of the action potential. In both neuron models, the illumination of the whole cell, in comparison to illumination of the soma, yielded bAP at the distal part of the apical tree that arrived almost 5 ms earlier and was 10 mV higher.

Fig. 6
figure 6

Back propagating action potential (bAP) modulated by illumination patterns. (a) bAP latency relative to time of somatic spike (top) and height (bottom), for somatic illumination (magenta) and whole cell illumination (blue). Marks: square for neuron model with both somatic and dendritic spiking (High: strong dendritic excitability), triangles: for the neuron model with only somatic Na\(^+\) spiking (Low: weak dendritic excitability). (b) Left: Morphological illustration of the cell showing a patterned illumination of a single dendritic branch (blue) together with the recording sites along the dendrite and a nearby reference branch (red). Right: AP latency measured from the onset of light and AP height of the induced bAPs (light pulse: 20 ms, 1 mW/mm\(^2\)) for the two models: High (top) and Low (bottom)

Full size image

4 Discussion and Conclusion

Traditionally, neural stimulation protocols have targeted predominantly the somatic compartment. Optogenetics offers for the first time the flexibility to localize the stimulation at almost arbitrary neuron sites. The effect of such patterned excitation on the spiking output is just beginning to be uncovered. This study aimed to theoretically investigate four, most commonly considered, such excitation patterns: focal somatic illumination, whole cell illumination, apical dendritic illumination and focal AIS/axon illumination.

  1. 1.

    Somatic Illumination In order to achieve 1 mW/mm\(^2\) threshold irradiance during somatic illumination, the cell needed to have a ChR2 conductance density of approximately 35 pS/\(\mu \)m\(^2\) (Table 3). Assuming the maximal conductance of a single channel is in the region 40–100 fS for the native ChR2 (Feldbauer et al. 2009; Bamberg et al 2008; Lin et al. 2009), then the neuron requires at least 350–900 ChR2 complexes per \(\mu \)m\(^2\) area of the membrane. In comparison, the conductance density of the Na\(_t\) channels in the soma is \(\sim \) \(20\,\)nS/\(\mu \)m\(^2\) (Kole et al. 2008) but the single channel conductance is more than three orders of magnitudes higher (\(\sim \) \(20\,\)pS (Tour et al. 1998)), so the density of the native Na\(_t\) channels was in fact of the same order of magnitude (\(\sim \)1000 complexes per \(\mu \)m\(^2\)). Figure 7 shows the impact of the ChR2 expression level (measured as the threshold irradiance) on the neural frequency-irradiance response curve. It illustrates that at low threshold irradiances (high expression of ChR2), the photocurrents can exceed the capacity of the cell to respond with generating APs (i.e. a neuron limited spiking) and the AP frequency can be damped, as in the case of neural excitation by strong current injection (also shown in Fig. 7). In contrast, at high threshold irradiances (low expression of ChR2), the photocurrents saturate before the maximal spiking frequency is achieved. Overall, the plots strengthen the need for a correct adjustment of the density (expression) level of ChR2 in order to achieve optogenetic experiments with a good dynamic range.

  2. 2.

    Dendritic Illumination A direct excitation of APs (i.e. suprathreshold depolarization) via illumination of the dendritic tree requires dendrites with a relatively strong excitability and a high density of the ChR2 (see Table 3). In the model used here, a dendritic illumination required almost double the density of ChR2 than the soma for a given threshold irradiance (Table 3). In the same way, for a given expression level, a dendritic illumination strategy requires almost double the light irradiance than a somatic illumination. This is despite the fact that illuminated area was approximately 18 times larger than the area of the soma. Thus, the dendritic illumination required in total 30 times more conductivity of ChR2 to achieve a similar threshold. The spikes for dendritic illumination show prolonged post spike depolarization.

  3. 3.

    Whole cell Illumination An illumination strategy that synchronizes the depolarization of the dendritic tree and the soma evokes spiking with particular characteristics including an increased bursting rate (in the case of high excitable dendrites, Fig. 3(a), a stronger after spike depolarization (in the case of low excitable dendrites, Fig. 4(a), and an enhanced back propagation of the APs, Fig. 6(b). Overall this type of illumination seems to be the most efficient. Our simulations show that a given threshold irradiance was achievable with only 6 % of ChR2 density needed in the case of somatic illumination (Table 3). In the same manner, for a given expression level of ChR2, whole cell illumination strategy requires an order of magnitude less light irradiance. However, as the effective illuminated area was much larger (27 times larger) than the area of the soma, the total amount of ChR2 ion channels that was required was in fact of a similar magnitude. Similarly, for a given ChR2 density, targeting the whole cell or jus the soma requires almost the same level of optical power. Interestingly, in both cases the threshold is not sensitivity to the excitability state of the dendritic tree. Our results imply that a stimulation strategy that combines a wide field illumination of the neural network using irradiances that are slightly below the whole cell threshold together with a somatic illumination of the targeted neuron using irradiances that are above the whole cell threshold will allow optogenetic stimulation with low light irradiance and high specificity.

  4. 4.

    Axonal Illumination Illuminations that target the AIS have lower voltage thresholds however they cannot sustain a prolonged illumination, as shown in Fig. 5(a) where only one or a few spikes are elicited. Overall, AIS illumination requires higher density of ChR2 to achieve a given threshold (in this study it required five times higher density of ChR2 than somatic illumination, see Tables 3 and 4. In the same manner, for a given expression level, AIS illumination has a higher illumination threshold. This observation allows us to predict that illuminations with irradiances at the somatic and the whole cell thresholds range will not co-excite nearby axon fibres of other cells.

    Comparisons between the ChR2 expression densities required to achieve a given threshold irradiance in the case of unmyelinated and myelinated axons (Table 4) showed that an unmyelinated AIS with an illumination length of 50 \(\mu \)m required an order of magnitude denser ChR2 than the neighbouring unmyelinated (naked) axon of 400 \(\mu \)m length. However, when similar size myelinated axon is illuminated about seventy times higher ChR2 density is required, but note that only ChR2 in the nodes of Ranvier will be effective. We also investigated how the level of the axon myelination affects the threshold of stimulation with ChR2. The variation of myelin was represented in the model by changes in the membrane capacitance. The insulating myelin sheaths had 25 times smaller membrane capacitance and the nodes had larger almost no active ion channels in the axon segments that are covered by myelin. In contrast, there was a relatively high (half of the AIS) density of Na\(^+\) channels in the nodes of Ranvier. ChR2 was expressed only in the nodes. We found that the threshold irradiance increases in proportion to the demyelination level of an axon see Table 4. Based on this observation we predict that measurements of the threshold irradiance (for example relative to the soma) could provide a simple (yet powerful) tool to remotely probe a loss of neural myelin sheath, which is a hallmark of several neurodegenerative diseases, including leukodystrophies and multiple sclerosis.

Fig. 7
figure 7

Impact of ChR2 density (which is inversely proportional the threshold irradiance \(E_t\)) on the \(F-E\) curve, where \(F\) is the initial spiking frequency (inverse time between the first two spikes) and \(E\) is the light irradiance. Top: Model with both perisomatic and dendritic action potentials (High). Bottom: Model with only perisomatic action potentials (Low). Simulated threshold irradiances were 0.1, 1 and 10 mW/mm\(^2\), which roughly correspond to experimentally observed high, normal and low expression levels

Full size image
Table 4 Minimum ChR2 expressions needed to achieve a threshold of 1 mW/mm\(^2\), 20 ms pulse by illuminating some axon parts
Full size table

Effects on Backpropagating Action Potentials

Our results show that co-excitation of the dendritic tree and the soma enhances the conductivity of the back propagating action potential. This observation brought us to speculate that a patterned illumination may be used to modulate the back propagating APs and hence spatially tailor the direction and amplitude of spike time dependent plasticity protocols. We simulated such a patterned illumination scenario in Fig. 6(b). It compares the latencies and heights of the bAPs along an illuminated dendritic branch and a nearby branch that is not illuminated. Overall, the results confirmed that the AP arrived to the tip of the illuminated dendritic branch at much shorter latencies and with higher amplitudes for both models of the dendritic excitability. This prediction can be experimentally tested by using patterned light sources, e.g. using micro-LED array (Grossman et al. 2010), phase modulation of light (Yang et al. 2011) or two-photon excitation (Papagiakoumou et al. 2010).

In summary, here we have described how ion channel mechanisms modulated by arbitrary external light pulses can be implemented in the NEURON simulation platform. This enables the use of computational modelling of optogenetic neurophysiology experiments, which we believe will substantially improve the interpretation of optogenetic results in neuroscience. Our methodology enables us to simulate the excitatory and inhibitory inputs to neurons at an arbitrary sub-cellular compartment. We used our models to derive illumination strategies that can improve the efficiency and selectivity of optogenetic protocols. Of particular interest was the examination of how different illumination patterns modulate dendritic activity, which affects activity-dependent synaptic plasticity and therefore influences information processing within neurons. We furthermore analysed the effect of a neurodegenerative process such as axon demyelination on the required stimulus, suggesting how this could be used in detecting loss of myelin sheaths. At the present time, our approach has been applied only to multi-compartment, morphologically realistic single neurons. We consider that it is now feasible to extend “computational optogenetics” to the level of neural networks, which will offer substantial insight into the perturbation and control of neural ensemble activity patterns, with applications to both understanding neural circuit function, and repairing brain disorders.