Analysis and Modeling of Subthreshold Neural Multi-Electrode Array Data by Statistical Field Theory

Multi-electrode arrays (MEA) are increasingly used to investigate spontaneous neuronal network activity. The recorded signals comprise several distinct components: Apart from artifacts without biological significance, one can distinguish between spikes (action potentials) and subthreshold fluctuations (local fields potentials). Here we aim to develop a theoretical model that allows for a compact and robust characterization of subthreshold fluctuations in terms of a Gaussian statistical field theory in two spatial and one temporal dimension. What is usually referred to as the driving noise in the context of statistical physics is here interpreted as a representation of the neural activity. Spatial and temporal correlations of this activity give valuable information about the connectivity in the neural tissue. We apply our methods on a dataset obtained from MEA-measurements in an acute hippocampal brain slice from a rat. Our main finding is that the empirical correlation functions indeed obey the logarithmic behavior that is a general feature of theoretical models of this kind. We also find a clear correlation between the activity and the occurrence of spikes. Another important insight is the importance of correctly separating out certain artifacts from the data before proceeding with the analysis.

1. Introduction

The multi-electrode arrays (MEA) system is becoming an increasingly important tool for investigations of neural activity, both in ex vivo brain tissue (e.g., a hippocampal slice preparation from rat or mouse Egert et al., 2002) and in in vitro neuronal cultures (e.g., from embryonic rodent brain tissue Illes et al., 2014 or human stem cells Heikkilä et al., 2009). This technology permits simultaneous long-term recordings from a fairly large number of extra-cellular electrodes. See e.g., Nam and Wheeler (2011) and Spira and Hai (2013) for general reviews of multi-electrode array technology.

Each electrode records alterations of the field potential elicited by spike activity of one or a few neurons in close vicinity of it. Extracellular spikes have an amplitude of 10–500 μV, and are considered as a manifestation of the intracellular action potential which has a much higher amplitude of 100 mV. Many methods have been developed for the detection and sorting of spike events (see e.g., Cotterill et al., 2016, for a recent review), and analysis of the statistical properties of spike trains is one of the major modes of investigating neural activity (see e.g., Rieke et al., 1997, for a pedagogical introduction to this field).

The focus of the present paper will however not be on the spikes, but rather on the subthreshold behavior of the potential during interspike intervals. This is often referred to as the local field potential (LFP). The extraction of these subthreshold fluctuations out of MEA-data streams which are contaminated with artifacts and spiking activity (see e.g., Waldert et al., 2013), as well as to develop the appropriate mathematical approaches applicable to describe their properties, are challenging issues in the research field of neuroscience.

The LFP is usually assumed to be a superposition of contributions from the neural activity in a fairly large neighborhood of an electrode and also gets modified by other types of cells than neurons. Typically it has an amplitude of a few μV, i.e., much less than the spikes. In contrast to the rather stereotyped spike waveforms from individual neurons, the LFP has a much more stochastic appearance, reflecting its origins from a large and rather heterogenous ensemble of neurons (Illes et al., 2014). Furthermore, spikes can only be recorded in relatively close vicinity to the recording electrode while the distance between recording electrode and LFP source can be several micrometers and millimeters (Destexhe and Bedard, 2013). Thereby, the LFP recorded in isolated brain structures actually carries important information about the transport properties of the intra-cellular medium (Bedard and Destexhe, 2015), the spatial and temporal structure of the neural activity (Destexhe and Bedard, 2013; Linden et al., 2014), and also the connectivity of neurons within brain tissues (Reichinnek et al., 2010). There are different functional aspects of neuronal circuits which can be revealed by modeling and analyzing LFP. Current-source density analysis is used to reveal the neural source of recorded LFP (Ness et al., 2015) which is still controversial (Riera et al., 2012). Since the pioneering work of Berger in the 1920s, LFP band-separations techniques has been devoted to analyzing the LFP in the frequency domain (see e.g., Buzsaki, 2011) with the aim of identifying physiologically and pathophysiologically relevant frequency bands. In addition, decomposing LFP into different frequency bands are used to correlate them to cognitive or motoric function as well as neuronal spiking activity. The aim here is to decipher brain activity in controlling perception, cognitive, motoric function and, in particular for the hippocampus, memory and learning abilities. Thus, the analysis of spike-LFP relationship represents another approach. Huge efforts are currently being done by model inversion approaches by creating artificial neuronal networks which produce realistic LFPs. In this approach, models of neural networks are combined with experimental data to identify the best-fit model. However, studies in which the applicability of mathematical or physical theories is evaluated by comparing the result of the model with experimental data are still rare but needed (Ness et al., 2015).

We aim to describe the spatiotemporal properties of subthreshold fluctuations in the rat hippocampal circuit by applying a mathematical description based on Gaussian statistical field theory to MEA data. Our study puts more emphasis on the spatio-temporal structure of correlations and less on oscillations. A generic feature of theoretical models in two spatial and one temporal dimension is a logarithmic behavior at short scales; as we will see this prediction is convincingly confirmed by our empirical material and is in a sense our main conceptual finding. From a perspective of practical electrophysiology, we would also like to emphasize the importance of correctly separating out the effects of certain artifacts without biological relevance before proceeding with the analysis.

2. Methods

2.1. Multi-Electrode Array Setup

Our dataset was acquired with a multi-electrode array system from Multi Channel Systems GmbH comprising 60 titanium/titanium nitride electrodes of 30 μm diameter arranged in a square grid pattern with 200 μm spacing on a non-conducting glass support. One of the electrodes served as a reference, and another one was not used, leaving 58 active electrodes. The voltage resolution, reflecting the binary representation of the data, was 2⁻¹⁶ × 10 mV ≃ 0.15μV. An 0.3 mm thick acute hippocampal slice from a 44 days old rat was fixed to the array field with a platina-nylon grid. See Figure 1 for a microscope image showing the positions of the electrodes and some of the relevant anatomical structures. Perfusion with a defined artificial cerebrospinal fluid (aCSF) provided the slice with glucose, a physiological salt concentration and osmolality. The layer of fluid above the slice had a thickness of several mm. The electrode potentials were sampled at 25 kHz during 600 s, yielding a total dataset of 870 × 10⁶ voltage measurements.

FIGURE 1

