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 three-compartment 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 (DFP) that contributes to the local field potential (LFP) of a neural population. This work aligns and satisfies the widespread dipole assumption that is motivated by the "open-field" configuration of the DFP around cortical pyramidal cells. Our reduced three-compartment scheme allows to derive networks of leaky integrate-and-fire (LIF) 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. (2008), and conclude that our biophysically motivated approach yields substantial improvement.

1. 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).

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 ID causes dendritic field potential (DFP).

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 and 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 (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, 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 “two-monopole,” 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).

2. Materials and methods

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.

2.1. Theory

We describe the ith cortical pyramidal neuron (Figure 1) from population 1 via the electronic equivalent (reduced) three-compartment 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 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.

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.

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 EE. These potentials drive EPSC IE, essentially consisting of sodium ions, through the cell plasma with resistance RE 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) EI provide a shortcut between the excitatory branch and the trigger zone, where inhibitory postsynaptic currents (IPSC) II (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 RI.

The cell membrane at the axon hillock itself is represented by the branch at the right hand side. Here, a capacitor C reflects the temporary storage capacity of the membrane. The serial circuit consisting of a battery EM and a resistor RM 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 EE, EI. Therefore, we refer to our model as to a “reduced compartment model” here.

The three compartments are coupled through longitudinal resistors, RA, RB, RC, and RD where RA, RB denote the resistivity of the cell plasma and RC, RD that of extracellular space (Holt and Koch, 1999).

Finally, the membrane voltage at the axon hillock U (the dynamical state variable) and the DFP V, which measures the drop in electrical potential along the extracellular resistor RD are indicated. For the aim of calculation, the mesh currents ID (the dendritic current), IB (the basal current), and IIF (the integrate-and-fire current) are indicated.

The circuit in Figure 2 obeys the following equations: 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 IIF equals the dendritic current ID. If, by contrast, also the inhibitory synapses are active, the dendritic current ID 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 ID flowing through the extracellular space of resistance RD, gives rise to a large DFP V.

In order to simplify the following derivations, we gauge the resting potential (Equation 4) to EM = 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 IIF through Equation (8):

Division by 1 + gE (RA + RD) 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 IIF 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 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 RD the DFP with parameters 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 and therefore the spike rate is enhanced by EPSPs but diminished by IPSPs. On the other hand, the dendritic shunting current ID 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 sE|I(t) are excitatory and inhibitory synaptic impulse response functions, respectively, and R is the spike train coming from presynaptic neuron j, when spikes were emitted at times tν. The additional time constant τ 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 θ(t)] from below 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 + τ after interrupting the dynamics for a refractory period τ.

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 wE > 0 for excitatory synapses and wI < 0 for inhibitory synapses, respectively.

Following Mazzoni et al. (2008) an individual postsynaptic current IE|I at a synapse between neurons i and j obeys where τE|I are decay time constants and τE|I are rise time constants of EPSC and IPSC, respectively. Auxiliary variables are denoted by xE|I, while FE|I prescribes presynaptic forcing with spike train Equation (30). Here, J = vwE|I 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 F0.

Derivating Equation (35) and eliminating xE|I 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 sE|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 and λ2 = −1/τE|I, entailing the Green's functions with the Heaviside step function Θ(t).

The constants AE|I, BE|I > 0 are obtained from the initial conditions sE|I(t) = 0, reflecting causality, and a suitable normalization

The initial condition yields AE|I = BE|I ≡ SE|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 (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.

2.2. 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.

Table 1

Parameters laid as in Mazzoni et al. (.

Synaptic efficacies/mV	On interneurons	On pyramidal neurons
GABA	2.7	1.7
Recurrent cortical AMPA	0.7	0.42
External thalamic AMPA	0.95	0.55

From these, we compute the synaptic weights through and

Next, we determine the factors r by virtue of Equation (23) through using the inhibitory synaptic conductivity , correspondingly, Equation (22) allows us to express α in terms of the excitatory synaptic weights through

From α we can determine the total excitatory synaptic conductivities gE 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 by choosing τ = 20ms (Mazzoni et al., 2008).

In order to also determine the DFP parameters Equations (26–28), we finally compute the ratios

The remaining electrotonic parameters RM, RA, RB, RC, and RD 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).

Table 2

Resistivities of cell membrane, cell plasma and extracellular space.

Medium	ρ/Ω cm
Cell membrane (at axon hillock)	5 × 107
Cell plasma (cytoplasm)	200
Extracellular space	333

Parameters from Rall (.

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, R) and extracellular (R, R) resistances, respectively, according to Equation (53), we consider a simple implementation 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 = πa2. 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 . 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

2.3. 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 (45–48). Primarily, we compare our LFP measure L4, proposed as the average of DFP, with the Mazzoni et al. LFP L1 which is defined as the sum of absolute values of GABA and AMPA currents (Equation 44). Additionally, we also compare all possible measures, namely, mean membrane potential , Mazzoni et al. LFP L1, average of Mazzoni et al. DFP L2, sum of DFP L3, and the average of DFP L4.

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) is considered static with θ = 18 mV and the reset potential E = 11 mV. The refractory period for excitatory neurons is τ = 2 ms while for inhibitory neurons it is τ = 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 τ=1 ms and the rise and decay times are given by Table 3.

