Inverse Current Source Density Method in Two Dimensions: Inferring Neural Activation from Multielectrode Recordings

To understand the brain at the systems level one needs precise information regarding the spatial and temporal activation of different neuronal populations. The recent developments in multielectrode construction have opened new possibilities for the electrophysiologist, providing the means to record extracellular potentials at tens to hundreds of closely spaced positions simultaneously (Csicsvari et al. 2003; Barthó et al. 2004; Buzsáki 2004; Blanche et al. 2005; Du et al. 2008). To take full advantage of this new technology, new data analysis methods must be developed to extract useful and precise information from the masses of data that we can now record relatively easily.

The high-frequency part of extracellular signals contains information about firing of action potentials in neurons within a distance of 0.1 mm or so around the individual recording contacts (Buzsáki 2004; Pettersen and Einevoll 2008). The low-frequency part, the local field potential (LFP), thought to be generated by synaptic inputs and their dendritic return currents, has a larger spatial spread due to volume conduction (Mitzdorf 1985; Pettersen et al. 2008; Katzner et al. 2009; Xing et al. 2009). The standard method of analysis for LFP has been to estimate the current-source density (CSD), i.e., the net volume density of current entering or leaving the extracellular tissue (Lorente de No 1947; Pitts 1952; Plonsey 1969; Nicholson and Freeman 1975; Freeman and Nicholson 1975; Mitzdorf 1985).

The Poisson equation provides the connection between the extracellular potential Φ and the current source density C under the assumption of passive spread in an homogeneous and isotropic medium: σ ΔΦ = –C, where Δ is the Laplace operator. To estimate C one may use the numerical second derivative in place of the Laplace operator (Pitts 1952; Nicholson and Freeman 1975; Freeman and Nicholson 1975). This approach has been commonplace in the analysis of recordings from linear (laminar) electrodes inserted perpendicular to cortical layers (e.g. Haberly and Shepherd 1973; Mitzdorf 1985; Schroeder et al. 1992; Ylinen et al. 1995; Lakatos et al. 2005; Lipton et al. 2006; Rajkai et al. 2008; de Solages et al. 2008). In this setting one has assumed that the extracellular potential (and by implication the current-source density) is constant in the lateral directions, an approximation that cannot always be justified (Nicholson and Freeman 1975; Pettersen et al. 2006; Einevoll et al. 2007; Pettersen et al. 2008; Katzner et al. 2009; Xing et al. 2009). Another problem with the standard numerical derivative approach is that it is impossible to estimate CSD at the outermost electrode contact positions unless one makes assumptions unjustified with complex electrode geometries (Vaknin et al. 1988; Łęski et al. 2007). Giving up the boundary becomes a particularly severe issue in two- and three-dimensional electrode geometries where the relative number of such boundary contacts becomes large (Łęski et al. 2007).

The estimation of current sources from recorded potentials has a long history in the interpretation of EEG (Guljarani 1998; He and Lian 2005; Nunez and Srinivasan 2006) and ECoG signals (e.g. Freeman 1980; Zhang et al. 2008). The inverse methods involve calculating a forward model of propagation of electric fields from the sources inside the brain to the recording electrodes on the scalp (EEG) or cortical surface (ECoG). The model is then inverted to estimate the sources from the recorded potentials. For these signals, which are recorded outside the source region, one common approach is to assume their sources to be a small number of mesoscopic dipoles located so far away from the electrode contacts that the far-field approximation can be evoked. However, more general source distributions have also been considered (for review see e.g. He and Lian 2005; Nunez and Srinivasan 2006).

For multielectrode contacts positioned inside neural tissue in the immediate vicinity of the neuronal sources, the far-field dipole approximation for calculating the forward solution is not applicable. Moreover, the approximation where one assumes the neuronal sources to be built up of pairs of two balanced current monopoles is unsuitable. In Lindén et al. (2010) (see Fig. 6) it was shown that both of these simplified source representations provided poor approximations of the local field potential (LFP) generated by neurons within distances of a millimeter or so. A better source representation for this situation is the continuous current-source density, and we have developed a forward-inverse scheme for CSD estimation called the inverse current source density (iCSD) method (Pettersen et al. 2006; Łęski et al. 2007; Wójcik and Łęski 2009). A main advantage of the iCSD method compared to the standard CSD analysis is that assumptions about the geometry of the CSD as well as electrical boundary conditions can be incorporated explicitly in the estimator.

The iCSD method assumes a specific form of the distribution of the current sources, for instance, linearly varying between the recording points. It connects parameters of the CSD with the potentials at the grid through forward calculation of the potential generated by the assumed distribution. While the idea is general and applicable to different geometries, until now it has been developed only for recordings with laminar (linear) electrodes (Pettersen et al. 2006) and recordings on a three-dimensional grid of electrode positions (Łęski et al. 2007).

Given the growing availability of two-dimensional electrode arrays, such as multi-shank laminar electrodes (Buzsáki 2004), it is presently important to provide efficient and general methods for CSD estimation from such recordings. Since the 3D iCSD cannot be applied directly to such recordings (as it requires at least a few contacts in all three spatial directions) in the present paper we develop a specific iCSD framework for the analysis of data from two-dimensional multielectrode arrays. This 2D iCSD method is validated on LFP data generated by several types of test sources where we know the true underlying source distribution and, thus, can make a quantitative assessment of the estimation accuracy. We first consider various types of Gaussian current sources, and then move on to analyze LFP generated by synaptically activated model populations consisting of about one thousand morphologically-reconstructed layer 5 neocortical pyramidal neurons (Mainen and Sejnowski 1996; Pettersen et al. 2008). For all tests we compare our iCSD estimates with estimates from the traditional CSD method. The proposed framework is further illustrated by application to a set of recordings with an 8 × 8 multi-shank electrode from the rat subiculum.

A key feature of the iCSD method is the explicit incorporation of assumptions regarding geometries of the underlying CSD sources. For the present situation with a 2D multi-electrode grid, a key parameter is the assumed spatial spread of CSD in the direction perpendicular to that grid (the standard CSD method neglects the variation of extracellular potential in this direction, an assumption that is physically unrealistic in that it does not correspond to any known CSD distribution). This is also the main difference between the 2D and the 3D iCSD; in the 3D case the behavior of CSD in the direction perpendicular to the main plane is calculated, whilst here it must be modeled. Usually, it is not clear a priori, which value of this parameter leads to optimal reconstruction. To circumvent this we found it useful to make iCSD estimates assuming different parameter values, and investigate how the salient features of the CSD pattern depend on the choice made (see Freeman 1980 for similar analysis of the influence the “focal depth” parameter on the analysis of ECoG data). To facilitate this approach we have developed a MATLAB toolbox with a simple graphical user interface to allow users to easily and rapidly investigate the dependency of the 2D iCSD estimates on both this parameter as well as the choice of boundary conditions. The toolbox bundled with three of the datasets used in this paper is provided under GNU General Public License v.3 or later and is available from the INCF software repository, http://software.incf.org/. The software and data can be used in published research provided this article is cited.

Inverse Current Source Density in Two Dimensions

Suppose that we measure the electric potential Φ at N points, x ₁, x ₂, …, x _N. Potentials are functions of the density of current sources C according to the Poisson equation (Nicholson and Freeman 1975):

