Weak electric fields promote resonance in neuronal spiking Analytical results from two-compartment cell and network models

Authors: Josef Ladenbauer, Klaus Obermayer

Abstract

Transcranial brain stimulation and evidence of ephaptic coupling have sparked strong interests in understanding the effects of weak electric fields on the dynamics of neuronal populations. While their influence on the subthreshold membrane voltage can be biophysically well explained using spatially extended neuron models, mechanistic analyses of neuronal spiking and network activity have remained a methodological challenge. More generally, this challenge applies to phenomena for which single-compartment (point) neuron models are oversimplified. Here we employ a pyramidal neuron model that comprises two compartments, allowing to distinguish basal-somatic from apical dendritic inputs and accounting for an extracellular field in a biophysically minimalistic way. Using an analytical approach we fit its parameters to reproduce the response properties of a canonical, spatial model neuron and dissect the stochastic spiking dynamics of single cells and large networks. We show that oscillatory weak fields effectively mimic anti-correlated inputs at the soma and dendrite and strongly modulate neuronal spiking activity in a rather narrow frequency band. This effect carries over to coupled populations of pyramidal cells and inhibitory interneurons, boosting network-induced resonance in the beta and gamma frequency bands. Our work contributes a useful theoretical framework for mechanistic analyses of population dynamics going beyond point neuron models, and provides insights on modulation effects of extracellular fields due to the morphology of pyramidal cells.

Introduction

The interaction between weak electric fields and neuronal activity in the brain has gained increased attention over the past decade [1–3]. These weak fields can be generated endogenously by populations of neurons [4–7] or through transcranial electrical stimulation [3, 8–10], and they can modify neural activity in various ways [4, 11–15]. Although the electric fields caused by this type of noninvasive intervention exhibit low magnitudes (≤ 1-2 V/m [8, 9, 16]) they can modulate neuronal spiking activity [3, 4, 17, 18] and lead to changes in cognitive processing [19, 20], offering a number of possible clinical benefits [21, 22]. The influence of extracellular fields on the subthreshold membrane voltage of single cells has been thoroughly studied and biophysically explained [17, 23–26]. How weak electric fields affect neuronal spiking activity and interact with network dynamics, however, is currently not well understood [3].

Multi-compartment models are useful tools to dissect effects in single neurons in the absence of input fluctuations [27–29], but they are not well suited to study neuronal and network spiking activity in noisy, in-vivo like conditions, because of their large complexity. Single-compartment (point) neuron models, on the other hand, allow for mechanistic analyses at the network level (see, e.g. [30–32]), but they lack the spatial structure that is required to biophysically describe the direct effect of an electric field on the membrane voltage [25]. Extracellular fields can only be effectively incorporated either in simplified, phenomenological ways [12, 14, 33] or by elaborate model extensions to accurately reproduce the effects on spatially elongated neurons [25] at the cost of an increased complexity of the system (and its analysis). Two-compartment neuron models feature a suitable compromise between biological rigor and analytical tractability in this regard, with minimal level of spatial detail necessary to biophysically take into account an extracellular electric field. Models of this type have proven useful to study the effects of constant fields on the activity of single neurons and synchronization of neuronal pairs [34–36]. Irrespective of extracellular fields simplified two-compartment neuron models have also been successfully applied, for example, to dissect neural coding of high-frequency signals [37, 38] or neurocircuit mechanisms [39]. However, powerful methods to study the spiking dynamics of these model neurons, at the levels of single cells and populations in a noisy, cortical setting, have been lacking.

Here we employ this model class to describe the activity of pyramidal (PY) neurons receiving in-vivo like fluctuating inputs in the presence and absence of an oscillatory extracellular field. The neuron model consists of a compartment for the soma and one for the apical dendrite, allowing to differentiate between synaptic inputs at the soma (including basal dendrite) and those at the distal (apical) dendrite. We first determine the model parameters semi-analytically to reproduce the behavior of a more sophisticated, spatially extended neuron model which involves the cable equation. We then develop an analytical method to effectively characterize the stochastic spiking dynamics of single neurons and sparsely coupled large populations of PY neurons and inhibitory (IN) interneurons. The method exploits the Fokker-Planck equation and a moment closure approximation technique.

We apply these tools to study (i) how a weak oscillatory field affects the spiking activity of neurons exposed to fluctuating synaptic inputs, (ii) how these effects compare to those of weak oscillatory inputs in the absence of an electric field, and (iii) how weak applied fields modulate network-induced oscillations. Our contribution sheds some light on the effects of extracellular fields at the population level. Furthermore, it provides useful methods for mechanistic studies on the dynamics of coupled compartmentalized spiking neurons that allow to broadly distinguish inputs according to the location of the synapses. Such a distinction is important, for example, in circuit models which involve different types of (inhibitory) neurons.

Results

Modeling approach

Pyramidal neuron model

The PY model neurons consist of two compartments, one for the soma and one for the (apical) dendrite, for which we consider trans-membrane capacitive currents, ionic leak currents, an approximation of the somatic Na^+^ current at spike initiation, an internal current, synaptic input currents and an extracellular electric field. The latter is defined as E(t)≔[Vse(t)-Vde(t)]/Δ, where Vse, Vde are the extracellular membrane potentials for the somatic and dendritic compartments, whose centers are separated with distance Δ. The dynamics of action potentials are simplified by a reset mechanism of the integrate-and-fire type at the soma. Fluctuating input currents at the soma and the dendrite, Is and Id, that mimic the compound effect of synaptic bombardment in vivo, are described by stochastic processes with means I¯s, I¯d and standard deviations σs, σd. In this model the dynamics of the somatic and the dendritic membrane voltage, Vs and Vd, respectively, are thus governed by two coupled differential equations together with a reset condition for spikes (for details see Methods section 1). A schematic circuit diagram of the model is depicted in Fig 1A.

Model parametrization

To determine the parameter values of the two-compartment model and assess whether its spatial complexity is adequate we semi-analytically fit the model to a biophysically more sophisticated, spatially elongated ball-and-stick neuron model (Fig 1B). That model involves the cable equation, an integrate-and-fire spike mechanism at the soma and an extracellular field which is assumed spatially homogeneous (see Methods section 4). Due to its mathematical complexity this model is not well suited for analyses of spiking dynamics and application in networks. The fitting method approximates the somatic responses of the ball-and-stick model in an efficient way without knowledge about the input or the electric field. In brief, for both models we analytically calculate the somatic subthreshold responses for small amplitude variations of the inputs and the electric field as well as the transient somatic voltage responses after a spike for threshold inputs, and subsequently apply a least-squares fit (Fig 1C–1E). Note that the membrane voltage dynamics at the soma directly affect spike timing and are therefore the most relevant (axon initial segments are absorbed by the somatic compartments and not separately included in the models). In this way we rapidly obtain an accurate reproduction of the relevant response properties of the ball-and-stick model (Fig 1C–1H). Specifically, the subthreshold responses, in terms of somatic impedances for inputs at the soma and dendrite (Fig 1C) and the somatic voltage response for an oscillatory applied field (Fig 1D), are well approximated, and spiking activity, in terms of spike coincidences (Fig 1F and 1G, see Methods section 5) and the spike rate (Fig 1H), are accurately predicted. Note that the accuracy in Fig 1F–1H results from the approximation quality for subthreshold responses (Fig 1C) and post-spike transients (Fig 1E). This quality of our fitting method is not restricted to the selected ball-and-stick morphology but extends to a range of plausible shapes (S1 Fig).

Characterization of spiking dynamics

