Revealing the distribution of transmembrane currents along the dendritic tree of a neuron from extracellular recordings
Abstract
Revealing the current source distribution along the neuronal membrane is a key step on the way to understanding neural computations; however, the experimental and theoretical tools to achieve sufficient spatiotemporal resolution for the estimation remain to be established. Here, we address this problem using extracellularly recorded potentials with arbitrarily distributed electrodes for a neuron of known morphology. We use simulations of models with varying complexity to validate the proposed method and to give recommendations for experimental applications. The method is applied to in vitro data from rat hippocampus.
https://doi.org/10.7554/eLife.29384.001Introduction
A variety of methods are used to investigate the electrophysiological properties of neurons. To date, patchclamp (Neher and Sakmann, 1976) is the most commonly used technique to monitor neuronal membrane potential. Despite its unquestionable utility, it remains challenging to monitor the activity of a cell at more than one or two sites. Extracellular recordings, on the other hand, deliver a more global picture of neural activity (Buzsáki et al., 2012; Einevoll et al., 2013a). With modern multielectrodes and microelectrode arrays, it is now possible to record neuronal activity from many thousands of channels (Buzsáki, 2004; Berdondini et al., 2005; Obien et al., 2014). However, this technique does not permit direct recording of membrane potentials but instead spiking activity [which may be of individual cells (singleunit activity, SUA) or multiple cells (multiunit activity or MUA, which is the mean firing rate of cell populations)] and components of postsynaptic activity visible at low frequencies (<300 Hz, socalled local field potential, LFP); see (Buzsáki et al., 2012; Einevoll et al., 2013a; Głąbska et al., 2017) for discussion.
So far, the main advantages of highdensity array recordings have been improved resolution of spike detection (Rey et al., 2015), as more cells can be identified in a single recording, improved stimulation precision (Hottowy et al., 2012; Chichilnisky, 2001), of particular importance for retinal neuroprosthetics, and new features observed in the profiles of slow fields (Ferrea et al., 2012). Recently, highdensity probes have also been used in studies of axon tracking (Bakkum et al., 2013; Lewandowska et al., 2016) and multisynaptic integration (Jäckel et al., 2017).
Multielectrode recordings have been traditionally used for improved spike sorting (Buzsáki, 2004; Berdondini et al., 2005; Obien et al., 2014; Muthmann et al., 2015) and for the reconstruction of current source densities (CSD) behind the recorded LFP (Pitts, 1952; Mitzdorf, 1985; Wójcik, 2015), although more specific methods have sinse been devised (Einevoll et al., 2013b; Głąbska et al., 2014; Głąbska et al., 2016). Several attempts have been made, using different approaches, to localize cells from multielectrode recordings. For example, accounting for the properties of electric field propagation in the tissue (Muthmann et al., 2015), that form the basis of CSD methods (Somogyvári et al., 2005, Somogyvári et al., 2012Somogyvári et al., 2012), or other triangulation approaches (Mechler et al., 2011; Mechler and Victor, 2012). We are not aware of any prior attempts, however, to reconstruct the CSD of individual neurons using their available morphologies, which we propose here.
This method assumes we have a set of extracellular recordings, coming from a specific neuron, whose morphology and location with respect to the electrode is known, collected with multiple contacts. This could be realized experimentally by patching a neuron close to the multielectrode and driving it through an intracellular injection or monitoring its activity to determine the contribution of this specific cell to the extracellular field. Computing spiketriggered average of the potential, which we do in our proofofconcept experiment, or driving the neuron with sinusoidal current and averaging the extracellular potential over periods of the driving current, are ways in which this could be achieved. When the recordings are complete, we inject dye into the cell and reconstruct its morphology. Thus, we obtain a set of synchronous multichannel extracellular recordings reflecting the activity of a single neuron whose morphology is also known, as well as the position of the neuron relative to the electrode contacts. Here, we show how to use this approach to infer the distribution of current sources based on cell morphology as they change in time. The data necessary to apply presented method have been available for some time (Henze et al., 2000; Gold et al., 2006), although recently have became much more comprehensive (Jäckel et al., 2017). While we believe the estimation of transmembrane currents along the cell morphology using this type of data has not been reported previously, similar questions have been posed by (Gold et al., 2006), who attempted to identify the biophysical properties of a neuron membrane based on the extracellular signature of the action potential. A similar strategy was used by Frey et al., 2009) in their studies of extracellular action potential shape observed with highdefinition multielectrode arrays.
The singlecell kernel Current Source Density method (skCSD) we introduce here is an application of the framework of the kernel Current Source Density method (Potworowski et al., 2012) to the data coming from a single cell. This is done by restricting current sources to cell morphology. This can be done efficiently for arbitrarily complex morphologies and arbitrary electrode configurations.
The importance of this work is that for the first time we show here how a collection of extracellular recordings in combination with cell morphology can be used to estimate how the current sources located on a studied cell contribute to the recorded field potential. Since it is feasible to acquire the relevant data, we believe that the method proposed here may be used to constrain the biophysical properties of the neuron membrane and facilitate consistency of the reconstructed morphology. Further, this method can help guide new discoveries by providing a more global picture of the current distribution based on neuronal morphology, leading to a coherent spatiotemporal view of synaptic drive and return currents of the observed neuron.
In the following Results section, we start with a highlevel overview of the skCSD method. We explain how it is applied and why it works. Then, we validate this method on several ground truth datasets obtained in simulations and apply it to data from a proofofconcept experiment. In the Discussion, we summarize our main findings and discuss the practical aspects and feasibility of experimental acquisition of the required data. Finally, in the Materials and methods section, we present the skCSD method in detail.
Results
The main result of this work is the introduction of the skCSD method, so we start here with a highlevel overview. The technical details are deferred to the Materials and methods section. Next, we study the properties of the skCSD reconstruction for three representative morphologies of increasing complexity and for different setups.
First, for a ballandstick neuron, we study the general quality of reconstruction of fine detail by considering CSD distributions in the form of standing waves of increasing spatial frequency which form the Fourier basis of any possible CSD profile. It is unlikely that standing waves would be naturally observed in a cell, therefore to better understand how the results for the Fourier space representation relate to a specific distribution which might arise in a physiological situation, we also consider reconstruction of sources for random synaptic activation of the ballandstick cell.
Secondly, we consider a Yshaped neuron with a single branching point, we check if skCSD can differentiate between synaptic activations located on the different branches close to the branching point. We also investigate the effects of random distribution of contacts on skCSD reconstruction. Finally, we investigate the possibility of skCSD reconstruction on a realistic model of a ganglion cell placed on a microelectrode array (MEA) as well as the sensitivity of the method to noise.
After establishing and validating skCSD on these fully controlled model datasets, to demonstrate neurophysiological viability, the CSD distribution was reconstructed for a pyramidal cell using the experimental spiketriggered averages of the recorded potentials.
Singlecell kernel current source density method, a highlevel overview
The goal of this section is to provide an overview of the proposed method for nonspecialists, with limited mathematical terminology and notation. A more formal discussion is provided in the Materials and methods section.
Assume we study a neuron, we have its morphology, we know how it is located with respect to a multielectrode used for extracellular recording, we also have a set of simultaneous recordings of extracellular potential generated by this cell collected with this multielectrode. In principle, the number and placement of the electrodes can be arbitrary. Also, the potential may be filtered, or we may consider the full spectrum of the signal, depending on whether we wish to focus more on synaptic contributions or consider the extracellular signatures of spiking. For now, we shall ignore the challenge of separating the part of the signal contributed by the studied neuron from background extracellular signals generated by nearby cells; we shall return to this issue in the Discussion.
We wish to reconstruct the distribution of current sources which generated the measured potentials. By assumption, we know the potential that comes from the studied cell, we wish to restrict the sources to lie on its morphology. To do this, we first represent the morphology by a closed line which we call the morphology loop. To construct it consider a onedimensional abstraction of the cell, where we ignore the thickness of the dendrites. Alternatively, you may imagine the graph constructed from the lines passing centrally through all the dendrites. Then, starting for example at the soma, we draw a line along this graph passing all sections along the dendrite, eventually reaching the starting point. The morphology loop is shown as the red line in Figure 1.
If we spread the morphology loop, we obtain a circle, which means all point of dendritic morphology have been mapped to a circle, and the opposite, any point on the circle has been mapped to the morphology. The mapping from the cell to the circle is not unique: we pass by most points on the dendrite twice, with the exception of the tips, which are visited once, and the branching points, which may be visited more than twice in our graph representation. So in most cases, a given point from the morphology corresponds to multiple points from the circle. The other mapping is unique: every point in the circle is mapped onto exactly one point on the morphology without skipping any. One way to think about the morphology loop is as a rubber band tightly wrapping around the neuron’s morphology.
We now want to consider the distribution of current sources on the morphology. We found it convenient, technically, to start with a distribution of sources along the morphology loop. Then, we wrap this distribution around the cell together with the loop. We do this by construction. We cover the morphology loop with a large number of identical but translated functions which we call the CSD basis functions denoted by ${\stackrel{~}{b}}_{j}(\mathbf{\mathbf{x}})$. There is a large flexibility here, but in practice we use Gaussians, so ${\stackrel{~}{b}}_{j}(\mathbf{\mathbf{x}})\propto \mathrm{exp}({(\mathbf{\mathbf{x}}{\mathbf{\mathbf{x}}}_{j})}^{2}/2{R}^{2})$. The number of basis functions we use, $M$, and the width of the basis function, $R$, are parameters of the method whose effect on results we discuss below.
We place these basis functions so that they uniformly cover the morphology loop. We require their centers to be uniformly spaced. When we wrap these functions around the morphology, passing a given dendrite twice will introduce more overlap. We only require that no two functions overlay: we want them to be independent which means each two must differ, here, in practice, be shifted with respect to each other. Once the CSD basis functions are in place, we compute the potential they generate in the whole space. The distribution of the potential in space coming from a single CSD basis function is called a basis function in the potential space, or potential basis function for short, and is denoted by ${b}_{j}(\mathbf{\mathbf{x}})$ (we drop the tilde).
We are now ready to start the estimation of the distribution of the current sources on the cell from the recording. This is a static procedure, in the sense that the estimated CSD at any moment in time depends only on the present value of the measurements. We are looking for the distribution of current sources in the form of a linear combination of the basis sources which we placed. This means a weighted sum of each CSD basis source
Then, the potential they generate, is the linear combination of the respective potential basis sources with the same weights,
Conceptually, we want to match the described potential with actual measurements and from this infer the weights. In practice, however, we cannot do it directly, because typically there will be many more basis functions, and therefore weights to be estimated, than the available measurements. To understand why this is a problem consider the simplest case, where you have two numbers, $x$ and $y$, to estimate from a single measurement which gives 1, and the physics of the problem gives equal contribution to the measurement, so that we must solve $x+y=1$. It is easy to see that without further constraints this has an infinite number of equally good solutions.
One way to solve this problem, which was our inspiration, was proposed by Vapnik (1998) (Appendix to Chapter 6: Estimating functions on the basis of indirect measurements), who effectively considered the problem of estimating one quantity (here: CSD) from measurements of its function (here: potential). Here, we combine it with the kernel trick (Schölkopf and Smola, 2002), which allows us to make indirect estimation in the highdimensional space of basis functions through computations in the space of measurements. We construct a kernel function which is a sum of products of the potential basis functions with themselves
This function, which takes two spatial arguments, can be understood as a similarity measure between the potentials at the two points. It is easy to see, that any model of potential we can construct from our potential basis functions can also be written as a linear combination of these kernel functions with one of the variables fixed, if we use sufficiently many kernels (large $L$)
The reason is that $K({\mathbf{\mathbf{x}}}_{l},\mathbf{\mathbf{x}})$ is a linear combination of all the basis functions spanning the potential space: $K({\mathbf{\mathbf{x}}}_{l},\mathbf{\mathbf{x}})={\sum}_{j=1}^{M}{b}_{j}({\mathbf{\mathbf{x}}}_{l}){b}_{j}(\mathbf{\mathbf{x}})$. Thus, if we take as many kernels as we have basis sources ($L=M$) and equate Equation 2 with Equation 1, we have to solve $M=L$ equations of the form
for ${\beta}_{l}$, or in other words, we have to find such points ${\mathbf{\mathbf{x}}}_{l}$ for which the above equation is solvable. This can be done, this is another way of saying that the functions ${b}_{j}(\mathbf{\mathbf{x}})$ form a basis. It is thus fair to say that the kernel functions $K({\mathbf{\mathbf{x}}}_{l},\mathbf{\mathbf{x}})$ and ${b}_{j}(\mathbf{\mathbf{x}})$ are equivalent basis. At this stage, it is not clear if one basis should be better than the other.
The rationale for using the kernels is provided by the Representer Theorem (Schölkopf and Smola, 2002), which shows that in the form of Equation 2 we can minimize the prediction error (sum of the squared differences between the predictions of our model and actual measurements) uniquely. Moreover, the solution obtained has as many parameters as there are measurements, and we take ${\mathbf{\mathbf{x}}}_{l}$ to be the measurement points. This is the advantage of the kernel approach over direct estimation of underlying model: here the number of parameters to be calculated is the same as the number of measurements $N$, much less than the original number of basis sources $M$, so it can be done.
Since we can expect our measurements to be noisy, a perfect fit will typically misrepresent the true potential, this is called overfitting. To counter this effect, we add a socalled regularization term to our error function to be minimized, which is the sum of squared parameters,
Thus, we want to simultaneously minimize the difference between our prediction and actual measurements while moderating the fluctuations of the interpolated potential in space. This only slightly changes the computation while significantly improving the stability of the solution. The result is
where $\mathbf{\mathbf{V}}$ is the vector of the measurements ${V}_{k}$, I the identity matrix, ${\mathbf{\mathbf{K}}}_{jk}=K({\mathbf{\mathbf{x}}}_{j},{\mathbf{\mathbf{x}}}_{k})$, and $\lambda $ the regularization parameter, which needs to be set.
We still need to obtain the CSD profile from the interpolated potential. In fact, we can do this without resorting directly to the basis functions. Replacing one of the potential basis functions in the definition of our kernel with the corresponding CSD basis function, we obtain what we call a crosskernel function,
Defining
the estimated CSD is
The skCSD is a modelbased analysis method since the specific model of CSD distribution we use, collecting of the CSD basis functions along with the model relating these functions to the measurement of potential, influences the result. This is advantageous, since all the assumptions are explicit and the user can see how they affect the result. All the estimation methods of any quantity contain assumptions, which in many cases are implicit and thus it is difficult to analyze how they affect the estimation. With all the parameters explicit we can study how their specific values affect the quality of solution. In particular, we wish to select parameters leading to the optimal solution. We do this using crossvalidation which we shall now explain.
We select a set of parameters: $R$, $N$, $\lambda $, which fixes the model. Then, going contact by contact, we ignore the signal recorded at that particular site and build models from the remaining signals. This model predicts the potential at the ignored contact. The difference between the prediction and the actual measurement is a measure of prediction quality for a given set of parameters. We then add squares of the differences between the actual measurements at every electrode and predictions from the respective models built from all signals except the reference. Scanning through a range of parameters we look for a minimum prediction error. We use the parameters minimizing the prediction error in the subsequent analysis.
Ballandstick neuron
Here, we consider the simplest neuron morphology, the socalled ballandstick model, which stands for the soma and a single dendrite. A virtual linear electrode was placed in parallel to the model cell 50 μm away, the electrodes were distributed evenly along the electrode extending for 600 μm.
Increasing the density and number of electrodes improves spatial resolution of the method
To study the spatial resolution of the skCSD method, we consider the ground truth membrane current source density distributions in the form of waves with increasing spatial frequencies
where $A=0.15nA/\mu m$ is the amplitude, $f\in \{0.5,1,1.5,\mathrm{\dots},12.5\}$ is the spatial frequency, $x$ is the position along the cell, $L$ is the length of the cell. Then, we compute the generated extracellular potential at the electrode locations. The laminar shank consisting of 8, 16 and 128 electrodes was placed 50 μm from the cell in parallel to the dendrite. Finite sampling of the extracellular space sets a limit to the spatial resolution of this method. Increasing the density of electrodes within the studied region leads to higher spatial precision. As shown in Figure 2, with 128 electrodes it is possible to reconstruct higher frequency distributions as compared to eight electrodes. This is reminiscent of the sampling theorem (Oppenheim et al., 1997), except here we measure the potential and reconstruct current sources, while in the sampling theorem we consider reconstruction of a continuous signal from discrete samples. What we observe is quite intuitive and typically observed in different discrete inverse methods (Hansen, 2010). Note that once we move to complex morphologies and random rather than regular electrode placement, the intuition we build here, that denser probing gives better spatial resolution, would hold true, even if the relation to the sampling theorem would be less apparent.
Reconstruction of random synaptic activation
Using the ballandstick neuron, we now place 100 synapses along the dendrite and stimulate them randomly in time. We simulate 70 ms of recordings from this synaptically activated cell. The stimulation is sufficiently strong to evoke spiking, see Materials and methods for details. The spiking is indicated by strong red spots in the lowest first two segments in Figure 3, which correspond to the soma. As can be seen, the reconstructed CSD distribution reflects the groundtruth, and the precision of reconstruction improves with an increasing number of contacts, which is reflected in the reduction of crossvalidation error. Notice how the reconstructed synaptic activity gets more precise with increased density of probing the potential. In particular observe how the width of the recovered synaptic activations and the somatic activations shrink with an increasing number of electrodes, which clearly shows improved resolution. This is consistent with our observations for the Fourier mode CSD profiles above. Not much change is seen in time, which is a consequence of the fact that skCSD, like all the CSD estimation methods, acts locally in time. That is, for every moment in time, the collection of potentials at this time, is analyzed. There is no direct relation to the past or future of the measured signals.
Simple branching morphology
Let us now study the effect of branching and breaking of rotational symmetry of the cell using the skCSD method. We consider here a simple $Y$shaped model neuron with one branching point (Figure 4B). We place two synapses, one on each branch (at segments 33 and 62, close to the branching point, see Figure 4D and Figure 5C). We consider both simultaneous and independent activation of these synapses, specifically, the first synapse was activated at 5, 45, 60 ms of the 70 ms long simulation, while the other was stimulated at 5, 25, 60 ms from the stimulation onset. Our goals were to determine if it was possible to separate the synaptic inputs located on two different branches, what happens at the branching point, how the arrangement of the electrodescell setup influences the reconstruction. We also wanted to determine if this method provides more detail about the current distribution on the cell than what is accessible from the interpolated potential and the CSD reconstructed with kCSD under the assumption of a smooth distribution of sources in space, which is the natural approach to try (Frey et al., 2009).
Differentiation of synaptic inputs located on different branches
To investigate the differentiation power of the proposed approach, we consider two placements of the cell with respect to the electrode grid. Plane $xy$, in which the cell is placed in parallel to the plane of electrodes 50 μm above (Figure 5A), and plane $xz$, where the cell is perpendicular to the grid, with the grid 50 $\mu m$ away from the dendritic shaft stemming from the soma, (Figure 5B,C). In Figure 5, each panel (A–C) shows the splineinterpolated extracellular potential (V), followed by standard kCSD reconstruction, both at the plane of the 4 × 16 electrode grid used for simulated measurement. Then, the ground truth and skCSD reconstruction are shown in the branching morphology representation in the plane containing the cell morphology. Each figure is superimposed with the morphology of the cell. The dark gray shapes are guides for the eye and are sums of circles placed along the morphology with a radius proportional to the amplitude of the sources located at the center of the circle. Panel A shows results for a synaptic input depolarizing one branch. Panel B shows the same current distribution as in the previous setup, but the cell is rotated by 90 degrees with respect to the grid. In panel C synaptic input is added to the other branch. Observe that in all three cases the interpolated potential and the standard CSD reconstruction, which can be drawn only in the plane of the electrode grid, do not differ significantly, hence they cannot distinguish between these three situations. On the other hand, the skCSD method correctly identified the synaptic inputs in all three cases.
Note that without the method proposed here, the most natural approach to analyze current sources is through use of the regular, population CSD. This approach was used, for example, to investigate the changing distribution of current sources during action potential generation using data from a highdefinition MEA (Frey et al., 2009). What we show in figures here and below, is that while CSD (kCSD) and skCSD are consistent, using the additional information about morphology renders significantly more detail about the activity studied.
The effect of electrode placement on skCSD reconstruction for Yshaped cell
In Figure 6, we show how the number and specific distribution of the electrodes affect the quality of the reconstruction in the case of simultaneous stimulation. Panel 5. A shows the ground truth data, that is the actual distribution of the transmembrane current sources, along the morphology. To visualize it simply, we used the interval representation, the soma is shown first, followed by one branch, followed by the other. Figure 6B shows the reconstruction results for regularly arranged 8 (4 × 2), 16 (4 × 4), 32 (4 × 8), and 64 (4 × 16) electrodes. In Figure 6C, we show reconstructions for five different random placements of the same number of electrodes as for the regular case. As expected, the skCSD method is able to recover the synaptic activations and the reconstruction resolution increases with the number of electrodes. Note that in certain cases, the random distribution is more efficient than the regular grid, which is probably due to more fortunate samplings of the area covered by the morphology.
The effect of basis on skCSD reconstruction for the Yshaped cell
To investigate reconstruction quality in the parameter space set by the number of basis functions ($M$), basis function width ($R$) and regularization parameter ($\lambda $), we used the simulation setup for the Yshaped morphology with 4 × 4 electrodes. Figure 7 shows the L1 reconstruction error for $M=32,128,512,1024$, $R=8,16,32,64,128$, and $\lambda ={10}^{5},{10}^{4},{10}^{3},{10}^{2},{10}^{1}$. As we can see, for the smallest basis ($M=32$) and small $\lambda $, the minimum error is obtained for wide basis sources, so that the basis functions have substantial overlap. This is necessary for the method to be able to reconstruct the family of test sources we considered. As the basis size increases, the reconstruction improves overall with minimum error obtained for narrow basis sources and small $\lambda $. The fact that we have two comparable minima for $M=32$, for small and large $\lambda $ (top and bottom right of the plot for $M=32$), means that the error we obtain by emphasizing the measurements (small $\lambda $) is comparable to the error we obtain by emphasizing the regularization term, which prevents overfitting (large $\lambda $) and in effect, reflects our doubt about precision of measurement. We interpret it here as the effect of insufficient basis size. This effect disappears with increasing basis size when a unique minimum appears for moderate values of $\lambda $ and for narrow basis sources, which can best resolve small details of the CSD to be reconstructed.
Reconstruction of current distribution on complex morphology
In this section, we consider the performance of the skCSD method in the case of a complicated, biologically realistic scenario. To achieve good spatial resolution, permitting detailed study of a cell with substantial extent, densely packed electrode arrays are required. In the present reconstruction we assumed a hexagonal grid arrangement with 17.5 μm interelectrode distance inspired by recent experiments on reconstructing axonal action potential propagation (Bakkum et al., 2013; Frey et al., 2009). We assumed a grid consisting of 936 contacts from which we used 128 for reconstruction to be consistent with the hardware of (Bakkum et al., 2013; Frey et al., 2009).
In the simulation we assumed an experimentally plausible scenario, where oscillatory current was injected to the soma of a neuron in a slice with other inputs impinging through a 100 excitatory synapses distributed on the dendritic tree. The simulated data consisted of two parts. During the first 400 ms, the cell was stimulated by the injected current as well as through the synapses. The amplitude of the injected current was 3.6 nA, the frequency of the current drive was around 6.5 Hz. During the second 400 ms the cell was stimulated only with the current. Figure 8 shows an example of the skCSD reconstruction at a time selected right after a spike was elicited by the cell. As we can see, neither the standard CSD reconstruction assuming smooth current distribution in space, nor the interpolated potential, give justice to the actual current distribution. At the same time, the skCSD reconstruction is quite a faithful reproduction of the ground truth. A movie comparing the ground truth with kCSD, interpolated potential, and skCSD reconstruction, in time, is provided as a supplementary material (Video 1).
Dependence of reconstruction on noise level
So far, we have assumed that the data are noisefree which is never true in an experiment. Both the measurement device and the neural tissue are potential sources of distorted data. To investigate how the performance of the method is influenced by noise, we added Gaussian white noise of differing amplitudes to the simulated extracellular recordings of Yshaped cell described above. Figure 9A shows the smoothed ground truth we used. The Yshaped neuron is placed on top of a MEA with a regular grid of 4 × 8 electrodes marked by asterisks. Figure 9B shows the noisefree reconstruction. Panel C–F of the figure show the reconstruction results for increasing measurement noise with signal to noise ratio, SNR$=16,4,1$. The signaltonoise ratio (SNR) here is the standard deviation of the simulated extracellular potentials normalized with the std of the added noise. The degradation of reconstruction visible in this figures is summarized in Figure 9C. As can be seen in the reconstruction plots (Figure 9D–F), the increasing noise actually does not seem to significantly alter the obtained reconstructions so the regularization is providing adequate correction, except for the noise on the order of signal (Figure 9F).
Dependence of reconstruction on the number and arrangement of recording electrodes
Reconstruction of the distribution of the current sources along the morphology with skCSD (just like the reconstructions of smooth population distributions with kCSD) formally can be attempted from an arbitrary set of recordings, even a single electrode. While we do not expect enlightening results at this extreme, it is natural to ask the following questions: (1) to what extent can we trust the reconstruction in a given case, (2) which of the reconstructed features are real and which are artifacts of the method, and (3) how can the optimal parameters be selected for this method. We will return to these issues in the Discussion. Here, we wish to investigate how the number of electrodes, the density of the grid, and the area covered by the MEA, affect the results.
To answer these questions, we selected a snapshot of simulation of the model of the ganglion cell described in the Materials and methods section, with the specific membrane current distribution shown in Figure 10A. In Figure 10B–H. we show seven different reconstructions assuming different experimental setups, with differing numbers of electrodes, covering different area.
In each case, we selected the width of basis functions and the regularization parameter for the method by minimizing the L1 error calculated for the first 1000 time steps of the simulation or crossvalidation error (L1T and CV in Figure 10I). To verify the quality of reconstruction we computed the L1 error between the ground truth and reconstruction for the remaining 5800 time steps of the simulation. We found that minimization of L1 error gave better results and L1V in Figure 10I shows the results for this case; however, the results obtained with minimization of CV error were often not much worse (not shown).
Given that L1 error can only be used where the ground truth is known, which is in simulations, we propose the following. Given the data necessary for application of the skCSD method, (neuronal morphology, positions of electrode contacts, and recorded signals) different CSD distributions should be assumed for the obtained morphology, reconstructions obtained for a range of parameters, then the L1 error could be used for optimization. Note that it is not necessary to actually simulate a model of the cell with proper membrane biophysics, which often is not known, although it might lead to more physiological test sources. It is sufficient to distribute the different sources along the morphology without making any assumptions concerning the biophysical properties of the neuronal membrane.
Once the parameters are obtained with this procedure, perform the analysis of actual experimental data with the obtained parameters. Performing the simulations and comparing the best reconstructions with the assumed ground truth has the further benefit of building intuition about which features of the real CSD survive in the reconstruction and which are distorted. This is another example of modelbased data analysis which we believe becomes inevitable with the growing complexity of experimental paradigms, such as the one considered here.
We feel that the above procedure is optimal, since it not only gives optimal parameters, but also allows one to investigate which features are recovered and which are misformed. However, if only parameters for estimation are needed, CV error could be used, which is simpler and the results are often comparable.
The results obtained in this study are consistent with our expectations: the quality of reconstruction improves with the coverage of the morphology by the electrodes, with increasing density of probing, and with increasing number of probes (Figure 10I). Interestingly, it seems, that it is difficult to improve the reconstruction beyond a certain level, in consequence, the setups with moderate densities (on the order of 200 μm IED) can easily compete with setups at the edge of current developments (40 μm IED, [(Berdondini et al., 2005). We believe that this is not a hard limit and that better results can be obtained here. This, however, requires further development of the methods.
Proofofconcept experiment: spatial current source distribution of spiketriggered averages
To examine the experimental feasibility of the skCSD method, we analyzed data from a patch clamp electrode and a linear probe with 14 working electrodes recording signals simultaneously from a hippocampal pyramidal cell in an in vitro slice preparation (see Materials and methods). As there is no ground truth data available in this case, the optimal width of the basis functions and the regularization parameter were selected using the L1 error and simulated data. To do this, we used the same simulation protocol as for the ganglion cell model. A snapshot of the reconstruction is shown in Figure 11 at the moment of firing. A 10 ms long video of the spike triggered average is shown in the supplementary materials (Video 2). At −0.05 ms the brief appearance of a sink (red) in the basal dendrites is visible which can be a consequence of the activation of voltage sensitive channels in the axon hillock, or the first axonal segment leading to the firing of the cell. Since there were no electrodes close to the axon initial segment, the skCSD method did not resolve it, instead it resolved to introduce the activity into the basal dendrite. This phenomenon is quickly replaced by a sink at the soma and in the proximal part of the apical dendritic tree, accompanied by sources (blue) in the basal and in the more distal apical dendrites. The extracellular potential on the second electrode reaches its minimum at 0.45 ms, which signals the peak of the spike. The deep red of the soma at this point signifies a strong sink, while the blue of the surrounding parts of the proximal apical and basal dendrites indicate the current sources set by the return currents. At 1.30 ms a source appears at the soma region, which indicates hyperpolarizing currents. Overall, the observed spatiotemporal CSD dynamics is dominated mostly by the somatic currents, responsible for the spike generation, and the corresponding counter currents. This example demonstrates the experimental feasibility of the skCSD method and may help in planning further experiments, aiming to reveal the spatial distribution and temporal dynamics of the synaptic input currents which evoke the firing of the neuron.
Discussion
Summary
In this work, we introduced a method to estimate the distribution of current sources (CSD) along the dendritic tree of a neuron given its known morphology and a set of simultaneous extracellular recordings of potential generated predominantly by this cell. First, assuming the ballandstick neuron model and a laminar probe parallel to the cell, we studied the basic viability of this method. We showed that introducing more electrodes to cover the same area leads to increased spatial resolution of the method allowing reconstruction of higher Fourier modes of the CSD generating the measured potentials (Figure 2). In a dynamic scenario of multiple synaptic inputs impinging on the cell, higher density of probes leads to higher reconstruction precision allowing us to distinguish individual inputs (Figure 3). Testing the reconstruction against the known CSD (the ground truth) shows a clear transition from poor to faithful reconstruction when the electrode distribution becomes dense enough to capture the fine detail of the CSD profile to be reconstructed (Figure 2.E). We also applied this method to more complex neuron morphologies, namely the Yshape and ganglion cell. As expected, the reconstructed CSD profiles became more detailed with increasing electrode number over a fixed area (Figures 6,10).
Using the Yshaped morphology we showed that (i) synaptic inputs activating different dendrites can be separated, Figure 5; (ii) skCSD provides meaningful information about the membrane CSD in cases, when the interpolated LFP and standard, population CSD analyses, are not informative, Figure 5; (iii) the reconstruction is not sensitive to a specific selection of electrode placement, Figures 6,9; and (iv) even significant additive noise (SNR = 1) is not prohibitive for the reconstruction, Figure 9.
Biologically, the most relevant example we considered was a ganglion cell model which we studied with virtual multielectrode arrays of different designs. The MEAs we considered differed with interelectrode distances for the simulated setups, as well as in the area they covered, ranging from an area close to the soma to roughly four times the size of the smallest square covering the whole morphology. The best results were obtained when we used the electrodes from the region which closely covered the cell (9.G and H); a reduction of interelectrode distance from 100 μm to 40 μm yielded less impressive results than selecting the electrodes from the smallest square covering the morphology. Our study, assumed realistic cell morphology of the ganglion cell and commercially available MEA designs, as well as realistic cell activity, showed that it is feasible to reconstruct the distribution of the current sources in realistic, noisy situations.
The skCSD method performed adequately for the proofofconcept experimental data, even if the nature of the experiment allowed only the reconstruction of the general features of the spiketriggered average spatiotemporal current source density distribution patterns.
While the traditional (population) CSD method was used to analyze membrane currents of single cells before (Buzsáki and Kandel, 1998; Bereshpolova et al., 2007; Frey et al., 2009), the first CSD method specific for investigating membrane currents on single cells was proposed by (Somogyvári et al., 2005). However, it assumed simplified, linear neuron morphologies. An important preprocessing step proposed there was separating the single neuron’s contributions to the extracellular potentials from the background activity. The novelty of the skCSD method proposed here is in its use of actual neuronal morphologies and in the underlying algorithmic solutions based on the kCSD method (Potworowski et al., 2012) which were initially used in the study of populations of neurons.
General comments
Observe that skCSD, like any other data analysis method, is not a magic wand. Technically, it can be applied to data coming from a single electrode just like the age profile of a human population can be estimated from a single specimen. Obviously, in both cases, the estimate would be a poor reflection of the distribution of interest. As we improve the sampling, the quality of the estimate improves, yet ultimately it is hard to judge a priori how many electrodes is enough and if our results obtain the required level of precision. We see two approaches to address this type of questions. One is through simulations, as we discussed. The other is analysis of the singular vectors arising in the decomposition of the matrix translating the measured potentials into the estimated CSD (Hansen, 2010); however, the necessary tools for kCSD and skCSD are still under development. We plan to investigate this further in the future.
Having obtained the distribution of currents it would be interesting to decompose it into physiologically meaningful components, such as synaptic currents, leak currents, voltagegated currents for different channels, etc. This seems rather challenging and we do not see a direct way of achieving this from experimental data. It is possible that an application of statistical decomposition methods will prove useful, as in the case of kCSD for population activity (Łęski et al., 2010; Głąbska et al., 2014). However, we find the contributions to the extracellular potential from individual currents highly counterintuitive (Głąbska et al., 2017).
Experimental recommendations
To attempt experimental application of skCSD we must have (1) an identified cell of known morphology, and (2) a set of simultaneous extracellular recordings of electric potential generated by this cell. Each aspect poses its challenges, some of which have been addressed here. Once we have the necessary data the natural question is how to select the parameters of the method in the specific context of a given setup, specific morphology, and recordings. Our investigations above give some indications. First, the electrodes selected for analysis should essentially uniformly cover the area spanned by the cell. This is illustrated in Figure 6, which shows that some degree of irregularity does not significantly affect the reconstruction. Secondly, the basis should be selected so that the basis sources could resolve the features on the membrane we are interested in (narrow basis functions) with sufficient multiplicity, that smooth coverage of the cell can be ensured (see Figure 7).
Clearly, as the irregularity of the setup grows we expect growing reconstruction errors. This can be studied with singular vectors of the operator transforming potential measurements into reconstructed sources, Equation 21, as discussed for example by (Hansen, 2010). We are convinced, however, that the most efficient approach to investigate the effects of these different parameters is through simulations. This is a natural place to apply the modelbased validation of data analysis (Denker et al., 2014). Our suggestion is to build a computational model of the cell. We believe that for the purpose of parameter selection assuming passive membrane in the dendrites should be sufficient, but of course, more realistic biophysical information may be included, especially if available. The model cell may be stimulated with synaptic input, with current injected, or even specific profiles of ground truth CSD may be placed along the cell. Then the extracellular potential must be computed at points where the actual electrodes are placed in the experiment. One can then investigate the effects of different parameter values on reconstruction and, for the analysis of actual experimental data, select those parameters minimizing prediction error on test data. The advantage of this procedure is twofold. First, we end up with a selection of parameters adapted for the specific problem at hand. Secondly, we build intuition regarding the interpretation of the results for our specific cell and setup. This approach is the only way to address arbitrary electrodecell configurations and to determine how much information we can extract in a given case.
We found that the best way to identify optimal parameters for reconstruction is by minimizing the L1 error between the reconstruction and the ground truth. Since we cannot have the ground truth in an experiment, but we can assume it in the modelbased validation, this is another argument for the modelbased validation approach. Obviously, to efficiently apply this technique, the appropriate computational tools must be available. We plan to develop and open framework facilitating such studies, meanwhile, the code used for the present study is available at https://github.com/csdori/skCSD (Cserpán, 2017; copy archived at https://github.com/elifesciencespublications/skCSD).
Challenges of recording extracellular potentials and obtaining the morphology from the same cell
Although recording extracellular potential with a MEA, filling a neuron with a dye, and reconstructing its morphology, are standard experimental techniques, using them simultaneously remains a challenge due to the size of the experimental devices which need to be arranged within a small volume. Cells in the vicinity of the MEA can be filled individually by intracellular or juxtacellular electrodes, or with bulk dying. Individual recording and dying with a glass electrode provides not only the morphology, but also unambiguous spike times, giving an opportunity to determine the extracellular potential fingerprint of the recorded cell on the MEA. Although these would be favorable data, intracellular recording less than 100 μm from the MEA is extremely challenging. Experimental setups featuring the necessary equipment already exist (Neto et al., 2016), but as far as we know, have not been used in this way. On the other hand, bulk dying techniques result in more filled neurons, although the quality of the dying, and thus the quality of the 3D morphology reconstructions, is considerably lower in these cases. Although there are methods for estimation of the cell position relative to the MEA (Somogyvári et al., 2005, Somogyvári et al., 2012), association of multiple optically labeled neurons with the recorded extracellular spike patterns is still unsolved.
Challenges of separating the activity of a single neuron from background
We propose two experimental scenarios one could apply to separate the activity of the studied neuron from the background. If we can sort the spikes elicited by the neuron of interest we can calculate the spiketriggered averages of the potentials reducing all uncorrelated contributions. Unfortunately, in live tissue, contributions from neighboring cells will have some correlations due to shared inputs. Separation of the contribution of the neuron of interest from the correlated background can be obtained in two ways. One is decomposition of the activity into meaningful components, for example, as shown by (Somogyvári et al., 2015), the high amplitude correlated oscillatory background of hippocampal theta activity can be extracted with independent component analysis, allowing the determination of celltypespecific time course of the synaptic input. Alternatively one could consider combining skCSD with population kCSD analysis, that is, consider basis sources covering the cell of interest as well as the space covering the whole population. This will be the subject of further study. The second experimental scenario to obtain the contributions to the extracellular potential from a specific cell is to drive the cell with intracellular current injection of known pattern, for example, with an oscillatory drive as we discussed (Figure 8), and by averaging over multiple periods (eventbased triggering). Again, further study is needed to establish the validity of this type of experimental procedure.
Challenges of using novel MEAs
Handling data from highdensity MEAs with thousands of electrodes will require further studies, as the large numbers of small singular values of the kernel matrix may introduce numerical sensitivity to the reconstruction. Also, optimal selection of electrodes in case of programmable MEAs merits further investigations. We believe it is best to address such issues when actual experiments are attempted.
Importance of this work
Traditional electrophysiology has focused on the electrical potential, which is relatively easy to access, from intracellular recordings, all kinds of patch clamp, juxtacellular, to extracellular and voltagesensitive dyes (Covey and Carter, 2015). While the relation of the actual measurement to the voltage at a point may significantly differ, often this is a reasonable interpretation, if needed, more realistic models of measurement can be considered, for example, averaging over the contact surface for extracellular electrodes, etc (Moulin et al., 2008; Ness et al., 2015).
Already in the middle of the 20th century, Walter H. Pitts observed that for recordings obtained with regular grids one can approximate the Poisson equation to estimate the distribution of current sources in the tissue (Pitts, 1952). His approach assumed recordings on a regular 3D grid, which was challenging to obtain for some 60 years (Łęski et al., 2007). However, with the work of Nicholson and Freeman (Nicholson and Freeman, 1975) 1D CSD analysis became attractive, as summarized by Ulla Mitzdorf (Mitzdorf, 1985). In 2012, we proposed how to overcome the restriction of regular grids with a kernel approach which both allows to use arbitrary distribution of contacts and corrects for noise (Potworowski et al., 2012). All the previous work, however, always assumed the contributions to the extracellular potential coming from the whole tissue and smooth in the estimation region.
In the present work, we show for the first time how a collection of extracellular recordings in combination with a cell morphology can be used to estimate the current sources located on the cell contributing to the recorded potential. Since it is now feasible experimentally to obtain the relevant data, we believe that the method proposed here may find its uses to constrain the biophysical properties of the neuron membrane, facilitate verification of morphological reconstructions, as well as guide new discoveries by offering a more global picture of the distribution of the currents along the cell morphology, giving a coherent view of the global synaptic bombardment and return currents within a cell.
Materials and methods
Overview of current source density reconstruction methods
Traditional CSD
Request a detailed protocolFor reader’s convenience, here we briefly present the basic ideas behind the traditional and recent approaches to reconstruction of current source density (CSD analysis). For a more complete review of CSD analysis see Wójcik (2015), for recent reviews of the relations between neural activity, current sources and the recordings see (Buzsáki et al., 2012; Einevoll et al., 2013a).
The relation between current sources in the tissue and the recording potentials is given by the Poisson equation
where $C$ stands for CSD and $V$ for the potential. While this can be studied numerically for nontrivial conductivity profiles (Ness et al., 2015), here we shall mostly assume a constant and homogeneous conductivity tensor, $\sigma $. In that case, the above equation simplifies to $C=\sigma \mathrm{\Delta}V$ and can be solved for $C$ given potential in the whole space. On the other hand, given the CSD in the whole space, the potential is given by
Walter Pitts observed that having recordings on a regular grid of electrodes we can estimate CSD by taking numerical second derivative of the potential (Pitts, 1952), we call this approach traditional CSD method. Pitt’s idea gained popularity only after Nicholson and Freeman popularized its use for laminar recordings (Nicholson and Freeman, 1975) in the cortex. In this setup, assuming the layers are infinite and homogeneous (Pettersen et al., 2006), the current source density at each layer can be estimated from
where ${z}_{j}$ is the position of the ${j}^{th}$ electrode and $h$ is the interelectrode distance.
Inverse CSD (iCSD)
Request a detailed protocolTo overcome limitations of the traditional approach, such as difficulty of handling the data at the boundary and hidden assumptions about the dimensions we do not probe, Pettersen et al. proposed a modelbased inverse CSD method (Pettersen et al., 2006). Initially proposed in 1D, the method was later generalized to other dimensionalities (Łęski et al., 2007; Łęski et al., 2011). Given a set of recordings ${V}_{1},\mathrm{\dots},{V}_{N}$ at regularly placed electrodes at ${\mathbf{\mathbf{x}}}_{1},\mathrm{\dots},{\mathbf{\mathbf{x}}}_{N}$ this method assumes a model of CSD parametrized with CSD values at the measurement points, $C(\mathbf{\mathbf{x}})={\sum}_{k=1}^{N}{C}_{k}{f}_{k}(\mathbf{\mathbf{x}})$, where ${f}_{k}(\mathbf{\mathbf{x}})$ are functions taking 1 at ${\mathbf{\mathbf{x}}}_{k}$, 0 at other measurement points, with the values at other points defined by the specific variant of the method, for example, spline interpolated in spline iCSD (Wójcik, 2015). Assuming the model $C(\mathbf{\mathbf{x}})$ one computes the potential at the electrode positions obtaining a relation between the model parameters, ${C}_{k}$, and the measured potential, ${V}_{k}$, which can be inverted leading to an estimate of the CSD in the region of interest.
Kernel CSD (kCSD)
Request a detailed protocolThe kernel Current Source Density method (Potworowski et al., 2012) can be considered a generalization of the inverse CSD. It is a nonparametric method which allows reconstructions from arbitrarily placed electrodes and facilitates dealing with the noise. Conceptually, the method proceeds in two steps. First, one does kernel interpolation of the measured potentials. Next, one applies a 'crosskernel' to shift the interpolated potential to the CSD. In 3D, in space of homogeneous and isotropic conductivity, this amounts to applying the Laplacian to the interpolated potential, Equation. (3). To handle all cases in a general way, including data of lower dimensionality or with nontrivial conductivity, we construct the interpolating kernel and crosskernel from a collection of basis functions. The idea is to consider current source density in the form of a linear combination of basis sources ${\stackrel{~}{b}}_{j}(\mathbf{\mathbf{x}})$, for example Gaussian,
where the number of basis sources $M\gg N$, the number of electrodes, and ${a}_{j}$ are the weights with which the basis sources are combined into the model CSD. Let ${b}_{j}(\mathbf{\mathbf{x}})$ be the contribution to the extracellular potential from ${\stackrel{~}{b}}_{j}(\mathbf{\mathbf{x}})$, which in 3D is
but in 1D or 2D we would need to take into account the directions we do not control in experiment (for example, along the slice thickness for a slice placed on a 2D MEA). Then, the potential will have a form
Since we cannot estimate $M$ coefficients ${a}_{j}$ from $N$ measurements for $N<M$, we construct a kernel for interpolation of the potential,
Then, any potential field $V(\mathbf{\mathbf{x}})$ spanned by ${b}_{i}(\mathbf{\mathbf{x}})$ can be written as
for some $L$, ${\mathbf{\mathbf{x}}}_{l}$, and ${\beta}_{l}$, but it minimizes the regularized prediction error
when $L=N$. Here, ${\mathbf{\mathbf{x}}}_{k}$ are the positions of the electrodes, ${V}_{k}$ are the corresponding measurements, $\lambda $ is the regularization constant. The minimizing solution is obtained for
where $\mathbf{\mathbf{V}}$ is the vector of the measurements ${V}_{k}$, and ${\mathbf{\mathbf{K}}}_{jk}=K({\mathbf{\mathbf{x}}}_{j},{\mathbf{\mathbf{x}}}_{k})$.
To estimate CSD we introduce a crosskernel
If we define
then the estimated CSD takes form of
where $\lambda $ is the regularization parameter and I the identity matrix; see (Potworowski et al., 2012) for derivation and discussion.
Spike CSD (sCSD)
Request a detailed protocolThe Spike CSD (Somogyvári et al., 2012) is the forerunner of the method presented here, as it aims to estimate the current source distribution of single neurons with unknown morphology. The sCSD method provides an estimation of the cellelectrode distance and uses a simplified model of the shape of the neuron to reach this. Separating potential patterns generated by different neurons is critical and it is obtained by clustering extracellular fingerprints of action potentials which are different for every neuron. The limitation of sCSD is the assumed simplified morphology of the model and low spatial resolution. Despite that, even with this simplified model, it was possible to demonstrate for the first time the EC observability of backpropagating action potentials in the basal dendrites of cortical neurons, the forward propagation preceding the action potential on the dendritic tree and the signs of the Ranviernodes (Somogyvári et al., 2012).
skCSD method
Request a detailed protocolThe singlecell kCSD method (skCSD), which we introduce in this work, is an application of the kCSD framework where we assume that the measured extracellular potential comes mainly from a cell of known morphology and known spatial relation to the MEA. To estimate the CSD in this case, we must cover the morphology of the cell with a collection of basis functions. To do this, a onedimensional parametrization of the cell morphology is needed. This could be done independently for each branch of the neuron or globally for the whole cell at once. While the first approach might seem easier, handling of the branching point is nontrivial. Instead, we decided to fit a closed curve on the morphology, which we call the morphology loop (Figure 1). This curve should cover all the segments of the cell, be as short as possible, and be aligned with the morphology. For example, in case of a ballandstick neuron, the curve starts at the soma, goes towards the tip of the dendrite, turns back, goes back to the soma, and closes there. One parameter $s$ is enough to unambiguously determine a position on this line, although most points on the morphology are mapped to two $s$ parameters. We also need a method to handle the branching points and guide the parametrization so that all the branches will be visited in an optimal way. This problem is a special case of the Chinese postman problem known from graph theory (Kwan, 1962). Given this information, we can distribute the basis functions ${\stackrel{~}{b}}_{j}(\mathbf{\mathbf{x}})$ along the morphology of the cell (Figure 1).
In practice, based on the morphology information we define an ordered sequence of all the segments such that the consecutive segments are always physically connected and preference is given to those neighbors which have not been visited yet. The process is continued until all the segments are covered and the last element in the sequence connects to the first element. Note that in the sequence the final segments of the branches are present once, the branching point multiple times and the intermediate ones twice. Then we fit a spline on the coordinates of the segments following the ordered sequence resulting in a morphology loop construction. The CSD basis functions are distributed along this loop uniformly. Any point $\mathbf{\mathbf{x}}\equiv (x,y,z)$ on the morphology can be parameterized with $s\in [0,l]$ on the loop:
where $l$ is twice the length of all the branches. Consider the following basis functions:
where ${s}_{i}$ is the location of the $i$th basis function on the morphology loop, $R$ its width.
The contribution to the extracellular potential from a basis source ${\stackrel{~}{b}}_{i}(s)$ is given by
As in kCSD, for CSD of the form
we obtain the extracellular potential as
As before, for estimation of potential we use kernel interpolation. Note that in this case the basis functions in the CSD space, $\stackrel{~}{b}(s)$, live on the morphology loop, while the basis functions in the potential space, ${b}_{i}(\mathbf{\mathbf{x}})$, live in the physical 3D space. To determine the current source density distribution along the fitted curve, we introduce the following kernel functions:
With these definitions the regularized solution for $C$ on the morphology loop is given by Equation 14:
To obtain the distribution of currents at a given point in space we need to sum the currents on the loop at points which are mapped to that physical position $\mathbf{\mathbf{x}}$:
Construction of ground truth data
Request a detailed protocolTo validate the method, we used simulated data which allows us to consider arbitrary cellelectrode setups and test various current patterns. The LFPy package (Lindén et al., 2013) was used to simulate the extracellular potential at arbitrarily placed virtual electrodes. We assumed the .swc morphology description format (Cannon et al., 1998) and the sections were further divided to segments. The coordinates of every segment’s ends were used to find the connections. Once the connection matrix was calculated, we used the Chinese postman algorithm to obtain the morphology loop. We calculated the potential using neuron models with various morphologies shown in Figure 4 and different input distributions, assuming one and twodimensional multielectrode arrays. We used toy models to better understand and characterize the method as well as a biologically realistic neuron model to estimate performance of skCSD in an experimentally realistic scenario.
The simplest setup we used was a ballandstick neuron recorded with a laminar probe. Various artificial CSD patterns and also biologically more realistic CSD distributions served as test distributions in order to quantify the spatial resolution and reconstruction errors. To generate the ground truth data we simulated a 500 μm long linear cell model of 52 segments in LFPy. The diameter of the two segments representing the soma was 20 μm, while the other segments were 4 μm wide. 100 synaptic excitation events were distributed randomly along this morphology in order to imitate a biologically realistic scenario.
To test the effect of branching on the results, a simple Yshaped morphology was used (Figure 4B). The synapses were placed at segments 33 and 62 on different branches close to the branching point. The first was stimulated at 5, 45, 60 ms, the other at 5, 25, 60 ms after the onset of the simulation. The idea here was to consider the inputs stimulated together and separately. The times of activation were randomly selected in such a way as to leave enough time for the membrane activity to settle down. This can be viewed as extreme cases of correlated and uncorrelated events. Note that the skCSD reconstruction is not affected directly by the temporal correlation of the synaptic inputs. Just like any CSD estimation method, skCSD is applied to the potentials recorded at a given point in time. Of course, it is affected indirectly, in the sense that slower or faster oscillating inputs lead to different spatial patterns due to filtering effects of the dendritic membrane: fast oscillations induce short dipoles, slow oscillations allow the current to spread along the cell leading to stronger dipoles (Lindén et al., 2010). We address these effects indirectly in Figure 2.
As a realistic example, we used a mouse retinal ganglion cell morphology (Kong et al., 2005) from NeuroMorpho.Org (Ascoli, 2006). In the simulations 608 segments were used. 100 synaptic excitation events were distributed randomly along this morphology within the first 400 ms of the simulation. The cell was also driven with an oscillatory current. In the dendrites, only passive ion channels were used.
Parameters of the simulations
Request a detailed protocolWe simulated three different model morphologies: ballandstick (BS), Yshaped (Y), and a ganglion cell (Gang). The Yshaped neuron was considered in two situations, when it was parallel (Y) or orthogonal (Yrot) to the MEA plane. The extracellular potential was computed at multiple points modeling different experimentally viable recording configurations (cell and setup). All combinations used are summarized in Table 1. The parameters describing the neuron membrane physiology are given in Table 2. The length of the simulation was 70 s in case of the ballandstick and Yshaped neurons, and 850 s for the ganglion cell model.
Parameters of synapses
Request a detailed protocolIn most simulations we modeled synaptic activity. We used synapses with discontinuous change in conductance at an event followed by an exponential decay with time constant $\tau $ (ExpSyn model as implemented in the NEURON simulator). When simulating the Yshaped neuron we placed two synapses with the following parameters: reversal potential: 0 mV , synaptic time constant: 2 ms, synaptic weight: 0.04 μS. The synapses were placed at segments 33 and 62 (see Figures 4,5). When simulating the other models (ballandstick and ganglion cell) we used the same type of synapse; however, the synaptic weights were a quarter of the above (0.01 μS) since they were more numerous (Table 2).
Measuring the quality of reconstruction
Request a detailed protocolTo validate the skCSD method, we need to consider two situations. When we know the ground truth — the actual distribution of sources which generated the measured potentials — we can compare the reconstruction with it. This is available directly only in simulations. In that case, we can measure the prediction error between the reconstruction and the original. However, the skCSD method by its nature gives smooth results. This is a consequence of kernel interpolation of the potential which occurs in the first step of the method. The same phenomenon occurs in regular CSD estimation (Wójcik, 2015). Thus, we can never recover the original CSD distribution but only a coarsegrained approximation. This is not a significant problem as the coarsegrained CSD should have equivalent physiological consequence. However, to compare the reconstructed density with the groundtruth, which is typically very irregular in consequence of multiple synaptic activations, we always smoothed the ground truth CSD with a Gaussian kernel. The width of the kernel was 15 μm for ballandstick model, while for the Yshaped and ganglion cell models we used 30 μm.
Thus, whenever ground truth was known, we computed L1 norm of the difference between the reconstruction ${C}^{*}$ and smoothed ground truth $C$ normalized by the L1 norm of $C$:
When analyzing experimental data we only have access to the noisy measurements and cannot apply the above strategy directly. Thus we consider two strategies. One is to use crossvalidation error (CV). In leaveoneout crossvalidation (Potworowski et al., 2012) we estimate CSD from all the measurements but one and compare estimated prediction with actual measurement on the removed electrode. Repeating this procedure for all the electrodes gives us a measure of prediction quality for a given set of parameters for this specific dataset. Scanning over some parameter range we identify optimal parameters as those giving minimum error. They are further used to analyze the complete data. The advantage of using crossvalidation error is that it does not require the knowledge of the ground truth current source density distribution and can still provide an estimation about the performance of the skCSD method. As this algorithm is quadratic in the number of electrodes, for large arrays one might prefer to use the leavepout crossvalidation instead. When we test how the quality of the reconstruction changes with the number of electrodes we use CV error normalized by the number of electrodes which can then be compared between different setups.
The other strategy we use and recommend in the experimental context, when we know the cell morphology and its geometric relation to the setup, as well as the measurements, is modelbased analysis. The idea is to simulate different current source distributions, either placing specific distribution by hand or by modeling activity of the cell assuming passive membrane and random or specific synaptic activations, both of which are relatively inexpensive both in computational time and coding complexity. This reduces the problem to the modeling case. We can use thus generated data (CSD and potentials) scanning for optimal reconstruction parameters to be used in analysis of actual experimental data from the setup.
To handle the effects of noise one should study its properties on electrodes, for example, assuming white measurement noise identify its variance, then tune the regularization parameter $\lambda $ on simulated sets with comparable simulated noise added.
Parameter selection
Request a detailed protocolTo apply the skCSD method, we need to decide upon the number of basis functions, set their width ($R$), and choose the regularization parameter $\lambda $. In this work, the number of basis function was set to 512 for all cases, which is at least twice the number of electrodes used. This is usually not a limitation, the more the better. For the basis width (Equation. (16)) we took the following values: 8, 16, 32, 64, 128 μm. Selection of the regularization parameter is not trivial (Potworowski et al., 2012; Hansen, 2010). Here, we tested the effect of the regularization parameter taking values of 0.00001, 0.0001, 0.001, 0.01, 0.1 The optimal parameters were identified by the lowest value of reconstruction error.
Visual representation of CSD on the morphology
Request a detailed protocolTo visualize the distribution of current sources and other quantities along a neuron morphology we use two representations of the cell:
Interval representation: we stack all the compartments consecutively along the yaxis so that the part of the dendrite stemming from the soma is shown first, followed by one branch, followed by the other. The order of the branches in the stack is taken from the morphology loop to make these representations consistent. The $x$axis either shows different time instants of the simulations or various distribution patterns.
Branching morphology representation: in this case a twodimensional projection of the cell is shown which is colored according to the amplitudes of the membrane current source densities at a time instant. To visually enhance the current events, gray circles proportional to the amplitude of CSD at a point are placed centered at the point to facilitate comprehension.
Experimental methods
In vitro experiment
Request a detailed protocolOne male Wistar rat (300 g) was used for the slice preparation procedure. The in vitro experiment was performed according to the EC Council Directive of November 24, 1986 (86/89/EEC) and all procedures were reviewed and approved by the local ethical committee and the Hungarian Central Government Office (license number: PEI/001/6959/2015). The animal was anesthetized with isoflurane (0.2 ml/100 g). Horizontal hippocampal slices of 500 µm thickness were cut with a vibratome (VT1200s; Leica, Nussloch, Germany). We followed our experimental procedures developed for human in vitro recordings (Kerekes et al., 2014), adapted to rodent tissue. Briefly, slices were transferred to a dual superfusion chamber perfused with artificial cerebrospinal fluid. Intracellular patchclamp recordings, cell filling, visualization and threedimensional reconstruction of the filled cell was performed as described in (Kerekes et al., 2014). For the extracellular local field potential recordings, we used a 16channel linear multielectrode (A16 × 1–2 mm50177A16, Neuronexus Technologies, Ann Arbor, MI), with an INTAN RHD2000 FPGAbased acquisition system (InTan Technologies, Los Angeles, CA). The system was connected to a laptop via USB 2.0. Wideband signals (0.1−7500 Hz) were recorded with a sampling frequency of 20 kHz and with 16bit resolution. The recorded neuron was held by a constant −40 nA current injection.
Data preprocessing
Request a detailed protocolOne hundred and fiftyfour spikes were detected on the 180 s long intracellular recording by 0 mV upward threshold crossing. A ± 5 ms wide time windows were cut around the moments of each spikes on each channels of the extracellular (EC) potential recordings and averaged, to access the fine details of the EC spatiotemporal potential pattern which accompanied the firing of the recorded neuron on all channels. Two channels were broken (2, 5); however, as the skCSD method allows retrieving CSD maps from arbitrarily distributed contacts, this has not prevented the analysis; the broken channels were excluded from further consideration. The averaged spatiotemporal potential maps were highpass filtered by subtracting a moving window average with 100 ms width. This filtering, together with the spiketriggered averaging procedure, ensured that the resulted EC potential map contains only the contribution from the actually recorded cell. The price we paid was filtering out EC signals of the spontaneous repetitive sharpwave like activity of the slice which was correlated by the firing of the recorded neuron and thus the presumptive synaptic inputs of the recorded neuron as well. An additional temporal smoothing by a moving average with 0.15 ms window was used to reduce the effect of noise.
Data availability

Neuron morphology Badea2011Fig2DuPublicly available at NeuroMorpho.Org (ID : NMO_10743).
References

Mobilizing the base of neuroscience data: the case of neuronal morphologiesNature Reviews Neuroscience 7:318–324.https://doi.org/10.1038/nrn1885

Highdensity electrode array for imaging in vitro electrophysiological activityBiosensors and Bioelectronics 21:167–174.https://doi.org/10.1016/j.bios.2004.08.011

Dendritic backpropagation and the state of the awake neocortexJournal of Neuroscience 27:9392–9399.https://doi.org/10.1523/JNEUROSCI.221807.2007

The origin of extracellular fields and currentsEEG, ECoG, LFP and spikesNature Reviews Neuroscience 13:407–420.https://doi.org/10.1038/nrn3241

Somadendritic backpropagation of action potentials in cortical pyramidal cells of the awake ratJournal of Neurophysiology 79:1587–1591.

Largescale recording of neuronal ensemblesNature Neuroscience 7:446–451.https://doi.org/10.1038/nn1233

An online archive of reconstructed hippocampal neuronsJournal of Neuroscience Methods 84:49–54.https://doi.org/10.1016/S01650270(98)000910

ReportReport From the 1st INCF Workshop on Validation of Analysis MethodsTechnical Report INCF.

Modelling and analysis of local field potentials for studying the function of cortical circuitsNature reviews. Neuroscience 14:770–785.https://doi.org/10.1038/nrn3599

Principles of Neural Coding37–59, Local field potential: biophysical origin and analysis, Principles of Neural Coding.

Microelectronic system for highresolution mapping of extracellular electric fields applied to brain slicesBiosensors and Bioelectronics 24:2191–2198.https://doi.org/10.1016/j.bios.2008.11.028

On the origin of the extracellular action potential waveform: A modeling studyJournal of Neurophysiology 95:3113–3128.https://doi.org/10.1152/jn.00979.2005

Discrete Inverse ProblemsSociety for Industrial and Applied Mathematics, Discrete Inverse Problems.

Intracellular features predicted by extracellular recordings in the hippocampus in vivoJournal of Neurophysiology 84:390–400.

Diversity of ganglion cells in the mouse retina: unsupervised morphological classification and its limitsThe Journal of comparative neurology 489:293–310.https://doi.org/10.1002/cne.20631

Extracting functional components of neural dynamics with Independent Component Analysis and inverse Current Source DensityJournal of Computational Neuroscience 29:459–473.https://doi.org/10.1007/s1082700902031

Intrinsic dendritic filtering gives lowpass power spectra of local field potentialsJournal of Computational Neuroscience 29:423–444.https://doi.org/10.1007/s1082701002454

Threedimensional localization of neurons in cortical tetrode recordingsJournal of Neurophysiology 106:828–848.https://doi.org/10.1152/jn.00515.2010

Dipole characterization of single neurons from their extracellular action potentialsJournal of Computational Neuroscience 32:73–100.https://doi.org/10.1007/s1082701103410

Current sourcedensity method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomenaPhysiological Reviews 65:37–100.

A new 3D finiteelement model based on thinfilm approximation for microelectrode array recording of extracellular action potentialIEEE Transactions on Biomedical Engineering 55:683–692.https://doi.org/10.1109/TBME.2007.903522

Spike Detection for Large Neural Populations Using High Density Multielectrode ArraysFrontiers in Neuroinformatics 9:28.https://doi.org/10.3389/fninf.2015.00028

Singlechannel currents recorded from membrane of denervated frog muscle fibresNature 260:799–802.

Validating silicon polytrodes with paired juxtacellular recordings: method and datasetJournal of Neurophysiology 116:892–903.https://doi.org/10.1152/jn.00103.2016

Theory of current sourcedensity analysis and determination of conductivity tensor for anuran cerebellumJournal of Neurophysiology 38:356–368.

Revealing neuronal function through microelectrode array recordingsFrontiers in Neuroscience 8:423.https://doi.org/10.3389/fnins.2014.00423

BookInvestigations on synaptic transmissionIn: Von Foerster H, editors. Cybernetics: Transactions of the 9th Conference. New York: Josiah Macy Foundation. pp. 159–166.

Kernel current source density methodNeural Computation 24:541–575.https://doi.org/10.1162/NECO_a_00236

Past, present and future of spike sorting techniquesBrain Research Bulletin 119:106–117.https://doi.org/10.1016/j.brainresbull.2015.04.007

BookLearning with Kernels: Support Vector Machines, Regularization, Optimization, and BeyondMIT press.

Determination of spatio temporal input current patterns of single hippocampal neurons based on extracellular potential measurementsProgram No. 267.02 2015 Neuroscience Meeting Planner. Washington, DC: Society for Neuroscience 2015. Online.

Localization of singlecell current sources based on extracellular potential patterns: the spike CSD methodEuropean Journal of Neuroscience 36:3299–3313.https://doi.org/10.1111/j.14609568.2012.08249.x

Modelbased source localization of extracellular action potentialsJournal of Neuroscience Methods 147:126–137.https://doi.org/10.1016/j.jneumeth.2005.04.002

Encyclopedia of Computational Neuroscience915–922, Current Source Density (CSD) Analysis, Encyclopedia of Computational Neuroscience.
Article and author information
Author details
Funding
Ministerstwo Nauki i Szkolnictwa Wyższego (2729/7.PR/2013/2)
 István Ulbert
Nemzeti Kutatási, Fejlesztesi és Innovacios Hivatal (K 113147)
 Zoltán Somogyvári
Nemzeti Agykutatasi Program (KTIA NAP 13120130001)
 Zoltán Somogyvári
Nemzeti Kutatasi, Fejilesztesi es Innovacios Hivatal (NN 118902)
 Zoltán Somogyvári
Nemzeti Agykutatasi Program (KTIA13NAPAIV/1 2 3 4 6)
 Daniel K Wójcik
Nemzeti Kutatasi, Fejlesztesi es Innovacios Hivatal (K119443)
 Lucia Wittner
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
Research supported by grants from the Polish Ministry of Science and Higher Education 2729/7.PR/2013/2, the Hungarian National Research, Development and Innovation Fund NKFIH K 113147, K 119443 and NN 118902, and the Hungarian National Brain Research Program KTIA NAP 13120130001 and KTIA13NAPAIV/1,2,3,4,6. We are grateful for Emese Pálfi and László Négyessy for the opportunity of using their Neurolucida setup at the Department of Anatomy, Histology and Embryology, Semmelweis University, to Chaitanya Chintaluri and Mark J. Hunt for critical reading of the manuscript, and to Joanna JędrzejewskaSzmek for help with testing and documenting the software for public exposition.
Ethics
Animal experimentation: The in vitro experiment was performed according to the EC Council Directive of November 24, 1986 (86/89/EEC) and all procedures were reviewed and approved by the local ethical committee and the Hungarian Central Government Office (license number: PEI/001/6959/2015).
Version history
 Received: June 7, 2017
 Accepted: November 16, 2017
 Accepted Manuscript published: November 17, 2017 (version 1)
 Version of Record published: December 5, 2017 (version 2)
 Version of Record updated: December 5, 2017 (version 3)
 Version of Record updated: February 20, 2018 (version 4)
Copyright
© 2017, Cserpán et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,034
 views

 379
 downloads

 12
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
Tissueclearing and labeling techniques have revolutionized brainwide imaging and analysis, yet their application to clinical formalinfixed paraffinembedded (FFPE) blocks remains challenging. We introduce HIFClear, a novel method for efficiently clearing and labeling centimeterthick FFPE specimens using elevated temperature and concentrated detergents. HIFClear with multiround immunolabeling reveals neuron circuitry regulating multiple neurotransmitter systems in a whole FFPE mouse brain and is able to be used as the evaluation of disease treatment efficiency. HIFClear also supports expansion microscopy and can be performed on a nonsectioned 15yearold FFPE specimen, as well as a 3month formalinfixed mouse brain. Thus, HIFClear represents a feasible approach for researching archived FFPE specimens for future neuroscientific and 3D neuropathological analyses.

 Neuroscience
Pavlovian fear conditioning has been extensively used to study the behavioral and neural basis of defensive systems. In a typical procedure, a cue is paired with foot shock, and subsequent cue presentation elicits freezing, a behavior theoretically linked to predator detection. Studies have since shown a fear conditioned cue can elicit locomotion, a behavior that  in addition to jumping, and rearing  is theoretically linked to imminent or occurring predation. A criticism of studies observing fear conditioned cueelicited locomotion is that responding is nonassociative. We gave rats Pavlovian fear discrimination over a baseline of reward seeking. TTLtriggered cameras captured 5 behavior frames/s around cue presentation. Experiment 1 examined the emergence of dangerspecific behaviors over fear acquisition. Experiment 2 examined the expression of dangerspecific behaviors in fear extinction. In total, we scored 112,000 frames for nine discrete behavior categories. Temporal ethograms show that during acquisition, a fear conditioned cue suppresses reward seeking and elicits freezing, but also elicits locomotion, jumping, and rearing  all of which are maximal when foot shock is imminent. During extinction, a fear conditioned cue most prominently suppresses reward seeking, and elicits locomotion that is timed to shock delivery. The independent expression of these behaviors in both experiments reveal a fear conditioned cue to orchestrate a temporally organized suite of behaviors.