where σ is the conductivity of the tissue, here assumed to be homogeneous and isotropic.

The main idea behind iCSD is to assume a specific distribution C(x,y,z) of the current sources from a class parameterized with parameters C = [C₁, C₂, …, C_N]. Once a model of CSD is assumed, it is a simple matter to evaluate the potentials, which would be measured at any point in space, by summing the contributions from every point source:

This leads to the relation

and if we choose the appropriate model of the CSD distribution, the operator F is linear and can be inverted. We thus obtain values of N parameters in terms of the measured potentials:

which gives the model CSD in its whole domain.

In this work we consider potentials measured on a two-dimensional regular grid of electrodes of N = n _x n _y points. Let us position the axes so that the plane of the grid corresponds to z = 0 and the electrode positions are

where Δx and Δy is the spacing in the x and y directions, m = 1,…,n _x , n = 1,…,n _y . We can uniquely decompose the CSD that generated the measured potentials into parts symmetric and anti-symmetric in z

where \( {{\text{C}}_{\text{S}}}\left( {x,y,z} \right) = \left( {C\left( {x,y,z} \right) + C\left( {x,y, - z} \right)} \right)/2,{C_A}\left( {x,y,z} \right) = \left( {C\left( {x,y,z} \right) - C\left( {x,y, - z} \right)} \right)/2 \). Substitution to Eq. 1b shows that the potentials measured in the plane of the grid come only from the symmetric part, so this is the only thing we can hope to reconstruct.

Our goal is to recover C _S (x,y,0). In the generic situation, we have no additional information about the distribution of CSD along the z axis, apart from the fact that it has finite extent. Therefore, we make the simplest possible assumption of a model CSD distribution which is a product of an a priori unknown two-dimensional profile c(x,y) and a specific symmetric profile in the perpendicular direction, i.e.

It is convenient to normalize the z-profile H(z) so that H(0) = 1.

The proposed approach would work best if anatomy of the probed volume suggests such a product structure of the CSD, for example, in case of a grid inserted perpendicularly to a layered brain structure. This is similar to the one-dimensional case where the most meaningful estimates of CSD can be obtained in laminar structures, such as cortex (Nicholson and Freeman 1975; Pettersen et al. 2006). If substantial extra knowledge of anatomy of the neighboring tissue or the processes taking place are known, one can develop specific variants of the method building more complex models of the assumed CSD. In general, however, we believe the assumptions we make are reasonable and, as our tests in the following pages show, the proposed approach does work adequately on a number of model sources.

For H(z) we take here either the step function, H(z) = 1 for  − h ≤ z ≤ h and H(z) = 0 otherwise, or a Gaussian: H(z) = exp(−z ² /2h ²). We choose the normalization H(0) = 1 so that c(x,y) = C(x,y,0). The function c(x,y) is in the simplest case an interpolation between the nodes. This could be either nearest neighbor interpolation (we set the value at any spatial point equal to the value at the nearest node), linear interpolation (the function describing values between the nodes is piecewise linear), or spline interpolation (which uses third-degree polynomials to produce distributions which vary more smoothly than in the linear case). In the simplest case we set c(x,y) to zero for positions outside the grid (we will refer to this case as “no boundary conditions”).

We showed for the three-dimensional case (Łęski et al. 2007) that the “no boundary conditions” assumption can lead to large reconstruction errors if some of the sources generating the fields are located outside the grid. This is because the reconstructed CSD tries to compensate for the effect of external sources with components within the grid, which leads to artifacts. One possible solution is to consider an extended layer of grid points, as shown in Fig. 1. Then at each of these additional points we can either (i) set the CSD to zero (which we denote by B boundary conditions), or (ii) copy the value from the nearest point of the original grid (denoted by D). For the latter case the values at additional points are not fixed. Such a procedure is well-defined for all points of the extended grid, including corners. We found that both variants of this approach improve the reconstruction fidelity within the grid (Łęski et al. 2007). In all cases the distribution of sources is completely described with N parameters.

Comparison of ‘no boundary treatment’ vs. b or d distribution of sources. In the distribution with no boundary treatment (a) the CSD is non-zero only inside the rectangular grid. To accommodate for sources lying outside of the grid, and to avoid artifacts in the reconstructed CSD, a grid with an additional layer of nodes is used (b), which results in an additional layer of non-zero CSD

Full size image

We now consider the situation where the CSD is assumed to be ‘step-wise’ constant (Pettersen et al. 2006), i.e., the CSD is assumed to be constant within the rectangle in the xy-plane assigned to each grid point. For this model the matrix F can be calculated as follows: Denote by (x _i , y _i ) and (x _j , y _j ) the positions of the i-th and j-th grid points, respectively, i and j go from 1 to N. The N parameters of the CSD distribution will be the values of C at the nodes, which we denote by C _i , i goes from 1 to N. The potential at the node i, Φ _i , is equal to

The matrix element F _ij is:

The construction of the matrix F for linear and spline interpolation is given in the Appendix. The procedure is analogous, but the calculations are more complex.

To compare the accuracy of the iCSD method quantitatively against the traditional CSD estimation method, we had to extend the traditional CSD method to predict spatially continuous CSD profiles. The following procedure was used: the grid of electrodes was extended (cf. Fig. 1) and at every point of this extension we copied the potential from the nearest point of the original boundary (Vaknin et al. 1988), which is non-ambiguous. The numerical second derivative was then calculated at the points of the initial grid and spline-interpolated in between.

Generation of Population Model Data

In order to test the new 2D iCSD method we generated model data for columnar populations of reconstructed layer-5 pyramidal cells receiving a combination of excitatory and inhibitory inputs resembling stimulus-evoked activation (Pettersen et al. 2008). In the simulation, we know the actual CSD generating the local-field potentials recorded (virtually) at the grid points of the multielectrode. We can use these to quantify the quality of CSD estimates from the recorded LFPs reliably. A similar procedure was used by Pettersen et al. (2008) to test the 1D iCSD method. We studied CSD and LFP generated by three such model columns spaced equally along a line as one “central” column surrounded by two “lateral” columns, and we assumed two positions of the virtual electrode grid with respect to these columns (see text below and Fig. 6a, b, for details). Two synaptic stimulation patterns were considered.

Neuron templates