We focus on spiking activity, in particular the dynamics of the (instantaneous) spike rate r. This quantity can be calculated exactly using the Fokker-Planck equation that governs the evolution of the joint probability density for the somatic and dendritic membrane voltage p(Vs, Vd, t), describing the stochastic dynamics in deterministic form (see Methods section 3). Since a numerical solution of this partial differential equation is very demanding in terms of implementation and especially computational effort we employ a moment closure approximation method. Specifically, we use p(Vs, Vd, t) = ps(Vs, t)pd(Vd|Vs, t), where ps is the marginal probability density for the somatic voltage and pd the probability density for the dendritic voltage conditioned on Vs, and approximate pd by a conditioned Gaussian probability density. The resulting system allows for convenient and efficient calculation of the spike rate responses for constant input statistics as well as weak sinusoidal variations of the input moments or the applied field. To this end only ordinary differential equations need to be solved. An evaluation of this method in comparison to numerical simulation for constant input moments is shown in Fig 1H. Spike rates from simulations are well matched for a range of plausible input moments. Moreover, this quality of agreement extends to different parametrizations of the two-compartment model that correspond to various ball-and-stick morphologies (S2A Fig).

Modulation of neuronal spiking activity

We first consider a PY neuron exposed to fluctuating inputs at the soma and the dendrite and a weak sinusoidal applied field. The noisy inputs drive the neuron to stochastic spiking activity that is influenced by the field. This effect can be quantified by the instantaneous spike rate across a large number of trials (obtained from numerical simulations with different realizations of the input), which is equivalent to the population-averaged spike rate for a large number of uncoupled PY neurons receiving independent inputs. The field leads to an oscillatory modulation of the spike rate that is accurately reproduced by our analytical calculation method (Fig 2A).

By measuring these spike rate responses over a range of field frequencies we observe a clear resonance in a biologically relevant, relatively narrow frequency band (Fig 2B). In other words, the spike rate oscillations are strongest for field oscillations in that frequency range. This resonance behavior is not shown by the subthreshold membrane voltage response to the field in the absence of suprathreshold fluctuating inputs. Interestingly, spike rate responses to weak sinusoidal modulations of the mean input at the soma or dendrite instead of an applied field do not exhibit such a resonance (Fig 2C and 2D). Response amplitudes typically decrease as the frequency of the input modulations increases, although spike rate responses remain elevated for somatic modulations of up to about 100 Hz. This behavior is robustly shown for various neuronal morphologies (S2B Fig). The overall response amplitudes, especially for input modulations at the dendrite, however, are strongly affected by the ratio between dendritic and somatic membrane surface. An increased ratio appears to cause increased spike rate responses. Notably, the analytical results are in good agreement with those of numerical simulations (Fig 2 and S2 Fig). This justifies the extensive application of our analytical method in the subsequent analyses.

Next, we assess the resonance behavior caused by the applied field for a range of biologically plausible input statistics (see Fig 3). Intriguingly, resonance robustly appears across these input conditions; that is, a resonance peak occurs for very different spiking statistics, including low rate, high variability (due to small input means I¯s, I¯d and large variances σs2, σd2) as well as high rate, low variability (due to large I¯s, I¯d and small σs2, σd2; cf. Figs 3 and 1H). Resonance frequency (Fig 3A) and strength (Fig 3B and 3C) increase with increasing input mean, but not necessarily with increasing input variance. For mean-dominated input (that is, large input mean and small input variance) the resonance frequency is similar to the baseline spike rate (compare the dashed curves for large input mean in Fig 3A with Fig 1H). In this case an increase in input variance interestingly leads to a decrease in resonance frequency and strength. For fluctuation-dominated input the resonance frequency is restricted to the the beta and low gamma frequency bands (∼15-40 Hz), and peak amplitudes vary between about 1-2 spikes/s.

How does the external field promote this resonance behavior? From the definition of the extracellular field and the circuit diagram in Fig 1A it becomes evident that the trans-membrane currents caused by the field at the soma and dendrite are opposed (using Kirchhoff’s law, that all incoming currents at a point of the circuit must sum to zero; see Eqs (1) and (2), where E(t) affects Vs and Vd with opposite sign). This indicates that the electric field effectively reflects anti-correlated inputs at the soma and dendrite. To examine the influence of input correlations more closely we applied sinusoidal modulations of the mean input at the soma and dendrite with a phase shift (Fig 4A). A resonance peak emerges as the phase shift increases from 0 (synchronized modulations / strong correlation) towards 180° (anti-synchronized modulations / negative correlation) and becomes most pronounced at that value. We further considered a time lag instead of a phase shift between the two mean input modulations (Fig 4B). The lag effectively reflects a difference in delays with which an oscillatory signal arrives at the two neuronal locations, for example, due to a relay population. Lags that are sufficiently large cause multiple resonance peaks. The frequency that corresponds to the most dominant peak (with positive frequency) decreases with increasing time lag, but is otherwise largely independent of the baseline input statistics (Fig 4C).

Modulation of network dynamics

To examine how an applied field interacts with network mechanisms to shape population dynamics we derived a mean-field network model from large populations of sparsely coupled PY neurons and IN neurons exposed to fluctuating background inputs and a spatially homogeneous (with respect to neuronal orientation) weak electric field. PY neurons are described by the two-compartment model and IN neurons by an established single-compartment spiking model (exponential integrate-and-fire), because IN neurons are spatially more compact and therefore the direct effect of the electric field on them is negligible [24]. Synaptic coupling is incorporated by delayed current pulses which cause post-synaptic potentials of reasonably small magnitude with distributed delays. In the derivation we apply a diffusion approximation and a moment closure method to transform the original network model into a manageable system of two Fokker-Planck equations—one for each population—coupled via synaptic input moments. This system describes the collective stochastic dynamics in a way that allows for convenient dissection of the network activity in terms of population-averaged instantaneous spike rates (see Methods section 6).

The network is parametrized to exhibit a baseline state of asynchronous (irregular) spiking activity that depends on the strength of the external drive. Network-induced resonance emerges for weak sinusoidal input modulations at the soma of PY neurons with strongest responses in the beta and gamma frequency bands (∼15-55 Hz, Fig 5). Depending on whether IN neurons target PY neurons only at the soma (Fig 5B, 5E and 5H) or only at the dendrite (Fig 5C, 5F and 5I) the resonance peak occurs at a higher or lower frequency compared to the mixed setup (Fig 5A, 5D, and 5G). This resonance behavior is more pronounced in a baseline state of increased activity (Fig 5A–5F) where the activity level of PY neurons seems more decisive than that of IN neurons (Fig 5G–5I). Interestingly, for input modulations at the dendrite of PY neurons such a clear resonance behavior does not appear.

When we consider a weak sinusoidal field instead of input modulations we observe a strongly amplified resonance of the population activity in the same frequency band. The prominent effect of an applied field in single cells thus carries over to PY-IN networks, boosting network-based oscillations that are mediated by an excitation-inhibition loop (as studied, for example, in [40]). Notably, the analytical results from the mean-field network are quantitatively validated by numerical simulations of networks with 10,000 neurons. The agreement is quite remarkable considering that the mean-field derivation requires additional approximations compared to the single cell setting. In sum these results indicate that an oscillatory electric field as generated, for example, by transcranial stimulation should be most effective to induce or maintain rhythmic network dynamics in the beta-gamma frequency range.

Discussion

Methodological aspects

Neuron model

