Analysis and Modelling of Subthreshold Neural Multi-electrode Array Data by Statistical Field Theory

TL;DR

This paper develops a Gaussian statistical field theory model to characterize subthreshold neural fluctuations in multi-electrode array data, revealing correlations related to neural connectivity and spike activity.

Contribution

The paper introduces a novel theoretical framework for analyzing subthreshold neural signals using Gaussian field theory, applied to real MEA data from rat hippocampal slices.

Findings

Empirical correlation functions follow logarithmic behavior predicted by the model.

Clear correlation between neural activity and spike occurrence.

Artefact removal is crucial for accurate analysis.

Abstract

Multi-electrode arrays (MEA) are increasingly used to investigate spontaneous neuronal network activity. The recorded signals comprise several distinct components: Apart from artefacts 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 behaviour that is a general feature of theoretical models of this kind. We also find a clear correlation between the activity and the occurence of spikes. Another important insight is the importance of correcly separating out certain artefacts from the data before proceeding with the analysis.

Equations

$$p_{\mathrm{linear}}(\mathbf{r},t^{\prime })=\frac{t^{\prime }-t+\Delta t_{\mathrm{spike}}}{2\Delta t_{\mathrm{spike}}}p(\mathbf{r},t+\Delta t_{\mathrm{spike}})-\frac{t^{\prime }-t-\Delta t_{\mathrm{spike}}}{2\Delta t_{\mathrm{spike}}}p(\mathbf{r},t-\Delta t_{\mathrm{spike}}).$$

$$⟨p(\mathbf{r},t)⟩=0.$$

$$⟨p({\mathbf{r}}_1,t_1)p({\mathbf{r}}_2,t_2)p({\mathbf{r}}_3,t_3)⟩=0$$

$$⟨p({\mathbf{r}}_1,t_1)p({\mathbf{r}}_2,t_2)p({\mathbf{r}}_3,t_3)p({\mathbf{r}}_4,t_4)⟩$$

$$⟨p({\mathbf{r}}_1,t_1)p({\mathbf{r}}_2,t_2)⟩=S(|{\mathbf{r}}_2-{\mathbf{r}}_1|,|t_2-t_2|)$$

$$\frac{\partial p(\mathbf{r},t)}{\partial t}=-\gamma \left(p(\mathbf{r},t)-\mu (t)\right)+\alpha {\nabla }^{2}p(\mathbf{r},t)+\xi (\mathbf{r},t).$$

$$p(\mathbf{r},t)={\int }_{-\infty }^{t}d{t}^{\prime }\int d^2{\mathbf{r}}^{\prime }G(\mathbf{r}-{\mathbf{r}}^{\prime },t-t^{\prime })\left(\gamma \mu (t^{\prime })+\xi ({\mathbf{r}}^{\prime },t^{\prime })\right),$$

$$G(\mathbf{r}-{\mathbf{r}}^{\prime },t-t^{\prime })=\frac{\exp \left(-\gamma (t-t^{\prime })-\frac{(\mathbf{r}-{\mathbf{r}}^{\prime }{)}^2}{4\alpha (t-t^{\prime })}\right)}{4\pi \alpha (t-t^{\prime })}$$

$$\frac{\partial }{\partial t}G(\mathbf{r}-{\mathbf{r}}^{\prime },t-t^{\prime })=\left(-\gamma +\alpha {\nabla }^2\right)G(\mathbf{r}-{\mathbf{r}}^{\prime },t-t^{\prime })$$

$$G(\mathbf{r}-{\mathbf{r}}^{\prime },0)={\delta }^{(2)}(\mathbf{r}-{\mathbf{r}}^{\prime }).$$

$$⟨\mu (t)⟩$$

$$⟨\xi (\mathbf{r},t)⟩$$

$$⟨\mu (t)\mu (t^{\prime })⟩$$

$$⟨\xi (\mathbf{r},t)\xi ({\mathbf{r}}^{\prime },t^{\prime })⟩$$

$$⟨\xi (\mathbf{r},t)\mu (t^{\prime })⟩$$

$$S(\rho ,\tau )=S_{\mathrm{slow}}(\tau )+S_{\mathrm{fast}}(\rho ,\tau ).$$

$$S_{\mathrm{slow}}(\tau )$$

$$=$$