A cell population resembling a layer-5 pyramidal-cell network in a neocortical column was built based on the compartmental model from Mainen and Sejnowski (1996). The compartmental model was downloaded from SenseLab’s ModelDB (http://senselab.med.yale.edu/; Hines et al. 2004; Migliore et al. 2003), and the simulation tool NEURON (http://www.neuron.yale.edu/; Carnevale and Hines 2006) was used to compute the neuronal dynamics. The neuron model has various active conductances in the axon segments, axon hillock, soma and dendrites. Similar to Pettersen et al. (2008), the original model of Mainen and Sejnowski (1996) was modified as follows: the active conductances in the dendrites were removed, the electrode was removed from soma, the whole neuron was rotated so that the primary dendrite was aligned to the positive y-axis, and the axon was then aligned straight downward along the negative y-axis.

Two different stimulus patterns were used for the pyramidal neuron model: one with combined apical excitation and basal inhibition, and one with combined basal excitation and basal inhibition. These two patterns correspond to stimulus patterns 2 (SIP2) and 4 (SIP4) in Pettersen et al. (2008). In both cases excitation and inhibition were tuned so that the soma potential was just below threshold value so that the neuron did not produce any action potentials.

Synaptic input was density-based, i.e., not based on point processes. Thus, the synaptic input was considered to be scattered throughout the whole branch to which it was applied, similar to Holt and Koch (1999). The excitatory synaptic input was conductance-based with an exponentially decaying temporal profile,

Here \( g_j^{\rm{e}} \) is the synaptic membrane conductance in branch j of the neuron, \( g_{{{ \max }}}^{\rm{e}} \) is the maximum conductance, τ _e is the time constant of the excitatory input, \( \Delta_j^{\rm{e}} \) is the onset time of synaptic input in branch j, and θ(t) is the Heaviside unit step function. \( \delta_j^{\rm{e}} \) is unity if branch j is set to receive excitatory input in the model, and zero if not. The synaptic input patterns, SIP2 and SIP4, both had the same excitatory time constant, τ _e  = 5 ms. The maximum conductance for the apically excited neurons (SIP2) was \( g_{{{ \max }}}^{\rm{e}} = 0.001 {\hbox{S/c}}{{\hbox{m}}^{{2}}} \) while the maximum conductance for the basally excited neurons (SIP4) was \( g_{{{ \max }}}^{\rm{e}} = 1 \times {10^{{ - 4}}} {\hbox{S/c}}{{\hbox{m}}^{{2}}} \).

The basal inhibitory input was similarly given by

An inhibitory time constant τ _i = 10 ms was used. The maximum conductance \( g_{{{ \max }}}^{\rm{i}} \) was adjusted so that the neuron did not produce action potentials. SIP2 had a maximal inhibitory membrane conductance of \( g_{{{ \max }}}^{\rm{i}} = 5 \times {10^4} {\hbox{S/c}}{{\hbox{m}}^2} \), while for SIP4 \( g_{{{ \max }}}^{\rm{i}} = 3 \times {10^4} {\hbox{S/c}}{{\hbox{m}}^2} \). For both stimulus patterns the inhibitory synaptic stimulus was applied to the soma and each branch of the basal dendrites, i.e., \( \delta_j^{\rm{i}} \) was unity only for the soma and these dendritic branches j. The basal excitation in SIP4, in contrast to the basal inhibition, did not include the soma, the first branch of the dominant apical dendrite, nor the side branches to this apical dendrite. The apical excitation of SIP2 was applied to all dendritic branches above the first branching point of the apical dendrite.

The total synaptic transmembrane current in segment k, being a part of branch j, is then given by

where E ^e = 0 is the reversal potential for the excitatory synapses, and E ^e = −70 mV is the reversal potential for the inhibitory synapses. V _k is the membrane potential in segment k of the neuron. The soma resting potential for the pyramidal neuron templates was −64.5 mV for both stimulus patterns. To include some temporal jitter in the onset of synaptic inputs, each branch’s synaptic input was started separately in time, i.e., the Δ _j 's were slightly different: the onset was chosen stochastically assuming a Gaussian distribution around the mean \( {\overline \Delta_{\rm{syn}}} \) with a standard deviation of \( {\sigma_{\rm{syn}}} = \sqrt {5} \) ms. The dynamical response of the pyramidal neuron to a particular synaptic-input pattern was computed in NEURON, and the transmembrane currents for all segments were written to file with a ten-digit precision. The extracellular potential for the neuronal templates was then computed using the line-source method (Holt and Koch 1999). The extracellular potential \( {\varphi_n}({{\bf r}},t) \) from a neuronal template n is then found by a sum over all segments of this neuron, i.e.,

where N _k is the number of segments in the pyramidal neuron, \( \Delta {s_k} \) is the length of line segment k of this neuron, r _nk is the radial distance (perpendicular to the segment) from segment number k, h _nk is the longitudinal distance (parallel to the dendritic segment) from the end of segment number k, \( {l_{{nk}}} = \Delta {s_k} + {h_{{nk}}} \) is the longitudinal distance from the start of the segment to the recording point, and I _nk (t) is the transmembrane segmental current (the ionic currents plus the capacitive current). The extracellular conductivity used in the simulations was σ = 0.3 S/m (Hamalainen et al. 1993). The extracellular population activity was calculated by linear superposition of the single pyramidal-neuron templates.

Structure of Population and Population Activity

The modeled layer-5 pyramidal populations have the typical sizes and spatial extensions of cortical columns, e.g., as seen for the layer-5 pyramidal population in rat barrel cortex (Beaulieu 1993; Feldmeyer and Sakmann 2000). They contained 1,000 pyramidal neurons, each randomly rotated around their primary dendrite and receiving the same spatial stimulus pattern. Their somas were positioned stochastically with uniform probability density within a cylinder oriented along the y-axis with height 0.4 mm and a diameter of 0.4 mm (cf. Fig. 5a). The mean synaptic onset times \( {\overline \Delta_{{{\rm{syn}},n}}} \) of the 1000 neurons in each population were chosen stochastically from a Gaussian distribution with a standard deviation of σ_sp = 5 ms. The probability distribution was truncated, i.e., set to zero, for times 2σ_sp = 10 ms smaller or larger than the mean value. This gives a maximum temporal translation of 20 ms between the neurons within the population.

Computation of Extracellular Potentials from Population Activity

Based on the set of stochastically chosen mean synaptic onset times (\( {\overline \Delta_{{{\rm{syn}},n}}} \)) for all 1,000 neurons in each population and the extracellular single-neuron templates in Eq. 4, the extracellular potential from the entire population was calculated as

where N = 1,000 is the number of neurons in the population, and ϕ _n (r, t) is the extracellular signature of neuron n with mean synaptic onset at time zero. The potential was computed at 1656 assumed recording positions defined by a 9 × 23 × 8 grid where the assumed recording positions in x-direction were from −0.7 mm to 0.7 mm with inter-contact distance of 0.2 mm, in the y-direction from −0.8 mm to 1.4 mm with an inter-contact distance of 0.1 mm, and in z-direction from −0.8 mm to 0.8 mm with inter-contact distance of 0.2 mm. In the y-direction the populations were centered so that the average soma position was in y = 0.

We tested several spatial placements of the populations of model neurons, which we call ‘central’ and ‘lateral’ (Fig. 6a, b). For the central population the horizontal centering was in (x,z) = (0,0), while the lateral populations were centered in (x,z) = (0,–0.6) mm and (x,z) = (0,0.6) mm. The lateral populations centered in (x,z) = (0,0.6) were produced by mirroring the populations centered in (x,z) = (0,–0.6) mm about the xy-plane, i.e. they shared the (mirrored) stochastic parameters discussed above (somatic placement, orientation, synaptic onsets and neuronal time-shifts). The center population was created with a new set of stochastic parameters.

To avoid singularities in Eq. 9 no dendritic segments were allowed to be closer to the assumed recording contacts than 5 μm; dendritic segments closer than this distance were given a radial distance of r _nk  = 5 μm when the extracellular potential was computed through Eq. 9.

Computation of the Actual Model CSD

When computing the potentials at assumed recording positions, all neural segments except for the somas were treated as linear current sources (Eq. 9). However, when computing the actual CSD for these populations, the current-sources were treated as point sources to improve computational efficiency. Each point source was then positioned at the center of its corresponding line segment.

The CSD is defined as the total current source within a small volume element, divided by the volume of this element. Because of the inhomogeneous nature of the CSD one cannot choose this volume element to be arbitrarily small without accepting a high degree of spatial noise due to the high variance in the number of sources in neighboring cubes. For this model study we computed the real CSD based on the total current-source within cubes with sides of 50 μm, which gave an acceptable high resolution without too much spatial noise.

Experimental Procedure

Rat Subicular Recordings

Adult male Wistar rats (200–400 g) were anaesthetized using urethane (1.5 g/kg; 30% w/v i.p.) and then placed in a stereotaxic frame with body temperature maintained at 37°C using a homoeothermic blanket (Harvard Apparatus, UK). A midline scalp incision was made to expose the skull and craniotomies were made above dorsal CA1 (AP: −4.5 mm, ML: 3.0 mm) and subiculum (AP: −8.4 mm, ML 3.5 mm) relative to Bregma (Gigg et al. 2000). The dura was then excised to allow insertion of electrodes. All procedures were in accordance with the Animals (Scientific Procedures) Act 1986, UK.

Electrodes and Stimulations

Dorsal subicular recording were made using a 64-contact probe (a8 × 8–5 mm 200-200-413, NeuroNexus Technologies, Michigan, USA). The configuration of the probe provided eight 413 μm² contacts at 200 μm spacing on each of the eight 5 mm long shanks. The probe was inserted at a 25° angle from vertical in the sagittal plane so that at target the shanks of the electrode were approximately perpendicular to the main axis of the subicular cell layer. A bipolar stimulating electrode consisting of two, twisted Teflon coated, tungsten microwires (125 μm diameter, insulated to the tips; Advent RM, UK) was placed in the alveus above dorsal CA1 (Fig. 2).

Histological verification of stimulation and recording electrodes (a). White circles show placement of electrode contacts (−8.2 from Bregma, 2.9 lateral to midline). The white square in the inset shows the position of the stimulating electrode during the recording (−4.5 mm from Bregma, 2.9 mm lateral to midline). The b panel shows superimposed mean evoked voltage responses at each electrode position. Maximal deflections occur mainly in the subicular region. Stimulating and recording electrode sites are identified on the nearest slice image to both points using Paxinos and Watson 1998. Black bars show 1 mm scale. Alv: alveus, CA1, CA3, DG: dentate gyrus, Post: postsubiculum

Full size image

Stimuli were triggered by 5 V analog pulses from a National Instruments card (PCI-6071E), controlled by programs written in LabVIEW (v8, National Instruments, USA). These pulses triggered a constant-current, isolated stimulator (DS3, Digitimer Ltd., UK). Stimuli were of fixed duration (0.2 ms). Once electrodes were at target an input–output curve was plotted of alvear stimulus intensity versus subicular response amplitude (e.g., Commins et al. 1998). All subsequent pulses were then set at half the intensity required to elicit the maximum response from the curve, in this case, 70 μA.

Single stimulation pulses were presented at a rate of 0.33 Hz for one minute. Perievent histograms of the mean fEPSP voltage response from these were calculated for every channel. As two channels had to be used for stimulation, the final result was an eight-by-eight contact profile minus two silent contacts. For the application of iCSD methods we filled these missing grid points with mean voltages from nearest neighbors. We have shown elsewhere that this procedure does not introduce much distortion for a small number of missing channels (Wójcik and Łęski 2009).

Data Acquisition

Signals were acquired using a Recorder64 (Plexon, USA) recording system, referenced to ground and amplified at source using a 20x gain AC-coupled headstage followed by pre-amplifier conditioning (total gain of 2500x). Other than the fixed system low-pass (6 kHz) no other filtering was applied. Local field potential signals were digitized at a sampling rate of 10 kHz per channel at 12 bit resolution and stored for offline analyses.

Histological Verification

At the end of the experiment, recording sites were determined by a combination of visual analyses of electrode tracks and lesions placed on the upper-most and lowest electrode contacts using a 30 μA positive current for 3 seconds (Townsend et al. 2002). Animals were subject to terminal anesthesia (2–3 ml of 30% urethane i.p.) and transcardially perfused with 100 ml of 0.9% saline followed by 150 ml of 10% formalin. The brain was removed and stored in 10% formalin for 24 hours, followed by immersion in 30% sucrose solution. Frozen 100 μm sections were made in the sagittal plane. and stained with Cresyl violet. Electrode placements were assigned with reference to the rat brain atlas of Paxinos and Watson (1998).

In this section we study the properties of the proposed method on several different datasets before we finally apply it to experimental data. We start with tests of reconstruction quality for surrogate data with a simple structure (Gaussian sources) and compare variants of the iCSD method with the traditional approach for CSD analysis. Then we analyze more complex data from simulation of columnar populations of layer-5 pyramidal neurons. Two different patterns of stimulation (basal and apical) are tested. We then observe and discuss the asymptotic independence of reconstructed error (large h limit). With the viability of our approach assured by tests on these model data, we then use the best variant of our method to analyze responses evoked by alvear stimulation recorded on 8 by 8 grid in the rat subiculum.

Tests on Gaussian Sources

Throughout the tests on Gaussian sources we assume that the potentials are measured on a grid of 8×8 electrodes, spaced by 0.2 mm in x and y directions. This choice was motivated by the experimental conditions used later.

Two-Dimensional Gaussian Sources

First we test the method on surrogate sources which have a product structure c(x,y)H(z), where c(x,y) is a sum of Gaussians, and H(z) is a step function with h = 0.5 mm (the exact formula for c(x,y) is given in the Appendix). We calculate the potentials generated by such sources on the assumed electrode grid. In the simplest case we assume that sources are non-zero only for x and y inside the electrode grid (Fig. 3a). For this case both linear and spline interpolation iCSD methods perform very well, provided we assume the correct value for h (Fig. 3c, d). The traditional (numerical second derivative) method also locates the sources, however, their shapes are distorted (Fig. 3b). To quantify the reconstruction fidelity we use normalized squared error:

Test of iCSD methods on sources based on two-dimensional Gaussian functions (arbitrary units). a original Gaussian sources (zero outside the box, product structure, h = 0.5 mm, potentials are measured on a grid of 8×8 electrodes), b reconstruction using traditional CSD method, that means numerical double derivative (smoothed), c reconstruction with linear iCSD method, d reconstruction with spline iCSD method e–g: the sources are the same as in a, but with h = 0.1 mm. Reconstruction of CSD (spline interpolation) with assumed h equal to 0.05 mm (e), 0.1 mm (f), 0.2 mm (g). Note different scales of the three plots: besides distortion of shapes of the sources there is also a global scaling due to different assumed values of h h–k: the sources inside the box are the same as in a, but now they are also non-zero outside the box. Reconstruction with spline iCSD method: h with no boundary treatment, i with b boundary conditions, j with d boundary conditions, k reconstruction using smoothed numerical double derivative

Full size image

where C(x, y, z = 0) are original sources and \( {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over C} }}(x,y) \)is the reconstructed CSD in the plane; M is a normalization constant. Since we are mostly interested in the spatial profile of the CSD, as the precise values are unavailable because the conductivity is usually unknown, we also studied the differences between the original profiles and scaled reconstructions. Specifically, as another measure of reconstruction error, we took:

where the integrals are over the area spanned by the electrode grid, and M is a normalization constant defined as previously.

For the traditional CSD method, the error e ₁ is 34%, while it is only 0.097% for the linear iCSD method and 0.019% for the spline iCSD method. Thus, if the assumed distribution is sufficiently smooth, one can reconstruct the original pattern faithfully from the limited information provided by the finite set of measurements. Some of the large error inherent in the traditional CSD method is due to incorrect estimation of the CSD amplitude and not the form of the spatial profile. However, even if the CSD amplitude is adjusted to minimize the error according to Eq. 13 one still finds an error of e ₂  = 32%, i.e., much larger than for the iCSD methods.

One of the reasons for using the inverse CSD method instead of the numerical second derivative, especially in two- or three-dimensional situations, is that the boundary data are better utilized. To illustrate this we calculated the errors (e ₁ ) of the reconstruction of the central part of the grid, i.e., the grid spanned by the electrodes 2–7 both in x and y. Not surprisingly, in all cases the errors are smaller: : e1 = 8.1% for the traditional approach, 0.069% for the linear iCSD and 0.0063% for the spline iCSD.

Usually we do not know the correct h and we must form an educated guess based on the available information (Fig. 3e, f, g). In the source data we used h = 0.1 mm and for reconstructions in Fig. 3e, f, g we took h equal to 0.05 mm, 0.1 mm and 0.2 mm, respectively. If the assumed h is different from the true h, then slight distortions appear in the shape of the reconstructed CSD. The errors are: e ₂  = 0.4% for assumed h = 0.05 mm, e ₂  = 0.019% for the assumed h equal to the true h (= 0.1 mm), and e ₂  = 2.1% if we assume h = 0.2 mm. Moreover, the amplitude of the reconstructed distribution varies with the assumed h (the CSD is scaled by a global factor which in this case is approximately equal to the ratio of true and assumed h).

