Skip to main content
Advertisement
  • Loading metrics

Extending Integrate-and-Fire Model Neurons to Account for the Effects of Weak Electric Fields and Input Filtering Mediated by the Dendrite

Abstract

Transcranial brain stimulation and evidence of ephaptic coupling have recently sparked strong interests in understanding the effects of weak electric fields on the dynamics of brain networks and of coupled populations of neurons. The collective dynamics of large neuronal populations can be efficiently studied using single-compartment (point) model neurons of the integrate-and-fire (IF) type as their elements. These models, however, lack the dendritic morphology required to biophysically describe the effect of an extracellular electric field on the neuronal membrane voltage. Here, we extend the IF point neuron models to accurately reflect morphology dependent electric field effects extracted from a canonical spatial “ball-and-stick” (BS) neuron model. Even in the absence of an extracellular field, neuronal morphology by itself strongly affects the cellular response properties. We, therefore, derive additional components for leaky and nonlinear IF neuron models to reproduce the subthreshold voltage and spiking dynamics of the BS model exposed to both fluctuating somatic and dendritic inputs and an extracellular electric field. We show that an oscillatory electric field causes spike rate resonance, or equivalently, pronounced spike to field coherence. Its resonance frequency depends on the location of the synaptic background inputs. For somatic inputs the resonance appears in the beta and gamma frequency range, whereas for distal dendritic inputs it is shifted to even higher frequencies. Irrespective of an external electric field, the presence of a dendritic cable attenuates the subthreshold response at the soma to slowly-varying somatic inputs while implementing a low-pass filter for distal dendritic inputs. Our point neuron model extension is straightforward to implement and is computationally much more efficient compared to the original BS model. It is well suited for studying the dynamics of large populations of neurons with heterogeneous dendritic morphology with (and without) the influence of weak external electric fields.

Author Summary

How extracellular electric fields—as generated endogenously or through transcranial brain stimulation—affect the dynamics of neuronal populations is of great interest but not well understood. To study neuronal activity at the network level single-compartment neuron models have been proven very successful, because of their computational efficiency and analytical tractability. Unfortunately, these models lack the dendritic morphology to biophysically account for the effects of electric fields, and for changes in synaptic integration due to morphology alone. Here, we consider a canonical, spatially extended model neuron and characterize its responses to fluctuating synaptic input as well as an oscillatory, weak electric field. In order to accurately reproduce these responses we analytically derive an extension for the popular integrate-and-fire point neuron models. We show that the dendritic cable acts as a filter for the synaptic input current, which depends on the input location, and that an electric field modulates the neuronal spike rate strongest at a certain (preferred) field frequency. These phenomena can be successfully reproduced using integrate-and-fire models, extended by a small number of components that are straightforward to implement. The extended point models are thus well suited for studying populations of coupled neurons with different morphology, exposed to extracellular electric fields.

Introduction

Extracellular electric fields in the brain and their impact on neural activity have gained a considerable amount of attention in neuroscience over the past decade. These electric fields can be generated endogenously [13] or through transcranial (alternating) current stimulation [46], and can modify the activity of neuronal populations in various ways [1, 79]. Although the fields generated by this type of noninvasive brain stimulation are rather weak (≤1 V/m [4, 5]) and do not directly elicit spikes, they can modulate spiking activity and lead to changes in cognitive processing, offering a range of possible clinical interventions [1012]. How external fields lead to changes of the membrane voltage in single cells has been studied in detail [1315]. However, their effects on population spike rate and the underlying mechanisms are largely unexplored.

Computational models of neurons exposed to electric fields offer a useful tool to gain a better understanding of these mechanisms. Multi-compartment models of neurons are well suited for corresponding investigations at the level of single cells and small circuits [16] but are too complex for a purposeful application in large populations. Single-compartment (point) neuron models of the integrate-and-fire (IF) type are well applicable to study the dynamics of large neuronal populations, due to their computational efficiency and analytical tractability [17]. However, typical IF model neurons lack the dendritic morphology required for a biophysical description of electric field effects. Furthermore, even in the absence of an extracellular field, the dendritic morphology strongly shapes neuronal response properties [18].

In this contribution, we extend the popular class of IF point neuron models to quantitatively account for morphology dependent modulations of neural activity due to: (i) dendritic influences on the integration of synaptic inputs and (ii) the effects of extracellular electric fields. Furthermore, we describe how oscillatory electric fields affect neuronal subthreshold and spiking activity and identify field-induced spike rate resonance. Specifically, we considered a canonical spatial pyramidal neuron model which consists of a somatic compartment and one (apical main) passive dendritic cable, and which is exposed to in-vivo like fluctuating synaptic input and an electric field. Based on that model we analytically derived an extension for the classical leaky and the refined exponential, [19], IF point neuron models in order to exactly reproduce the subthreshold dynamics of the spatial model for arbitrary parametrizations. We then evaluated the extended IF models by quantitatively comparing their spiking activity with the spiking activity of the corresponding spatial model. Finally, we used these models to study the effects of an oscillating electric field (due to the presence of the dendritic cable) on the spike rate dynamics.

Results

Our derivation of the extended point neuron model consists of two steps. We first calculate the somatic membrane voltage of a ball-and-stick (BS) model in response to subthreshold synaptic inputs at the soma and the distal dendrite and to a time-varying, spatially homogeneous, extracellular electric field. This involves solving a generalized cable equation [20]. Second, we seek to exactly reproduce this voltage response in the point neuron model by deriving additional model components (see Fig 1): two linear temporal filters, one for each input location, to be applied to the “raw” synaptic input and one additional input current to describe the field effect. The model components are given in analytical form and depend on the parameters of the BS model and the electric field. We refer to the model equipped with the new components as the extended point (eP) neuron model. We first derive this extension for the well-known leaky IF (LIF) neuron model, and present the extension adapted for the exponential IF (EIF) neuron model in a separate section.

thumbnail
Fig 1. Diagram of the extended point neuron model.

Top: Visualization of the ball-and-stick, BS, (left) and the extended point, eP, (right) neuron models. Both models receive synaptic input currents at the soma and the distal dendrite, Is(t) and Id(t), and are exposed to an external electric field E(t). Ls(t) and Ld(t) denote the additional input filters describing the dendritic effects. IE(t) denotes the additional input current describing the field effect. VeP(t) and VBS(0, t) denote the corresponding membrane voltages (at the soma). Bottom: Electrical circuit diagram for the subthreshold dynamics of the BS model. For a description of the parameters and their values see Table 1.

https://doi.org/10.1371/journal.pcbi.1005206.g001

Models

The BS neuron model consists of a lumped soma attached to a passive dendritic cable of length L. The dynamics of its membrane voltage, when receiving synaptic inputs at the soma, Is(t), and the distal dendrite, Id(t), and when exposed to a spatially homogeneous external electric field, E(t), are governed by the cable equation: (1) subject to the boundary conditions: (2) (3) at the soma (x = 0) and the end of the dendrite (x = L). VBS denotes the deviation of the membrane voltage from rest, Vrest, VBS(x, t) ≔ VBS,i(x, t) − VBS,e(x, t) − Vrest, where VBS,i and VBS,e are the intra- and extracellular potentials. The effects of a spike are described by the IF-type reset condition for the soma: (4) and by a short refractory period of length Tref during which VBS(0, t) is clamped at the reset value Vr. Spike times are defined by the times at which the somatic membrane voltage VBS(0, t) crosses the spike voltage value Vs from below. cm denotes the membrane capacitance, gm the membrane conductance, and gi the internal (axial) conductance of a dendritic cable segment of unit length. Cs and Gs are the somatic membrane capacitance and leak conductance. The exponential term with threshold slope factor ΔT and effective threshold voltage VT approximates the somatic sodium current at spike initiation [19]. For details see Methods.