Model-based investigations on how extracellular electric fields impact neuronal spiking activity face a methodological the neuronal spatial extent plays an important biophysical role, however, morphologically detailed neuron models are exceedingly complex for use in noisy, in-vivo like conditions, especially in large networks. Two-compartment neuron models feature the minimal level of spatial detail required to biophysically account for the direct polarization effects caused by the field. At the same, simplified models of this type retain the computational efficiency and, in part, the analytical tractability of point neurons that are widely used to study network dynamics. Our approach is based on a two-compartment spiking PY neuron model which accounts for an external electric field and includes fluctuating synaptic input at the soma and/or (apical) dendrite. Using semi-analytical techniques this model was rapidly calibrated to accurately approximate the (subthreshold somatic and spiking) response properties of a morphologically more plausible ball-and-stick model neuron.

Based on that spatially elongated model we have recently developed an extension for point neuron models to match the subthreshold somatic responses to synaptic inputs and an applied weak electric field [25]. The methodology presented here entails several advantages compared to our previous the fitted two-compartment model involves a low-dimensional differential equation in contrast to an integro-differential equation (in case of the extended point model to account for dendritic filtering effects), allowing for faster numerical simulation, straightforward implementation of networks, and a natural way to account for an extracellular field. Notably, the last feature highlights a benefit over networks of point model neurons for which field effects were previously implemented in a simplified, phenomenological way [12, 14, 33]. Importantly, the two-compartment model further allowed for the application of analytical methods to study spiking dynamics.

Spike rate dynamics

We devised a method to efficiently study the spike rate dynamics of these two-compartment neurons and the population activity of sparsely coupled large networks. It employs the Fokker-Planck partial differential equation and a moment closure technique for dimension reduction which allows to numerically solve the resulting system in reasonable time. These solutions yielded an accurate and efficient approximation of the instantaneous (population-averaged) spike rate (cf. Figs 1H, 2 and 5 and S2 Fig).

Similar moment closure methods were previously applied in different contexts, such as integrate-and-fire point model neurons with synaptic dynamics (see, e.g. [41, 42] and references therein). In Ref. [41] the issue was raised that moment closure may not be applicable for certain ranges of parameter values. For the setting and regions of parameter space considered throughout the present study the approximation method was well suited. This might not be the case for very small or large values of the input parameters (causing unphysiological behavior).

We approximated the conditional probability density pd(Vd|Vs, t) by a conditioned Gaussian, thereby closing the system of lower dimensional equations at the 3^rd^ central moment of Vd (assuming the 3^rd^ and higher cumulants are zero). Closure at lower order moments lead to markedly less accurate results due to similar timescales for the somatic and dendritic voltage and strong coupling between those variables. On the other hand, the room for improvement is vanishingly small such that closure at higher order moments does not justify its increased implementation complexity and computational demands.

Notably, moment closure at order 1 is equivalent to the adiabatic approximation frequently applied for adaptive integrate-and-fire neurons in the presence of noise [32, 43–45]. In that case the conditional first centered moment (of the adaptation variable given membrane voltage) is approximated by the corresponding unconditioned moment. This usually leads to an accurate reproduction of the spike rate dynamics due to the circumstance that the timescales of the two variables are separated. Such an adiabatic approximation has also been applied for two-compartment Purkinje model neurons [37] which possess large dendritic trees. A substantial difference in somatic and overall dendritic capacitance justifies the assumption of separated timescales for those models. Our approach, in contrast, is also valid for model parametrizations that yield rather similar timescales for soma and dendrite, and therefore suitable for pyramidal cells.

It is worth noting that our analytical techniques also allow for correlated somatic and dendritic input fluctuations (cf. Methods section 3); an investigation into such input correlations, however, is beyond the scope of the present study. The methodological results presented here may further be used to derive a simple spike rate model in form of a low-dimensional differential equation from two-compartment model neurons. This may be achieved by adapting available approximation methods [32] to our reduced description based on the Fokker-Planck equation.

Described phenomena

In our modeling study we focused on pyramidal neurons from cortex in a spontaneous state, in which the neural membrane potentials typically fluctuate close to the threshold for most of the time and spikes are emitted at irregular intervals. For a single neuron an oscillatory electric field with an amplitude of 1 V/m caused a small polarization effect that leads to a modulation of the spiking probability. At the level of a local population of similarly oriented neurons these small effects became apparent in the collective activity (population spike rate) which exhibited a modulation of a few spikes/s per neuron on average. This type of stochastic resonance is a hypothesized mechanism, perhaps the principle one, by which transcranial electrical stimulation causes immediate functional effects [3].

We found that spiking activity is most strongly modulated in a rather narrow frequency band. This resonance robustly emerged for a range of plausible input statistics. The resonance frequency varied depending on the inputs; for the biologically relevant regime of strong input fluctuations, however, the resonance peak was constrained to the beta and gamma frequency bands (∼15-45 Hz for single neurons and ∼15-55 Hz for networks). The phenomenon clearly differs from the spike rate resonance for oscillatory input modulations shown in integrate-and-fire point neurons (see, e.g. [46] or Fig 7C in [32]), which only occurs for mean-dominated inputs (hence, fairly regular spiking) and where the resonance frequency rapidly increases as the input increases.

Recent modeling results suggest that this resonance frequency may be largely determined by the location of background input (soma vs. distal dendrite) [25]. Our results confirm that resonance frequencies are higher for mainly dendritic background inputs compared to mainly somatic background inputs (cf. Fig 3A). However, unless the input at one location is completely extinguished (as in [25]) the statistics of background inputs at either location appear to play the dominant role in determining this frequency. It may also be noted at this point that spike rate resonance frequencies are lower for model neurons that include an exponential nonlinearity (at the soma) compared to purely leaky integrate-and-fire neurons [25].

To date, experimental evidence for frequency-dependent modulation of neuronal activity by extracellular fields is very sparse (see [47] for a review). Weak electric fields alternating at 30 Hz have been shown to increase spiking coherence of pyramidal cells in rat hippocampal slices [18], and fields with high-frequency components have been evidenced to entrain spiking activity in ferret primary visual cortex more effectively than fields that only contain low-frequency components [4]. To the best of our knowledge, the effects of multiple field frequencies have not yet been experimentally assessed. Therefore, most of our results on spike rate resonance are currently not confirmed and may be regarded as predictions.

Interestingly, these resonance effects at the suprathreshold level were not shown by the subthreshold membrane voltage, whose amplitude monotonically decreased with increasing field frequency (cf. Figs 1D and 2B). The latter phenomenon is in accordance with electrophysiological the subthreshold response amplitude (which scales linearly with the field amplitude [23]) is of the same order of magnitude as that measured in pyramidal cells and decreases with increasing field frequency [17]. Note that the parameter values of our model were not optimized to reproduce this behavior quantitatively but can be adjusted accordingly in a straightforward way (cf. S1 Fig) [25].

The two-compartment model naturally accounts for input filtering caused by the dendrite (for example, as described in [25]). Notably, the somatic impedance and spike rate response for input modulations at the soma exhibited rather distinct dependencies on input frequency (compare Figs 1C and 2C), which may be attributed to the nonlinearity caused by spiking. The presence of the dendrite lead to increased spike rate responses to weak oscillatory input modulations with high frequency at the soma, also for strongly fluctuating background inputs. This is consistent with previous findings for Purkinje neurons in vitro and using a two-compartment model similar to the one used here [37] despite marked differences in the morphology between Purkinje and pyramidal cells. Specifically, the ratio between dendritic and somatic membrane surface (hence, capacitance) are quite distinct for these two types of neurons, which explains certain differences in the results.

Methods

Python code for the models and methods presented below is freely available at GitHub: https://github.com/neuromethods/two-compartment-models-and-weak-electric-fields.