In the experimental setting the sources often extend beyond the electrode grid. Hence, we tested the iCSD methods on potentials calculated using the complete spatial extent of the Gaussian sources. The situation now changes dramatically and the variants of the iCSD method assuming non-zero distribution only inside the grid work very poorly (Fig. 3h): the reconstruction error is almost 500%. However, we can overcome this by using B or D boundary conditions (Fig. 3i, j), which leads to reconstruction errors e ₁ of 8.4% and 2.4%, respectively. The traditional method (Fig. 3k) gives a reconstruction error of e ₁  = 24% (e ₂  = 19%), which is much better than iCSD with no boundary conditions, but substantially worse than iCSD with proper boundary treatment. If we consider only the central grid, then the e ₁ error for iCSD assuming no sources on the outside is 19%, for B and D boundary conditions e ₁ is 1.3% and 0.29%, respectively, and for the traditional CSD we get e ₁  = 11%.

Three-Dimensional Gaussian Sources

The most interesting test case of Gaussian data consists of truly three-dimensional sources. The set of test sources C we use here is a sum of four, three-dimensional Gaussians whose centers are at z ranging from −0.3 mm to 0.6 mm. The sources were chosen in such a way that in the plane of the grid, z = 0, the distribution was the same as in the product case studied before (see Fig. 3a, the exact formula is given in the Appendix). We also used the same two-dimensional 8×8 grid of electrodes. The questions of interest are now what h and which boundary conditions would work best. Our tests on product sources (Fig. 3e, f, g) suggest that the reconstructed CSD will be defined up to a multiplicative constant depending on the assumed h. Therefore, to compare the quality of reconstruction for different h we have to scale the reconstructed CSD by a constant and so in this section we use the error e ₂ as defined by Eq. 13.

We performed the reconstructions assuming h = (0.05 × 2ⁿ) mm, n = 0,…,6. The reconstruction errors for the method with D boundary conditions are shown in Fig. 4d where the two curves are for H(z) being either a step function or a Gaussian (see figure caption). To understand the meaning of different values of reconstruction error we plot the actual CSD in the plane of the (virtual) recording grid (Fig. 4a) and two examples of reconstructed CSD (with α set to the values found to minimize the error for the two different choices of h). Fig. 4b shows results for h = 1.6 mm which corresponds to the smallest error of reconstruction in the studied regime (e ₂  ~ 10%). Figure 4c shows the reconstructed CSD for h = 0.1 mm which corresponds to large reconstruction error (e ₂  = 20%). From this test one may conclude that these two different choices for the H(z) profile (step and Gaussian profiles) give similar results. The reconstruction error is high for very small h and reaches minimum for h ~ 1 mm, which is roughly the extent of the test Gaussian sources. For comparison, the error for the traditional CSD method is in this case e ₂  = 26%.

Reconstruction error (Eq. 13) for inferring two-dimensional current source-density distribution on a 2D section through a set of three-dimensional Gaussian sources. The spline iCSD method with D boundary conditions was used. a Original CSD in the plane of the measurement grid (arbitrary units). b Reconstructed CSD for h = 1.6 mm which corresponds to the smallest error of reconstruction in the studied regime (e ₂  ~ 10%). c Reconstructed CSD for h = 0.1 mm which corresponds to large error of reconstruction (e ₂  ~ 20%). d Dependence of reconstruction error (e) on the assumed width of the sources (h). The dotted curve is for H(z) being a step function, for the solid curve H(z) is a Gaussian. The reconstructions in b and c are scaled with a parameter α to minimize the reconstruction error, see text for details. In b and c the function h(z) is a step function

Full size image

Tests on Population Model Data

For tests on our population model data, i.e., the data from forward modeling on a columnar population of layer-5 pyramidal neurons, we used potentials calculated on an array of 8 × 23 positions, spaced by 0.2 mm in x direction and by 0.1 mm in y direction. Here, as in the case of three-dimensional Gaussian sources, we can reconstruct the CSD up to a multiplicative constant (at least if we are not in the large-h limit, see the next section). However, unlike the previous case, we can now reconstruct the whole time-course of the activity. Therefore, as the error measure we can take

where the double integral is over the area spanned by the electrode grid. To calculate this integral we need the actual model sources in the grid plane, i.e., C(x,y,z = 0), and for this we used the actual model sources averaged locally in cubes of edge length 0.05 mm (see Materials and Methods). To obtain \( \hat{C}(x,y) \) we simply used the value of the reconstructed sources in the middle of the cube. The normalization constant M was

The (time-independent) constant α in Eq. 14 was chosen to minimize the error e ₃ . In Fig. 5 we present reconstructions of activity of a simulated single column of layer-5 pyramidal neurons (Fig. 5a) in the planar section defined by the virtual recording electrodes. The electrode grid is placed along the axis of the apical dendrites of the simulated population and passes through the center of the cylinder in which the somas are located.

Placement of the simulated population of 1000 layer-5 pyramidal neurons a. All somas are contained within the cylinder. Three example neurons are shown. Panels b and d show the actual CSD generated by the population stimulated apically (b) or basally (d) (excitatory synaptic input). c, e: reconstruction of the CSD with spline iCSD method, boundary conditions D, h = 0.2 mm

Full size image