$$S_{{\mu }^2}(\tau )=\left(1-\frac{1}{{\gamma }^2}\frac{{\partial }^2}{\partial {\tau }^2}\right)S_{\mathrm{slow}}(\tau ).$$

$$S_{\mathrm{fast}}(\rho ,\tau )$$

$$=$$

$$S_{\mathrm{fast}}(0,\tau )$$

$$=$$

$$S_{\mathrm{fast}}(\rho ,0)$$

$$=$$

$${\nu }_{\mathrm{total}}\simeq 9.5\mathrm{Hz}.$$

$$\alpha$$

$$\gamma$$

$${\sigma }^2$$

$$1/\gamma$$

$$\sqrt{\alpha /\gamma }$$

$$\sqrt{{\sigma }^2/\alpha }$$

$$T=\frac{1}{\gamma +\alpha (2\pi /\lambda {)}^2}$$

Full text

Analysis and Modelling of Subthreshold Neural Multi-electrode Array Data by Statistical Field Theory

Måns Henningson

Department of Physics, Chalmers University of Technology

S-412 96 Göteborg, Sweden

mans@chalmers.se

Sebastian Illes

Institute of Neuroscience and Physiology, University of Gothenburg

Box 430, S-405 30 Göteborg, Sweden

sebastian.illes@neuro.gu.se

Abstract:

Multi-electrode arrays (MEA) are increasingly used to investigate spontaneous neuronal network activity. The recorded signals comprise several distinct components: Apart from artefacts 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 behaviour that is a general feature of theoretical models of this kind. We also find a clear correlation between the activity and the occurence of spikes. Another important insight is the importance of correcly separating out certain artefacts from the data before proceeding with the analysis.

1 Introduction

The multi-electrode array system (MEA) 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 [1]) and in in vitro neuronal cultures (e.g. from embryonic rodent brain tissue [2] or human stem cells [3]). This technology permits simultaneous long-term recordings from a fairly large number of extra-cellular electrodes. See e.g. [4][5] 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 $\mu 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. [6] 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. [7] 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 behaviour 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. [8]), as well as to develop the appropriate mathematical approaches applicable to describe their properties, are challenging issues in the research field of neuroscience.