In the proposed IF point neuron extension, that is, the eP model, the deviation of the membrane voltage, VeP, from rest is governed by (5) and by the reset condition: (6) where VeP is clamped to for the duration of the refractory period Tref after every spike. CeP and GeP are the membrane capacitance and leak conductance. The scaling factor α ensures an equal membrane voltage response to the depolarizing current described by the exponential terms in both models (BS and eP). We consider two versions of these models separately. First, we treat the LIF versions in detail, where we omit the exponential terms in Eqs 2 and 5; specifically, by taking the limit ΔT → 0 (and setting Vs = VT). In the subsequent part we then consider the (full) EIF versions of the BS and eP models. Below we explain in detail how the components of the point model extension are derived: the linear input filters Ls(t), Ld(t), the additional input current equivalent to the field effect, IE(t), and, in case of the (full) EIF type models, the scaling factor α. The analytical expressions of these model components are given in Eqs 10, 13, 20 and 21 (for the LIF case), and in Eqs 2226 (for the EIF case). To mimic the remaining depolarization along the dendritic cable after each spike, we choose an elevated reset voltage for all point neuron models: .

For comparison we also use a point neuron model (of LIF and EIF type, respectively) without the extension, that is, Ls(t) = Ld(t) = δ(t) and α = 1, and we fit the parameters of that model to best reproduce the activity of the BS model for equal synaptic inputs (details see below). We refer to this model as the P model.

The somatic input filter for the LIF model

We first consider the BS and eP model neurons of the LIF type (i.e, ΔT → 0, Vs = VT) receiving subthreshold synaptic input at the soma in the absence of an electric field (E(t) = 0, IE(t) = 0, Id(t) = 0). To avoid ambiguity we use the superscript Is for the membrane voltage variables in this case. The somatic membrane voltage response of the BS model (Eqs 13) can be calculated as (see Methods) (7) (8) where indicates the temporal Fourier transform and ω = 2πf denotes angular frequency. The somatic membrane voltage response of the eP model (Eq 5) is given by (9) The dendritic filter Ls required to exactly reproduce the somatic membrane voltage response of the BS model, i.e., , must then be equal to ratio of the impedances of both models: (10) where z(ω) is given by Eq 8. In the following, we choose the membrane capacitance and conductance of the eP model to be equal to the corresponding somatic quantities of the BS model, CeP = Cs, GeP = Gs. To see the necessity of the filter, let us consider the P model (no dendritic filter, ) whose subthreshold response is given by (11) Because of the additional frequency-dependent term in the denominator of Eq 7 compared to Eq 11, it is not possible to adjust the parameters CP and GP of the P model such that for all frequencies ω. The somatic response of the BS model can only be approximated in this case.

Fig 2A shows the impedances, , m ∈ {BS|x = 0, eP, P}, of the three neuron models for an example set of parameter values for the BS model. The two parameters of the P model (CP and GP) were determined by matching the steady-state somatic voltage, , and minimizing the mean squared distance between and over the visualized range of input frequencies. The impedance of the eP model matches the impedance of the BS model exactly while the impedance of the P model deviates substantially, in particular for larger frequencies.

thumbnail
Fig 2. Impedance and filters for somatic inputs.

A: Impedances , , and of the three neuron models as a function of input frequency. B: Gain and phase of the input filter as a function of input frequency. The neuronal morphology varied as indicated, in terms of dendritic cable length (350 μm, 700 μm, 1050 μm), cable diameter (0.6 μm, 1.2 μm, 1.8 μm) and soma diameter (5 μm, 10 μm, 15 μm). * indicates the default parameter values. For all other parameter values used see Table 1.

https://doi.org/10.1371/journal.pcbi.1005206.g002

Fig 2B–2D show the amplitudes and phases of the input filter for various sets of parameters for the BS morphology. is always a high-pass filter, which attenuates the somatic inputs at lower and amplifies them at higher frequencies. This effect is more pronounced for a larger dendritic and a smaller somatic compartment. It becomes stronger with increasing ratio of dendritic over somatic size. Nevertheless, the filter does not differ qualitatively for changes in neuron morphology.

We next compare how well the point neuron models eP and P reproduce the spiking activity of the BS model neuron. For this purpose we consider an in vivo-like fluctuating synaptic input current Is(t) described by an Ornstein-Uhlenbeck process (see Methods). The model outputs are compared over a range of values for the input mean and standard deviation σs. The parameter values of the P model were adjusted to best reproduce the spike train of the BS model (see Methods for details). Fig 3A displays the time series of the somatic membrane voltage of the three models in response to the same input currents—a weak (small , σs) and a strong current (large , σs). For both input currents, the eP model well reproduces the somatic voltage dynamics of the BS model. Consequently, the spike times are also well reproduced. There is, however, a mismatch between the voltage traces during short periods (of less than approximately 10 ms duration) after spikes have occurred. This discrepancy is a result of the remaining dendritic depolarization after a spike has occurred in the BS model, which is only approximated by the elevated reset voltage (see section Models above) in the point neuron models. In comparison, the P model performs worse in reproducing the BS membrane voltage dynamics, particularly the fast fluctuations are poorly recovered. This is expected from the mismatch in the impedance for high frequencies (cf. Fig 2A).

thumbnail
Fig 3. Reproduction of spiking activity for somatic inputs using LIF type models.

A: Membrane voltage traces of the BS (blue), eP (green) and P (red) neuron models in response to a weak ( pA, σs = 11.94 pA, top) and a strong input current ( pA, σs = 33.34 pA, bottom). The parameter values of the P model were tuned independently to maximize the coincidence factor ΓBS,P for each set of input parameters. B: Coincidence factor for the BS and eP model spike trains, ΓBS,eP (left), and for the BS and P model spike trains, ΓBS,P (right) as a function of input mean and standard deviation σs. C: Difference ΓBS,eP − ΓBS,P between the coincidence factors shown in B. D: Spike rate difference of the BS and eP models (left) and of the BS and P models (right) as a function of and σs. E: Spike rate of the BS neuron model. The input parameters used in A are indicated in B. Results presented in B-E are averages over 6 noise realizations. The parameter values of the BS model are listed in Table 1.

https://doi.org/10.1371/journal.pcbi.1005206.g003

In Fig 3B–3E we compare spiking activity in terms of spike coincidences and spike rates for a wide range of input parameters. We used the spike coincidence measure Γ which quantifies the similarity between two spike trains for a given precision of 3 ms (see Methods). The maximum value of 1 indicates an optimal match, i.e., spike times always coincide, a value of 0 corresponds to pure chance, i.e., the degree of coincidences for two Poisson spike trains with equal rates. The P model was fitted to the BS model for each input (in terms of , σs) separately. The parameters of the eP model, on the other hand, are constant and do not depend on the input at all. The eP model very accurately reproduces the BS spike times for small spike rates (Γ ≥ 0.9 for small and σs), see Fig 3B and 3E. This performance decreases only slightly for increasing σs (noise dominated input) and somewhat stronger for increasing (mean dominated input). Generally, Γ decreases with increasing spike rates. This can be attributed to the transient periods after spikes during which the dendritic cable is still loaded and the membrane voltages of both neuron models deviate. Those periods do not depend on the spike rate and therefore have a stronger deteriorating effect when the interspike intervals are smaller. In addition, when σs is small the model neurons spike repetitively in a rather clock-like manner, with comparable rate but most likely out of phase due to mismatches caused by the membrane voltage resets. This helps understand the rather low values of Γ for mean dominated inputs. The spike rate of the BS model is also reproduced quite well by the eP model, which underestimates it only slightly (Fig 3D). Spike coincidence and spiking rate reproduction of the eP model can be improved even further by additionally tuning the reset voltage using Γ or the spike rate distance as a cost function. The P model, in comparison, is substantially worse in reproducing the spike times at small spike rates and only slightly better than the eP model for large spike rates (Fig 3B and 3C). The spike rate of the BS model is slightly overestimated by the P model (Fig 3D). Even though the parameters of the P model were optimized in an input-dependent manner the eP model leads to an improved reproduction of the BS spiking activity overall.

In summary, the dendritic cable implements a high pass filter for inputs at the soma. Due to the derived filter for somatic inputs, the eP model—without having fitted any of its parameters—well reproduces the BS model dynamics for subthreshold and suprathreshold inputs. Notably, the computation time required for the BS model was at least 25 times that of the eP model, using measurements on a single core of a desktop computer.

The distal input filter for the LIF model