Figure 1. Microscope image of a hippocampal slice on the multi-electrode array. The inter-electrode distance is 200 μm. The reference electrode (just below the middle of the leftmost column) and the unused electrode (just to the left of the upper right corner) are indicated. We also give the approximate positions of the regions CA1, CA3, and the Dentate Gyrus as can be determined by usual anatomical considerations.

Conventionally, various filters are applied to the measured signals, but several studies demonstrate that such procedures do not remove spike components in LFPs (subthreshold activity) (Ray et al., 2008; Quilichini et al., 2010; Ray and Maunsell, 2011), and we will not use this approach. See however Figure 2 for the high and low pass filtered (above or below 50 Hz respectively) raw data recorded on the 59 electrodes (including the reference electrode just below the middle of the leftmost column).

FIGURE 2

Figure 2. Left: The high-pass (> 50 Hz) filtered signals on the 59 electrodes (including the reference electrode). Right: The low-pass (< 50 Hz) filtered signals. These figures are produced with a 2nd order Butterworth IIR filter. Each square comprises the entire 10 min registration and a ±50μV voltage interval.

2.2. Spikes

Since the spikes are not our primary interest, but rather obscure the analysis of the much smaller subthreshold fluctuations, they must be detected and removed from the dataset. We do this by a rather simple algorithm, which certainly leaves much room for improvements but is sufficient for our purposes: To detect spike events on an electrode at spatial point r, we consider the difference d(r, t) between the potential p(r, t) at time t and its average during a preceding time interval of some length Δt_average. We consider a spike to be fired at time t if the magnitude |d(r, t)| of this deviation then attains its maximum in the time window of length 2Δt_window centered at t and exceeds a threshold value d_threshold. In view of the typical ms timescale of the action potential dynamics, we used Δt_average = 10 ms and Δt_window = 2 ms. Concerning the threshold value, the value d_threshold = 20μV certainly misses many true but smaller spike events, but since our goal here is merely to remove large events that would interfere with the subsequent analysis, this is not a matter of great concern to us. On the other hand, picking a too low threshold value would give many false positives, and would lead us to remove large time intervals of intense neural activity, i.e., precisely the data that is our prime interest. In any case, the precise values of these parameters are not critical for our discussion.

A spike at time t can now be removed by replacing the true potential p(r, t′) in the interval $t-\Delta t_{spike}<t^{\prime }<t+\Delta t_{spike}$ for some time Δt_spike by the linear interpolating function:

$\begin{array}{l}{p}_{linear}\left(\text{r},t^{\prime }\right)=\frac{t^{\prime }-t+\Delta t_{spike}}{2\Delta t_{spike}}p\left(\text{r},t+\Delta t_{spike}\right)\\ -\frac{t^{\prime }-t-\Delta t_{spike}}{2\Delta t_{spike}}p\left(\text{r},t-\Delta t_{spike}\right).& (1)\end{array}$

We used Δt_spike = 2 ms, which in view of the observed spiking frequency leads to an almost negligible loss of subthreshold data, while still cutting out all large potential deviations. Henceforth, p(r, t) will always refer to the potential with all the spikes removed in this way.

2.3. The Stochastic Field Theory

Once the spikes have been removed, our aim is to describe the dynamics of the remaining subthreshold fluctuations. Our approach is to construct a simple model of this as a stochastic process which reproduces the main features of our dataset. Viewing the potential as the sum of a very large number of independent small contributions from different sources indicates (by the central limit theorem of statistics) that it should be normally distributed. This agrees well with the properties of our dataset, and it is thus a reasonable first approximation to limit ourselves to Gaussian models. We choose the reference potential so that the expectation value of the potential p(r, t) vanishes at all spatial points r and times t:

$\begin{array}{ll}\langle p\left(\text{r},t\right)\rangle =0.& (2)\end{array}$

All information is now contained in the two-point function 〈p(r₁, t₁)p(r₂, t₂)〉, and the higher-point functions can be expressed in terms of this by the Isserlis' theorem (in statistical physics mostly known as Wick's theorem), e.g.,

$\begin{array}{ll}\langle p\left({\text{r}}_1,t_1\right)p\left({\text{r}}_2,t_2\right)p\left({\text{r}}_3,t_3\right)\rangle =0& (3)\end{array}$

and

$\begin{array}{l}\langle p\left({\text{r}}_1,t_1\right)p\left({\text{r}}_2,t_2\right)p\left({\text{r}}_3,t_3\right)p\left({\text{r}}_4,t_4\right)\rangle \\ =\langle p\left({\text{r}}_1,t_1\right)p\left({\text{r}}_2,t_2\right)\rangle \langle p\left({\text{r}}_3,t_3\right)p\left({\text{r}}_4,t_4\right)\rangle \\ +\langle p\left({\text{r}}_1,t_1\right)p\left({\text{r}}_3,t_3\right)\rangle \langle p\left({\text{r}}_2,t_2\right)p\left({\text{r}}_4,t_4\right)\rangle \\ +\langle p\left({\text{r}}_1,t_1\right)p\left({\text{r}}_4,t_4\right)\rangle \langle p\left({\text{r}}_2,t_2\right)p\left({\text{r}}_3,t_3\right)\rangle .& (4)\end{array}$

Clearly, average neural activity depends both on the spatial location (related to different anatomical structures) and on time (reflecting the appearance of specific events during the course of the registration). However, in particular in view of the finite amount of data available, a natural first step of the analysis is to disregard these aspects. To begin with, we will thus make the assumption that the stochastic process is stationary in time as well as homogeneous and isotropic in space. We then have:

$\begin{array}{ll}\langle p\left({\text{r}}_1,t_1\right)p\left({\text{r}}_2,t_2\right)\rangle =S\left(|{\text{r}}_2-{\text{r}}_1|,|t_2-t_2|\right)& (5)\end{array}$

for some covariance function S(ρ, τ) which will be our primary object of study. Both of these assumptions certainly represent important oversimplifications, and later in the paper we will consider more general models.

On short time-scales (up to about 100 ms or so), the potential p(r, t) fluctuates around a slowly varying equilibrium potential μ(t) that is more or less independent of the spatial position r. We propose to describe this by a local, Gaussian, Markovian stochastic model of the form:

$\begin{array}{ll}\frac{\partial p\left(\text{r},t\right)}{\partial t}=-\gamma \left(p\left(\text{r},t\right)-\mu \left(t\right)\right)+\alpha {\nabla }^{2}p\left(\text{r},t\right)+\xi \left(\text{r},t\right).& (6)\end{array}$