In the figure we present the reconstruction at one particular point in time for each type of stimulation (apical (5 B, C) or basal (5 D, E) excitation), the complete time course of the simulated activity compared with the best reconstruction is available as a video file in the supplementary material, also available at http://www.youtube.com/icsd2d. We can see that for both types of stimulation the estimation of CSD using splines effectively smoothes the real activity. This is very pronounced in the case of basal excitation (Fig. 5d), as here the actual CSD activity is more sparse. The reconstruction errors are 3% for apical stimulation (Fig. 5c vs. b), and 58% for basal stimulation (Fig. 5e vs. d). Note that error on the order of 50% indicates inadequacy of the chosen electrode grid to resolve the structures arising at the very small spatial scales in the second case. As seen in Fig. 5d and e the iCSD method is nevertheless able to reconstruct the gross features of the actual CSD distribution.

We compared the quality of reconstruction for different numbers of simulated barrel columns, different types of stimulation, different placement of the electrode grid with respect to the sources, and different assumed h (Fig. 6). Panels A and B in Fig. 6 show the placement of the electrode grid in the xz-plane. The 8 × 23 grid of electrodes extends for 1.4 mm in the x (Fig. 6a) or z (Fig. 6b) direction and for 2.2 mm in the y-direction. Panels C and D correspond to apical excitatory stimulation, while panels E and F show the results for basal excitatory stimulation. The spatial setup shown in panels A and B is the same for both stimulation schemes (panels C and E are for setup shown in panel A, panels D and F correspond to the setup depicted in panel B). The results show that for each configuration there appears to be an optimal value of h, as in the case of Gaussian sources (Fig. 6c, d; Fig. 6e, f). The large difference in CSD reconstruction errors in cases of apical (Fig. 6c, d) and basal (Fig. 6e, f) stimulation is the consequence of very different spatial structure of CSD in these two cases. For apical stimulation, we observe CSD activity of high absolute values over a large spatial extent, slowly varying in space (Fig. 5b), therefore, the spline approximation between points on the scale set by the electrode grid is reasonable. On the other hand, for basal stimulation, we observe fine-grained activity, i.e., activity at spatial scales smaller than the grid spacing (Fig. 5c), and the spline-based approximation is, as expected, less able to account for the actual activity. In the case with one central column the optimal h is 0.1–0.3 mm; this translates to an assumed thickness of CSD distribution 2 h of 0.2–0.6 mm, which is close to the diameter of the cylinder containing the somas used in the simulation.

Reconstruction fidelity of model CSD data for different placements of sources and recording grid. a, b: top-view of the grid of electrodes (x′s) and cylinders containing somas of the simulated layer-5 pyramidal cells (circles). The two setups a and b differ in the placement of the electrode grid with respect to the axis going through centers of three columns. Note that in the setup b the electrode grid is slightly off-center, which allows the impact of such placement on the reconstruction error to be estimated. c–f: reconstruction error of the whole time-course of CSD for varying number of active columns and h. c and e correspond to setup a, d and f correspond to setup b. c and d: apical excitatory stimulation, e and f: basal excitatory stimulation

Full size image

If we now consider the configuration with lateral sources, the situation depends on whether the grid is placed in the orthogonal plane (Fig. 6a) or aligns with the new columns (Fig. 6b). In the first case (Fig. 6a), the optimal h grows slightly, which is the effect of increasing size of the sources perpendicularly to the grid In the latter case (Fig. 6b), the optimal h changes very little with the addition of new columns. This is consistent with the constant extent of the sources in the direction perpendicular to the grid. The dependence of reconstruction error on h is more pronounced in the case of apical excitation, but for both types of stimulation the minimum is obtained for values close to the size of the source generators. This confirms our conclusions from the tests on Gaussian sources that the choice of h should be based on the expected extent of the activity along the axis orthogonal to the grid. Note also that for the setup shown in Fig. 6b, where the electrode grid is slightly off-center, the errors for the single central column are higher than for the setup shown in Fig. 6a. This is because the CSD distribution assumed in our iCSD scheme is symmetrical in the axis perpendicular to the grid. Indeed, the LFP generated by any CSD distribution is equal to that generated by a distribution symmetrized with respect to the plane of the electrodes. The recorded LFP will be the same whether all the CSD is ‘behind’ or ‘in front of’ the electrode grid, or shared in any fraction between the two sides. Therefore, for any data set, we can only hope to reconstruct the symmetrized part of the CSD, unless there is additional anatomical or physiological information indicating specific asymmetry, which can then be included in the method explicitly.

We also calculated the reconstruction errors for the CSD estimated using the traditional method for the situations depicted in Fig. 6. Typically, the errors of the traditional method are slightly higher than in the inverse method for the best choice of h. For example, the errors for the configuration shown in Fig. 6a and apical excitation (Fig. 6c) are 13.0%, 11.7%, and 11.2% for a single column, two columns, and three columns, respectively. One situation where the traditional method gives significantly higher errors is when there is activity at the boundary of the grid (Fig. 6b, lowermost of the three columns). For three active columns (circles in Fig. 6d) the error of the reconstruction obtained with the traditional method is as high as 19.6%. These findings are compatible with the results for Gaussian sources, where there was always significant activity at the boundary. Thus, in this case we see that the main benefit of using iCSD instead of traditional approach is significant gain of precision at the boundaries.

Selection of h

The parameter h is the main a priori unknown parameter specifying the 2D iCSD method, and the choice of h is, therefore, an important issue in practical applications. Optimally, h should be chosen based on the expected extent of neural activity, known anatomy, or the results of forward modeling studies according to the expected size of a typical active population.

As seen in Figs. 4 and 6 the assumed activity depth influences the estimation error. However, in these examples the estimation errors seem to converge towards a constant value above h = 1 mm, which corresponds to an assumed activity depth 2 h of 2 mm. This is particularly visible for the case of basal stimulation. This reflects that the 2D iCSD estimator itself approaches a fixed ‘large-h’ form when h becomes large enough. It is shown in the Appendix that the typical large-h transition value depends on the geometry of the multielectrode, in particular the distance between recording grid points. This large-h transition value is potentially of practical importance: in a system where neural activity is more widespread than this value in the direction orthogonal to the 2D grid, the ‘large-h’ 2D iCSD method can be used and the uncertainty due to lack of knowledge of the true effective value of h will be minimal. The large-h limit of the 2D iCSD method is studied in the Appendix.

Analysis of Experimental Data

To test the utility of the proposed method in analysis of real data we applied it to the set of simultaneous recordings obtained in the rat subiculum using an 8 by 8 multielectrode, as outlined in the “Materials and Methods” section. We assumed h = 0.5 mm and Gaussian profile of H(z). Figure 7a–h compares the reconstructed current-sources using the iCSD 2D with D boundary conditions (marked 2D iCSD) with the interpolated averaged potentials in time from the stimulation. The interpolated data were superimposed on top of anatomical borders (Fig. 7i) according to the histology (Fig. 2).