Table 3

Synaptic rise (τ.

Synaptic times	τr/ms	τd/ms
GABA	0.25	5
AMPA on interneurons	0.2	1
AMPA on pyramidal neurons	0.4	2

Parameters laid as in Mazzoni et al. (2008)

Moreover, synaptic efficacies, JE|I, 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, c0 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 c0 can range between 1.2 and 2.6 spikes/ms. The parameters of the Ornstein–Uhlenbeck process are τ = 16 ms and the standard deviation σ=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 weights . This implies evaluating the second order process (Equations 35, 36) with a different forcing term. Specifically, starting from IE|I(t) ≡ wE|I EE|I(t) = sE|I(t) * FE|I and pre-multiplying both sides with and subsequently re-arranging we obtain the desired forcing term . Note further that by expanding the term FE|I with Equation (37) and using Relation (49) we finally obtain .

3. 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.

Figure 3

Dynamics of the network and LFP comparisons: the three columns represent different runs of the network for three different rates, 1.2, 1.6, and 2.4 spikes/ms. In each column, all panels show the same 250 ms (extracted from 2 s simulations). The first panels (A–C) represent thalamic inputs with the different rates. The second panels (D–F) corresponds to a raster plot of the activity of 200 pyramidal neurons. The third panels (G–I) depict average instantaneous firing rate (computed on a 1 ms bin) of interneurons (blue) and fourth panels (J–L) correspond to average instantaneous firing rate of pyramidal neurons. The fifth panels (M–O) show the Mazzoni et al. LFP L1 from Equation (45). Finally, the last panels (P–R) depict our proposed LFP measure L4, which is the average of dendritic field potential (DFP) (Equation 48).

The focus is to compare our proposed measure L4, defined as mean of the DFP (Equation 48), with the Mazzoni et al. LFP L1 from Equation (45). In Figure 3 one sees two main striking differences between the two measures, namely in frequency and in amplitude. Specifically, L1 responds instantaneously to the spiking network activity by means of high frequency oscillations. Moreover, L1 also exhibits a large amplitude. In contrast, our mean DFP L4 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 L2 (Equation 48) and additionally the mean membrane potential (the standard considered in the neuroscientific literature). These are shown in Figure 4.

Figure 4

Comparison of different LFP measures when the network receives constant signal with three different rates (1.2, 1.6, and 2.4 spikes/ms). Again, only 250 ms is represented (extracted from 2 s simulation). The first panels (A–C) corresponding to the different rates shows the most widespread LFP measure used in the literature, namely average membrane potential . The second panels (D–F) shows the Mazzoni et al. LFP L1 from Equation (45). The third panels (G–I) displays the average of the Mazzoni et al. DFP L2 (Equation 46). Similarly, the fourth panels (J–L) shows the total, L3, (Equation 47) and the last panels (M–O) depicts the averaged, L4, (Equation 48) LFP measure. Note the different amplitude scales between measures.

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 L1 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 L1 or L2 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.

In Figure 5 we finally contrast the power spectra of the different LFP measures.

Figure 5

Comparison of power spectra of the various LFP measures when the network receives constant signal with three different rates (1.2, 1.6, and 2.4 spikes/ms). The first panels (A–C) corresponding to the different rates shows the power spectrum of the average membrane potential . The second panels (D–F) and third panels (G–I) show power spectra of the total and average of L1 and L2 corresponding to Mazzoni et al. (2008), respectively. The fourth panels (J–L) and the last panels (M–O) display power spectra of the L3 and L4 measures from our model, respectively. Note we show the full spectrum up to 5 kHz only for convenience due to the fine sample rate.

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, L1 and L2. In contrast, our LFP measures L3 and L4 fare better, and in particular, L4 decays at an approximately similar rate as the average membrane potential.

4. 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-to-spike 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, 2010, 2011) we aligned the parameters of our model with the model of Mazzoni et al. (2008) who approximated DFP as the sum of moduli of excitatory and inhibitory synaptic currents (Equation 44). From both approaches, we computed four different LFP estimates: L1, the sum of Mazzoni et al. DFP, L2, the population average of Mazzoni et al. DFP, L3 the sum of our dipole DFP, and L4 the population average of our dipole DFP (Equations 45–48).

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 ID 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 Ca2+ and a Ca2+-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 non-linear equations, with Hodgkin and Huxley (1952) type activation, but it will fall short from explicit and analytical observation equations.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.