Here ${\nabla }^2=\sum _{i=1}^2{\partial }_i{\partial }_i$ is the Laplacian operator in two spatial dimensions. The relaxation constant γ represents the tendency of the potential to return to its equilibrium value μ(t), and the diffusion constant α represents the tendency of spatial inhomogeneities to be smoothed out. The last term, which in stochastic modeling is usually referred to as “noise,” represents the contributions from the neural activity of a large number of neurons in the vicinity, much as molecular impacts drive Brownian motion. The usefulness of this description is related to the time scale of changes in the equilibrium potential μ(t) being larger than about 100 ms. For more background on statistical field theory (see e.g., Itzykson and Drouffe, 1991).

With initial data given in the far past so that its influence can be neglected, the solution to this equation is:

$\begin{array}{l}p(\text{r},t)={\displaystyle {\int }_{-\infty }^{t}d}{t}^{\prime }{\displaystyle \int d^2}{\text{r}}^{\prime }G(\text{r}-{\text{r}}^{\prime },t-t^{\prime })\left(\gamma \mu \right(t^{\prime })\\ +\xi ({\text{r}}^{\prime },t^{\prime })),& (7)\end{array}$

where the Green's function

$\begin{array}{ll}G\left(\text{r}-{\text{r}}^{\prime },t-t^{\prime }\right)=\frac{exp\left(-\gamma \left(t-t^{\prime }\right)-\frac{{\left(\text{r}-{\text{r}}^{\prime }\right)}^2}{4\alpha \left(t-t^{\prime }\right)}\right)}{4\pi \alpha \left(t-t^{\prime }\right)}& (8)\end{array}$

obeys the differential Equation

$\begin{array}{ll}\frac{\partial }{\partial t}G\left(\text{r}-{\text{r}}^{\prime },t-t^{\prime }\right)=\left(-\gamma +\alpha {\nabla }^2\right)G\left(\text{r}-{\text{r}}^{\prime },t-t^{\prime }\right)& (9)\end{array}$

and the initial condition

$\begin{array}{ll}G\left(\text{r}-{\text{r}}^{\prime },0\right)={\delta }^{\left(2\right)}\left(\text{r}-{\text{r}}^{\prime }\right).& (10)\end{array}$

Here and in the sequel, spatial integrals ∫ d²r are always taken over the infinitely extended plane. Boundary conditions at infinity (provided by the decay of the Green's function) are such that partial integrations do not generate any boundary contributions. The idea of Equation (7) and similar equations below is that because of the linearity of the model (6), there is a linear relationship between the driving input (represented by the last factor of the integrand) and the potential. The properties of the Green function ensure that this is indeed a solution to Equation (6). (For a further discussion on Green's function techniques for solving linear partial differential Equations, see e.g., Arfken et al., 2012).

The equilibrium potential μ(t) and the driving term ξ(r, t) are both assumed to have vanishing expectations values

$\begin{array}{l}\langle \mu \left(t\right)\rangle =0\\ \langle \xi \left(\text{r},t\right)\rangle =0& (11)\end{array}$

leading indeed to a vanishing expectation value for the potential p(r, t). We furthermore assume the covariance function of μ(t) to be given by some slowly varying function $S_{{\mu }^2}(\tau )$, whereas the driving term is assumed to be white both in space and time and uncorrelated with μ(t):

$\begin{array}{l}\langle \mu \left(t\right)\mu \left(t^{\prime }\right)\rangle =S_{{\mu }^2}\left(|t-t^{\prime }|\right)\\ \langle \xi \left(\text{r},t\right)\xi \left({\text{r}}^{\prime },t^{\prime }\right)\rangle ={\sigma }^2{\delta }^{\left(2\right)}\left(\text{r}-{\text{r}}^{\prime }\right)\delta \left(t-t^{\prime }\right)\\ \langle \xi \left(\text{r},t\right)\mu \left(t^{\prime }\right)\rangle =0.& (12)\end{array}$

Here the constant σ² represents the intensity of the neural activity. Accordingly, we can now decompose the covariance function appearing in Equation (5) as:

$\begin{array}{ll}S\left(\rho ,\tau \right)=S_{slow}\left(\tau \right)+S_{fast}\left(\rho ,\tau \right).& (13)\end{array}$

The first term in Equation (13) represents the contributions from the slow oscillations and can be expressed in terms of the covariance function $S_{{\mu }^2}(\tau )$ of the equilibrium potential. More precisely:

$\begin{array}{l}{S}_{slow}\left(\tau \right)={\gamma }^2{\displaystyle {\int }_{-\infty }^0}d{t}^{\prime }\int d^2{\text{r}}^{\prime }{\displaystyle {\int }_{-\infty }^{\tau }}d{t}^{\prime \prime }\int d^2{\text{r}}^{\prime \prime }\\ G\left(-{\text{r}}^{\prime },-t^{\prime }\right)G\left(-{\text{r}}^{\prime \prime },\tau -t^{\prime \prime }\right)S_{{\mu }^2}\left(|t^{\prime }-t^{\prime \prime }|\right)\\ \\ ={\gamma }^2{\displaystyle {\int }_{-\infty }^0}d{t}^{\prime }{\displaystyle {\int }_{-\infty }^{\tau }}d{t}^{\prime \prime }exp\left(-\gamma \left(\tau -t^{\prime }-t^{\prime \prime }\right)\right)\\ S_{{\mu }^2}\left(|t^{\prime }-t^{\prime \prime }|\right)\\ \\ =\gamma exp\left(-\gamma \tau \right){\displaystyle {\int }_0^{\tau }}d{\tau }^{\prime }cosh\left(-\gamma {\tau }^{\prime }\right)S_{{\mu }^2}\left({\tau }^{\prime }\right)\\ +\gamma cosh\left(-\gamma \tau \right){\displaystyle {\int }_{\tau }^{\infty }}d{\tau }^{\prime }exp\left(-\gamma {\tau }^{\prime }\right)S_{{\mu }^2}\left({\tau }^{\prime }\right).& (14)\end{array}$

In principle, this may be inverted to express $S_{{\mu }^2}(\tau )$ in terms of S_slow(τ):

$\begin{array}{ll}{S}_{{\mu }^2}\left(\tau \right)=\left(1-\frac{1}{{\gamma }^2}\frac{{\partial }^2}{\partial {\tau }^2}\right)S_{slow}\left(\tau \right).& (15)\end{array}$