We next consider subthreshold synaptic input at the distal dendrite instead of somatic input, but otherwise the same setup as in the previous section. Here we use superscipt Id for the membrane voltage variables to better distinguish from the previous scenario. The somatic membrane voltage response of the BS model can be expressed as (see Methods) (12) where z(ω) is given by Eq 8. In order to reproduce that voltage response using the eP model, for which (cf. Eq 9), we obtain (13) As in the previous section, we choose CeP = Cs, GeP = Gs. In contrast to the somatic input filter Ls the filter Ld for distal inputs exhibits low pass properties for various BS morphologies, see Fig 4A. The shape of this filter is largely independent of the soma size. Compared to the attenuation of low frequency in case of somatic input, the filter gain for high frequency dendritic input is much lower. This results in a stronger filtering effect for dendritic inputs than for somatic inputs.

thumbnail
Fig 4. Distal input filter and reproduction of spiking activity using LIF type models.

A: Gain and phase of the input filter as a function of frequency. The neuronal morphology varied as indicated, in terms of dendritic cable length (350 μm, 700 μm, 1050 μm), cable diameter (0.6 μm, 1.2 μm, 1.8 μm) and soma diameter (5 μm, 10 μm, 15 μm). * indicates the default parameter values. B: Coincidence factor for the BS and eP model spike trains, ΓBS,eP (left), and for the BS and P model spike trains, ΓBS,P (right) as a function of input mean and standard deviation σd. C: Difference ΓBS,eP − ΓBS,P between the coincidence factors shown in B. D: Spike rate difference of the BS and eP models (left) and of the BS and P models (right) as a function of and σd. E: Spike rate of the BS neuron model. Results presented in B-E are averages over 6 noise realizations. The default parameters values of the BS model are listed in Table 1.

https://doi.org/10.1371/journal.pcbi.1005206.g004

An evaluation of the distal input filter in terms of reproduction of BS spiking activity (Γ and rates) is shown in Fig 4B–4E for a range of input mean and standard deviation σd values. For comparison we used the P model (without filter) whose parameters were tuned to best reproduce the spike train of the BS model for each input (i.e., (, σd)-pair) separately. The eP model very accurately reproduces the BS spike times for small spike rates (Γ ≥ 0.9 for small and σd). The accuracy drops somewhat as increases, which can be explained as in the previous section. Interestingly, the performance does not deteriorate with increasing spike rate in general; it remains high if the noise intensity σd is sufficiently strong (Γ ≥ 0.8 for σd ≥ 80 pA, independent of in the considered range). The spike rate of the BS model is somewhat underestimated by the eP model (Fig 4D). It should be noted that the spike rate reproduction could be substantially improved by an increased reset voltage value , as the remaining dendritic depolarization after spikes is more pronounced in case of distal input compared to somatic input. The computational speed-up of the eP model here is the same as in the previous section. The P model, in comparison, is less accurate across all inputs (Fig 4B–4D), even though its parameters depend on the input.

In summary, the dendritic cable implements a low pass filter for inputs at the distal dendrite, and due to the corresponding derived filter the eP model reproduces the BS model dynamics for subthreshold and suprathreshold inputs much better than the P model.

Effect of an extracellular electric field on the neuronal dynamics

We now consider an extracellular electric field—in addition to the synaptic inputs—to which the neuron is exposed to. We characterize the effects of that field on the subthreshold somatic membrane voltage and spiking dynamics of the BS neuron and we determine an explicit expression for the additional input current of the extended point neuron model to reproduce these effects. The electric fields we are interested in are oscillatory, spatially uniform on the neuronal scale and weak such as induced by transcranial brain stimulation [6]. In the following, we consider a field with amplitude E1 and angular frequency φ, (14) Recall that VBS,e(x, t) is the extracellular potential. The BS subthreshold somatic membrane voltage response to this field, , is determined by Eqs 13. Using the temporal Fourier transform the solution can be expressed analytically as (15) where z(ω) is given by Eq 8 (see Methods). Note, that we again neglect the exponential current in this section (LIF case, but see next section for the EIF case). In the time domain this yields (16) (17) where arg(x) denotes the argument of the complex number x. The overall subthreshold response in presence of the electric field and synaptic input can be decomposed as (18) with , and given by Eqs 7, 12 and 15. For the eP model, on the other hand, we have (19) To guarantee an equal subthreshold response in both models, i.e., , we obtain the following expression for the additional input current, (20) (21) where we set CeP = Cs and GeP = Gs (as in the previous sections). It should be noted that these results are not restricted to sinusoidal field variations, as considered here, and can be easily adjusted for any time-varying or constant description of the electric field using its Fourier transform.

The equivalent input current IE(t) as well as the somatic subthreshold sensitivity to the field, |A(φ)|/E1 and the phase shift between oscillating membrane voltage and field, arg(A(φ)), with A(φ) from Eq 17, are shown in Fig 5. Interestingly, the amplitude of IE(t) increases with increasing field frequency (Fig 5A), while the sensitivity decreases (Fig 5B). The sensitivity curve changes quantitatively, but not qualitatively, with varying neuronal morphology (Fig 5B). Specifically, its dependence on the field frequency becomes more pronounced with increasing ratio of dendritic size over somatic one. The cable length has the strongest impact in this respect. Notably, the morphology parameters can be adjusted such that the sensitivity curve well matches with empirical results obtained from rat hippocampal pyramidal cells in vitro. The phase shift between the somatic membrane voltage and field oscillations also depends on the field frequency. It exhibits an anti-phase relation for slow oscillations, and decreases with increasing frequency (Fig 5B).

thumbnail
Fig 5. Input current equivalent to the field effect and somatic sensitivity.

A: Input current IE to reproduce the effect of a 1 V/m field in the eP model. Its amplitude and phase shift relative to the field as a function of field frequency. B: Neuron sensitivity to the field, i.e., the ratio between its somatic membrane voltage amplitude and the field amplitude, and phase shift between the oscillatory membrane voltage and the field. The neuronal morphology varied as indicated, in terms of dendritic cable length (350 μm, 700 μm, 1050 μm), diameter (0.6 μm, 1.2 μm, 1.8 μm) and soma diameter (5 μm, 10 μm, 15 μm). ▲ indicate values obtained from electrophysiological recordings of rat hippocampal pyramidal cells [15]. * indicates the default parameter values. For all other parameter values used see Table 1.

https://doi.org/10.1371/journal.pcbi.1005206.g005

We next assess how the electric field affects spiking activity for a range of field frequencies using the BS and eP models. For that purpose, we simulated both model neurons subject to the field and noisy synaptic input at the soma or at the distal dendrite. The synaptic drive alone is strong enough to cause stochastic spiking with rate r0. The oscillatory field leads to an oscillatory spike rate modulation quantified as r1(φ)sin(φt + ψ(φ)) around the constant baseline spike rate r0 (see Methods for details). Note that this spike rate modulation measure is related to the frequently used spike field coherence measure.

The amplitude r1 and phase shift ψ of the spike rate modulation for various somatic inputs (in terms of and σs), a range of field oscillation frequencies and two field strengths are shown in Fig 6. The eP model well reproduces the spike rate dynamics of the BS model exposed to the field for all considered field and input parameter values. The amplitude r1 increases linearly with increasing field magnitude E1. In contrast to the subthreshold sensitivity to the field (cf. Fig 5B), the spike rate modulation exhibits a clear resonance in the beta and gamma frequency bands across the different inputs. In other words, the spike rate oscillations are strongest for field oscillations of that frequency range. The amplitude peak is more pronounced for stronger inputs and most prominent when the input is dominated by its mean (large , small σs). This resonance amplitude rapidly increases with increasing baseline spike rate—by increasing both, mean and standard deviation of the background input from small values—and saturates at about r0 = 30 Hz (Fig 6, center). The resonance frequency shifts rather gradually from the beta to the gamma range as the baseline spike rate increases from a few spikes per second to about 60 Hz. The phase shift ψ varies around π, depending on the input and field frequency. Note that ψ = π implies that the probability of spiking is largest at the trough of the field oscillation. This results from the orientation of the field, which, in case of E(t) = E0 > 0, induces a (hyper-)polarized somatic membrane voltage.

thumbnail
Fig 6. Spike rate modulation due to an electric field for somatic inputs.

