A biophysical observation model for field potentials of networks of leaky integrate-and-fire neurons

We present a biophysical approach for the coupling of neural network activity as resulting from proper dipole currents of cortical pyramidal neurons to the electric field in extracellular fluid. Starting from a reduced threecompartment model of a single pyramidal neuron, we derive an observation model for dendritic dipole currents in extracellular space and thereby for the dendritic field potential that contributes to the local field potential of a neural population. This work aligns and satisfies the widespread dipole assumption that is motivated by the"open-field"configuration of the dendritic field potential around cortical pyramidal cells. Our reduced three-compartment scheme allows to derive networks of leaky integrate-and-fire models, which facilitates comparison with existing neural network and observation models. In particular, by means of numerical simulations we compare our approach with an ad hoc model by Mazzoni et al. [Mazzoni, A., S. Panzeri, N. K. Logothetis, and N. Brunel (2008). Encoding of naturalistic stimuli by local field potential spectra in networks of excitatory and inhibitory neurons. PLoS Computational Biology 4 (12), e1000239], and conclude that our biophysically motivated approach yields substantial improvement.


INTRODUCTION
Since Hans Berger's 1924 discovery of the human electroencephalogram (EEG) (Berger, 1929), neuroscientists achieved much progress in clarifying its neural generators (Creutzfeldt et al., 1966a,b;Nunez and Srinivasan, 2006;Schomer and Lopes da Silva, 2011). These are the cortical pyramidal neurons, as sketched in Figure 1, that possess a long dendritic trunk separating mainly excitatory synapses at the apical dendritic tree from mainly inhibitory synapses at the soma and at the perisomatic basal dendritic tree (Creutzfeldt et al., 1966a;Spruston, 2008). In addition, they exhibit an axial symmetry and are aligned in parallel to each other, perpendicular to the cortex' surface, thus forming a palisade of cell bodies and dendritic trunks. When both kinds of synapses are simultaneously active, inhibitory synapses generate current sources and excitatory synapses current sinks in extracellular space, hence causing the pyramidal cell to behave as a microscopic dipole surrounded by its characteristic electrical field, the dendritic field potential (DFP). The densely packed pyramidal cells form then a dipole layer whose superimposed currents give rise to the local field potential (LFP) of neural masses and eventually to the EEG (Nunez and Srinivasan, 2006;Lindén et al., 2010;Lindén et al., 2011;Schomer and Lopes da Silva, 2011).
Despite of the progress from experimental neuroscience, theoretically understanding the coupling of complex neural network dynamics to the electromagnetic field in the extracellular space poses challenging problems; some of them have been addressed to some extent by Bédard et al. (2004); Bédard and Destexhe (2009), and Bédard and Destexhe (2012).
In computer simulation studies, neural mass potentials, such as LFP and EEG are most realistically simulated by means of multicompartmental models (Protopapas et al., 1998;Sargsyan et al., 2001;Lindén et al., 2010;Lindén et al., 2011). Lindén et al. (2010) calculated the current dipole momentum of the DFP for single pyramidal and stellate cells, based on several hundreds compartments of the dendritic trees. Their results were in compliance with the standard dipole approximation of the electrostatic multipole expansion in the far-field (more than 1 mm remote from the dendritic trunk), but they found rather poor agreement with that approximation in the vicinity of the cell body. For comparison they also computed a "two-monopole" model of one synaptic current and its counterpart, the somatic return current, estimated from the current dipole momentum of the whole dendritic tree. This "two-monopole" model, which corresponds to an electrically equivalent single dipole model, obtained from the decomposition of the dendrite into two compartments, better approximates the true current dipole momentum in the vicinity of the pyramidal neuron. By superimposing the DFPs of pyramidal cells to the ensemble LFP, Lindén et al. (2011) found that LFP properties cannot be attributed to the far-field dipole approximation.
However, realistic multicompartmental models are computationally too expensive for large-scale neural network simulations. Therefore, various techniques have been proposed and employed to overcome computational complexity. These include networks of point models (i.e., devoid from any spatial representation), based on conductance models (Hodgkin and Huxley, 1952;Mazzoni et al., 2008), population density models FIGURE 1 | Sketch of a cortical pyramidal neuron with extracellular current dipole between spatially separated excitatory (open bullet) and inhibitory synapses (filled bullet). Neural in-and outputs are indicated by the jagged arrows. Dendritic current I D causes dendritic field potential (DFP). (Omurtag et al., 2000), or firing rate models (Wilson and Cowan, 1972), which can be seen as a sub class of population density models, with uniform density distribution (Chizhov et al., 2007). In these kinds of models, mass potentials such as LFP or EEG are conventionally described as averaged membrane potential. A different class of models are neural mass models (Jansen and Rit, 1995;Wendling et al., 2000;David and Friston, 2003;Rodrigues et al., 2010), where mass potentials are estimated either through sums (or actually differences) of excitatory postsynaptic potentials (EPSP) (David and Friston, 2003) or of excitatory postsynaptic currents (EPSC) (Mazzoni et al., 2008).
In particular, the model of Mazzoni et al. (2008) which is based on Brunel and Wang (2003), recently led to a series of follow-up studies (Mazzoni et al., 2010(Mazzoni et al., , 2011 addressing the correlations between numerically simulated and experimentally measured LFP/EEG with spike rates by means of statistical modeling and information theoretic measures. In all of the above point models and their extension to population models, it is assumed that the extracellular space is iso-potential and the majority of studies thereby neglect the effect of extracellular resistance. That is, the extracellular space constitutes a different and isolated domain with no effect on neuronal dynamics. In this article we extend the ad hoc model of Mazzoni et al. (2008) toward a biophysically better justified approach, taking the dipole character of extracellular currents and fields into account. Basically, our model corresponds to the "twomonopole," or, equivalent dipole model of Lindén et al. (2010) which gave a good fit of the DFP close to the cell body of a cortical pyramidal neuron. However, we aim to keep the simplicity of the Mazzoni et al. (2008) model in terms of computational complexity, by endowing the extracellular space with resistance and by keeping point-like neuronal circuits. That is, in our case we do not quite consider point neurons, nor spatially extended models with detailed compartmental morphology, yet an intermediate level of description is achieved. To this end we propose a reduced three-compartmental model of a single pyramidal neuron (Destexhe, 2001;Wang et al., 2004;beim Graben, 2008), and derive an observation model for the dendritic dipole currents in the extracellular space and thereby for the DFP that contributes to the LFP of a neural population. Interestingly, our reduced three-compartmental model enables us to derive a leaky integrate-and-fire (LIF) mechanism [as for a point model (Mazzoni et al., 2008)], with additional observation equations for the DFP, which all together allows to study the relationship between spike rates and LFP. Our derivations also nicely map realistic electrotonic parameters to phenomenological parameters considered in Mazzoni et al. (2008). Mazzoni et al. (2008) consider three populations of neurons, namely excitatory cortical pyramidal cells (population 1), inhibitory cortical interneurons (population 2), and excitatory thalamic relay neurons (population 3), passing sensory input to the cortex that is simulated by a random (Erdős-Rényi) graph of K = 4000 pyramidal and L = 1000 interneurons with connection probability P = 0.2.

THEORY
We describe the ith cortical pyramidal neuron (Figure 1) from population 1 via the electronic equivalent (reduced) threecompartment model (Figure 2) (Destexhe, 2001;Wang et al., 2004;beim Graben, 2008), which is parsimonious to derive our observation model: one compartment for the apical dendritic FIGURE 2 | Proposed electronic equivalent circuit for a pyramidal neuron (reduced three-compartmental model). Note that the apical and basal dendrites are not true compartments since capacitors are not explicitly represented, rather, these are implicitly taken into account via EPSP and IPSP static functions, thus keeping computational complexity low.
tree, another one for soma and perisomatic basal dendritic tree (Lindén et al., 2010), and the third-actually a LIF unit-for the axon hillock where membrane potential is converted into spike trains by means of an integrate-and-fire mechanism.
Excitatory synapses are represented by the left-most branch, where EPSP at a synapse between a neuron j from population 1 or 3 and neuron i act as electromotoric forces E E ij . These potentials drive EPSC I E ij , essentially consisting of sodium ions, through the cell plasma with resistance R E ij from the synapse toward the axon hillock.
The middle branch describes the inhibitory synapses between a neuron k from population 2 and neuron i. Here, inhibitory postsynaptic potentials (IPSP) E I ik provide a shortcut between the excitatory branch and the trigger zone, where inhibitory postsynaptic currents (IPSC) I I ik (essentially chloride ions) close the loop between the apical and perisomatic dendritic trees. The resistivity of the current paths along the cell plasma is given by R I ik . The cell membrane at the axon hillock itself is represented by the branch at the right hand side. Here, a capacitor C i reflects the temporary storage capacity of the membrane. The serial circuit consisting of a battery E M and a resistor R M denotes the Nernst resting potential and the leakage conductance of the membrane, respectively (Johnston and Wu, 1997). Finally, a spike generator (Hodgkin and Huxley, 1952;Mazzoni et al., 2008) (indicated by a "black box") is regarded of having infinite input impedance. Both, EPSP and IPSP result from the interaction of postsynaptic receptor kinetics with dendritic low-pass filtering in compartments one and two, respectively (Destexhe et al., 1998;Lindén et al., 2010). Hence the required capacitances, omitted in Figure 2, are already taken into account by E E ij , E I ik . Therefore, we refer to our model as to a "reduced compartment model" here.
The three compartments are coupled through longitudinal resistors, i denote the resistivity of the cell plasma and R C i , R D i that of extracellular space (Holt and Koch, 1999).
Finally, the membrane voltage at the axon hillock U i (the dynamical state variable) and the DFP V i , which measures the drop in electrical potential along the extracellular resistor R D i are indicated. For the aim of calculation, the mesh currents I D i (the dendritic current), I B i (the basal current), and I IF i (the integrate-and-fire current) are indicated.
The circuit in Figure 2 obeys the following equations: (1) Here, p is the number of excitatory and q is the number of inhibitory synapses connected to neuron i.

The circuit described by Equations
(1-7) shows that the neuron i is likely to fire when the excitatory synapses are activated. Then, the integrate-and-fire current I IF i equals the dendritic current I D i . If, by contrast, also the inhibitory synapses are active, the dendritic current I D i is shunted between the apical and perisomatic basal dendritic trees and only a portion could evoke spikes at the trigger zone (Equation 4). On the other hand, the large dendritic current I D i flowing through the extracellular space of resistance R D i , gives rise to a large DFP V i . In order to simplify the following derivations, we gauge the resting potential (Equation 4) to E M = 0, yielding From Equation (5) we obtain the individual EPSC's as And accordingly, the individual IPSC's from Equation (6) Inserting Equation (9) into Equation (1) yields the excitatory dendritic current where we have introduced the excitatory dendritic conductivity Likewise we obtain the inhibitory dendritic currents from Equations (2) and (10) as with the inhibitory dendritic conductivity With these results, we obtain an interface equation for an observation model as follows. Rearranging Equation (11) yields Next, we eliminate I IF i through Equation (8): gives the desired expression for the extracellular dendritic dipole current: with the following electrotonic parameters In order to derive the evolution equation we consider the integrate-and-fire current I IF i that is given through Equation (3). The individual EPSCs and IPSCs have already been obtained in Equations (9) and (10), respectively. Inserting Equation (13) into Equation (3) yields Next we insert our interface equation Equation (16) and also Equation (8): and obtain after some rearrangements and after multiplication with Observation models for field potentials the dynamical law for the membrane potential at axon hillock: where we have introduced the following parameters: • time constants • excitatory synaptic weights • inhibitory synaptic weights Using the result Equation (20), we can also eliminate the temporal derivative in the interface equation Equation (16) through which yields And eventually, by virtue of Equation (7) after multiplication with with parametersw The change in sign of the inhibitory contribution from Equation (20) to Equation (25) has an obvious physical interpretation: In Equation (20), the change of membrane potential U i and therefore the spike rate is enhanced by EPSPs but diminished by IPSPs. On the other hand, the dendritic shunting current I D i in Equation (25) is large for both, large EPSPs and large IPSPs.
From Equation (20) we eventually obtain the neural network's dynamics by taking into account that postsynaptic potentials are obtained from presynaptic spike trains through temporal convolution with postsynaptic impulse response functions, i.e., where s E|I i (t) are excitatory and inhibitory synaptic impulse response functions, respectively, and R j is the spike train coming from presynaptic neuron j, when spikes were emitted at times t ν . The additional time constant τ L is attributed to synaptic transmission delay (Mazzoni et al., 2008). These events are obtained by integrating Equation (20) with initial condition where E is some steady-state potential (Mazzoni et al., 2008). If at time t = t ν the membrane reaches a threshold [with possibly a time-dependent activation threshold θ i (t)] from below dU i (t) dt > 0 then an output spike δ(t − t ν ) is generated, which is then followed by a potential resetting as follows Additionally, the integration of the dynamical law is restarted at time t = t ν+1 + τ rp after interrupting the dynamics for a refractory period τ rp . Inserting Equation (29) into Equation (20) entails the evolution equation of the neural network where the signs had been absorbed by the synaptic weights, such that w E ij > 0 for excitatory synapses and w I ik < 0 for inhibitory synapses, respectively.
Following Mazzoni et al. (2008) an individual postsynaptic current I E|I ij at a synapse between neurons i and j obeys with spike train Equation (30). Here, J ij = vw E|I ij denotes synaptic gain with v = 1 mV as voltage unit.
Note that Equation (37) is essentially a weighted sum of delta functions, such that a single spike can be assumed as particular forcing with some constant F 0 . Derivating Equation (35) and eliminating x E|I ij transforms Equations (35, 36) into a linear second-order differential equation with constant coefficients Equation (39) with the particular forcing Equation (38) is solved by a Green's function s E|I i (t) such that the general solution of Equation (39) is obtained as the temporal convolution For t = 0, Equation (39) assumes its homogeneous form and is easily solved by means of the associated characteristic polynomial with roots λ 1 = −1/τ E|I d and λ 2 = −1/τ E|I r , entailing the Green's functions with the Heaviside step function (t).
The constants A E|I , B E|I > 0 are obtained from the initial conditions s E|I i (t) = 0, reflecting causality, and a suitable normalization The initial condition yields A E|I = B E|I ≡ S E|I , while the remaining constant due to normalization. Therefore, the normalized Green's functions are those of Brunel and Wang (2003) Now, we are able to compare our DFP V i (Equation 25) with the estimate of Mazzoni et al. (2008) which is given (in our notation) as the sums of the moduli of excitatory and inhibitory synaptic currents, i.e., where "MPLB" refers to the authors Mazzoni et al. (2008). From Equations (25) and (44), respectively, we compute two models of the LFP. First, by summing DFP across all pyramidal neurons (beim Graben and Kurths, 2008;Mazzoni et al., 2008), and, second by taking the DFP average (Nunez and Srinivasan, 2006), which yields where K is number of pyramidal neurons.

PARAMETER ESTIMATION
Next, we relate the electrotonic parameters of our model to the phenomenological parameters of Mazzoni et al. (2008). To this end, we first report their synaptic efficacies in Table 1. From these, we compute the synaptic weights through and Next, we determine the factors r i by virtue of Equation (23) through using the inhibitory synaptic conductivityḡ GABA = 1 nS, correspondingly, Equation (22) allows us to express α ij in terms of the From α ij we can determine the total excitatory synaptic conductivities g E i according to Equation (17) through and hence Inserting next Equation (18) into Equation (21) yields Equation (52) could constraint the choice of the membrane capacitance C i by choosing τ i = 20 ms (Mazzoni et al., 2008).
In order to also determine the DFP parameters Equations (26-28), we finally compute the ratios

The remaining electrotonic parameters R
i are estimated from cell geometries as follows. The resistance R of a volume conductor is proportional to its length and reciprocally proportional to its cross-section A, i.e., where ρ is the (specific) resistivity of the medium. Table 2 shows the resistivities of the three kinds of interest which then allows to evaluate the various volume conductor resistances according to Equation (53). We consider a total dendritic length of 2 = 20 μm and a dendritic radius of a = 7 μm, that are generally subjected to variation. Equally, parameters that were allowed to vary are the length and radius of the axon hillock, yet herein we consider a length of 2 = 20 μm and radius of a = 0.5 μm (Mainen et al., 1995;Destexhe, 2001;Kole and Stuart, 2012). To evaluate the intracellular (R A , R B ) and extracellular (R D , R C ) resistances, respectively, according to Equation (53), we consider a simple implementation

Extracellular space 333
Parameters from Rall (1977); Mainen et al. (1995); Kole and Stuart (2012), and Gold et al. (2007). Note that the resistivity of the cell membrane has to be related to the constant membrane thickness (≈10 nm).
where the length is half of the dendritic length (i.e., basal and apical length are symmetrical, but this can be broken). However, the cross sectional area for the cytoplasm is simply A = πa 2 . Finally, the area of the axon hillock is simply the surface area of a cylinder.
In order to also determine the cross-section of extracellular space between dendritic trunks we make the following approximations. We assume that dendritic trunks are parallel aligned cylinders of radius a and length that are hexagonally dense packed. Then the centers of three adjacent trunks form an equilateral triangle with side length 2a and hence area 2 √ 3a 2 . The enclosed space is then given by the difference of the triangle area and the area of three sixth circle sectors, therefore Hence, the cross-section of extracellular space surrounding one trunk is

SIMULATIONS
Subsequently, we implement an identical network to the one considered by Mazzoni et al. (2008) with Brian Simulator, that is a Python-based environment (Goodman and Brette, 2009). However, the derivations from the previous section enables the possibility of setting a dipole observable that measures the local DFP on each pyramidal neurons, given by Equation (25). This allows then to define a mesoscopic LFP observable, which can be equated either as averaged DFP or simply given as the sum of DFP, given by Equations ( For completeness, we briefly summarize the description of the network [we refer the reader to Mazzoni et al. (2008) for details]. The network models a cortical tissue with LIF neurons, composed of 1000 inhibitory interneurons and 4000 pyramidal neurons, which are described by the evolution Equation (34). The threshold crossings given by Equation (32)  θ i = 18 mV and the reset potential E = 11 mV. The refractory period for excitatory neurons is τ rp = 2 ms while for inhibitory neurons it is τ rp = 1 ms. The network connectivity is random and sparse with a 0.2 probability of directed connection between any pair of neurons. The evolution of synaptic currents, fast GABA (inhibitory) and AMPA (excitatory) are described via the second order evolution Equations (35, 36), which are activated by incoming presynaptic spikes represented by Equation (30). The latency of the postsynaptic currents is set to τ L = 1 ms and the rise and decay times are given by Table 3. Moreover, synaptic efficacies, J E|I ij , for simulation were presented in Table 1. Note that Relation (49) then allows to determine the synaptic weights. Additionally, all neurons receive external thalamic excitatory inputs, that is, via AMPA-type synapses, which are activated by random Poisson spike trains, with a time varying rate that is identical for all neurons. Specifically, the thalamic inputs are the only source of noise, which attempts to account for both cortical heterogeneity and spontaneous activity. This is achieved by modeling a two level noise, where the first level is an Ornstein-Uhlenbeck process superimposed with a constant signal and the second level is a time varying inhomogeneous Poisson process. Thus, we have the following time varying rate, λ(t), that feeds into inhomogeneous Poisson process: where η(t) represents Gaussian white noise, c 0 represents a constant signal (but equally could be periodic or other), and the operation [·] is the threshold-linear function, [x] + = x if x > 0, [x] + = 0 otherwise, which circumvents negative rates. The constant signal c 0 can range between 1.2 and 2.6 spikes/ms. The parameters of the Ornstein-Uhlenbeck process are τ n = 16 ms and the standard deviation σ n = 0.4 spikes/ms. For complete exposition, we note that from an implementation viewpoint (within the Brian simulator), a copy of the postsynaptic impulse response function (Equation 29) has to be evaluated to calculate the DFP (Equation 25) with weightsw E|I ij . This implies evaluating the second order process (Equations 35, 36) with a different forcing term. Specifically, and premultiplying both sides withw E|I ij and subsequently re-arranging we obtain the desired forcing termF Parameters laid as in Mazzoni et al. (2008). further that by expanding the term F E|I ij with Equation (37) and using Relation (49) we finally obtainF E|I ij =w E|I ij τ i vR j (t).

RESULTS
Following Mazzoni et al. (2008), the network simulations are run for 2 s with three different noise levels, specifically, receiving a constant signal with three different rates 1.2, 1.6, and 2.4 spikes/ms as depicted in Figure 3. Note that these input rates do not mean that a single neuron fires at these high rates. Rather, it can be obtained from multiple neurons that jointly fire with slower, yet desynchronized, rates converging at the same postsynaptic cell. The Poisson process ensures that this is well represented.
The focus is to compare our proposed measure L 4 , defined as mean of the DFP (Equation 48), with the Mazzoni et al. LFP L 1 from Equation (45). In Figure 3 one sees two main striking differences between the two measures, namely in frequency and in amplitude. Specifically, L 1 responds instantaneously to the spiking network activity by means of high frequency oscillations. Moreover, L 1 also exhibits a large amplitude. In contrast, our mean DFP L 4 measures comparably to experimental LFP, that is, in the order of millivolts, and although it responds to population activity, it has a relatively smoother response. Actually one can realize that our LFP estimate represents low-pass filtered thalamic input.
The physiological relevance of this is not yet clear in our work. However, recent work (Poulet et al., 2012) shows that desynchronized cortical state during active behavior is driven by a centrally generated increase in thalamic action potential firing (i.e., thalamic firing controls cortical states). Thus, it seems that cortical synchronous activity is suppressed when thalamic input increases, thereby suggesting that cortical desynchronized states to be related to sensory processing. This work further quantifies these observations by applying Fast Fourier Transform (FFT) to cortical EEG and subsequently comparing with thalamic firing rate by means of Pearson correlation coefficient. Unfortunately they do not quantify the amount of thalamic oscillations contained within the cortical EEG.
Yet, to keep a comparable comparison between measures, we also compute the average of the Mazzoni et al. DFP L 2 (Equation 48) and additionally the mean membrane potential (the standard considered in the neuroscientific literature). These are shown in Figure 4.
Clearly, in terms of time profile, the summed and averaged observables are similar within the same class of LFP measures. However, in all cases the Mazzoni et al. LFP L 1 exhibits a significantly larger order of magnitude, which diverges substantially from experimental LFP amplitudes, typically varying between 0.5 and 2 mV (Lakatos et al., 2005;Niedermeyer, 2005). In contrast, although the mean DFP is not contained within the interval from 0.5 to 2 mV it arguably performs better. However, we do concede further work is required. Some gains in improving the different LFP measures can be achieved by applying for example a weighted average, which would mimic the distance of an electrode to a particular neuron by means of a lead field kernel (Nunez and Srinivasan, 2006). For example, a convolution of either L 1 or L 2 with a Gaussian kernel (representing the distance to a neuron), would yield a measure that captures better the LFP or better the DFP of the nearest neurons. However, further work will be required to properly quantify the gain when space is taken into account.

Frontiers in Computational
In Figure 5 we finally contrast the power spectra of the different LFP measures.
One interesting feature is that the power spectrum of the Mazzoni et al. LFP measures decays much more slowly that the average membrane potential for higher frequencies. This observation is true for both, L 1 and L 2 . In contrast, our LFP measures L 3 and L 4 fare better, and in particular, L 4 decays at an approximately similar rate as the average membrane potential.

DISCUSSION
In this article we derived a model for cortical dipole fields, such as DFP/LFP from biophysical principles. To that aim we decomposed a cortical pyramidal cell, the putative generator of those potentials, into three compartments: the apical dendritic tree as the place of mainly excitatory (AMPA) synapses, the soma and the perisomatic dendritic tree as the place of mainly inhibitory (GABA) synapses, and the axon hillock as the place of wave-tospike conversion by means of an integrate-and-fire mechanism. From Kirchhoff 's laws governing an electronic equivalent circuit of our model, we were then able to derive the evolution equation for neural network activity (Equation 34) and, in addition, an observation equation (25) for the dendritic dipole potential contributing to the LFP of a cortical population.
In order to compare our approach with another model discussed in the recent literature (Mazzoni et al., 2008(Mazzoni et al., , 2010(Mazzoni et al., , 2011 we aligned the parameters of our model with the model of Mazzoni et al. (2008)  Our results indicate two main effects between our dipole LFP measures and those of Mazzoni et al. Firstly, the measures based on Mazzoni et al. (2008) systematically overestimate LFP amplitude by almost one order of magnitude. One reason for that could be attributed to the direct conversion of synaptic current into voltage without taking extracellular conductivity into account, as properly done in our approach. Yet, another, even more crucial reason is disclosed by our equivalent circuit (Figure 2). In our approach there is just one extracellular current I D flowing from the perisomatic to the apical dendritic tree. In the model of Mazzoni et al. (2008), however, two synaptic currents that might be of the same order of magnitude are superimposed to the DFP. Secondly, the measures based on Mazzoni et al. (2008) also systematically overestimate LFP frequencies. This could probably be attributed partly to spurious higher harmonics introduced by computing absolute values. Moreover, taking the power spectrum shows that the Mazzoni et al. (2008) measure decays much more slowly than the average membrane potential, which is at variance with experimental data.
However, at the current stage, both models, that of Mazzoni et al. (2008) and our own, agree with respect to the polarity of DFP and LFP. The measures based on Mazzoni et al. (2008) have positive polarity simply due to the moduli. On the other hand, also the direction of current dipoles in our model is constrained by the construction of the equivalent circuit (Figure 2) where current sources are situated at the perisomatic and current sinks are situated at apical dendritic tree. Taking this polarity as positive also entails positive DFP and LFP that could only change in strength. However, it is well known from brain anatomy that pyramidal cells appear in at least two layers, III and VI, of neocortex. This is reflected in experiments when an electrode traverses different layers by LFP polarity reversals, and, of course, by the fact that LFP and EEG oscillate between positive and negative polarity. Adapting our model to this situation could be straightforwardly accomplished in the framework of neural field theory by fully representing space and simulating layered neural fields (Amari, 1977;Jirsa and Haken, 1996;beim Graben, 2008). By contrast such a generalization is impossible at all with the model of Mazzoni et al. (2008) due to the presence of absolute values.
On theses grounds we have good indication that our measure is an improvement to the Mazzoni et al. LFP measures, and, quite importantly, it is biophysically better motivated than the ad hoc model of Mazzoni et al. (2008). However, much considerable effort is still required to underpin all the relevant LFP mechanisms and to better represent experimental LFP/EEG dynamics.
Finally, our work provides a new framework where DFPs and the relationship between firing rates and local fields can be explored without the extreme demand on computational complexity involved in multicompartmental modeling (Protopapas et al., 1998;Sargsyan et al., 2001;Lindén et al., 2010;Lindén et al., 2011) by adopting reduced compartment circuits. For example, we envisage to extend our recent work which maps firing rate model (derived from LIF models) to population density models (Chizhov et al., 2007), but now incorporating our observational DFP model. In addition, our framework is analytically amenable and thus can be applied to any linear differential equation, for instance, GIF (Gif-sur-Yvette Integrate Fire) models, which are improvements to the LIF models and compute more accurately spike activations (Rudolph-Lilith et al., 2012). Also resonant membranes (mediated by Ca 2+ and a Ca 2+ -activated K + ionic currents) that describe sub-threshold oscillations and which can be easily expressed by linear equations (Mauro et al., 1970) can be incorporated in our derivations. We note, however, that our framework can be applied to nonlinear equations, with Hodgkin and Huxley (1952) type activation, but it will fall short from explicit and analytical observation equations.