ORIGINAL RESEARCH article
A biophysical observation model for field potentials of networks of leaky integrate-and-fire neurons
- 1Bernstein Center for Computational Neuroscience Berlin, Berlin, Germany
- 2Department of German Language and Linguistics, Humboldt-Universität zu Berlin, Berlin, Germany
- 3Centre for Robotics and Neural Systems, School of Computing and Mathematics, University of Plymouth, Plymouth, UK
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.
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.
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 EEij. These potentials drive EPSC IEij, essentially consisting of sodium ions, through the cell plasma with resistance REij 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) EIik provide a shortcut between the excitatory branch and the trigger zone, where inhibitory postsynaptic currents (IPSC) IIik (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 RIik.
The cell membrane at the axon hillock itself is represented by the branch at the right hand side. Here, a capacitor Ci 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 EEij, EIik. Therefore, we refer to our model as to a “reduced compartment model” here.
The three compartments are coupled through longitudinal resistors, RAi, RBi, RCi, and RDi where RAi, RBi denote the resistivity of the cell plasma and RCi, RDi that of extracellular space (Holt and Koch, 1999).
Finally, the membrane voltage at the axon hillock Ui (the dynamical state variable) and the DFP Vi, which measures the drop in electrical potential along the extracellular resistor RDi are indicated. For the aim of calculation, the mesh currents IDi (the dendritic current), IBi (the basal current), and IIFi (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 IIFi equals the dendritic current IDi. If, by contrast, also the inhibitory synapses are active, the dendritic current IDi 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 IDi flowing through the extracellular space of resistance RDi, gives rise to a large DFP Vi.
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 IIFi through Equation (8):
Division by 1 + gEi (RAi + RDi) 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 IIFi 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
And eventually, by virtue of Equation (7) after multiplication with RDi the DFP
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 Ui and therefore the spike rate is enhanced by EPSPs but diminished by IPSPs. On the other hand, the dendritic shunting current IDi 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|Ii(t) are excitatory and inhibitory synaptic impulse response functions, respectively, and Rj 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 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 wEij > 0 for excitatory synapses and wIik < 0 for inhibitory synapses, respectively.
Following Mazzoni et al. (2008) an individual postsynaptic current IE|Iij at a synapse between neurons i and j obeys
where τE|Id are decay time constants and τE|Ir are rise time constants of EPSC and IPSC, respectively. Auxiliary variables are denoted by xE|Iij, while FE|Iij prescribes presynaptic forcing
with spike train Equation (30). Here, Jij = vwE|Iij 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|Iij 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|Ii(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|Id and λ2 = −1/τE|Ir, 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|Ii(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 Vi (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
Table 1. Parameters laid as in Mazzoni et al. (2008).
From these, we compute the synaptic weights through
Next, we determine the factors ri by virtue of Equation (23) through
using the inhibitory synaptic conductivity , correspondingly, Equation (22) allows us to express αij in terms of the excitatory synaptic weights through
From αij we can determine the total excitatory synaptic conductivities gEi according to Equation (17) through
Inserting next Equation (18) into Equation (21) yields
Equation (52) could constraint the choice of the membrane capacitance Ci by choosing τi = 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 RMi, RAi, RBi, RCi, and RDi 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 (RA, RB) and extracellular (RD, RC) 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
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 θ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, JE|Iij, 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 τ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 weights . This implies evaluating the second order process (Equations 35, 36) with a different forcing term. Specifically, starting from IE|Iij(t) ≡ wE|Iij EE|Iij(t) = sE|Ii(t) * FE|Iij and pre-multiplying both sides with and subsequently re-arranging we obtain the desired forcing term . Note further that by expanding the term FE|Iij with Equation (37) and using Relation (49) we finally obtain .
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.
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.
We thank Michelle Lilith, Claude Bédard, Alain Destexhe, and Jürgen Kurths for fruitful discussion. In addition, we would like to thank Samantha Adams for providing help with Brian Simulator installations and initial discussions of Brian usage. This research was supported by a DFG Heisenberg grant awarded to PbG (GR 3711/1-1).
Bédard, C., and Destexhe, A. (2012). “Modeling local field potentials and their interaction with the extracellular medium,” in Handbook of Neural Activity Measurement, eds R. Brette and A. Destexhe (Cambridge, MA: Cambridge University Press), 136–191.
beim Graben, P. (2008). “Foundations of neurophysics,” in Lectures in Supercomputational Neuroscience: Dynamics in Complex Brain Networks, Springer Complexity Series, eds P. b. Graben, C. Zhou, M. Thiel, and J. Kurths (Berlin: Springer), 3–48.
Brunel, N., and Wang, X.-J. (2003). What determines the frequency of fast network oscillations with irregular neural discharges? I. synaptic dynamics and excitation-inhibition balance. J. Neurophysiol. 90, 415–430.
Creutzfeldt, O. D., Watanabe, S., and Lux, H. D. (1966a). Relations between EEG phenomena and potentials of single cortical cells. I. evoked responses after thalamic and epicortical stimulation. Electroencephalogr. Clin. Neurophysiol. 20, 1–18.
Creutzfeldt, O. D., Watanabe, S., and Lux, H. D. (1966b). Relations between EEG phenomena and potentials of single cortical cells. II. spontaneous and convulsoid activity. Electroencephalogr. Clin. Neurophysiol. 20, 19–37.
Destexhe, A., Mainen, F., and Sejnowski, T. J. (1998). “Kinetic models of synaptic transmission,” in Methods in Neuronal Modelling. From Ions to Networks, eds C. Koch and I. Segev (Cambridge, MA: MIT Press), 1–25.
Lakatos, P., Shah, A., Knuth, K., Ulbert, I., Karmos, G., and Schroeder, C. (2005). An oscillatory hierarchy controlling neuronal excitability and stimulus processing in the auditory cortex. J. Neurophysiol. 94, 1904–1911.
Mazzoni, A., Panzeri, S., Logothetis, N. K., and Brunel, N. (2008). Encoding of naturalistic stimuli by local field potential spectra in networks of excitatory and inhibitory neurons. PLoS Comput. Biol. 4:e1000239. doi: 10.1371/journal.pcbi.1000239
Mazzoni, A., Whittingstall, K., Brunel, N., Logothetis, N. K., and Panzeri, S. (2010). Understanding the relationships between spike rate and delta/gamma frequency bands of LFPs and EEGs using a local cortical network model. Neuroimage 52, 956–972.
Niedermeyer, E. (2005). “The normal EEG of the waking adult,” in Electroencephalography: Basic Principles, Clinical Applications, and Related Fields, 5th Edn, eds E. Niedermeyer and F. L. D. Silva (Philadelphia, PA: Lippincott Williams & Wilkins), 167–192.
Protopapas, A., Vanier, M., and Bower, J. M. (1998). “Simulating large networks of neurons,” in Methods in Neuronal Modelling. From Ions to Networks, eds C. Koch and I. Segev (Cambridge, MA: MIT Press), 461–498.
Rall, W. (1977). “Core conductor theory and cable properties of neurons,” in Handbook of Physiology – The Nervous System, Cellular Biology of Neurons, Vol. 1, ed E. R. Kandel (Bethesda, MD: American Physiological Society), 39–97.
Rudolph-Lilith, M., Dubois, M., and Destexhe, A. (2012). Analytical integrate-and-fire neuron models with conductance-based dynamics and realistic postsynaptic potential time course for event-driven simulation strategies. Neural Comput. 34, 1426–1461.
Sargsyan, A. R., Papatheodoropoulos, C., and Kostopoulos, G. K. (2001). Modeling of evoked field potentials in hippocampal CA1 area describes their dependence on NMDA and GABA receptors. J. Neurosci. Methods 104, 143–153.
Schomer, D. L., and Lopes da Silva, F. H. (eds.). (2011). Niedermayer's Electroencephalography. Basic Principles, Clinical Applications, and Related Fields. 6th Edn. Philadelphia, PA: Lippincott Williams & Wilkins.
Wang, X.-J., Tegnér, J., Constantinidis, C., and Goldman-Rakic, P. S. (2004). Division of labor among distinct subtypes of inhibitory neurons in a cortical microcircuit of working memory. Proc. Natl. Acad. Sci. U.S.A. 101, 1368–1373.
Keywords: biophysics, neural networks, leaky integrate-and-fire neuron, current dipoles, extracellular medium, field potentials
Citation: beim Graben P and Rodrigues S (2013) A biophysical observation model for field potentials of networks of leaky integrate-and-fire neurons. Front. Comput. Neurosci. 6:100. doi: 10.3389/fncom.2012.00100
Received: 03 September 2012; Accepted: 14 December 2012;
Published online: 04 January 2013.
Edited by:Dimitris Pinotsis, University College London, UK
Reviewed by:Gaute T. Einevoll, Norwegian University of Life Sciences, Norway
Julien Modolo, Lawson Health Research Institute, Canada
Copyright © 2013 beim Graben and Rodrigues. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.
*Correspondence: Peter beim Graben, Department of German Language and Linguistics, Humboldt-Universität zu Berlin, Unter den Linden 6, D–10099, Berlin, Germany. e-mail: email@example.com