1. Two-compartment neuron model

The two-compartment neuron model consists of two differential equations for the dynamics of the somatic and the dendritic membrane voltage, Vs and Vd, respectively, together with a reset condition of the integrate-and-fire type (for similar models with and without an extracellular field see [35, 37, 38, 48]), CsdVsdt+Iion,s(Vs)=Gi(Vd-Vs-ΔE(t))+Is(t),(1) CddVddt+GdVd=Gi(Vs-Vd+ΔE(t))+Id(t),(2) Iion,s(Vs)≔GsVs−GeΔTeVs−VTΔT,(3) ifVs(t)≥VththenVs(t)←Vr,(4) where Vs, Vd are defined by the difference between the deviations Vsi, Vdi of the intracellular membrane potentials from the leak reversal potential (assumed identical for soma and dendrite) and the extracellular membrane potentials Vse, Vde for the somatic and dendritic compartment, respectively, Vs(t)≔Vsi(t)-Vse(t),Vd(t)≔Vdi(t)-Vde(t).(5)

Cs, Cd and Gs, Gd denote the capacitances and leak conductances of the somatic and dendritic membranes. The exponential term with conductance Ge, threshold slope factor ΔT and effective threshold voltage VT approximates the rapidly increasing Na^+^ current at spike initiation [49]. Gi is the internal conductance between the somatic and the dendritic compartment, Δ is the spatial distance between their centers, and E denotes the extracellular electric field, defined by E(t)≔Vse(t)-Vde(t)Δ.(6)

Is and Id are the synaptic input currents at the soma and dendrite, respectively. When Vs increases beyond VT, it diverges to infinity in finite time due to the exponential term, which defines a spike. In practice, however, the spike is said to occur when Vs reaches a given threshold value Vth > VT. The downswing of the spike is not explicitly modeled; instead, when Vs passes Vth (from below), the somatic membrane voltage is instantaneously reset to a lower value Vr, cf. (4).

Eqs (1) and (2) are the current balance equations for the center points of the two compartments according to the electrical circuit diagram in Fig 1A. This can be seen by using Kirchhoff’s law, that all incoming currents at a circuit point must sum to zero, and the definitions (5) and (6) which imply Vdi-Vsi=Vd-Vs-ΔE. We consider an applied weak sinusoidal field, E(t)=E0+E1sin(ωt),(7) with offset E0 = 0, amplitude E1 = 1 V/m and angular frequency ω, unless stated otherwise.

The synaptic inputs are fluctuating currents that mimic the compound effect of synaptic bombardment in-vivo, described by Is(t)≔I¯s(t)+σs(t)ξs(t),(8) Id(t)≔I¯d(t)+σd(t)ξd(t)(9) with time-varying moments I¯s, I¯d and σs, σd, and uncorrelated unit Gaussian white noise processes ξs, ξd, i.e., 〈ξs(t)ξd(t + τ)〉 = δ(τ)δsd, where 〈 ⋅ 〉 denotes expectation (with respect to the ensemble of noise realizations at times t and t + τ) and δsd is the Kronecker delta.

2. Calculation of subthreshold responses

We analytically calculate the somatic membrane voltage response for small amplitude variations of the synaptic inputs Is(t), Id(t) and a weak oscillatory field E(t), which do not elicit spikes. Considering that the somatic voltage evolves sufficiently below the effective threshold value VT allows us to neglect the exponential term in Eq (1) (i.e., the somatic membrane is purely leaky and capacitive). Using the Fourier transform of Eqs (1) and (2) we then V^s(ω)=I^s(ω)Z^2CIs(ω)+I^d(ω)Z^2CId(ω)+GiΔE^(ω)[Z^2CId(ω)-Z^2CIs(ω)],(10) where Z^2CIs(ω)=1/[(Csiω+Gs+Gi)-Gi2/(Cdiω+Gd+Gi)],(11) Z^2CId(ω)=Z^2CIs(ω)Gi/(Cdiω+Gd+Gi)(12) are the somatic impedances for inputs at the soma and at the dendrite, respectively. .^ indicates a Fourier transformed variable and ω denotes angular frequency.

3. Calculation of spike rate responses

To assess spiking activity we solved the stochastic differential equations Eqs (1)–(4), (7)–(9) numerically using the Euler-Maruyama integration scheme with time steps between 0.01 ms and 0.05 ms. In addition, we employed analytical calculations as described in the following.

Fokker-Planck system

For improved readability we rewrite Eqs (1)–(3) with extracellular field according to Eq (7) and synaptic input given by Eqs (8) and (9) in compact dVsdt=f(Vs)+aVd+μs(t)+σss(t)ξs(t)+σsd(t)ξd(t),(13) dVddt=bVs+cVd+μd(t)+σds(t)ξs(t)+σdd(t)ξd(t),(14) where the coefficients on the right hand side depend on the parameters of the system described in Methods section 1. Note that here σsd = σds = 0, since the input fluctuations at the soma and dendrite are uncorrelated; however, the methods in this section may also be applied in scenarios where any of these parameters is nonzero and varies over time.

The dynamics of the joint membrane voltage probability density p(Vs, Vd, t) for this system plus reset condition (4) are governed by the Fokker-Planck equation (see, e.g. [32, 44, 50]) ∂p∂t+∂qs∂Vs+∂qd∂Vd=0,(15) where qs and qd are the probability fluxes for the somatic and dendritic membrane voltage, respectively, given by qs≔[f(Vs)+aVd+μs(t)]p-σss2(t)2∂p∂Vs-σ˜(t)2∂p∂Vd,(16) qd≔[bVs+cVd+μd(t)]p-σdd2(t)2∂p∂Vd-σ˜(t)2∂p∂Vs(17) with σ˜(t)≔σss(t)σds(t)+σsd(t)σdd(t), subject to the boundary p(Vth,Vd,t)=0(absorbingboundary)(18) limVs→-∞qs(Vs,Vd,t)=0(reflectingboundary)(19) limVd→±∞qd(Vs,Vd,t)=0(reflectingboundaries)(20) limVs↘Vrqs(Vs,Vd,t)-limVs↗Vrqs(Vs,Vd,t)=qs(Vth,Vd,t)(re-injection).(21)

The last condition, a re-injection of probability flux, accounts for the voltage reset in the neuron model. We obtain the instantaneous spike rate as r(t)=∫-∞∞qs(Vth,Vd,t)dVd.(22)

Dimension reduction

Solving the 2+1 dimensional Fokker-Planck partial differential equation (PDE) system (15)–(21) numerically is possible in principle, but computationally demanding. Here we reduce the dimension of this PDE system, which greatly reduces the effort to calculate a) the steady-state spike rate, and b) spike rate responses to weak sinusoidal variations of the input moments or the applied field. For that purpose only ordinary differential equation (ODE) systems need to be solved (see further below). This advantage in terms of implementation and particularly computation becomes crucial for derived mean-field models of multiple neuronal populations (see Methods section 6).