Because of the second derivative, it is however difficult to achieve an accurate estimate of $S_{{\mu }^2}(\tau )$ with the available data, and we will not develop this approach further.

The second term in Equation (13) represents the contributions from the driving term and can be expressed in terms of the intensity σ². A short computation gives

$\begin{array}{l}{S}_{\text{fast}}(\rho ,\tau )={\displaystyle {\int }_{-\infty }^{0}d}{t}^{\prime }{\displaystyle \int d^2}{\text{r}}^{\prime }{\displaystyle {\int }_{-\infty }^{\tau }d}{t}^{\prime \prime }{\displaystyle \int d^2}{\text{r}}^{\prime \prime }\\ G(-{\text{r}}^{\prime },-t^{\prime })G(\rho -{\text{r}}^{\prime \prime },\tau -t^{\prime \prime }){\sigma }^2{\delta }^{\left(2\right)}({\text{r}}^{\prime }-{\text{r}}^{\prime \prime })\\ \delta (t^{\prime }-t^{\prime \prime })\\ ={\sigma }^2{\displaystyle {\int }_{-\infty }^{0}d}{t}^{\prime }{\displaystyle \int d^2}{\text{r}}^{\prime }\\ \frac{\exp \left(-\gamma (\tau -2t^{\prime })-\frac{{(-{\text{r}}^{\prime })}^2}{4\alpha (-t^{\prime })}-\frac{{(\rho -{\text{r}}^{\prime })}^2}{4\alpha (\tau -t^{\prime })}\right)}{{\left(4\pi \alpha \right)}^2(-t^{\prime })(\tau -t^{\prime })}\\ ={\sigma }^2{\displaystyle {\int }_{-\infty }^{0}d}{t}^{\prime }{\displaystyle \int d^2}{\text{r}}^{\prime }\\ \frac{\begin{array}{l}\exp (-\gamma (\tau -2t^{\prime })-\frac{\tau -2t^{\prime }}{4\alpha (-t^{\prime })(\tau -t^{\prime })}\hfill \\ {\left({\text{r}}^{\prime }-\frac{(-t^{\prime })\rho }{\tau -2t^{\prime }}\right)}^2-\frac{{\rho }^2}{4\alpha (\tau -2t^{\prime })})\hfill \end{array}}{{\left(4\pi \alpha \right)}^2(-t^{\prime })(\tau -t^{\prime })}\\ ={\sigma }^2{\displaystyle {\int }_{-\infty }^{0}d}{t}^{\prime }\frac{\exp \left(-\gamma (\tau -2t^{\prime })-\frac{{\rho }^2}{4\alpha (\tau -2t^{\prime })}\right)}{4\pi \alpha (\tau -2t^{\prime })}.& (16)\end{array}$

For fixed ρ or τ, this is a monotonously decreasing function of τ or ρ respectively. In general it cannot be expressed in terms of any well-known elementary or special functions. However, for vanishing spatial separation, i.e., ρ = 0, it is given by:

$\begin{array}{l}{S}_{fast}\left(0,\tau \right)={\sigma }^2{\displaystyle {\int }_{-\infty }^0}d{t}^{\prime }\frac{exp\left(-\gamma \left(\tau -2t^{\prime }\right)\right)}{4\pi \alpha \left(\tau -2t^{\prime }\right)}\\ \\ =\frac{{\sigma }^2}{8\pi \alpha }\Gamma \left(0,\gamma \tau \right)\\ \\ =\frac{{\sigma }^2}{8\pi \alpha }\left(-log\left(\gamma \tau \right)-{\gamma }_{EM}+{O}\left(\gamma \tau \right)\right),& (17)\end{array}$

where Γ is the (upper) incomplete Gamma-function and γ_EM = 0.5772… is the Euler-Mascheroni constant. Similarly, for vanishing temporal separation, i.e., τ = 0, we instead have:

$\begin{array}{l}{S}_{\text{fast}}(\rho ,0)={\sigma }^2{\displaystyle {\int }_{-\infty }^{0}d}{t}^{\prime }\frac{\exp \left(-\gamma (-2t^{\prime })-\frac{{\rho }^2}{4\alpha (-2t^{\prime })}\right)}{4\pi \alpha (-2t^{\prime })}\\ =\frac{{\sigma }^2}{4\pi \alpha }{K}_0\left(\sqrt{\gamma /\alpha }\rho \right)\\ =\frac{{\sigma }^2}{4\pi \alpha }(-\log \left(\sqrt{\gamma /\alpha }\rho \right)-{\gamma }_{\text{EM}}\\ +\log 2+O\left(\sqrt{\gamma /\alpha }\rho \right)),& (18)\end{array}$

where K₀ is a modified Bessel-function.

The most important aspects of the results Equations (17) and (18) are that they exhibit the logarithmic dependence of the covariance function for short temporal and spatial separations respectively with coefficients that are directly related to the parameters of the model. Such logarithmic behavior is a generic feature of field theories in two spatial dimensions regardless of the details of the model, but does not hold in other dimensions.

3. Results

3.1. Spikes

With our choices Δt_average = 10 ms, Δt_refractory = 2 ms, and d_threshold = 20μV, the dataset had a total spike firing frequency of:

$\begin{array}{ll}{\nu }_{total}\simeq 9.5Hz.& (19)\end{array}$

The spikes were rather unevenly distributed, both in time over the 600 s registration and over the 58 electrodes: About 49% of all spikes were fired in the Dentate Gyrus, where they were mostly negative and tended to occur in short burst of less than 1 s, and 48% were fired the CA3 region, where they were mostly positive and the spiking frequency fluctuated on times scales of about 100 s. (The remaining 3% tended to occur in the DG/CA-3 intermediate area.) Although the total number of spikes in these two areas were very nearly equal, their temporal distributions were quite different and give no evidence for any causal connection. See Figures 3, 4 for the spatial and temporal distribution of the spikes. See Figure 5 for examples of negative and positive spikes and their removal by linear interpolation.

FIGURE 3

Figure 3. Individual spiking frequencies detected on the different electrodes represented by the radii of the dots. The most spiking electrode (in the Dentate gyrus) had a spiking frequency of about 1.7 Hz.

FIGURE 4

Figure 4. Left: Histogram of the temporal distribution of spikes during the 600 s registration in the Dentate Gyrus in 10 s bins. Right: Histogram of the temporal distribution of spikes during the 600 s registration in the CA3 region in 10 s bins.

