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

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.

Multi-compartment models are useful tools to dissect effects in single neurons in the absence of input fluctuations [27][28][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][31][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][35][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)  Table 1). C: amplitude of subthreshold somatic impedances for inputs at the somaẐ I s and dendriteẐ I d as a function of input frequency (using Eqs (11), (12) and (67)-(69)). D: amplitude of subthreshold somatic voltage responses to a sinusoidal electric field with amplitude E 1 = 1 V/m (and E 0 = 0) as a function of field frequency (using Eqs (10)- (12) and (66)-(69)). E: somatic voltage transients after a spike for constant threshold inputs at the soma (left) and dendrite (right) of the BS model (black) and 2C model (green); for details see Methods section 4. F: example time series of the somatic voltage in response to fluctuating inputs, with spike times indicated, as well as voltage histograms of both models (right) and voltage density p s (green line; analytically calculated). Note that although p s appears Gaussian this is not an assumption of our method. G: spike coincidence measure Γ between the spike trains of both models for � I d ¼ 3 pA, s d ¼ 5 pA ffi ffi ffi ffi ffi ffi ms p ðleftÞ and � I s ¼ 3 pA, s s ¼ 15 pA ffi ffi ffi ffi ffi ffi ms p ðrightÞ. Γ = 1 indicates an optimal match with precision Δ c = 5 ms, Γ = 0 indicates chance (for details see Methods section 5). H: spike rates from numerical simulation (symbols) and analytical calculation (green lines; cf. Methods section 3) for input parameter values as in G. The grey square marks the parameter values from F. All curves in C-E and H represent analytical results.
https://doi.org/10.1371/journal.pcbi.1006974.g001 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(V s , V d , 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 where p s is the marginal probability density for the somatic voltage and p d the probability density for the dendritic voltage conditioned on V s , and approximate p d 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 balland-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 s 2 s , s 2 d ) as well as high rate, low variability (due to large � I s , � I d and small s 2 s , s 2 d ; 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 (green circles), as well as amplitude of normalized subthreshold somatic voltage responses (magenta dashed line). C and D: amplitude of spike rate responses (green) and of normalized subthreshold somatic voltage responses (impedances; magenta dashed lines) to sinusoidal modulations of the mean input at the soma, � I s ðtÞ ¼ � I 0 s þ � I 1 s sinðotÞ, with amplitude � I 1 s ¼ 0:5 pA (C) or at the dendrite, as a function of modulation frequency in the absence of an extracellular field. Baseline (constant) input statistics in C and D correspond to those in B. Remark: the spike rate responses to the field and those to modulations of the mean input are mathematically linked (by Eq (57)) such that the former can be calculated from responses to modulations at the soma and dendrite, respectively, using response amplitude and phase information. 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 V s and V d 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 Resonance frequency (argmax ω r 1 (ω)/(2π); A), amplitude (max r 1 (ω); B) and strength (max r 1 (ω)/r 1 (0); C) of the oscillatory spike rate with amplitude r 1 (cf. Methods section 3) in response to a sinusoidal applied field with E 1 = 1 V/m and angular frequency ω for a range of input statistics. Left: 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.

Methodological aspects
Neuron model. Model-based investigations on how extracellular electric fields impact neuronal spiking activity face a methodological challenge: 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 twocompartment 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 work: 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 p d (V d |V s , t) by a conditioned Gaussian, thereby closing the system of lower dimensional equations at the 3 rd central moment of V d (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][44][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 integrateand-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 observations: 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.

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, V s and V d , 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]), where V s , V d are defined by the difference between the deviations V i s , V i d of the intracellular membrane potentials from the leak reversal potential (assumed identical for soma and dendrite) and the extracellular membrane potentials V e s , V e d for the somatic and dendritic compartment, respectively, C s , C d and G s , G d denote the capacitances and leak conductances of the somatic and dendritic membranes. The exponential term with conductance G e , threshold slope factor Δ T and effective threshold voltage V T approximates the rapidly increasing Na + current at spike initiation [49]. G i 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Þ ≔ V e s ðtÞ À V e d ðtÞ D : ð6Þ I s and I d are the synaptic input currents at the soma and dendrite, respectively. When V s increases beyond V T , 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 V s reaches a given threshold value V th > V T . The downswing of the spike is not explicitly modeled; instead, when V s passes V th (from below), the somatic membrane voltage is instantaneously reset to a lower value V r , 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) We consider an applied weak sinusoidal field, with offset E 0 = 0, amplitude E 1 = 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