To reduce the dimension of the PDE system we utilize a moment closure approximation method. The (full) probability density p can be expressed in terms of the marginal probability density for the somatic voltage, ps, and the conditional probability density for the dendritic voltage, pd, as p(Vs, Vd, t) = ps(Vs, t) pd(Vd|Vs, t). Note that pd is characterized by a (potentially) infinite number of conditioned moments {η~~d,1~~(Vs, t), η~~d,2~~(Vs, t), …}. The method approximates pd by considering only the first k moments as described below (see [41] for the application of such a method in a different setting). We transform the PDE system (15)–(21) into a system of 1+1-dimensional PDEs ∂ps∂t=L0(ηd,1)[ps]∂psηd,1∂t=L1(ηd,1,ηd,2)[ps]∂psηd,2∂t=L2(ηd,1,ηd,2,ηd,3)[ps]…(23) together with the associated boundary conditions by multiplying Eqs (15)–(21) with Vdl for l ∈ {0, 1, 2, …} and integrating over Vd assuming that p and qd tend sufficiently fast to zero for Vd → ±∞, i.e., limVd→±∞Vdlp(Vs,Vd,t)=0, limVd→±∞Vdlqd(Vs,Vd,t)=0. Each linear operator Ll in (23) depends on the next higher conditioned moment η~~d,l+1~~ as indicated, hence the system is (potentially) infinitely large. Note that we have omitted the obvious arguments Vs, t for ps, η~~d,l~~ for improved readability, the linear operators are specified further below.

We close the system (23) at k = 3 by setting the 3^rd^ central moment of Vd (as well as higher cumulants) to zero, such that ηd,3=3ηd,1ηd,2-2ηd,13, thereby assuming that pd can be sufficiently well approximated by a conditioned Gaussian probability density, pd(Vd|Vs,t)∝exp{-[Vd-ηd,1(Vs,t)]22[ηd,2(Vs,t)-ηd,12(Vs,t)]}.(24)

For a motivation of this assumption see the remark below Eq (38). In the following we specify the system of 3 coupled 1+1-dimensional PDEs we obtain in this way. Using the definitions p~~s,1~~ ≔ psη~~d,1~~ and p~~s,2~~ ≔ psη~~d,2~~ the system is given ∂ps∂t+∂us∂Vs=0(25) ∂ps,1∂t+∂us,1∂Vs-[bVs+μd(t)]ps-cps,1+σ˜(t)2∂ps∂Vs=0(26) ∂ps,2∂t+∂us,2∂Vs-2[bVs+μd(t)]ps,1-2cps,2-[σds2(t)+σdd2(t)]ps+σ˜(t)∂ps,1∂Vs=0(27) with us=∫-∞∞qsdVd=[f(Vs)+μs(t)]ps+aps,1-σss2(t)+σsd2(t)2∂ps∂Vs(28) us,1=∫-∞∞VdqsdVd=[f(Vs)+μs(t)]ps,1+aps,2-σss2(t)+σsd2(t)2∂ps,1∂Vs+σ˜(t)2ps(29) us,2=∫-∞∞Vd2qsdVd=[f(Vs)+μs(t)]ps,2+ahs-σss2(t)+σsd2(t)2∂ps,2∂Vs+σ˜(t)ps,1(30) and hs=3ps,1ps,2ps-2ps,13ps(31) subject to the conditions ps(Vth,t)=0(32) limVs→-∞us(Vs,t)=limVs→-∞us,1(Vs,t)=limVs→-∞us,2(Vs,t)=0(33) limVs↘Vrus(Vs,t)-limVs↗Vrus(Vs,t)=us(Vth,t)(34) limVs↘Vrus,1(Vs,t)-limVs↗Vrus,1(Vs,t)=us,1(Vth,t)(35) limVs↘Vrus,2(Vs,t)-limVs↗Vrus,2(Vs,t)=us,2(Vth,t)(36) and requiring that ps is initialized such that ∫-∞Vthps(Vs,0)dVs=1. The spike rate is then obtained by r(t)=us(Vth,t).(37)

Note that ps(Vth, t) = 0 implies p~~s,1~~(Vth, t) = p~~s,2~~(Vth, t) = 0. The conditions (33) follow from condition (19) and a self-consistency requirement with respect to the dynamics of the unconditioned moments ηd,1(t)=∫-∞Vthps,1(Vs,t)dVs and ηd,2(t)=∫-∞Vthps,2(Vs,t)dVs. The latter can be seen by integration of Eqs (26) and (27) over Vs and comparison with the moment equations obtained by successive integration of Eq (15) over Vs, multiplication by Vd or Vd2, respectively, and integration over Vd. Conditions (34)–(36) follow from (21). Note also that us,1(Vth,t)=r(t)ηd,1(Vth,t),us,2(Vth,t)=r(t)ηd,2(Vth,t).(38)

Remark: The assumption that pd can be sufficiently well approximated by a conditioned Gaussian is supported by the circumstance that for subthreshold inputs and an electric field which keep the somatic voltage (sufficiently) below the spike threshold the approximation is excellent. In that case pd(Vd|Vs, t) is indeed a conditioned Gaussian probability density (because the exponential term in Eq (1) as well as conditions (18) and (21) are negligible). Since this is not the case for stronger inputs (that cause spiking) the reproduction performance of this approximation needs to be evaluated (cf. Figs 1E, 1G and 2 and Discussion).

Steady state

For the steady state (in case of constant parameters μs, μd, σss, σds, σsd, σdd) we obtain, by setting the time derivatives in Eqs (25)–(27) to zero, the 6-dimensional ODE system dusdVs=0dpsdVs=2[(f(Vs)+μs)ps+aps,1-us]σss2+σsd2(39) dus,1dVs=(bVs+μd)ps+cps,1-σ˜2dpsdVs(40) dps,1dVs=2[(f(Vs)+μs)ps,1+aps,2+σ˜ps/2-us,1]σss2+σsd2(41) dus,2dVs=2(bVs+μd)ps,1+2cps,2+(σds2+σdd2)ps-σ˜dps,1dVs(42) dps,2dVs=2[(f(Vs)+μs)ps,2+ahs+σ˜ps,1-us,2]σss2+σsd2(43) subject to the conditions (32)–(36) (with time dependence omitted).

We solve this nonlinear ODE system (nonlinearity due to hs, cf. Eq (31)) with variable coefficients numerically by integrating Eqs (39)–(43) backwards from Vth with ps(Vth) = p~~s,1~~(Vth) = p~~s,2~~(Vth) = 0, us(Vth) = 1, u~~s,1~~(Vth) = η~~d,1~~(Vth), u~~s,2~~(Vth) = η~~d,2~~(Vth) to a sufficiently small (lower bound) voltage value Vlb, taking into account the “jump” conditions (34)–(36), and determine η~~d,1~~(Vth) and η~~d,2~~(Vth) such that u~~s,1~~(Vlb) = u~~s,2~~(Vlb) = 0, cf. condition (33), is fulfilled. Then, the scaling factor r (the spike rate) of the obtained solution is determined such that ∫VlbVthps(Vs)dVs=1 holds. The solution was achieved by means of a Python implementation using a root finding algorithm provided by the package Scipy [51] (optimize.root, modification of the Powell hybrid method [52]) and low-level virtual machine acceleration through the package Numba [53].

Response to modulations

In order to characterize the spike rate dynamics we calculate the response to small-amplitude sinusoidal variations of the mean input around a baseline, I¯s(t)=I¯s0+I¯s1sin(ωt), I¯d(t)=I¯d0+I¯d1sin(ωt), or to a weak sinusoidal field, cf. Eq (7). These modulations translate to sinusoidal modulations of μs(t) and μd(t) in the system (25)–(36). The parameters σss, σds, σsd and σdd remain constant.