FIGURE 5

Figure 5. Left: Example of a negative spike (potential as a function of time) in the Dentate Gyrus and its removal by linear interpolation. Right: Example of a positive spike in the CA3 region and its removal by linear interpolation.

3.2. Artifacts

After the spikes had been removed, we computed the covariance function S(ρ, τ) by using the entire dataset sampled at 25 kHz. The magnitude of the covariance function S(0, 0) at vanishing spatial and temporal separation, i.e., the variance of the signal, had a magnitude of about 15 μV². Two features of the covariance function S(ρ, τ) appeared to be artifacts without biological significance:

• There was an almost perfectly periodic component with a period of about 145 ms (corresponding to 6.9 Hz with some overtones) and a maximal amplitude of about 0.12 μV² that persisted essentially undamped until τ = 10 s or more. The extreme and persistent regularity of this phenomenon makes it clear that it originated within the electronics of the multi-electrode array system. Although the magnitude was quite modest, we still found it advantageous (and straightforward) to subtract this component from the covariance function, since its time scale was so close to those of biological relevance.

• For ρ = 0, i.e., at vanishing spatial separation, there was a component with a pronounced peak in the interval 0 < τ < 0.2 ms (i.e., during 5 sampling intervals) with a maximal amplitude of about 4 μV². The short spatial range (less than the electrode spacing) and the short time-scale involved strongly suggested that this phenomenon was due to essentially independent errors in the individual voltage measurements (about 2 μV) with an extremely short correlation time (about 0.2 ms). Because of the large magnitude, it was necessary to take this component properly taken into account, although its time scale of course was much shorter than those of biological phenomena.

See Figure 6 for the appearance of these two artifact components in the covariance function. Henceforth S(ρ, τ) will always refer to the covariance function after these artifacts had been removed by subtracting the two temporal profiles exhibited in the figure from the raw-data covariance function.

FIGURE 6

Figure 6. Left: Two periods of the long-term artifact with 145 ms periodicity. This was obtained by subtracting a 145 ms moving average from the raw-data covariance function. For time lags exceeding about 1 s, the difference had an almost perfectly periodic appearance, which was then extend down to zero time lag. Right: The rapidly decaying artifact due to errors in the individual voltage measurements.

3.3. Slow Dynamics

For large values of τ (> 100 ms), the main feature of the covariance function S(ρ, τ) was a damped oscillatory behavior mainly in the 0 - 2 Hz frequency range, which was largely independent of the spatial separation ρ. This slow oscillation is possibly of biological relevance, but we will not attempt to analyze or model it in the present paper. (See e.g., Buzsáki and Draguhn, 2004, for a review of LFPs with different frequency bands.) See Figure 7 for this long-time behavior of the covariance function. We took S_slow(τ) in Equation (13) to be given by S(ρ_large, τ) for ρ_large = 1.7 mm, i.e., the largest spatial separation available to us. As can be seen from Figure 8 (lowest curve in the left panel), for such a large spatial separation the covariance was essentially independent of the time lag up to about 100 ms, so S_fast(ρ_large, τ) is negligible.

FIGURE 7

Figure 7. Long-time covariance at vanishing spatial separation as a function of time lag.

FIGURE 8

Figure 8. Left: Short-time covariance at spatial separations 0, 0.2, and 1.7 mm (top, middle and lower curve) as a function of time lag. Right: Covariance at vanishing time lag as a function of spatial separation. In all cases, we exhibit an average over all pairs of electrodes with the indicated spatial separation.

3.4. Fast Dynamics

For small values of τ (< 100 ms), the covariance function S(0, τ) at vanishing spatial separation indeed increased logarithmically as τ approaches zero. For ρ > 0, this increase was cut off so that the equal time covariance S(ρ, 0) has a finite value that increases logarithmically as ρ approaches zero. See Figure 8 for the temporal and spatial dependence of the short-time covariance function. Note the logarithmic abscissa axis in these figures!

The measured values of S_fast(ρ, τ) = S(ρ, τ) − S(ρ_large, τ) were fitted to the the theoretical prediction Equation (16). This is shown in Figure 9 with the parameter values

$\begin{array}{l}\alpha \simeq 0.0025{mm}^2{ms}^{-1}\\ \gamma \simeq 0.0030{ms}^{-1}\\ {\sigma }^2\simeq 0.035{\mu V}^2{mm}^2{ms}^{-1}.& (20)\end{array}$

As can be seen from the figure, the agreement between theory and experiment was excellent, providing a convincing and a priori falsifiable confirmation of the validity of our approach and simplifying assumptions. (In view of our still rather restricted dataset, we refrain from quoting any specific uncertainty range of these parameters.)

FIGURE 9

Figure 9. Left: Fitting S_fast(ρ, 0) to Equation (18). (Equal time covariance as a function of spatial separation). Right: Fitting S_fast(0, τ) to Equation (17). (Covariance as a function of time lag at vanishing spatial separation.)

The definition and determination of the three quantities α, γ and σ² constitute the main results of the present work. An equivalent, but in many respects more illuminating presentation of the results is to combine these parameters into characteristic time, length and voltage scales:

$\begin{array}{lll}1/\gamma & \simeq & 330ms\\ \sqrt{\alpha /\gamma }& \simeq & 0.91mm\\ \sqrt{{\sigma }^2/\alpha }& \simeq & 3.7\mu V.& (21)\end{array}$

The 330 ms time scale of these “fast” fluctuations may seem uncomfortably close to the time scale of the “slow” fluctuations of the equilibrium potential μ(t) (which seems to be around 1 s). In this context, we remark that the time scale can be generalized to:

$\begin{array}{ll}T=\frac{1}{\gamma +\alpha {\left(2\pi /\lambda \right)}^2}& (22)\end{array}$

for fluctuations of some finite wavelength λ. In the long wave-length limit λ → ∞ we recover 1/γ, whereas for the inter-electrode distance λ = 0.2 mm we instead get T ≃ 0.4 ms. So for wave-lengths relevant for investigating the local dynamics, there is no problem with the time scale. We also remark that the 25 kHz sampling frequency is clearly high enough.

The length scale 0.91 mm is comparable to the extent of the entire multi-electrode array. However, the dimensions of the slice of neural tissue are considerably larger, so there is no need to worry about finite size effects. More importantly, the length scale is sufficiently large compared to the 0.2 mm inter-electrode distance to assure the validity of this experimental approach to the study of subthreshold fluctuations.