Calculation of subthreshold responses
We analytically calculate the somatic membrane voltage response for small amplitude variations of the synaptic inputs I s (t), I d (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 V T 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) are the somatic impedances for inputs at the soma and at the dendrite, respectively.: indicates a Fourier transformed variable and ω denotes angular frequency.

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 form: 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(V s , V d , t) for this system plus reset condition (4) are governed by the Fokker-Planck equation (see, e.g. [32,44,50]) where q s and q d are the probability fluxes for the somatic and dendritic membrane voltage, respectively, given by 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 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, p s , and the conditional probability density for the dendritic voltage, p d , as p(V s , V d , t) = p s (V s , t) p d (V d |V s , t). Note that p d is characterized by a (potentially) infinite number of conditioned moments {η d,1 (V s , t), η d,2 (V s , t), . . .}. The method approximates p d 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 together with the associated boundary conditions by multiplying Eqs (15)- (21) with V l d for l 2 {0, 1, 2, . . .} and integrating over V d assuming that p and q d tend sufficiently fast to zero for Each linear operator L l 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 V s , t for p s , η 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 V d (as well as higher cumulants) to zero, such that Z d;3 ¼ 3Z d;1 Z d;2 À 2Z 3 d;1 , thereby assuming that p d can be sufficiently well approximated by a conditioned Gaussian probability density, 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 ≔ p s η d,1 and p s,2 ≔ p s η d,2 the system is given by: and requiring that p s is initialized such that R V th À 1 p s ðV s ; 0ÞdV s ¼ 1. The spike rate is then obtained by Note that p s (V th , t) = 0 implies p s,1 (V th , t) = p s,2 (V th , t) = 0. The conditions (33) follow from condition (19) and a self-consistency requirement with respect to the dynamics of the unconditioned moments Z d;1 ðtÞ ¼ R V th À 1 p s;1 ðV s ; tÞdV s and Z d;2 ðtÞ ¼ R V th À 1 p s;2 ðV s ; tÞdV s . The latter can be seen by integration of Eqs (26) and (27) over V s and comparison with the moment equations obtained by successive integration of Eq (15) over V s , multiplication by V d or V 2 d , respectively, and integration over V d . Conditions (34)-(36) follow from (21). Note also that u s;1 ðV th ; tÞ ¼ rðtÞZ d;1 ðV th ; tÞ; u s;2 ðV th ; tÞ ¼ rðtÞZ d;2 ðV th ; tÞ: Remark: The assumption that p d 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 p d (V d |V s , 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 subject to the conditions (32)-(36) (with time dependence omitted). We solve this nonlinear ODE system (nonlinearity due to h s , cf. Eq (31)) with variable coefficients numerically by integrating Eqs (39)-(43) backwards from V th with p s (V th ) = p s,1 (V th ) = p s,2 (V th ) = 0, u s (V th ) = 1, u s,1 (V th ) = η d,1 (V th ), u s,2 (V th ) = η d,2 (V th ) to a sufficiently small (lower bound) voltage value V lb , taking into account the "jump" conditions (34)- (36), and determine η d,1 (V th ) and η d,2 (V th ) such that u s,1 (V lb ) = u s,2 (V lb ) = 0, cf. condition (33), is fulfilled. Then, the scaling factor r (the spike rate) of the obtained solution is determined such that R V th V lb p s ðV s ÞdV s ¼ 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 0 s þ � I 1 s sin ðotÞ, � I d ðtÞ ¼ � I 0 d þ � I 1 d sin ðotÞ, 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 with small m 1 s > 0 and m 1 d > 0 (thereby introducing a companion system) and approximate the solution to first order, rðtÞ ¼ r 0 þr 1 ðoÞe iot , wherer 1 is a complex variable from which the response amplitude r 1 and phase shift ψ of the oscillatory spike rate r 0 + r 1 (ω)sin(ωt+ ψ(ω)) can be extracted in a straightforward way: r 1 ¼ jr 1 j, c ¼ 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 p s ðV s ; tÞ ¼ p 0 s ðV s Þ þp 1 s ðV s ; oÞe iot (analogously for p s,1 , p s,2 , u s , u s,1 , u s,2 ). For fixed (angular) frequency ω we obtain the following ODE system (neglecting terms of second and higher order in m 1 s ; m 1 d ): where x α , x β , x γ solve the homogeneous part (m 1 g;2 ðV th Þ ¼ r 0 , and "jump" conditions (53)- (55). Note that r 0 , Z 0 d;1 ðV th Þ and Z 0 d;2 ðV th Þ are known from the solution for the steady state system. x d solves the inhomogeneous system (45)-(50) with m 1 s ¼ 1 (m 1 d ¼ 0) and condition x d (V th ) = 0. These solutions are obtained numerically by backward integration from V th to V lb .r m s 1 ðoÞ together withẐ 1 d;1 ðV th Þ;Ẑ 1 d;2 ðV th Þ are then calculated by solving the linear equation system that arises to satisfy the condition u 1 s ðV lb Þ ¼û 1 s;1 ðV lb Þ ¼û 1 s;2 ðV lb Þ ¼ 0. The solution method forr m d 1 ðoÞ is completely analogous with the difference thatr m d 1 appears instead ofr m s 1 in Eq (56) and x d solves the inhomogeneous system (45)- (50) . We obtain the amplitude of the spike rate response to an applied field from whereas the response modulations to sinusoidal mean input at the soma and dendrite in the absence of an oscillatory field are given byr

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 together with a reset condition if Vð0; tÞ � V th then Vð0; tÞ where c m = cD d π is the membrane capacitance, g m = % m D d π is the membrane conductance and g i = % i (D d /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 D d is the cable diameter. c s ¼ cD 2 s p and g s ¼ % m D 2 s p are the somatic membrane capacitance and leak conductance, respectively, with soma diameter D s . The exponential term with threshold slope factor Δ T and effective threshold voltage V T 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 V th from below (cf. spike mechanism of the 2C model). I s (t), I d (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 I s (t), I d (t) and a weak oscillatory field E(t), which do not elicit spikes. We consider that the somatic voltage evolves sufficiently below V T , 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; oÞ ¼ WðxÞÛ ðoÞ ¼ WðxÞ We obtain the system of differential equations where: indicates a (temporally) Fourier transformed variable. The solution of this system can In addition we calculate the somatic voltage time series in response to subthreshold input for a given initial condition V(x, 0) = V 0 (x). Using separation of variables and the Laplace transformṼ ðx; sÞ ¼ WðxÞŨ ðsÞ ¼ WðxÞ with complex variable s we obtain where: indicates a (temporally) Laplace transformed variable. We solve this system and obtaiñ where ±z(s) are the roots of the characteristic polynomial g i y 2 = c m s + g m of Eq (72), given by zðsÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The somatic voltage time series V(0, t) is then computed by inverse transformingṼ ð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 fitV s ðoÞ toV ð0; oÞ using Eqs (10) and (66) over a range of angular frequencies ω 2 [0, ω max ], requiring that the voltage values for ω = 0 (i.e. the steady states) match exactly. This constraint determines three parameters, where we have introduced the electrotonic length constant l ≔ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi g i =g m p . The remaining three subthreshold parameters C s , C d and G s are obtained using the method of least squares, with ω max /(2π) = 10 kHz. In the second step we determine the reset voltage V r 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 V s (t) to V(0, t) across the time interval t 2 [0, τ s ] with initial conditions V s (0) = V r , V d (0) = (G i V th + I d )/(G d + G i ) and Vð0; 0Þ ¼ V 0 r , Vðx; 0Þ ¼ I s coshððL À xÞ=lÞ þ I d ½coshðx=lÞ þ g s =ðlg m Þ� g s coshðL=lÞ þ lg m sinhðL=lÞ 0 < x < L ð82Þ for threshold somatic input I s = V th [g s + λg m tanh(L/λ)], I d = 0 and for threshold dendritic input I d = V th [g s cosh(L/λ) + λg m sinh(L/λ)], I s = 0 simultaneously using the method of least squares. τ s = C s /(G s + G i ) is the somatic membrane time constant of the 2C model. Note that we consider threshold input as constant current that yields V = V T (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 V s (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 G e = C s g s /c s (cf. Eqs (1) and (59)). The values for Δ T , V T and V th 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.

Spike coincidence measure
To quantify the similarity between the spike trains of the two-compartment and the ball-andstick model neurons we used the coincidence factor Γ defined by [62] where N c is the number of coincident spikes with precision (i.e., maximal temporal separation) Δ c , N BS and N 2C are the number of spikes in the spike trains of the ball-and-stick and the two-compartment models, respectively. hN c i = 2rΔ c N BS is the expected number of coincidences generated by a homogeneous Poisson process with spike rate r = N 2C /T as exhibited by the two-compartment model, where T is the duration of the spike train. The factor N ¼ 1 À 2rD 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.

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, because IN neurons do not exhibit an elongated spatial morphology compared with PY neurons. We used C IN = 0.2 nF and G IN = 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 x 2 {s, d, IN}. We consider large numbers of connections K x and reasonably small coupling strengths J x . 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 meanfield limit N ! 1) the overall synaptic input can be approximated by a mean part with