For mathematical convenience we write the modulations in complex form μs(t)=μs0+μs1eiωtμd(t)=μd0+μd1eiωt,(44) with small μs1>0 and μd1>0 (thereby introducing a companion system) and approximate the solution to first order, r(t)=r0+r^1(ω)eiωt, where r^1 is a complex variable from which the response amplitude r1 and phase shift ψ of the oscillatory spike rate r0 + r1(ω)sin(ωt+ ψ(ω)) can be extracted in a straightforward r1=|r^1|, ψ=arg(r^1) (see, e.g. [46] for a similar type of analysis in a different setting). Note that also the state variables of this solution take the form ps(Vs,t)=ps0(Vs)+p^s1(Vs;ω)eiωt (analogously for p~~s,1~~, p~~s,2~~, us, u~~s,1~~, u~~s,2~~). For fixed (angular) frequency ω we obtain the following ODE system (neglecting terms of second and higher order in μs1,μd1): du^s1dVs=-iωp^s1dp^s1dVs=2[(f(Vs)+μs0)p^s1+μs1ps0+ap^s,11-u^s1]σss2+σsd2(45) du^s,11dVs=-iωp^s,11+(bVs+μd0)p^s1+μd1ps0+cp^s,11-σ˜2dp^s1dVs(46) dp^s,11dVs=2[(f(Vs)+μs0)p^s,11+μs1ps,10+ap^s,21+σ˜p^s1/2-u^s,11]σss2+σsd2(47) du^s,21dVs=2(bVs+μd0)p^s,11+2μd1ps,10+2cp^s,21+(σds2+σdd2)p^s1-σ˜dp^s,11dVs(48) dp^s,21dVs=2[(f(Vs)+μs0)p^s,21+μs1ps,20+ah^s1+σ˜p^s,11-u^s,21]σss2+σsd2(49) with h^s1=[4(ps,10)3(ps0)3-3ps,10ps,20(ps0)2]ps1+[3ps,20ps0-6(ps,10)2(ps0)2]ps,11+3ps,10ps0ps,21(50) subject to the conditions p^s1(Vth;ω)=0(51) limVs→-∞u^s1(Vs;ω)=limVs→-∞u^s,11(Vs;ω)=limVs→-∞u^s,21(Vs;ω)=0(52) limVs↘Vru^s1(Vs;ω)-limVs↗Vru^s1(Vs;ω)=u^s1(Vth;ω)(53) limVs↘Vru^s,11(Vs;ω)-limVs↗Vru^s,11(Vs;ω)=u^s,11(Vth;ω)(54) limVs↘Vru^s,21(Vs;ω)-limVs↗Vru^s,21(Vs;ω)=u^s,21(Vth;ω),(55) where .^ indicates a complex valued variable that depends on ω. Note that this ODE system depends on the (steady state) solution of the system (39)–(43) through ps0, ps,10 and ps,20.

The linear (complex valued) ODE system with variable coefficients (45)–(55) can be conveniently solved in the following way. The desired spike rate response solution can be written as r^1(ω)=μs1r^1μs(ω)+μd1r^1μd(ω) where the solution component r^1μs(ω) corresponds to μd1=0 and, vice versa, r^1μd(ω) corresponds to μs1=0. We first describe how we obtain r^1μs(ω) (r^1μd(ω) is calculated in an analogous way, see further below). The solution xs(Vs;ω)≔(p^s1,p^s,11,p^s,21,u^s1,u^s,11,u^s,21) associated with r^1μs(ω) can be decomposed into (omitting the argument ω below for improved readability) xs(Vs)=r^1μsxα(Vs)+η^d,11(Vth)xβ(Vs)+η^d,21(Vth)xγ(Vs)+xδ(Vs),(56) where xα, xβ, xγ solve the homogeneous part (μs1=0, μd1=0) of the ODE system (45)–(50) with p^x1(Vth)=p^x,11(Vth)=p^x,21(Vth)=0 for x ∈ {α, β, γ}, u^α1(Vth)=1, u^β1(Vth)=u^γ1(Vth)=0, u^α,11(Vth)=ηd,10(Vth), u^α,21(Vth)=ηd,20(Vth), u^β,11(Vth)=r0, u^β,21(Vth)=0, u^γ,11(Vth)=0, u^γ,21(Vth)=r0, and “jump” conditions (53)–(55). Note that r0, ηd,10(Vth) and ηd,20(Vth) are known from the solution for the steady state system. xδ solves the inhomogeneous system (45)–(50) with μs1=1 (μd1=0) and condition xδ(Vth) = 0. These solutions are obtained numerically by backward integration from Vth to Vlb. r^1μs(ω) together with η^d,11(Vth),η^d,21(Vth) are then calculated by solving the linear equation system that arises to satisfy the condition u^s1(Vlb)=u^s,11(Vlb)=u^s,21(Vlb)=0. The solution method for r^1μd(ω) is completely analogous with the difference that r^1μd appears instead of r^1μs in Eq (56) and xδ solves the inhomogeneous system (45)–(50) with μd1=1 (μs1=0).

We obtain the amplitude of the spike rate response to an applied field from r^1(ω)=GiΔE1Cdr^1μd(ω)-GiΔE1Csr^1μs(ω),(57) whereas the response modulations to sinusoidal mean input at the soma and dendrite in the absence of an oscillatory field are given by r^1(ω)=I¯s1r^1μs(ω)/Cs and r^1(ω)=I¯d1r^1μd(ω)/Cd, respectively.

4. Parametrization via ball-and-stick model

In the following we describe a semi-analytical technique to fit the two-compartment (2C) model to a biophysically more detailed, spatially extended ball-and-stick (BS) model. In particular, the parameter values of the 2C model are determined to best approximate the somatic voltage dynamics of the BS model. This is done in an efficient way using analytical results for the voltage dynamics of both models, and it does not depend on a particular choice of parameter values for the input or the extracellular field. This part may thus be regarded as a reduction of the BS model.

Ball-and-stick model

The BS neuron model consists of a finite passive dendritic cable of length L with lumped somatic compartment at the proximal end, x = 0, and a sealed end boundary condition at the distal extremity, x = L. It includes capacitive and leak currents across the membrane, an approximation of the spike-generating sodium current at the soma, an internal current (along the cable) and synaptic input currents at the soma and distal dendrite as well as a spatially homogeneous but time-varying external electric field (for details on the derivation of this model see [25]). The dynamics of the model are governed by cm∂V∂t-gi∂2V∂x2+gmV=0,0<x<L,(58) cs∂V∂t-gi∂V∂x+gsV-gsΔTeV-VTΔT=Is(t)-giE(t),x=0,(59) ∂V∂x=Id(t)gi+E(t),x=L(60) together with a reset condition ifV(0,t)≥VththenV(0,t)≔Vr′,(61) where cm = cDd π is the membrane capacitance, gm = ϱm Ddπ is the membrane conductance and gi = ϱi(Dd/2)^2^π is the internal (axial) conductance of a dendritic cable segment of unit length. c is the specific membrane capacitance (in F/m^2^), ϱi is the specific internal conductance (in S/m), ϱm is the specific membrane conductance (in S/m^2^) and Dd is the cable diameter. cs=cDs2π and gs=ϱmDs2π are the somatic membrane capacitance and leak conductance, respectively, with soma diameter Ds. The exponential term with threshold slope factor ΔT and effective threshold voltage VT approximates the rapidly increasing Na^+^ current at spike initiation [49]. Spike times are defined by the times at which the somatic membrane voltage V(0, t) crosses the threshold voltage value Vth from below (cf. spike mechanism of the 2C model). Is(t), Id(t) and E(t) are described by Eqs (8), (9) and (7), respectively. The parameter values are provided in Table 1.

To generate spike trains we simulated the BS neuron model using a semi-implicit numerical scheme (Crank-Nicolson method; see, e.g. appendix C of [59]) that was extended for stochasticity as proposed in [60], and by applying the tridiagonal matrix algorithm. Discretization steps were 5 μs for time and L/200 for space (along the dendrite).