Finally the voltage scale 3.7μV is safely smaller than the spikes (which we have cut off at 20 μV). However, it is quite comparable both to the errors in the individual voltage measurements (about 2 μV) and the amplitude of the slow fluctuations of the equilibrium potential μ(t), so it is important to carefully separate these three phenomena.

3.5. Activity

Sofar we have considered the parameters α and γ as well as the activity σ² to be constants. This is reasonable for α and γ, at least if we view these constants as reflecting only the passive electric transport properties of the intracellular medium and not the propagation of signals along the axons. But the activity should rather be described by a function σ²(r, t) of space and time, reflecting the characteristics of the neuronal populations in the different anatomical regions as well as the time course of the neural processing. The value σ² ≃ 0.035μV²mm² ms⁻¹ that we have determined should thus be regarded as a spatial and temporal average.

Retracing the steps leading to Equation (16), we find that with a non-constant activity σ²(r, t), this expression is no longer valid. However, the leading logarithmic divergence of Equation (17), which originates from the short distance behavior of the model, still holds. Since S_slow(t) is regular for small t, we thus have:

$\begin{array}{ll}\langle p\left(\text{r},t\right)p\left(\text{r},t+\delta t\right)\rangle =\frac{{\sigma }^2\left(\text{r},t\right)}{8\pi \alpha }\left(-log\left(\gamma \delta t\right)+{O}\left(1\right)\right).& (23)\end{array}$

The diffusion constant α is of course already known. The values of δt can e.g., be chosen in the interval 0.2 to 10 ms. Since we have only a single measurement of the potential p(r, t) for each value of r and t, we can only estimate such expectation values by averaging over a rather large time interval (at least about 100 ms) around t, which limits the temporal resolution of the method.

Averaging over the entire 600 s registration, we found that the temporal mean $\overline{{\sigma }^2(\text{r},t)}$ of the activity was concentrated in the Dentate Gyrus and the CA3 region just like the spikes, but much more spread out. There was however also substantial activity in the area intermediate between these two regions (where essentially no spikes occur), whereas the CA1 region showed very little activity. See Figure 10 for this spatial distribution of activity. Comparison can be made with the spatial distribution of spikes in Figure 3.

FIGURE 10

Figure 10. Temporal mean over the 600 s registration of the activity on the different electrodes. The most active electrode (in the Dentate gyrus) has an activity of about 0.11 μV²mm²ms⁻¹.

Averaging over time intervals of 1 s instead, we could investigate the temporal dependance of the activity in the different regions. We found a clear correlation with the spiking in the Dentate Gyrus and the CA3 region. In the intermediate non-spiking area, the pattern was more reminiscent of the CA3 region than the Dentate Gyrus. The CA1 region showed a rather constant lower activity. See Figure 11, which should be compared with the corresponding temporal distributions of spikes in Figure 4.

FIGURE 11

Figure 11. Mean activity in each of the four quadrants of the multi-electrode array (roughly corresponding to the Dentate Gyrus, the DG/CA3-intermediate area, the CA3 region and the CA1 region clockwise from the lower right corner) as a function of time during the 600 s registration. The temporal resolution in these graphs is 1 s.

One sees that activity and spiking are indeed different phenomena, although there seems to exist some connection between them. Again, we take the view that the spiking frequency registered on the different electrodes reflects not only what is going at that location in the tissue but also on how a few individual neurons happen to be in more or less close contact with the electrodes.

3.6. Connectivity

In contrast to the potential p(r, t), which is a stochastic variable, we have considered the activity σ²(r, t) to be a given function of space and time. Ultimately, one would of course like to formulate some (deterministic or stochastical) dynamical model for it, but we will not pursue this here and instead content ourselves with a purely descriptive treatment. While the activity directly influences the variance of the signal p(r, t), it shows essentially no correlation with the mean of p(r, t). This is another indication that the equilibrium potential μ(t), while serving as a common voltage reference for the entire network, may not be of immediate biological relevance.

A very useful quantity for characterizing the activity σ²(r, t) is its covariance function between separate points r₁ and r₂ at some time lag Δt:

$\begin{array}{l}Cov\left({\sigma }^2\left({\text{r}}_1,t\right),{\sigma }^2\left({\text{r}}_2,t+\Delta t\right)\right)=\\ \overline{{\sigma }^2\left({\text{r}}_1,t\right){\sigma }^2\left({\text{r}}_2,t+\Delta t\right)}-\overline{{\sigma }^2\left({\text{r}}_1,t\right)}\overline{{\sigma }^2\left({\text{r}}_2,t\right)}.& (24)\end{array}$

Here and in the sequel, an overline denotes an average over the time t. In particular, we have the temporal autocovariance function

$\begin{array}{l}Cov\left({\sigma }_{mean}^2\left(t\right),{\sigma }_{mean}^2\left(t+\Delta t\right)\right)=\\ \overline{{\sigma }_{mean}^2\left(t\right){\sigma }_{mean}^2\left(t+\Delta t\right)}-\overline{{\sigma }_{mean}^2\left(t\right)}\overline{{\sigma }_{mean}^2\left(t+\Delta t\right)},& (25)\end{array}$

where

$\begin{array}{ll}{\sigma }_{mean}^2\left(t\right)=\frac{1}{{Vol}_{\Omega }}{\displaystyle {\int }_{\Omega }}{d}^2\text{r}{\sigma }^2\left(\text{r},t\right)& (26)\end{array}$

is the spatial mean of the activity. (We take the domain Ω to cover the entire multi-electrode array.) Empirically, we find that

$\begin{array}{ll}Cov\left({\sigma }_{mean}^2\left(t\right),{\sigma }_{mean}^2\left(t+\Delta t\right)\right)~exp\left(-\beta \Delta t\right),& (27)\end{array}$

with decay constant

$\begin{array}{ll}\beta \simeq 0.1s^{-1}.& (28)\end{array}$

(In these last formulas, the expressions are in fact independent of the time t appearing in the left hand sides.)

Similarly, we can investigate the spatial autocovariance function:

$\begin{array}{l}Cov\left({\sigma }^2\left(\text{r},t\right),{\sigma }^2\left(\text{r}+\Delta \text{r},t\right)\right)=\\ \overline{{\sigma }^2\left(\text{r},t\right){\sigma }^2\left(\text{r}+\Delta \text{r},t\right)}-\overline{{\sigma }^2\left(\text{r},t\right)}\overline{{\sigma }^2\left(\text{r}+\Delta \text{r},t\right)}.& (29)\end{array}$