Top/bottom, left/right: Spike rate modulation of the BS (blue) and the eP (green) models due to an oscillating electric field as a function of its frequency, for different field amplitudes (E1 = 1 V/m, solid lines; E1 = 10 V/m, dashed lines) and somatic inputs: pA, σs = 11.94 pA (top left), pA, σs = 33.34 pA (top right), pA, σs = 11.94 pA (bottom left), and pA, σs = 33.34 pA (bottom right). Magenta lines show the spike rate modulation of the eP model for which IE was given by IE(t) = I1 sin(φt + ϕ) with constant amplitude I1 = |B(0.5/(2π))|, phase shift ϕ = arg(B(0.5/(2π))), B from Eq 21 and E1 = 10 V/m. Note the different amplitude scales in the top panel. Results for larger field amplitude (E1 = 10 V/m) are not displayed for the mean driven regime (top right), because spike rate modulation amplitudes exceeded the baseline rate in that case, which impedes the modulation quantification procedure (see Methods). Center: Resonance frequency argmax(r1) and amplitude max(r1) of the spike rate modulation of the eP model as a function of baseline spike rate r0, which was changed by simultaneously increasing (, σs) from (4.25 pA, 8.89 pA) to (8.12 pA, 36.40 pA).

https://doi.org/10.1371/journal.pcbi.1005206.g006

To examine the importance of the specific shape of IE(t), we also considered an alternative sinusoidal input current IE(t) = I1 sin(φt + ϕ) for the eP model. Note that the amplitude and phase shift of that current are constant across different field frequencies. Using that current, the typical resonance of the spike rate modulation due to the field cannot even roughly be reproduced (Fig 6).

Let us now inspect spike rate modulation due to the field in presence of distal dendritic inputs instead of somatic ones. In Fig 7 the results are shown for various distal inputs (in terms of and σd). Interestingly, for all considered distal dendritic inputs, spike rate modulation amplitudes increase monotonically with the field frequency over the whole considered range (up to 1 kHz, see Discussion for an explanation). Similarly as for somatic inputs, modulation is strongest for mean dominated (large , small σd) distal inputs, and the phase shift ψ varies around π. Overall, the eP model well reproduces the modulation observed in the BS model.

thumbnail
Fig 7. Spike rate modulation due to an electric field for distal dendritic inputs.

Spike rate modulation of the BS (blue) and the eP (green) models due to an oscillating electric field as a function of its frequency, for different distal dendritic inputs: pA, σd = 33.04 pA (top left), pA, σd = 111.2 pA (top right), pA, σd = 33.04 pA (bottom left), and pA, σd = 111.2 pA (bottom right). Magenta lines show the spike rate modulation of the eP model for which IE was given by IE(t) = I1 sin(φt + ϕ) with constant amplitude I1 = |B(0.5/(2π))| and phase shift ϕ = arg(B(0.5/(2π))) with B from Eq 21 and E1 = 10 V/m. Note the different amplitude scales in the upper panel.

https://doi.org/10.1371/journal.pcbi.1005206.g007

Extension for EIF model neurons

In the previous sections, we considered only capacitive and leak currents through the neuronal membrane; the model extension presented there applies to the LIF type model neurons. Here, we consider the BS and eP models described by Eqs 13, 5 without neglecting the exponential term, that approximates the voltage dependent sodium current at spike initiation. That is, we derive and evaluate the model extension for model neurons of the EIF type.

To derive the required model components Ls(t), Ld(t), α and IE(t) we linearize the exponential terms in Eqs 2 and 5 around a baseline voltage value V0 and then proceed similarly as above. Specifically, we calculate the subthreshold somatic membrane voltage response of the BS model, using the (temporal) Fourier transform, and obtain four response components: , where , and denote the voltage response components to Is, Id and E, respectively, and the additional term is due to the (linearized) exponential term. These four voltage response components are given by the explicit expressions Eqs 4143 in the Methods section. For the eP model, on the other hand, we can also calculate the subthreshold membrane voltage response in the Fourier domain, , given by Eq 48. By requiring equal subthreshold responses, , we obtain the following explicit expressions for the components Ls, Ld, α and IE, considering the electric field defined in Eq 14: (22) (23) (24) (25) (26) where z(ω) is given by Eq 8. The scaling factor α guarantees that the voltage response component caused by the exponential term, , is reproduced. In other words, α ensures that the spike initiation current, described by the exponential term, leads to an equal steady state in both models. Note that the two filters for EIF neurons and those for LIF neurons depend on input frequency in qualitatively the same way (by comparing Eqs 22 and 23 with Eqs 10 and 13).

We assessed the reproduction of BS spiking activity by the extended EIF model for somatic inputs using the spike coincidence factor Γ and estimated spike rates (Fig 8). Here again the parameter values of the P model were adjusted to maximize ΓBS,P for each input separately. The range of input parameter values was chosen to obtain similar spike rates as in Fig 3. Despite the linearization in the derivation, the eP model achieves a correct reproduction of the BS spike trains (Γ ≥ 0.7 for a wide range of input parameters). In particular, ΓBS,eP is large for small spike rates (small and σs) and decreases for increasing (towards mean dominated input), see Fig 8A and 8D. The eP model tends to underestimate the firing rate of the BS model (Fig 8C). This discrepancy in the rate could be reduced by optimizing the point model reset voltage, , to better account for the remaining dendritic cable depolarization in the BS model. Similarly, an improved performance of the eP model in terms of spike train reproduction could be achieved by tuning this reset voltage. The P model, on the other hand, rather poorly reproduces the BS spiking dynamics for small input noise intensity (Γ ≤ 0.6 for σs ≤ 30 pA, see Fig 8A). Overall, also in presence of the exponential term the eP model clearly outperforms the simpler P model for small spike rates (ΓBS,eP − ΓBS,P ≥ 0.3 for small and σs) and achieves similar performance for high spiking activity (Fig 8B).

thumbnail
Fig 8. Reproduction of spiking activity for somatic inputs using EIF type models.

A: Coincidence factor for the BS and eP model spike trains, ΓBS,eP (left), and for the BS and P model spike trains, ΓBS,P (right) as a function of input mean and standard deviation σs. The parameter values of the P model were optimized to maximize ΓBS,P for each input (i.e., (, σs)-pair) independently. B: Difference ΓBS,eP − ΓBS,P between the coincidence factors shown in B. C: Spike rate difference of the BS and eP models (left) and of the BS and P models (right) as a function of and σs. D: Spike rate of the BS neuron model. Results presented in A-D are averages over 6 noise realizations. The parameter values of the BS model are listed in Table 1.

https://doi.org/10.1371/journal.pcbi.1005206.g008

The reproduction of spiking activity of the BS model was also assessed for distal dendritic inputs. The range of input parameters ( and σd) was adjusted to obtain similar BS spike rates as for the LIF case. The eP model performs well, in particular for small spike rates or sufficiently strong noise intensity; its performance decreases in the mean driven regime (Fig 9A). On the contrary the P model fails to reproduce the BS spiking activitiy (see Fig A in S1 Text for more details).

thumbnail
Fig 9. Reproduction of spiking activity for dendritic inputs and spike rate modulation due to an electric field using EIF type models.

A: Coincidence factor for the BS and eP model spike trains, ΓBS,eP (left) as a function of input mean and standard deviation σs. Difference ΓBS,eP − ΓBS,P between the coincidence factors obtained with the eP and the P models (center). The parameter values of the P model were optimized to maximize ΓBS,P for each input (i.e., (, σs)-pair) independently. Spike rate of the BS neuron model (right). Results are averages over 6 noise realizations. The parameter values of the BS model are listed in Table 1. B: Spike rate modulation of the BS (blue) and the eP (green) models due to an oscillating electric field as a function of its frequency, for different distal somatic inputs: pA, σs = 24.08 pA (left), pA, σs = 68.21 pA (right). Magenta lines show the spike rate modulation of the eP model for which IE was given by IE(t) = I1 sin(φt + ϕ) with constant amplitude I1 = |B(0.5/(2π))| and phase shift ϕ = arg(B(0.5/(2π))) with B from Eq 21 and E1 = 10 V/m. C: Same as B for dendritic synaptic input instead of somatic one: pA, σd = 57.73 pA (left), pA, σd = 203.41 pA (right).

https://doi.org/10.1371/journal.pcbi.1005206.g009