Calculation of subthreshold responses

We analytically calculate the somatic membrane voltage response of the BS model for small variations of the synaptic inputs Is(t), Id(t) and a weak oscillatory field E(t), which do not elicit spikes. We consider that the somatic voltage evolves sufficiently below VT, which allows us to neglect the exponential term in Eq (59) (cf. Methods section 2). The linear PDE (58) together with the boundary conditions (59) and (60) is then solved using separation of variables V(x, t) = W(x)U(t) and the Fourier transform V^(x,ω)=W(x)U^(ω)=W(x)∫-∞∞U(t)e-iωtdt.(62)

We obtain the system of differential equations cmiωV^-gi∂2V^∂x2+gmV^=00<x<L,(63) csiωV^-gi∂V^∂x+gsV^=I^s(ω)-giE^(ω)x=0,(64) ∂V^∂x=I^d(ω)gi+E^(ω)x=L,(65) where .^ indicates a (temporally) Fourier transformed variable. The solution of this system can be expressed as V^(0,ω)=I^s(ω)Z^BSIs(ω)+I^d(ω)Z^BSId(ω)+giE^(ω)Z^BSId(ω)-Z^BSIs(ω) with Z^BSIs(ω)=1/[csiω+gs+z(ω)gitanh(z(ω)L)],(67) Z^BSId(ω)=Z^BSIs(ω)/cosh(z(ω)L),(68) where ±z(ω) are the roots of the characteristic polynomial giy^2^ = gm + cmiω of Eq (63), z(ω)=gm+gm2+ω2cm22gi+sgn(ω)i-gm+gm2+ω2cm22gi.(69)

Note that Z^BSIs and Z^BSId are the somatic impedances for inputs at the soma and the distal dendrite, respectively.

The response to a sinusoidal field variation, E(t) = E1 sin(φt), with constant Is and Id can be expressed in the time domain as V(0,t)=|E1[Z^BSId(φ)-Z^BSIs(φ)]|sin(φt+arg[Z^BSId(φ)-Z^BSIs(φ)]).(70)

In addition we calculate the somatic voltage time series in response to subthreshold input for a given initial condition V(x, 0) = V0(x). Using separation of variables and the Laplace transform V˜(x,s)=W(x)U˜(s)=W(x)∫0∞U(t)e-stdt,(71) with complex variable s we obtain cm(sV˜−V0)−gi∂2V˜∂x2+gmV˜=00<x<L,(72) cs(sV˜−V0)−gi∂V˜∂x+gsV˜=I˜s(s)−giE˜(s)x=0,(73) ∂V˜∂x=I˜d(s)gi+E˜(s)x=L,(74) where .˜ indicates a (temporally) Laplace transformed variable. We solve this system and obtain V˜(0,s)=[I˜s(s)+csV0(0)]Z˜BSIs(s)+[I˜d(s)+cm∫0Lcosh(z(s)(L−y))V0(y)dy]Z˜BSId(s)+giE˜(s)Z˜BSId(s)−Z˜BSIs(s) with Z˜BSIε(s)=1/[csis+gs+z(s)gitanh(z(s)L)],(76) Z˜BSId(s)=Z˜BSIs(s)/cosh(z(s)L),(77) where ±z(s) are the roots of the characteristic polynomial giy^2^ = cms + gm of Eq (72), given by z(s)=cms+gmgi.(78)

The somatic voltage time series V(0, t) is then computed by inverse transforming V˜(0,s) using an efficient numerical method [61].

Parameter fitting

We approximate the somatic voltage dynamics of the BS model by the 2C model in two steps. First, we fit V^s(ω) to V^(0,ω) using Eqs (10) and (66) over a range of angular frequencies ω ∈ [0, ω~~max~~], requiring that the voltage values for ω = 0 (i.e. the steady states) match exactly. This constraint determines three parameters, Gd=(gs-Gs)cosh(L/λ)+λgmsinh(L/λ),(79) Gi=Gd/[cosh(L/λ)-1],(80) Δ=gi/Gi.(81) where we have introduced the electrotonic length constant λ≔gi/gm. The remaining three subthreshold parameters Cs, Cd and Gs are obtained using the method of least squares, with ω~~max~~/(2π) = 10 kHz. In the second step we determine the reset voltage Vr by approximating the transient BS somatic voltage time series immediately after a spike elicited by threshold somatic input and threshold dendritic input for E(t) = 0. Specifically, we fit the post-spike voltage time series Vs(t) to V(0, t) across the time interval t ∈ [0, τs] with initial conditions Vs(0) = Vr, Vd(0) = (GiVth + Id)/(Gd + Gi) and V(0,0)=Vr′, V(x,0)=Iscosh((L-x)/λ)+Id[cosh(x/λ)+gs/(λgm)]gscosh(L/λ)+λgmsinh(L/λ)0<x<L(82) for threshold somatic input Is = Vth[gs + λgm tanh(L/λ)], Id = 0 and for threshold dendritic input Id = Vth[gs cosh(L/λ) + λgm sinh(L/λ)], Is = 0 simultaneously using the method of least squares. τs = Cs/(Gs + Gi) is the somatic membrane time constant of the 2C model. Note that we consider threshold input as constant current that yields V = VT (calculated from the linear subthreshold model systems). The voltage time series V(0, t) of the BS model is rapidly computed using the Laplace transform, Eq (75), and the voltage time series Vs(t) of the 2C model is calculated analytically in a straightforward way (linear ODE system). To guarantee an equal effectiveness of the exponential term on the somatic membrane voltage dynamics in the 2C model compared to the BS model we set Ge = Cs gs/cs (cf. Eqs (1) and (59)). The values for ΔT, VT and Vth are set equal to those of the BS model. Notably, this fitting method is very efficient, since it involves analytical results and does not depend on specific realizations of time series for the neuronal input or extracellular field.

5. Spike coincidence measure

To quantify the similarity between the spike trains of the two-compartment and the ball-and-stick model neurons we used the coincidence factor Γ defined by [62] ΓBS,2C=Nc-⟨Nc⟩(NBS+N2C)/21N,(83) where Nc is the number of coincident spikes with precision (i.e., maximal temporal separation) Δc, NBS and N2C are the number of spikes in the spike trains of the ball-and-stick and the two-compartment models, respectively. 〈Nc〉 = 2rΔcNBS is the expected number of coincidences generated by a homogeneous Poisson process with spike rate r = N2C/T as exhibited by the two-compartment model, where T is the duration of the spike train. The factor N=1-2rΔc normalizes Γ~~BS,2C~~ to a maximum value of 1, which is attained if the spike trains match optimally (with precision Δc). Γ~~BS,2C~~ = 0 on the other hand would result from a homogeneous Poisson process with rate that corresponds to the spike train of the two-compartment model, and therefore indicates pure chance. Here we used Δc = 5 ms.

6. Two-population network