Taking the spatial mean, we here find

$\begin{array}{ll}\frac{1}{{Vol}_{\Omega }}{\displaystyle {\int }_{\Omega }}{d}^2\text{r}Cov\left({\sigma }^2\left(\text{r},t\right),{\sigma }^2\left(\text{r}+\Delta \text{r},t\right)\right)~exp\left(-\kappa \left|\Delta \text{r}\right|\right),& (30)\end{array}$

with decay constant

$\begin{array}{ll}\kappa \simeq 1.4{mm}^{-1}.& (31)\end{array}$

See Figure 12 for the corresponding autocorrelation functions (normalized to 1 for Δt = 0 and Δr = 0 respectively). Note the logarithmic scales! The deviations from exponential decay for small time lags and distances can be attributed to the measurement errors, which in the temporal case are smoothed out over 1 s by our data analysis.

FIGURE 12

Figure 12. Left: Fitting the spatial autocorrelation function of the activity at equal time to the exponential expression 30. Right: Fitting the temporal autocorrelation function of the spatial mean activity to the exponential expression 27.

The behavior of these covariance functions should be relevant for the understanding of neural connectivity and communication. The exponential decay is qualitatively rather different from the logarithmic behavior characteristic of passive transport in two spatial dimensions as we have investigated for the potential p(r, t). This possibly indicates that the activity is propagated by some more “active” mechanism, for which we do not have any specific proposal. However, the temporal scale of about 10 s is long enough that the biological significance of these slow changes may be questioned, in which case they should probably be attributed to drifting conditions during the registration. Indeed, the biologically relevant information transfer in the neural tissue should be encoded in fluctuations of the activity at much shorter time-scales, which we are however unable to probe with our present methods. On the other hand, the spatial scale of about 0.7 mm is quite similar to the scale $\sqrt{\alpha /\gamma }\simeq 0.91mm$ set by the diffusion process, and again indicates that the multi-electrode array is adequate for investigating these phenomena.

Finally, we considered the Pearson correlation coefficient of the activity ${\sigma }^2(r_1,t)$ and ${\sigma }^2(r_2,t)$ separately for all pairs of adjacent electrodes, which gives a way of investigating the local neural connectivity. A priori, such a correlation can be weak or strong regardless of the mean and variances of the two activities under consideration. With |Δr| = 0.2 mm the average correlation coefficient was about 0.50, but varied considerably between 0.1 and 0.9 for the different pairs. With sufficiently strong inhibitory connections, one could in principle also imagine negative correlation coefficients in the interval -1 to 0, but these did not occur in our dataset. Highly correlated pairs indicated a path of information flow from the Dentate Gyrus to the region CA3 with a hint of a continuation toward CA1, in agreement with the expectations from anatomical considerations. (Actually, our methods cannot determine the direction of this information flow, since the correlation is invariant under the exchange of two electrodes.) See Figure 13 for an attempt at a graphical rendering of this connectivity pattern.

FIGURE 13

Figure 13. Left: Histogram of the distribution of correlation coefficients for the activity on all pairs of adjacent electrodes. Right: The connectivity between adjacent electrodes. The thickness of the lines is proportional to the fourth power of the correlation coefficient for the corresponding pair of electrodes.

4. Discussion

Our main finding is that the LFP can be remarkably well-described by a Gaussian statistical field theory in two space and one time dimension. Depending on the electrical properties of the perfusion liquid above the tissue sample, one may argue that this should be modeled as a three-dimensional rather than a two-dimensional system. This would give a qualitatively rather different model, in which correlations decay as the inverse of the distance rather than logarithmically in both space and time. However, our two-dimensional model fits the data excellently, whereas such a three-dimensional model would be in clear disagreement. So we have been able to make a clear and falsifiable theoretical prediction and verify it experimentally in what we think is a convincing manner. From a perspective of practical electrophysiology, we would like to emphasize that this analysis must be preceded by a correct elimination of certain artifacts of no biological significance.

Thus, we have provided a proof of concept of a new approach for studying neural circuit function and applied it to the dataset described above. By our approach, we have described the mean activity (Figure 11) and correlation (Figure 13) of subthreshold fluctuations within specific hippocampal sub-regions. It appears that the connection between CA3 and CA1 in this particular isolated hippocampal ex vivo slice preparation is not preserved. Even though this represents a drawback of our used MEA data set, the presented connection between DG and CA3 demonstrates that our approach allows for the identification and visualization of connected sub-regions in isolated brain-slice preparations.

We also computed specific mean values for the parameters characterizing the duration, spatial distribution and amplitude of subthreshold fluctuations in all hippocampal sub-regions. Since we aim to present a proof-of-concept of our approach by using data sets collected only in only one hippocampal slice preparation, we did not perform a hippocampal sub-region specific classification of subthreshold fluctuations. Of course, such parametric description of sub-neuronal network properties within the hippocampus is quite interesting and will be addressed in future studies.

Our method to extract subthreshold fluctuations out of MEA data sets can be used to uncover spatial and temporal correlations of sub-hippocampal neuronal circuits within brain slice preparations, which can not be achieved by analyzing the localization of spike activity. Indeed, while the possibility to detect spikes is largely determined by the accidental proximity of a neuron to an electrode, the activity as we define it should be a robust concept. Thus, describing the spatial and temporal properties of subthreshold fluctuations in ex vivo or in vitro neuronal circuits may represent a better approach to uncover functional connectivity within neuronal circuits than analysis of synchronous bursting. However, of course also the activity as we have defined it is to a large extent determined by the population of nearby neurons, so it is not obvious to separate these aspects from each other. Indeed, although we have defined the activity without any reference to detected spike events, it still shows a clear correlation with these, and its spatial distribution and correlations agree well with expectations from anatomical considerations. See e.g., Giugliano et al. (2004) for a discussion of the relationship between collective network phenomena and the spiking of individual neurons.

It would be interesting to try to get a better understanding of the laws underlying the slow dynamics, rather than just describing the resulting equilibrium potential. This can be done with a larger dataset, provided that the longterm stability of the preparation can be assured.