In summary, the somatic and distal dendritic input filters obtained for EIF neurons are qualitatively similar to the ones obtained for LIF neurons. The eP model, in contrast to the P model, well reproduces the BS model dynamics for subthreshold and suprathreshold inputs—also for the EIF case.

Spike rate modulations due to an oscillatory electric field using EIF type model neurons for synaptic background input at the soma or distal dendrite are displayed in Fig 9 (see also Fig B and Fig C in S1 Text for additional parameter values of the background input). Similarly to the LIF case, spike rate modulation amplitudes do not decrease monotonically with the field frequency. For somatic background input, we find spike rate resonance in the beta and gamma frequency range, similarly as shown by LIF type models. However, in case of distal dendritic input, EIF neurons exhibit resonance peaks in the high gamma frequency band, in contrast to LIF neurons, whose resonance frequency is substantially larger (see Discussion for an explanation). For both input locations the spike rate modulations shown by the BS model are well reproduced by the eP model and resonance amplitudes are stronger for large spike rates (i.e., large , σs and large , σd, respectively).

Discussion

In this contribution we presented an extension for IF point model neurons to accurately reflect the filtering of synaptic inputs caused by the presence of a dendrite and the effects of weak, oscillatory electric fields on neuronal activity. Based on a canonical BS neuron model, we analytically derived additional components for LIF and EIF point neuron models to exactly reproduce the subthreshold voltage dynamics of the spatially extended BS neuron.

These new components consist of (i) two linear filters applied to synaptic inputs depending on their location (soma or distal dendrite) and (ii) an additional input current quantifying the field effect on the membrane voltage. The EIF point model requires an additional scaling parameter to accurately match the BS voltage dynamics. Exhaustive evaluations for suprathreshold in-vivo like fluctuating inputs demonstrated that the BS spiking activity is well reproduced by the extended point neuron model in both cases (LIF and EIF). Optimizing the parameters of the standard LIF and EIF models without the derived extension components, however, does not suffice to adequately reproduce the BS model dynamics.

Due to their computational efficiency the extensions of the point neuron models are well suited for application in large networks to investigate, for example, the effects of neuronal morphology and electrical fields on neuronal spiking activity at the population level. Additionally, our methodological results serve as a building block to derive mean-field descriptions for the collective (spike rate) dynamics of large coupled populations [17, 21, 22, 23, chapter 4.2]. An implementation of the presented models using Python (for the eP model) and NEURON (for the BS model) is freely available at https://github.com/nigroup/IF_extension.

Below, we summarize our results on the obtained input filters and the field effects on neuronal dynamics.

Synaptic input filtering due to the dendrite

We have demonstrated that synaptic input is integrated at the soma in distinct ways due to the presence of the dendrite, depending on the input site. Distal dendritic input is low-pass filtered (cf. Fig 4A), in accordance with previous results [24], whereas somatic input is high-pass filtered (cf. Fig 3B). The latter effect is consistent with recent measurements from Purkinje cells and with theoretical results [18] which show a similar change in somatic impedance due to the presence of a dendritic tree (Fig. 4 in [18], in comparison with Fig 2A here). Consequently, the presence of a dendrite can lead to an enhanced neuronal spiking response to high-frequency somatic inputs [18], which may be further amplified by the dendritic effect on the sharpness of spikes at the axon initial segment [25]. The derived IF model extension enables efficient analyses of the BS spike rate response to modulations of the input current—which are, however, not within the scope of this paper.

There are two different strategies for taking into account complex neuron morphologies in models while keeping numerical simulation computationally efficient. One option is to reduce the number of compartments while retaining important properties of the dendritic tree [26]. Alternatively, one can extend point neuron models with temporal kernels which are calibrated to reproduce the somatic membrane voltage response to synaptic inputs as observed in complex morphological cells [27, 28]. Our approach is of the latter type, with the advantage that the temporal kernels (filters) are analytically derived from the underlying morphological BS model.

A similar extension for point model neurons to reproduce dendritic input integration of model cells with complex morphology has been recently proposed in [29]. Using the Green’s function formalism a synapse model was developed, whose computational complexity practically allows for only a small number of synaptic input locations. Based on the BS model we were able to derive input filters for point model neurons using only the Fourier transform (without having to rely on the Green’s function) and these filters are simple to implement.

We have demonstrated that our extended model outperforms the simpler point neuron model in terms of spike train reproduction. Overall, it performs well for suprathreshold inputs, particularly in case of distal inputs and for somatic inputs that are not too strong. That performance could be further enhanced by optimizing the reset voltage to better reflect the remaining dendritic membrane depolarization in the BS model after each spike, as was mentioned previously.

In our study we have considered passive dendrites. Nonlinear (spike-generating) currents along the dendrite, which cause nonlinear synaptic input integration [3032], could be incorporated using our approach in a “quasi-active” framework [24]. This would involve solving the cable equation with linearized nonlinear components, similarly as for the exponential terms used here (EIF case).

Effects of weak electric fields on neuronal activity

We investigated in detail the effects of a spatially homogeneous, oscillating, weak electric field, as induced by transcranial electrical stimulation, on the activity of the BS neuron. Such a one-dimensional spatial (cable plus soma) model provides a good approximation for neurons with elongated (apical) dendrites exposed to a uniform extracellular electrical field as long as the dendritic (apical main) cable is not substantially smaller than its electrotonic length [33, chapter 2.5]. Following the somatic doctrine [6], we focused on the effects of the field that are due to the polarization of the membrane voltage at the soma. We analytically calculated the subthreshold voltage response, whose properties are in accordance with electrophysiological observations: the response magnitude scales linearly with the field amplitude [13], as shown by the sensitivity in Fig 5. This sensitivity is of the same order of magnitude as that measured in pyramidal cells [15], i.e., around 0.30 mm for low frequency fields, and decreases with increasing field frequency in a morphology dependent manner [14]. For non-uniform electric fields, e.g., as generated by point source stimulation, however, the sensitivity can be roughly constant for frequencies up to at least 100 Hz [8]. Interestingly, such a behavior can also be observed for a uniform field in case of a rather short dendritic cable (cf. Fig 5B).

While polarization effects due to direct current fields have been extensively studied [3436], the effects of time-varying fields are less well understood. The response of the subthreshold membrane voltage to time-varying fields has been calculated in [37] for a finite dendritic cable with leaky currents at one end, and in [38] for a spatially non-uniform field. Using a one-dimensional cable model [33] showed that the electrotonic length is a key quantity that determines the neuronal subthreshold response to an electric field. Specifically, elongated neurons are less sensitive to high frequency fields than compact ones. How the voltage response to an input current at a particular location along the cable depends on input frequency is largely determined by the membrane time constant. In case of an electrical field, however, which corresponds to symmetrical stimulation at both ends of the cable, the voltage response is also strongly influenced by currents flowing through the low-resistant intracellular medium. This results in an enhanced high frequency response to an extracellular field when compared to an input current [33, chapter 5].

Nevertheless, a somatic compartment was not considered in these studies. Using the BS model we have shown that the relative size of the soma compared to the dendritic cable substantially affects the neuronal sensitivity to the field.

Further, we found frequency-dependent spike rate modulation (and hence, spike field coherence) caused by the electric field. Unlike neuronal subthreshold sensitivity, spike rate modulation amplitude did not decrease with the field frequency and its precise relationship to field frequency depended on the synaptic input location. Spike rate modulation exhibited a clear resonance in the beta and gamma frequency bands in presence of only somatic inputs (cf. Fig 6 and B in S1 Text), whereas for only distal dendritic inputs, spike rate modulation amplitudes are strongest at much larger frequencies (cf. Fig 7 and C in S1 Text). This can be linked to a theoretical result showing that the response of single-compartment model neurons to high frequency inputs is stronger for larger autocorrelation times of a fluctuating synaptic input current [39]. Since fluctuating synaptic inputs arriving at the distal dendrite are low-pass filtered, the autocorrelation time of the corresponding input current felt by the soma is increased (or rather limited from below). Spike rate resonance frequencies were lower for EIF neurons as compared to LIF neurons, in particular for background inputs only at the distal dendrite. This may be explained by the fact that the presence of the exponential term, describing the spike initiating sodium current, decreases the rate response to high frequency inputs [19] (see also the analytical results in [18]). In all cases, the amplitude of the modulation also depended on the input strength (input mean and noise intensity), but its relationship to field frequency was not strongly affected by the input parameters. Recently it has been shown that Purkinje neurons, due to their large dendritic trees, exhibit spike rate resonance at rather high frequencies in response to somatic input modulations and in the presence of noisy dendritic input [18, Fig 5], which is qualitatively similar to the field-induced resonance effects described here (cf. Figs 7 and 9C). It should be noted, however, that an oscillatory (spatially uniform) external field corresponds to oscillatory input currents with opposite sign at the soma and the distal dendrite, respectively (cf. Eqs 2 and 3). The effects of the field can thus not be easily anticipated from those of an input current modulation at the soma alone. Furthermore, the dendritic membrane surface compared to the somatic one for Purkinje cells [18] is substantially larger than that of pyramidal neurons as considered here, which additionally impedes to directly relate the results.