In this section we describe the network model of a large number N of sparsely and randomly coupled pyramidal (PY) neurons and inhibitory (IN) interneurons, and derive a mean-field system from it. Each PY neuron is described by the 2C model, Eqs (1)–(4), (7)–(9), and each IN neuron is described by an exponential integrate-and-fire (point) neuron model, CINdVdt=-GINV+GINΔTeV-VTΔT+IIN(t),(84) ifV≥VththenV←Vr′,(85) because IN neurons do not exhibit an elongated spatial morphology compared with PY neurons. We used CIN = 0.2 nF and GIN = 10 nS. The model neurons receive fluctuating external and recurrent synaptic input and are exposed to an applied weak electric field E(t) that is spatially uniform. For fields induced by transcranial brain stimulation [10] this is a valid assumption. Each PY neuron receives inputs from Ks IN neurons at the soma and Kd IN neurons at the dendrite and each IN neuron receives inputs from KIN PY neurons. Synaptic coupling is described by delayed current pulses that produce postsynaptic potentials of size Js, Jd or JIN (depending on the location of the synapse). Specifically, the input currents for neuron k, with somatic and dendritic membrane voltage V~~s,k~~, V~~d,k~~ (in case of a PY neuron) or overall membrane voltage Vk (in case of an IN neuron), are given by Ix,k(t)≔I¯x,ext(t)+σx,ext(t)ξx,k(t)+Cx∑lJk,lx∑mδ(t-tlm-dk,lx),(86) x ∈ {s, d, IN}, where Jk,lx∈{Jx,0}, tlm is the m-th spike time of neuron l and dk,lx the coupling delay between neuron l and k. The delays are independently sampled according to the probability distribution with density pdx(d)≔Γ(d)τdx-τrx(e-d/τdx-e-d/τrx)(87) x ∈ {s, d, IN}. We consider large numbers of connections Kx and reasonably small coupling strengths Jx.

We simulated networks of this type in the presence of a sinusoidal weak field and, alternatively, sinusoidal modulations of the mean inputs at the soma or dendrite of PY neurons. Network simulations were performed using the Python-based Brian2 simulator [63]. In addition, we derived a mean-field description and analytically computed network activity in terms of population-averaged spike rates, as described in the following.

Derived mean-field model and resonance analysis

For large networks (in the mean-field limit N → ∞) the overall synaptic input can be approximated by a mean part with additive fluctuations, Ix,k(t)≈I¯x(t,rdx)+σx(t,rdx)ζx,k(t),(88) I¯x(t,rdx)≔I¯x,ext(t)+CxJxKxrdx(t),(89) σx2(t,rdx)≔σx,ext2(t)+Cx2Jx2Kxrdx(t),(90) x ∈ {s, d, IN} with delayed spike rates rds=pdsrIN, rdd=pddrIN, rdIN=pdIN*rPY, and unit white Gaussian noise process ζ~~x,k~~ that is uncorrelated to that of any other neuron (see, e.g. [32]). This step is valid under the assumptions of (i) sufficient presynaptic activity, (ii) that neuronal spike trains can be approximated by Poisson processes and (iii) that the correlations between the fluctuations of synaptic inputs for different neurons vanish. Note that the latter assumption is supported by sparse and random synaptic connectivity.

The previous approximation allows us to express the collective spiking dynamics in terms of a coupled system of Fokker-Planck PDEs, one for the PY population (cf. Methods section 3) and one for the IN population, where the coupling is mediated through the synaptic input moments I¯x and σx2 (cf. Eqs (89) and (90)). Specifically, the system for the PY population consists of Eq (15) with qs≔[Gi(Vd-Vs-ΔE(t))-Iion,s(Vs)+I¯s(t,rds)]pCs-σs2(t,rds)2Cs2∂p∂Vs,(91) qd≔[Gi(Vs-Vd+ΔE(t))-GdVd+I¯d(t,rdd)]pCd-σd2(t,rdd)2Cd2∂p∂Vd(92) subject to the conditions (18)–(21) and rPY given by Eq (22). The system for the IN population is specified by ∂pV∂t+∂qV∂V=0,(93) qV≔[GIN(ΔTeV-VTΔT-V)+I¯IN(t,rdIN)]pCIN-σIN2(t,rdIN)2CIN2∂p∂V(94) subject to the conditions pV(Vth,t)=0,limV→-∞qV(V,t)=0,(95) limV↘Vr′qV(V,t)-limV↗Vr′qV(V,t)=qV(Vth,t)(96) and spike rate given by rIN = qV(Vth, t).

To analyze the resonance properties of this mean-field network model we proceed similarly as in Methods section we consider either a weak sinusoidal electric field, Eq (7), or weak sinusoidal modulations of the external mean input at the soma or dendrite of PY neurons, I¯x,ext(t)=I¯x,ext0+I¯x,ext1eiωt, x ∈ {s, d}. We assume a parametrization of the network such that without these modulations (E1=I¯x,ext1=0) it exhibits asynchronous activity represented by a fixed point solution of the mean-field system. We write the modulations in complex form and express the first order population spike rate responses as rPY(t)=rPY0+r^PY1(ω)eiωt and rIN(t)=rIN0+r^IN1(ω)eiωt. To obtain the stationary (steady state) components rPY0 and rIN0 we solve the system (39)–(43) with μx0=I¯x,ext0/Cx+JxKxrIN0 and σx2=σx,ext2/Cx2+Jx2KxrIN0, x ∈ {s, d} (and with the corresponding boundary conditions) together with the respective system for the IN population with μIN0=I¯IN,ext0/CIN+JINKINrPY0 and σIN2=σIN,ext2/CIN2+JIN2KINrPY0 via fixed point iteration of the form (rPY,n+10,rIN,n+10)=F(rPY,n0,rIN,n0), where n denotes the iteration number. Recall that σsd = σds = 0 in Eqs (13) and (14), as noted above. The response components r^PY1 and r^IN1 are obtained from r^PY1(ω)=I¯s,ext1−GiΔE1Csr^1μs(ω)+JsKsp^ds(ω)r^IN1(ω)[r^1μs(ω)+Jsr^1σs2(ω)]+I¯d,ext1+GiΔE1Cdr^1μd(ω)+JdKdp^dd(ω)r^IN1(ω)[r^1μd(ω)+Jdr^1σd2(ω)],(97) r^IN1(ω)=JINKINp^dIN(ω)r^PY1(ω)[r^1μIN(ω)+JINr^1σIN2(ω)],(98) where p^dx(ω)=(τrxiω+1)-1(τdxiω+1)-1, x ∈ {s, d, IN} is the Fourier transformed delay density (cf. Eq (87)) and r^1μx, r^1σx2 are the population spike rate response components for sinusoidal modulations of μx, σx2 of unit amplitude. Note that Eqs (97) and (98) can be jointly solved in a straightforward way. In addition to first order responses to modulations of μs and μd we therefore need to calculate the responses for weak sinusoidal modulations of σs2 and σd2. That is, the rate response solution of the system (25)–(36) for σss2(t)=σs,02+σs,12eiωt and σdd2(t)=σd,02+σd,12eiωt is required. This is done in an analogous way as explained above for sinusoidal modulations of μs and μd (see Methods section 3). In particular, we solve (45)–(55) with σss2=σs,02, σdd2=σd,02, σsd2=σds2=0, and where the inhomogeneous terms μs1ps0, μs1ps,10, μs1ps,20 and 2μd1ps,10 in Eqs (45), (47), (49) and (48) are replaced by -σs,122dps0dVs, -σs,122dps,10dVs, -σs,122dps,20dVs and σd,12ps0, respectively. The resulting system can be numerically solved as explained in Methods section 3, using Eq (56) where r^1μs is replaced by r^1σs2 and r^1σd2, respectively, and xδ solves the adjusted inhomogeneous system (with σs,12=1, μs1=μd1=σd,12=0 and σd,12=1, μs1=μd1=σs,12=0, respectively). For the simpler case of the exponential integrate-and-fire point model for the IN population (steady state and response modulations) we refer to Refs. [32, 46].

Weak electric fields promote resonance in neuronal spiking activity: Analytical results from two-compartment cell and network models