Our model has a small number free parameters, the values of which can be readily determined by fitting the experimentally measured correlations. We expect that these parameters will provide robust and reproducible quantities suitable for comparative studies between brain tissue samples from different anatomical regions and developmental stages under various physiological and patho-physiological conditions. It can also be valuable to study patient-specific neuronal circuits obtained from induced pluripotent stem cell technology. In the future, it would thus be very interesting to apply these methods to more datasets.

In a different direction, the understanding of the passive transport properties of the neural preparation developed in this article should be useful also for analyzing spiking events. Indeed, such an analysis is complicated by the fact that signals spread between the different electrodes, so a natural approach is to begin by reconstructing the local sources of these events by inverse methods based on our model. We plan to return to these issues in forthcoming publications.

Author Contributions

SI performed the experimental work, read and improved the manuscript. MH conceived the theoretical analysis, wrote most of the manuscript.

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.

Acknowledgments

We have benefitted from discussions with Eric Hanse. The research of MH was supported by a grant from the Kristina Stenborg foundation. The research of SI was supported by a grant from the Alzheimerfonden (AF-556051).

References

Arfken, G. B., Weber, H. J., and Harris, F. E. (2012). Mathematical Methods for Physicists. Waltham, MA: Academic Press.

Google Scholar

Bedard, C., and Destexhe, A. (2015). Local field potential interaction with the extracellular medium. Encycloped. Comput. Neurosci. 1540–1547. doi: 10.1007/978-1-4614-7320-6_720-1

CrossRef Full Text | Google Scholar

Buzsaki, G. (2011). Rhythms of the Brain. New York, NY: Oxford University Press.

Google Scholar

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Cotterill, E., Charlesworth, P., Thomas, C. W., Paulsen, O., and Eglen, S. J. (2016). A comparison of computational methods for detecting bursts in neuronal spike trains and their application to human stem cell-derived neuronal networks. J. Neurophysiol. 116, 306–321. doi: 10.1152/jn.00093.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., and Bedard, C. (2013). Local field potential. Scholarpedia 8:10713. doi: 10.4249/scholarpedia.10713

CrossRef Full Text | Google Scholar

Egert, U., Heck, D., and Aertsen, A. (2002). Two-dimensional monitoring of spiking networks in acute brain slices. Exp. Brain Res. 142, 268–274. doi: 10.1007/s00221-001-0932-5

PubMed Abstract | CrossRef Full Text

Giugliano, M., Darbon, P., Arsiero, M., Luscher, H. R., Streit, J. (2004). Single-neuron discharge properties and network activity in dissociated cultures of neocortex. J. Neurophysiol. 92, 977–996. doi: 10.1152/jn.00067.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Heikkilä, T. J., Ylä-Outinen, L., Tanskanen, J. M., Lappalainen, R. S., Skottman, H., Suuronen, R., et al. (2009). Human embryonic stem cell-derived neuronal cells form spontaneously active neuronal networks in vitro. Exp. Neurol. 218, 109–116. doi: 10.1016/j.expneurol.2009.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Illes, S., Jakab, M., Beyer, F., Gelfert, R., Couillard-Despres, S., Schnitzler, A., et al. (2014). Intrinsically active and pacemaker neurons in pluripotent stem cell-derived neuronal populations. Stem Cell Rep. 2, 323–336. doi: 10.1016/j.stemcr.2014.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Itzykson, C., and Drouffe, J.-M. (1991). Statistical Field Theory. Vols I and II. Cambridge, UK: Cambridge University Press.

Google Scholar

Linden, H., Tetzlaff, T., Potjans, T. C., Pettersen, K. H., Grun, S., Diesmann, M., et al. (2014). Modeling the spatial reach of the LFP. Neuron 72, 859–872. doi: 10.1016/j.neuron.2011.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Nam, Y., and Wheeler, B. C. (2011). In vitro microelectrode array technology and neural recordings. Crit. Rev. Biomed. Eng. 39, 45–62. doi: 10.1615/CritRevBiomedEng.v39.i1.40

PubMed Abstract | CrossRef Full Text | Google Scholar

Ness, T. V., Chintaluri, C., Potworowski, J., Leski, S., Glabska, H., Wojcik, D. K., et al. (2015). Modelling and analysis of electrical potentials recorded in microelectrode arrays (MEAs). Neuroinformatics 13, 403–426. doi: 10.1007/s12021-015-9265-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Quilichini, P., Sirota, A., and Buzsaki, G. (2010). Intrinsic circuit organization and theta-gamma oscillation dynamics in the entorhinal cortex of the rat. J. Neurosci. 30, 11128–11142. doi: 10.1523/JNEUROSCI.1327-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ray, S., and Maunsell, J. H. (2011). Different origins of gamma rhythm and high-gamma activity in macaque visual cortex. PLoS Biol. 9:e1000610. doi: 10.1371/journal.pbio.1000610

PubMed Abstract | CrossRef Full Text | Google Scholar

Ray, S., Hsiao, S. S., Crone, N. E., Franaszczuk, P. J., and Niebur, E. (2008). Effect of stimulus intensity on the spike-local field potential relationship in the secondary somatosensory cortex. J. Neurosci. 28, 7334–7343. doi: 10.1523/JNEUROSCI.1588-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Reichinnek, S., Kunsting, T., Draguhn, A., and Both, M. (2010). Field potential signature of distinct multicellular activity patterns in the mouse hippocampus. J. Neurosci. 30, 15441–15449. doi: 10.1523/JNEUROSCI.2535-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieke, F., Warland, D., de Ruyter van Steveninck, R., and Bialek, W. (1997). Spikes. Exploring the Neural Code. Cambridge, MA: MIT Press.

Riera, J. J., Ogawa, T., Goto, T., Sumiyoshi, A., Nonaka, H., Evans, A., et al. (2012). Pitfalls in the dipolar model for the neocortical EEG sources. J. Neurophysiol. 108, 956–975. doi: 10.1152/jn.00098.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Spira, M. E., and Hai, A. (2013). Multi-electrode array technologies for neurscience and cardiology. Nat. Nanotechnol. 8, 83–94. doi: 10.1038/nnano.2012.265

PubMed Abstract | CrossRef Full Text | Google Scholar

Waldert, S., Lemon, R. N., and Kraskov, A. (2013). Influence of spiking activity on cortical local field potentials. J. Physiol. 591, 5291–5303. doi: 10.1113/jphysiol.2013.258228

PubMed Abstract | CrossRef Full Text | Google Scholar