Existing experimental studies on the modulation of neuronal activity by extracellular fields have considered a small number of field frequencies (see [40] for a review). Therefore, our results on spike rate resonance are currently not completely confirmed and may be regarded as predictions. In accordance with our findings weak alternating electric fields (of 30 Hz) have been shown to increase the spiking coherence of pyramidal cells in rat hippocampal slices [41], where this increase was proportional to the subthreshold membrane polarization. Moreover, spatially uniform extracellular fields with high-frequency components entrained spiking activity in ferret primary visual cortex more effectively than fields that only contain low-frequency components [1, Fig. S6]. Our predictions on spike rate modulation by an oscillating electric field are thus in agreement with current knowledge and are informative for future experimental studies. Those results may further be of potential interest for the design of transcranial electrical stimulation protocols.

Regarding the point model extension, we analytically derived an expression for an input current to reproduce the effect of the field as extracted from the biophysically grounded BS model. The amplitude and phase of this input current depend on the parameters of the BS neuron and the electric field. Previously, simple phenomenologically obtained input currents have been used for point neuron network simulations, with either constant amplitudes (across frequencies) [1, 9] or amplitudes fitted to electrophysiological data [7]. Interestingly, the latter study used an input current whose magnitude decreases with increasing frequency, in contrast to the equivalent current we obtained (whose magnitude increases with frequency up to 10 kHz). The neuronal subthreshold sensitivity in that study and the ones shown here, however, are similar. This apparent discrepancy in the currents describing the field effect may be explained by the impedance of the applied model neurons, which naturally influences the equivalent input current. In [7] the model parameters (and thus the impedance) were not fitted to real cells; hence it is unlikely that the model impedance matched with the impedance of the cells from which the current amplitudes were estimated [15]. The successful reproduction of the BS spike rate modulation due to the field by the eP model presented here supports the high-pass properties of the equivalent input current.

In the present study, we derived an extension for point neuron models of the LIF and EIF types. Additional model variables with slow dynamics [42] may also be included in this framework, in order to reflect, for example, effects of slowly deactivating potassium channels that mediate spike rate adaptation and associated characteristic neuronal response properties [43, 44]. In that case, a separation of timescales argument could be used to derive the model extension.

The results we extracted from a canonical spatial neuron model provide insight into the effects of cellular morphology on synaptic input integration and the impact of extracellular electric fields on neuronal activity. In particular, the presented point model extension, which is straightforward to implement and efficient to simulate, shall give rise to comprehensive computational investigations of neuronal population activity entrainment due to transcranial stimulation.

Methods

The ball-and-stick (BS) neuron model

Model derivation.

The BS neuron model consists of a finite passive dendritic cable of length L with a lumped somatic compartment at one extremity x = 0 and a sealed-end at the other. We consider this neuron model exposed to synaptic inputs at the soma Is(t) and the distal dendritic end Id(t) and to an electric field E(t) (see Fig 1). The electric field is spatially uniform at the scale of the neuron, which is a valid assumption for fields induced by transcranial brain stimulation [6]. Assuming a homogeneous, purely ohmic medium (see [45] for the cable equation in a non-ohmic medium), the subthreshold dynamics of the membrane voltage along the dendritic cable are governed by [20] (27) (28) where VBS(x, t) ≔ VBS,i(x, t) − VBS,e(x, t) − Vrest, with intra- and extracellular potentials VBS,i and VBS,e, respectively. cm = cDdπ is the membrane capacitance per unit length, gi = ϱi(Dd/2)2 π is the internal (axial) conductance per unit length and gm = ϱm Ddπ is the membrane conductance per unit length. c is the specific membrane capacitance (in F/m2), ϱi is the specific internal conductance (in S/m), ϱm is the specific membrane conductance (in S/m2) and Dd is the cable diameter. Note, that the rightmost equality in Eq 27 is due to our assumption of a spatially uniform electric field E(x, t) ≡ E(t).

At the proximal end of the dendritic cable, x = 0, we consider a lumped soma, assuming that the somatic diameter Ds is small compared to the cable length L. The corresponding boundary condition is given by [46] (29) and thus (30) where and are the somatic membrane capacitance and leak conductance, respectively. At the distal end of the dendritic cable, x = L, we have [46] (31) due to the synaptic input Id(t), and therefore, (32) The subthreshold voltage dynamics of the BS model are thus determined by Eqs 27, 30 and 32. The spiking mechanism is implemented by the reset condition 4 with refractory period (see Models in the section Results).

Calculation of the subthreshold somatic response.

To analytically calculate the somatic membrane voltage response of the BS model, we consider small variations of the synaptic inputs Is(t), Id(t) and a weak oscillatory electric field E(t). This allows us to linearize the exponential term in Eq 30 around a baseline voltage value V0 to obtain (33) for x = 0. Note that in case of a purely leaky and capacitive neuronal membrane (i.e., without the exponential term, in the limit ΔT → 0) the linearization above is not required and the response calculated below is also exact for larger (subthreshold) synaptic inputs and electric field magnitudes. The linear partial differential Eq 27 together with the boundary conditions 33 and 32 can be solved using separation of variables VBS(x, t) = W(x)U(t) and the temporal Fourier transform (34) where ω = 2πf denotes angular frequency. We obtain the system of differential equations (35) (36) (37) where indicates the (temporal) Fourier transform and δ(.) the Dirac delta function. The solution of the second order linear differential Eq 35 is given by (38) where ±z(ω) are the roots of the characteristic polynomial gi λ2 = gm + cm iω of Eq 35: (39) (same as Eq 8). The coefficients a1(ω) and a2(ω) are calculated by inserting from Eq 38 in Eqs 36 and 37 to finally obtain (40) with (41) (42) and (43) The function sech(x) = cosh(x)−1 refers to the hyperbolic secant. Here, , and denote respectively the voltage response components to Is, Id and E. is the voltage “response” caused by the (linearized) exponential term. Since E(t) = E1 sin(φt), the response to the field can be expressed in the time domain as (44) (45)

The extended point (eP) neuron model

The subthreshold voltage dynamics of the eP model is specified by Eq 5 which is complemented by the reset condition 6 together with a refractory period (see Models in the section Results).

Calculation of the subthreshold response.

We consider again small variations of the synaptic inputs Is(t), Id(t) and the current IE(t) that corresponds to a weak oscillatory electric field E(t). Linearizing the exponential term in Eq 5 around the steady-state somatic voltage value, V0, we obtain (46) Note that here again in case of a purely leaky and capacitive membrane (ΔT → 0) the linearization above is not required and the response calculated below is also exact for larger (subthreshold) inputs. Using the Fourier transform on Eq 46 yields (47) which can be easily solved to obtain (48)

Numerical simulation

Synaptic input and electric field.

To generate spike trains we considered in-vivo like noisy synaptic inputs Is(t), Id(t). Specifically, Ix(t), x ∈ {s, d} was described by an Ornstein-Uhlenbeck process (49) with mean , correlation time τ and standard deviation σx. ξx(t) is a Gaussian white noise process, i.e. with zero mean 〈ξx(t)〉 = 0 and delta autocorrelation 〈ξx(t)ξx(t+t′)〉 = δ(t′), where 〈.〉 denotes the expectation operator. Eq 49 was numerically solved using the method described in [53].

The electric field was described by (50) (same as Eq 14) with amplitude E1 and angular frequency φ = 2πf. The values for all parameters are provided in Table 1.