Laminar profile of network activity in rat subiculum. Panels a–h: interpolated data from 8 × 8 MEA recording grid (1.6 × 1.6 mm) in dorsal subiculum (iCSD at left, potential at right). Insets in the left panels show post-stimulus time for each frame. Black bar in a shows 1 mm scale. Anatomical borders are indicated in i with matched average potentials (red box) shown in j (scale bar: 1 mV/15 ms; w.m. white matter; data shown 15 ms post-stimulus). In a–d single-pulse alvear stimulation (* in i) produces a peak negative voltage (red) across the subicular cell layer (SUBp; panel a voltage) with bordering negativity (molecular layer (SUBm) and basal dendrites) that moves over time to the SUBp/SUBm border. There are also smaller positive–negative responses in presubiculum (PrS) and in the border region between PrS and Sub. This pattern is also reflected in the 2D iCSD patterns as sinks (red) and sources (blue). In c–h the main SUB pattern appears to ‘split’ around a central quiescent zone, fade (d–e) and then strengthen again (f–g) before fading slowly

Full size image

A paper centred on a large 2D data set will be submitted elsewhere, hence, we provide here a basic interpretation of the present responses. Alvear activation evokes an antidromic subicular population spike (SUBp, peaks at 1.7 ms) across the full extent of the subicular cell layer (Menendez de la Prida 2003). The main SUBp sink back-propagates across the cell layer towards the molecular layer (1.7–2.1 ms sequence), then ‘splits’ and finally fades away slowly (from 2.5–3.8 ms; along SUBp/SUBm). We interpret this ‘split’ as feed-forward activation from antidromically-activated subicular pyramids that project laterally (but perhaps not to the middle of SUB) with the likely recruitment of local inhibitory cells (Harris et al. 2001). A strengthening of the main split pattern and a reversal of the activity at SUB/PrS border (and PrS) likely reflects feed-back inhibition of SUBp, producing an ‘inhibitory’ source and ‘passive’ sink. After back-propagation the region lying along the ‘split’ has very little synaptic activity, suggesting that the proposed feedback mechanisms on either side of this region are not active here.

Figure 8 shows the CSD reconstructed using different h in the iCSD method (0.05 mm, 0.5 mm, and 3.2 mm respectively in Fig. 8a–c) and using the traditional method (Fig. 8d). To obtain the values at the boundary layer in the traditional CSD reconstruction we used the 2D analog of the Vaknin procedure. The shapes recovered using these methods are very similar. If we use h = 0.5 mm as reference the errors e₂ of the four plots are 3.3%, 0%, 0.8%, and 2.4%. The largest difference is in the amplitude of the CSD, which is 263, 133, 126, and 77 for Fig. 8a–d, respectively (arbitrary units). This example shows that in some cases traditional CSD and iCSD may lead to equivalent results. We expect this to happen especially when there is little activity close to the edges, or when the inter-electrode distance is large compared to the spatial scale of activity (here we have sinks and sources separated by just 0.2 mm, that is, one inter-contact distance). The latter claim is supported by tests similar to the ones shown in Fig. 3, but with longer and more narrow sources, resembling the experimental activity shown in Fig. 8 (test data not shown). In this case we also obtained reconstruction errors on the order of 9% for the spline iCSD variants (B and D) only mildly better than for the traditional method (e₂ =10%) with the difference between the iCSD and traditional method around 6%.

Reconstruction of CSD from experimental data (t = 1.7 ms after stimulation). Three reconstructions with the iCSD method (spline, Gaussian z-profile) with different assumed h (a–c) and the traditional method with Vaknin procedure (d) give similar shapes of the CSD. The predicted amplitudes are different (see text), each plot is rescaled using the amplitude of the CSD to facilitate comparison of the shape

Full size image

Among the many techniques available for the study of information processing in the brain, electrophysiology stands out when it comes to temporal resolution of the signal. Some obstacles in its use are (a) the technical problems of simultaneous signal recording at multiple sites and (b) the large spatial range of electric fields measurable with these electrodes. It is now feasible to implant around a hundred electrodes within a relatively small brain volume to record signals simultaneously, allowing for remarkable spatial and temporal resolution. The amount of information coming from such experiments, matched with adequate analytical tools, are used for precise identification of single units (Buzsáki 2004), brain machine interfacing (Nicolelis 2001) and in studies of LFPs for a precise description of population activity of neural structures and spatial localization of synaptic connections.

Although a remarkable spatial resolution in recorded potentials is achievable with modern multielectrodes, this does not automatically afford a correspondingly high spatial resolution in the estimated neural activity due to the inherent long-range properties of LFPs. The Inverse Current Source Density method proposed in Pettersen et al. (2006) for laminar 1D recordings, developed by Łęski et al. (2007) for 3D data and here for 2D grids, allows more robust reconstruction of the sources and sinks generating the measured LFPs than previous methods. One specific advantage, particularly important in 2D (and in 3D), is that the iCSD method appears to recover the CSD close to the boundary of the electrode-contact grid more accurately than the traditional CSD method. With the increasing availability of multi-shank, multi-contact electrodes and microelectrode arrays there is now a growing need for data analysis methods to match the sophistication of the sensor hardware. In this regard, we believe that the 2D variant proposed here, with its freely available implementation, will find immediate use in such electrophysiological studies.

Previous Studies of Two-dimensional CSD

The links between CSD and LFPs were discussed in full three-dimensional generality previously by Nicholson and co-workers (Nicholson 1973; Nicholson and Freeman 1975; Nicholson and Llinás 1975). Their approach has been used most frequently, however, for analyses of laminar recordings in one dimension (e.g. Haberly and Shepherd 1973; Mitzdorf 1985; Schroeder et al. 1992; Ylinen et al. 1995; Lakatos et al. 2005; Lipton et al. 2006; Rajkai et al. 2008; de Solages et al. 2008). There are several studies that involved the estimation of CSD in two dimensions (Novak and Wheeler 1989; Shimono et al. 2000; Shimono et al. 2002; Lin et al. 2002; Phongphanphanee et al. 2008). However, all of these applied the traditional approach through estimation of the second numerical derivative, and most applied the specific technique (smoothing followed by differentiation) proposed by Shimono et al. (2000). The problem with boundary values was observed by Novak and Wheeler (1989), who refrained from the analysis of such edge data. On the other hand, in the papers by Shimono et al. (2000, 2002), Lin et al. (2002), and Phongphanphanee et al. (2008), there is no explicit discussion of boundary treatment.

A notable difference between the standard (numerical second derivative) and the inverse CSD approaches is that in iCSD we can specify the parameter h, which translates to the assumed thickness of the region with sources. The standard CSD method has an implicit large-h assumption: it implies that no field escapes in the third dimension (the double derivative in the Poisson equation is assumed to be zero in the third dimension). Although the 2D iCSD methods require an explicit assumption of h which is not always easy to justify, it is better to have this assumption explicit than forcing a large-h assumption, which is done implicitly in the standard CSD method.

We should stress that we do not want to imply here that the traditional CSD approach is invalid or should be abandoned, as it is simple to apply and in many cases may lead to reasonable estimates of CSD (see our experimental example). However, from the general viewpoint presented here, we believe 2D iCSD is better grounded. Indeed, we are convinced that the framework of 2D iCSD offers a viable and meaningful alternative, which can be further extended to incorporate additional physiological knowledge as it becomes available.

Main Results

In this article we have developed several variants of the iCSD 2D. The basic difference between them comes from the assumed structure of the CSD: either constant within boxes or interpolated (linear or spline) between the points of an electrode grid. To accommodate sources located outside the probed region one can expand the grid spanning the CSD (Fig. 1). The CSD at these extra points can be set to duplicate the values from neighboring points or alternatively set to zero. Our tests on several surrogate data sets, including Gaussian sources of different spatial structure (Fig. 3) and extracellular potentials generated from a model population of 1000 pyramidal cells (Fig. 5), indicate that the optimal approach is the spline interpolated iCSD on an extended grid with duplicated boundaries.