The local field potential 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 $\mu V$, i.e. much less than the spikes. In contrast to the rather stereotyped spike waveforms from individual neurons, the local field potential has a much more stochastic appearance, reflecting its origins from a large and rather heterogenous ensemble of neurons [11]. Furthermore, spikes can only be recorded in relatively close vicinity to the recording electrode while the distance between recording electrode and local field potentials source can be several micrometers and millimeters [10]. Thereby, the local field potential recorded in isolated brain structures actually carries important information about the transport properties of the intra-cellular medium [9], the spatial and temporal structure of the neural activity [12][10], and also the connectivity of neurons within brain tissues [14]. There are different functional aspects of neuronal circuits which can be revealed by modelling and analyzing local field potentials. Current-source density analysis is used to reveal the neural source of recorded local field potentials [13] which is still controversial [15]. Since the pioneering work of Berger in the 1920s, local field potential band-separations techniques has been devoted to analyzing the local field potential in the frequency domain (see e.g. [16]) with the aim of identifying physiologically and pathophysiologically relevant frequency bands. In addition, decomposing local field potentials 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 local field potentials. 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 [13].

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 behaviour 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 artefacts 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 $\mathrm{\mu m}$ diameter arranged in a square grid pattern with 200 $\mathrm{\mu 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^{-16}\times 10\,\mathrm{mV}\simeq 0.15\,\mathrm{\mu V}$. An 0.3 mm thick acute hippocampal slice from a 44 day 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\times 10^{6}$ voltage measurements.

Conventionally, various filters are applied to the measured signals, but several studies demonstrate that such procedures do not remove spike components in local field potentials (subthreshould activity) [17][18][19], 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).

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 algoritm, which certainly leaves much room for improvements but is sufficient for our purposes: To detect spike events on an electrode at spatial point ${\bf r}$, we consider the difference $d({\bf r},t)$ between the potential $p({\bf r},t)$ at time $t$ and its average during a preceeding time interval of some length $\Delta t_{\mathrm{average}}$. We consider a spike to be fired at time $t$ if the magnitude $|d({\bf r},t)|$ of this deviation then attains its maximum in the time window of length $2\Delta t_{\mathrm{window}}$ centered at $t$ and exceeds a threshold value $d_{\mathrm{threshold}}$. In view of the typical ms timescale of the action potential dynamics, we used $\Delta t_{\mathrm{average}}=10$ ms and $\Delta t_{\mathrm{window}}=2$ ms. Concerning the threshold value, the value $d_{\mathrm{threshold}}=20\,\mathrm{\mu 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({\bf r},t^{\prime})$ in the interval $t-\Delta t_{\mathrm{spike}}<t^{\prime}<t+\Delta t_{\mathrm{spike}}$ for some time $\Delta t_{\mathrm{spike}}$ by the linear interpolating function

[TABLE]

We used $\Delta t_{\mathrm{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({\bf 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({\bf r},t)$ vanishes at all spatial points ${\bf r}$ and times $t$:

[TABLE]

All information is now contained in the two-point function $\langle p({\bf r}_{1},t_{1})p({\bf r}_{2},t_{2})\rangle$, 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.

[TABLE]

and

[TABLE]

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

[TABLE]

for some covariance function $S(\rho,\tau)$ 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({\bf r},t)$ fluctuates around a slowly varying equilibrium potential $\mu(t)$ that is more or less independent of the spatial position ${\bf r}$. We propose to describe this by a local, Gaussian, Markovian stochastic model of the form

[TABLE]

Here $\nabla^{2}=\sum_{i=1}^{2}\partial_{i}\partial_{i}$ is the Laplacian operator in two spatial dimensions. The relaxation constant $\gamma$ represents the tendency of the potential to return to its equilibrium value $\mu(t)$, and the diffusion constant $\alpha$ represents the tendency of spatial inhomogeneities to be smoothed out. The last term, which in stochastic modelling 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 $\mu(t)$ being larger than about 100 ms. For more background on statistical field theory, see e.g. [21].

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

[TABLE]

where the Green’s function

[TABLE]

obeys the differential equation

[TABLE]

and the initial condition

[TABLE]

Here and in the sequel, spatial integrals $\int d^{2}{\bf 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 10 and similiar equations below is that because of the linearity of the model (9), 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 9. (For a further discussion on Green’s function techniques for solving linear partial differential equations, see e.g. [22].)

The equilibrium potential $\mu(t)$ and the driving term $\xi({\bf r},t)$ are both assumed to have vanishing expectations values

[TABLE]

leading indeed to a vanishing expectation value for the potential $p({\bf r},t)$. We furthermore assume the covariance function of $\mu(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 $\mu(t)$:

[TABLE]

Here the constant $\sigma^{2}$ represents the intensity of the neural activity. Accordingly, we can now decompose the covariance function appearing in (8) as

[TABLE]

The first term in (19) 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

[TABLE]

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

[TABLE]

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 (19) represents the contributions from the driving term and can be expressed in terms of the intensity $\sigma^{2}$. A short computation gives

[TABLE]

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

[TABLE]

where $\Gamma$ is the (upper) incomplete Gamma-function and $\gamma_{\mathrm{EM}}=0.5772\ldots$ is the Euler-Mascheroni constant. Similarly, for vanishing temporal separation, i.e. $\tau=0$, we instead have

[TABLE]

where $K_{0}$ is a modified Bessel-function.

The most important aspects of the results (38) and (43) 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 behaviour 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 $\Delta t_{\mathrm{average}}=10$ ms, $\Delta t_{\mathrm{refractory}}=2$ ms, and $d_{\mathrm{threshold}}=20\,\mathrm{\mu V}$, the dataset had a total spike firing frequency of

[TABLE]

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

3.2 Artefacts

After the spikes had been removed, we computed the covariance function $S(\rho,\tau)$ 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 $\mathrm{\mu V}^{2}$. Two features of the covariance function $S(\rho,\tau)$ appeared to be artefacts 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 $\mu V^{2}$ that persisted essentially undamped until $\tau=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 $\rho=0$, i.e. at vanishing spatial separation, there was a component with a pronounced peak in the interval $0<\tau<0.2$ ms (i.e. during 5 sampling intervals) with a maximal amplitude of about 4 $\mu V^{2}$. 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 $\mu 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 artefact components in the covariance function. Henceforth $S(\rho,\tau)$ will always refer to the covariance function after these artefacts had been removed by subtracting the two temporal profiles exhibited in the figure from the raw-data covariance function.

3.3 Slow dynamics

For large values of $\tau$ ($>100$ ms), the main feature of the covariance function $S(\rho,\tau)$ was a damped oscillatory behaviour mainly in the 0 - 2 Hz frequency range, which was largely independent of the spatial separation $\rho$. 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. [20] for a review of local field potentials with different frequency bands.) See figure 7 for this long-time behaviour of the covariance function. We took $S_{\mathrm{slow}}(\tau)$ in (19) to be given by $S(\rho_{\mathrm{large}},\tau)$ for $\rho_{\mathrm{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_{\mathrm{fast}}(\rho_{\mathrm{large}},\tau)$ is negligible.

3.4 Fast dynamics

For small values of $\tau$ ($<100$ ms), the covariance function $S(0,\tau)$ at vanishing spatial separation indeed increased logarithmically as $\tau$ approaches zero. For $\rho>0$, this increase was cut off so that the equal time covariance $S(\rho,0)$ has a finite value that increases logarithmically as $\rho$ 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_{\mathrm{fast}}(\rho,\tau)=S(\rho,\tau)-S(\rho_{\mathrm{large}},\tau)$ were fitted to the the theoretical prediction (29). This is shown in figure 9 with the parameter values

[TABLE]

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

The definition and determination of the three quantities $\alpha$, $\gamma$ and $\sigma^{2}$ 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:

[TABLE]

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 $\mu(t)$ (which seems to be around 1 s). In this context, we remark that the time scale can be generalized to

[TABLE]

for fluctuations of some finite wavelength $\lambda$. In the long wave-length limit $\lambda\rightarrow\infty$ we recover $1/\gamma$, whereas for the inter-electrode distance $\lambda=0.2$ mm we instead get $T\simeq 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\mathrm{\mu V}$ is safely smaller than the spikes (which we have cut off at 20 $\mathrm{\mu V}$). However, it is quite comparable both to the errors in the individual voltage measurements (about 2 $\mathrm{\mu V}$) and the amplitude of the slow fluctuations of the equilibrium potential $\mu(t)$, so it is important to carefully separate these three phenomena.

3.5 Activity

Sofar we have considered the parameters $\alpha$ and $\gamma$ as well as the activity $\sigma^{2}$ to be constants. This is reasonable for $\alpha$ and $\gamma$, 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 $\sigma^{2}({\bf 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 $\sigma^{2}\simeq 0.035\,\mathrm{\mu V}^{2}\mathrm{mm}^{2}\,\mathrm{ms}^{-1}$ that we have determined should thus be regarded as a spatial and temporal average.

Retracing the steps leading to (29), we find that with a non-constant activity $\sigma^{2}({\bf r},t)$, this expression is no longer valid. However, the leading logarithmic divergence of (38), which originates from the short distance behaviour of the model, still holds. Since $S_{\mathrm{slow}}(t)$ is regular for small $t$, we thus have

[TABLE]

The diffusion constant $\alpha$ is of course already known. The values of $\delta t$ can e.g. be chosen in the interval 0.2 ms to 10 ms. Since we have only a single measurement of the potential $p({\bf r},t)$ for each value of ${\bf 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}({\bf 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.

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

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({\bf r},t)$, which is a stochastic variable, we have considered the activity $\sigma^{2}({\bf 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({\bf r},t)$, it shows essentially no correlation with the mean of $p({\bf r},t)$. This is another indication that the equilibrium potential $\mu(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 $\sigma^{2}({\bf r},t)$ is its covariance function between separate points ${\bf r}_{1}$ and ${\bf r}_{2}$ at some time lag $\Delta t$:

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

with decay constant

[TABLE]

(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

[TABLE]

Taking the spatial mean, we here find

[TABLE]

with decay constant

[TABLE]

See figure 12 for the corresponding autocorrelation functions (normalized to $1$ for $\Delta t=0$ and $\Delta{\bf 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.

The behaviour 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 behaviour characteristic of passive transport in two spatial dimensions as we have investigated for the potential $p({\bf 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.91\mathrm{mm}$ 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 $|\Delta{\bf 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 towards 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.

4 Discussion

Our main finding is that the local field potential 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 modelled 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 preceeded by a correct elimination of certain artecfacts 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. [23] 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.

Acknowledgments

We have benefitted from discussions with Eric Hanse.

The research of M.H. was supported by a grant from the Kristina Stenborg foundation.

The research of S.I. was supported by a grant from the Alzheimerfonden (AF-556051).