Ball-and-stick neuron model.

The BS neuron model was numerically solved using the NEURON simulation environment [54]. We applied a finite difference space discretization scheme with 50 segments for the dendritic cable and the implicit (or backward) Euler time discretization scheme. The integration time step was fixed to 0.05 ms when the exponential term was omitted (ΔT → 0) and 0.025 ms otherwise. Decreasing the time step size and increasing the number of segments did not lead to noticeable changes in the membrane voltage time series. The time-varying extracellular potential was included using the built-in “extracellular” mechanism in NEURON.

Point neuron model.

The point neuron models were numerically solved using the forward Euler time discretization scheme and the same time step as used for the BS model. The solution method was implemented in Python using the library “Numba” for fast computation. Both point model neurons received the same realization of synaptic input Ix(t), x ∈ {s, d} as the BS model neuron. The linear filters Lx(t), x ∈ {s, d} in the eP model were implemented using the fast Fourier transform (FFT) of Ix(t) and the inverse FFT of . The membrane capacitance and conductance of the eP model were chosen as equal to the corresponding somatic quantities of the BS model, CeP = Cs, GeP = Gs. Except stated otherwise, the corresponding parameter values of the P model were determined by fitting the BS model spiking activity in terms of spike coincidences. Specifically, GP was chosen such that the steady-state (somatic) membrane voltage of the BS model is matched exactly, i.e. , x ∈ {s, d}, with impedance , m ∈ {BS|x = 0, P} using Eqs 7, 12 and 11. The value for CP was then selected to maximize the coincidence factor Γ (defined below) between 52 s lasting spike trains of the BS and P model neurons for each shown pair (, σx) of input parameter values. In presence of the exponential term in the models (ΔT > 0) we used V0 = Vr, which is the steady-state (somatic) voltage in the absence of synaptic input. Several other values, within the range [Vr, Vs], were tested for V0. This did not lead to a substantial improvement of the reproduction performance.

Analysis methods for the spike trains

Spike coincidence measure.

To quantify the similarity between the spike trains of the different model neurons we used the coincidence factor Γ defined by [55] (51) where Ncoinc is the number of coincident spikes with precision (i.e., maximal temporal separation) Δ, Nref and Ncomp are the number of spikes in the reference spike train and in the one being compared to it, respectively. 〈Ncoinc〉 = 2rΔNref is the expected number of coincidences generated by a homogeneous Poisson process with the same spike rate r = Ncomp/T as that shown by the compared spike train, where T is the spike train duration. The factor normalizes Γref,comp to a maximum value of 1 which is reached if the spike trains match optimally (with precision Δ). Γref,comp = 0 on the other hand would result from a homogeneous Poisson process with the same rate as for the reference spike train and thus indicates pure chance.

Spike rate resonance and phase shift measure.

To examine and compare their suprathreshold responses to an oscillatory electric field E(t) (Eq 50), we simulated the neuron models subject to both a field and a noisy synaptic background current. The latter was located either at the soma or at the distal dendrite Ix(t), x ∈ {s, d} (Eq 49). We considered regimes (in terms of (, σx)-pairs) where the synaptic drive is sufficiently strong to cause the neuron to spike stochastically with rate r0. The sinusoidal field then causes a modulation of the spike rate that becomes apparent over many trials (i.e., independent realizations of Ix(t)). This quantity can also be thought of as the spike rate averaged over a population of neurons which individually receive a noisy synaptic drive but collectively respond to the same oscillatory field. This population, or trial-averaged, instantaneous spike rate can be expressed as (52) where r1 and ψ denote respectively the amplitude and phase shift, both depending on the angular frequency φ of the field. Note that in the eP model the field effect is described by the oscillatory current IE(t). To estimate the spike rate modulation at a given field frequency we first extracted and collected the field phase ϕs ∈ [0, 2π) for each spike time ts, such that E(ts) = E0 sin(ϕs). These phases were calculated from 944 trials of at least 26 s duration each, for which the first 2 s were disregarded to avoid transients, and only complete field cycles were considered. We then computed a spike rate histogram from the set {ϕs} using 20 equally sized bins and finally applied a sinusoid of the form F(ϕ) = r0 + r1 sin(ϕ + ψ) with ϕ ∈ [0, 2π) to fit that histogram using the method of least squares, where r0 was given by the histogram mean value. In this way we obtained r1 and ψ.

Supporting Information

S1 Text. Supplementary Figures.

Fig A: Reproduction of spiking activity for dendritic inputs using EIF type models. Fig B: Spike rate modulation due to an electric field for somatic inputs using neuron models of the EIF type. Fig C: Spike rate modulation due to an electric field for distal dendritic inputs using neuron models of the EIF type.

https://doi.org/10.1371/journal.pcbi.1005206.s001

(PDF)

Author Contributions

  1. Conceptualization: JL FA.
  2. Funding acquisition: JL KO.
  3. Investigation: FA JL.
  4. Methodology: FA JL.
  5. Project administration: JL KO.
  6. Software: FA JL.
  7. Supervision: JL KO.
  8. Validation: FA JL.
  9. Visualization: FA JL.
  10. Writing – original draft: FA JL KO.
  11. Writing – review & editing: FA JL KO.

