The goal of this section is to provide an overview of the proposed method for non-specialists, 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 one-dimensional 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.

Figure 1

Download asset Open asset

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 ${\tilde{b}}_j(\mathbf{𝐱})$. There is a large flexibility here, but in practice we use Gaussians, so ${\tilde{b}}_j(\mathbf{𝐱})\propto \exp (-{(\mathbf{𝐱}-{\mathbf{𝐱}}_j)}^2/2R^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{𝐱})$ (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 high-dimensional 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{𝐱}}_l,\mathbf{𝐱})$ is a linear combination of all the basis functions spanning the potential space: $K({\mathbf{𝐱}}_l,\mathbf{𝐱})={\sum }_{j=1}^M{b}_j({\mathbf{𝐱}}_l)b_j(\mathbf{𝐱})$. 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{𝐱}}_l$ for which the above equation is solvable. This can be done, this is another way of saying that the functions $b_j(\mathbf{𝐱})$ form a basis. It is thus fair to say that the kernel functions $K({\mathbf{𝐱}}_l,\mathbf{𝐱})$ and $b_j(\mathbf{𝐱})$ 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{𝐱}}_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 so-called 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{𝐕}$ is the vector of the measurements $V_k$, I the identity matrix, ${\mathbf{𝐊}}_{jk}=K({\mathbf{𝐱}}_j,{\mathbf{𝐱}}_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 cross-kernel function,

Defining

the estimated CSD is

The skCSD is a model-based 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 cross-validation 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.

Here, we consider the simplest neuron morphology, the so-called ball-and-stick 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.

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.

Figure 2

Download asset Open asset

Using the ball-and-stick 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 ground-truth, and the precision of reconstruction improves with an increasing number of contacts, which is reflected in the reduction of cross-validation 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.

Figure 3

Download asset Open asset

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

Figure 4

Download asset Open asset

Figure 5

Download asset Open asset

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 spline-interpolated 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 high-definition 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.

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.

Figure 6

Download asset Open asset

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 Y-shaped 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 ${\displaystyle \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 over-fitting (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.

Figure 7

Download asset Open asset

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

Figure 8

Download asset Open asset

Video 1

Download asset

So far, we have assumed that the data are noise-free 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 Y-shaped cell described above. Figure 9A shows the smoothed ground truth we used. The Y-shaped neuron is placed on top of a MEA with a regular grid of 4 × 8 electrodes marked by asterisks. Figure 9B shows the noise-free 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 signal-to-noise 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).

Figure 9

Download asset Open asset

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.

Figure 10

Download asset Open asset

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 cross-validation error (L1-T 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 L1-V 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 model-based 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.

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 spatio-temporal 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.

Figure 11

Download asset Open asset

Video 2

Download asset

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 ball-and-stick 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 Y-shape 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 Y-shaped 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 multi-electrode arrays of different designs. The MEAs we considered differed with inter-electrode 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 inter-electrode 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 proof-of-concept experimental data, even if the nature of the experiment allowed only the reconstruction of the general features of the spike-triggered average spatio-temporal 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.

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, voltage-gated 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 counter-intuitive (Głąbska et al., 2017).

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 model-based 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 two-fold. 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 electrode-cell 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 model-based validation, this is another argument for the model-based 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/elifesciences-publications/skCSD).

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.

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 spike-triggered 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 cell-type-specific 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 (event-based triggering). Again, further study is needed to establish the validity of this type of experimental procedure.

Handling data from high-density 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.

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 voltage-sensitive 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.

The single-cell 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 one-dimensional 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 non-trivial. 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 ball-and-stick 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 ${\tilde{b}}_j(\mathbf{𝐱})$ 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{𝐱}\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 ${\tilde{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, $\tilde{b}(s)$, live on the morphology loop, while the basis functions in the potential space, $b_i(\mathbf{𝐱})$, 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{𝐱}$:

To validate the method, we used simulated data which allows us to consider arbitrary cell-electrode 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 two-dimensional 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 ball-and-stick 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 Y-shaped 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.

We simulated three different model morphologies: ball-and-stick (BS), Y-shaped (Y), and a ganglion cell (Gang). The Y-shaped neuron was considered in two situations, when it was parallel (Y) or orthogonal (Y-rot) 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 ball-and-stick and Y-shaped neurons, and 850 s for the ganglion cell model.

In 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 Y-shaped 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 (ball-and-stick 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).

To 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 coarse-grained approximation. This is not a significant problem as the coarse-grained CSD should have equivalent physiological consequence. However, to compare the reconstructed density with the ground-truth, 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 ball-and-stick model, while for the Y-shaped 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 cross-validation error (CV). In leave-one-out cross-validation (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 cross-validation 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 leave-p-out cross-validation 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 model-based 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.

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

To visualize the distribution of current sources and other quantities along a neuron morphology we use two representations of the cell:

1. Interval representation: we stack all the compartments consecutively along the y-axis 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.

2. Branching morphology representation: in this case a two-dimensional 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.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

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 13-1-2013-0001 and KTIA-13-NAP-A-IV/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ędrzejewska-Szmek for help with testing and documenting the software for public exposition.

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/695-9/2015).

1. Received: June 7, 2017
2. Accepted: November 16, 2017
3. Accepted Manuscript published: November 17, 2017 (version 1)
4. Version of Record published: December 5, 2017 (version 2)
5. Version of Record updated: December 5, 2017 (version 3)
6. Version of Record updated: February 20, 2018 (version 4)

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