The main free parameter is the assumed width of the sources perpendicular to the multielectrode grid, h. Tests on model data indicate that the performance is best when h is in the order of the actual size of the current sources. Thus, it is beneficial to estimate h from anatomical studies of the target region. However, reasonable deviations from the best choice of h deteriorate the profile of the reconstructed sources only very slightly. We also observed stabilization of the error values for the normalized CSD with growing h, which corresponds to assuming infinite extent of sources in the direction perpendicular to the grid. For the cortical pyramid model (Fig. 5) such stability was reached at about h = 1 mm (i.e., assumed width of 2 mm) for our multielectrode with a grid-spacing of 0.2 mm. This error stabilization reflects convergence of the 2D iCSD method in the large-h limit. This convergence is independent of the neural sources and only depends on the geometry of the grid-points of the multielectrode. Thus, for a given spatial extension of neural activity in the direction orthogonal to the 2D grid, the large-h limit can, in principle, be reached by reducing the distance between grid points in the 2D multielectrode. In many real biological applications, one may already be at this large-h limit. In such situations, one may simply choose a sufficiently large h when constructing the iCSD estimator. In case of very limited information on the possible extent and distribution of the sources, we recommend experimenting with several different values of h and looking for features appearing stably across the different iCSD reconstructions.

To test our method in practice we applied it to a set of recordings from rat subiculum (Figs. 2, 7). Compared with the interpolated potentials, iCSD plots provide more precise localization of neural activity in subiculum following alvear stimulation. The observed activity pattern points to the heterogeneous structure of connections in the subiculum and supports the columnar anatomical model proposed by Harris et al. (2001).

Graphical User Interface Tool for CSD Analysis

As supplementary material we provide a MATLAB toolbox containing the scripts used in our analysis together with a simple GUI (Fig. 9), allowing easy calculation of different variants of iCSD from voltage data recorded on 2D regular grids. The viewer is bundled with three of the datasets used in the paper: two model datasets (apical and basal excitatory input to a population of pyramidal cells), and one experimental (data recorded in the rat subiculum). It is also possible to import and analyze user data provided as a workspace variable or an array stored in a MAT file. The viewer makes it straightforward to switch between traditional CSD and different variants of inverse CSD, to adjust boundary conditions, and to test different values of the h parameter in case of iCSD. The resulting figures can be exported as PNG files.

Graphical tool for CSD analysis. The main window shows details of the dataset, the CSD method, two panels with voltage and CSD data, and controls for browsing through the dataset

Full size image

The software and the data can be used freely in research provided this paper is cited in any material using the results of their application. The toolbox is provided under GNU General Public License version 3 or later and is available from the INCF software repository, http://software.incf.org/, where a toolbox for iCSD analysis of 1D linear recordings, CSDplotter (Pettersen et al. 2006), can also be found.

Future Developments

Inverse Current Source Density developed in one (Pettersen et al. 2006), two (here), and three dimensions (Łęski et al. 2007) is a flexible framework that allows us to incorporate different assumptions about the distribution of sources or the geometry of probed structures. So far, however, it has been developed only on regular grids with the assumption of uniform and homogeneous tissue conductivity. Conductivity in any structure is neither completely homogeneous nor uniform (see for example recent measurements reported in Logothetis et al. 2007 or a study of conductivity in the rat barrel cortex, Goto et al. 2010), therefore, one way of developing the method would be to consider inhomogeneous or anisotropic media. Equation 1a then becomes \( \nabla (\sigma \nabla \,\Phi ) = - C \), where σ is a tensor which can take different values at different spatial positions. The solution to this equation would then be used to construct the forward model matrix F. The generic situation can only be treated by solving the equation numerically, but in several important cases the solution can be given explicitly. For instance, if σ is constant in space then one can always rotate the Cartesian coordinate system so that σ becomes a diagonal matrix: σ = diag(σ _x ,σ _y ,σ _z ). Then the solution Φ for the unit source located at the origin is

Also some important cases of inhomogeneity can be solved analytically. One such case is conductivity that is constant within planar layers but has different value in each layer (Gold et al. 2006; López-Aguado et al. 2001), then the method of images can be used to express the solution Φ as a series.

Another way of generalizing the method would be to consider arbitrary distributions of recording points. This seems to be of immediate use, as it is becoming increasingly easy to insert large numbers of independent electrodes probing large volumes of the brain, but not necessarily on regular grids. This will be an object of our further study. Until now we have assumed that the grid spanning the CSD is the same as the electrode grid. This need not necessarily be so and one could consider electrode arrays with differing geometries, e.g., electrode contacts not forming an exact straight line (Barthó et al. 2004).

An important topic we have only briefly covered in the present work is the stability of the results with respect to noise (see the Appendix). The current-source density analysis requires us to calculate—explicitly in ‘classical’ CSD and implicitly in iCSD—the spatial derivative of the potential. This leads to amplification of noise present in the signals, see detailed analysis in (Freeman 1980) in the context of ECoG. The thorough study of the influence of different sources of noise (measurement noise, uncertainty in the electrodes’ positions, etc.) is beyond the scope of this paper. However, the analysis of the condition number of the matrix F for different iCSD variants (step, linear, spline) suggests that there is a tradeoff between how well a given method resolves the shape of the sources and how stable it is with respect to noise.

It would be useful to establish general range of applicability of the proposed methods. Unfortunately, to make general statements regarding the optimality or stability of the iCSD method one would have to know the class of possible distributions of fields in the brain, which is not available. This is why we tested the quality of different methods used for CSD estimation on both artificial test data and simulated data from synaptically activated populations of biophysically detailed neuron models with reconstructed morphologies. As illustrated by Fig. 5b and d the ground truth (“true CSD”) is quite “noisy” in the sense that it varies from pixel to pixel. Nevertheless, for all situations encountered we found both the spline-iCSD and linear-iCSD methods to be stable in the sense that it gave results in good agreement with the know ground truth (see, e.g., Fig. 5c and e). In addition, we found no signs of instability when applying our method to the experimental test data. We only encountered instability in predictions in test cases where the current sources were positioned far away from the recording grid, which resemble situation in the EEG or ECoG source estimation (not shown). Such situations might require use of more involved approaches but they are not our goal: we are interested in estimating sources close the recording grid within the same plane. It would be interesting to see how more involved approaches from EEG or ECoG source localization would perform in this context (He and Lian 2005; Nunez and Srinivasan 2006). Also, it would be interesting what one could achieve in the context of CSD estimation using Bayesian techniques, which allow to incorporate the prior knowledge (e.g. anatomical) in a natural way (Baillet and Garnero 1997; Schmidt et al. 1999).

The toolbox for CSD analysis in 2D is publicly available from the INCF software repository, http://software.incf.org/ (project name iCSD 2D). The toolbox consists of MATLAB scripts, a GUI viewer, and three example data sets. The software and the data can be used freely in research provided this paper is cited in any material using the results of their application. The supplementary video is available at http://www.youtube.com/icsd2d and can be used freely in non-commercial education. For other usage, contact the authors.