References

  1. 1. Fröhlich F, McCormick DA. Endogenous electric fields may guide neocortical network activity. Neuron. 2010;67(1):129–143. pmid:20624597
  2. 2. Buzsáki G, Anastassiou CA, Koch C. The origin of extracellular fields and currents—EEG, ECoG, LFP and spikes. Nat Rev Neurosci. 2012;13(6):407–420. pmid:22595786
  3. 3. Einevoll GT, Kayser C, Logothetis NK, Panzeri S. Modelling and analysis of local field potentials for studying the function of cortical circuits. Nat Rev Neurosci. 2013;14(11):770–785. pmid:24135696
  4. 4. Neuling T, Wagner S, Wolters CH, Zaehle T, Herrmann CS. Finite-element model predicts current density distribution for clinical applications of tDCS and tACS. Front Psychiatry. 2012;3:83. pmid:23015792
  5. 5. Datta A, Bansal V, Diaz J, Patel J, Reato D, Bikson M. Gyri-precise head model of transcranial direct current stimulation: improved spatial focality using a ring electrode versus conventional rectangular pad. Brain Stimul. 2009;2(4):201–207. pmid:20648973
  6. 6. Bikson M, Reato D, Rahman A. Cellular and network effects of transcranial direct current stimulation. In: Miniussi C, Paulus W, Rossini PM, editors. Transcranial Brain Stimul. CRC Press; 2012. p. 55–92. https://doi.org/10.1201/b14174-5
  7. 7. Reato D, Rahman A, Bikson M, Parra LC. Low-intensity electrical stimulation affects network dynamics by modulating population rate and spike timing. J Neurosci. 2010;30(45):15067–15079. pmid:21068312
  8. 8. Anastassiou CA, Perin R, Markram H, Koch C. Ephaptic coupling of cortical neurons. Nat Neurosci. 2011;14(2):217–223. pmid:21240273
  9. 9. Ali MM, Sellers KK, Fröhlich F. Transcranial alternating current stimulation modulates large-scale cortical network activity by network resonance. J Neurosci. 2013;33(27):11262–11275. pmid:23825429
  10. 10. Marshall L, Helgadóttir H, Mölle M, Born J. Boosting slow oscillations during sleep potentiates memory. Nature. 2006;444(7119):610–613. pmid:17086200
  11. 11. Berenyi A, Belluscio M, Mao D, Buzsaki G. Closed-loop control of epilepsy by transcranial electrical stimulation. Science. 2012;337(6095):735–737. pmid:22879515
  12. 12. Herrmann CS, Rach S, Neuling T, Strüber D. Transcranial alternating current stimulation: a review of the underlying mechanisms and modulation of cognitive processes. Front Hum Neurosci. 2013;7:279. pmid:23785325
  13. 13. Bikson M, Inoue M, Akiyama H, Deans JK, Fox JE, Miyakawa H, et al. Effects of uniform extracellular DC electric fields on excitability in rat hippocampal slices in vitro. J Physiol. 2004;557(Pt 1):175–90. pmid:14978199
  14. 14. Radman T, Ramos RL, Brumberg JC, Bikson M. Role of cortical cell type and morphology in subthreshold and suprathreshold uniform electric field stimulation in vitro. Brain Stimul. 2009;2(4):215–228. pmid:20161507
  15. 15. Deans JK, Powell AD, Jefferys JGR. Sensitivity of coherent oscillations in rat hippocampus to AC electric fields. J Physiol. 2007;583(Pt 2):555–565. pmid:17599962
  16. 16. Tiganj Z, Chevallier S, Monacelli E. Influence of extracellular oscillations on neural communication: a computational perspective. Front Comput Neurosci. 2014;8:9. pmid:24570661
  17. 17. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J Comput Neurosci. 2000;8(3):183–208. pmid:10809012
  18. 18. Ostojic S, Szapiro G, Schwartz E, Barbour B, Brunel N, Hakim V. Neuronal morphology generates high-frequency firing resonance. J Neurosci. 2015;35(18):7056–7068. pmid:25948257
  19. 19. Fourcaud-Trocmé N, Hansel D, van Vreeswijk C, Brunel N. How spike generation mechanisms determine the neuronal response to fluctuating inputs. J Neurosci. 2003;23(37):11628–11640. pmid:14684865
  20. 20. Rattay F. Analysis of models for external stimulation of axons. IEEE Trans Biomed Eng. 1986;BME-33(10):974–977. pmid:3770787
  21. 21. Richardson MJE. Dynamics of populations and networks of neurons with voltage-activated and calcium-activated currents. Phys Rev E. 2009;80(2):1–16. pmid:19792172
  22. 22. Augustin M, Ladenbauer J, Obermayer K. How adaptation shapes spike rate oscillations in recurrent neuronal networks. Front Comput Neurosci. 2013;7:1–11. pmid:23450654
  23. 23. Ladenbauer J. The collective dynamics of adaptive neurons: insights from single cell and network models. PhD thesis, Technische Universität Berlin; 2015.
  24. 24. Koch C. Biophysics of computation: information processing in single neurons. Oxford University Press; 2004.
  25. 25. Eyal G, Mansvelder HD, de Kock CPJ, Segev I. Dendrites impact the encoding capabilities of the axon. J Neurosci. 2014;34(24):8063–71. pmid:24920612
  26. 26. Pinsky PF, Rinzel J. Intrinsic and network rhythmogenesis in a reduced Traub model for CA3 neurons. J Comput Neurosci. 1994;1(1–2):39–60. pmid:8792224
  27. 27. Dayan P, Abbott LF. Theoretical neuroscience. vol. 806. Cambridge, MA: MIT Press; 2001.
  28. 28. Jolivet R, Lewis TJ, Gerstner W. Generalized integrate-and-fire models of neuronal activity approximate spike trains of a detailed model to a high degree of accuracy. J Neurophysiol. 2004;92(2):959–976. pmid:15277599
  29. 29. Wybo WAM, Stiefel KM, Torben-Nielsen B. The Green’s function formalism as a bridge between single- and multi-compartmental modeling. Biol Cybern. 2013;107(6):685–694. pmid:24037222
  30. 30. Migliore M, Shepherd GM. Emerging rules for the distributions of active dendritic conductances. Nat Rev Neurosci. 2002;3(5):362–370. pmid:11988775
  31. 31. Zhou D, Li S, hui Zhang X, Cai D. Phenomenological incorporation of nonlinear dendritic integration using integrate-and-fire neuronal frameworks. PLOS One. 2013;8(1):e53508. pmid:23308241
  32. 32. Zhang D, Li Y, Rasch MJ, Wu S. Nonlinear multiplicative dendritic integration in neuron and network models. Front Comput Neurosci. 2013;7:56. pmid:23658543
  33. 33. Malik NA. Frequency-dependent response of neurons to oscillating electric fields. PhD thesis, University of Warwick; 2011.
  34. 34. Cartee LA, Plonsey R. The transient subthreshold response of spherical and cylindrical cell models to extracellular stimulation. IEEE Trans Biomed Eng. 1992;39(1):76–85. pmid:1572684
  35. 35. Plonsey R, Altman KW. Electrical stimulation of excitable cells-a model approach. P IEEE. 1988;76(9):1122–1129.
  36. 36. Tranchina D, Nicholson C. A model for the polarization of neurons by extrinsically applied electric fields. Biophys J. 1986;50(6):1139–56. pmid:3801574
  37. 37. Monai H, Omori T, Okada M, Inoue M, Miyakawa H, Aonishi T. An analytic solution of the cable equation predicts frequency preference of a passive shunt-end cylindrical cable in response to extracellular oscillating electric fields. Biophys J. 2010;98(4):524–533. pmid:20159148
  38. 38. Anastassiou CA, Montgomery SM, Barahona M, Buzsaki G, Koch C. The effect of spatially inhomogeneous extracellular electric fields on neurons. J Neurosci. 2010;30(5):1925–1936. pmid:20130201
  39. 39. Brunel N, Chance FS, Fourcaud N, Abbott LF. Effects of synaptic noise and filtering on the frequency response of spiking neurons. Phys Rev Lett. 2001;86(10):2186–2189. pmid:11289886
  40. 40. Reato D, Rahman A, Bikson M, Parra LC. Effects of weak transcranial alternating current stimulation on brain activity—a review of known mechanisms from animal studies. Frontiers in human neuroscience. 2013;7:687. pmid:24167483
  41. 41. Radman T, Su Y, An JH, Parra LC, Bikson M. Spike timing amplifies the effect of electric fields on neurons: implications for endogenous field effects. J Neurosci. 2007;27(11):3030–3036. pmid:17360926
  42. 42. Brette R, Gerstner W. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J Neurophysiol. 2005;94(5):3637–3642. pmid:16014787
  43. 43. Ladenbauer J, Augustin M, Shiau L, Obermayer K. Impact of adaptation currents on synchronization of coupled exponential integrate-and-fire neurons. PLOS Comput Biol. 2012;8(4):e1002478. pmid:22511861
  44. 44. Ladenbauer J, Augustin M, Obermayer K. How adaptation currents change threshold, gain and variability of neuronal spiking. J Neurophysiol. 2014;111(5):939–953. pmid:24174646
  45. 45. Bédard C, Destexhe A. Generalized cable theory for neurons in complex and heterogeneous media. Phys Rev E. 2013;88(2):022709. pmid:24032866
  46. 46. Tuckwell HC. Introduction to theoretical neurobiology. vol. 1. Cambridge University Press; 1988. https://doi.org/10.1017/CBO9780511623271
  47. 47. Migliore M, Ferrante M, Ascoli GA, Signal GAA. Signal propagation in oblique dendrites of CA1 pyramidal cells. J Neurophysiol. 2005;94(6):4145–4155. pmid:16293591
  48. 48. Mainen ZF, Sejnowski TJ. Influence of dendritic structure on firing pattern in model neocortical neurons. Nature. 1996;382(6589):363–366. pmid:8684467
  49. 49. Rattay F. The basic mechanism for the electrical stimulation of the nervous system. Neuroscience. 1999;89(2):335–346. pmid:10077317
  50. 50. Major G, Larkman AU, Jonas P, Sakmann B, Jack J. Detailed passive cable models of whole-cell recorded CA3 pyramidal neurons in rat hippocampal slices. J Neurosci. 1994;14(8):4613–4638. pmid:8046439
  51. 51. Spruston N. Pyramidal neuron. Scholarpedia. 2009;4(5):6130.
  52. 52. Badel L, Lefort S, Brette R, Petersen CCH, Gerstner W, Richardson MJE. Dynamic I-V curves are reliable predictors of naturalistic pyramidal-neuron voltage traces. J Neurophysiol. 2008;99(2):656–666. pmid:18057107
  53. 53. Gillespie D. Exact numerical simulation of the Ornstein-Uhlenbeck process and its integral. Phys Rev E. 1996;54(2):2084–2091. pmid:9965289
  54. 54. Hines ML, Carnevale NT. NEURON: A tool for neuroscientists. Neurosci. 2001;7(2):123–135.
  55. 55. Gerstner W, Kistler WM. Spiking neuron models: Single neurons, populations, plasticity. Cambridge University Press; 2002. https://doi.org/10.1017/CBO9780511815706