kCSD-python, reliable current source density estimation with quality control

Interpretation of the extracellular recordings can be difficult due to the long range of electric field but can be facilitated by estimating the density of current sources (CSD). Here we introduce kCSD-python, an open Python package implementing Kernel Current Source Density (kCSD) method, and introduce several new techniques to facilitate CSD analysis of experimental data and interpretation of the results. We investigate the limitations imposed by noise and assumptions in the method itself. kCSD-python allows CSD estimation for arbitrary distribution of electrodes in 1D, 2D, and 3D, assuming distributions of sources in tissue, a slice, or in a single cell, and includes a range of diagnostic aids. We demonstrate its features in a Jupyter notebook tutorial to facilitate uptake by the community.

1 Introduction 24 Extracellular potential recordings are a mainstay of neurophysiology. However, the long 25 range of electric field still makes their interpretation challenging despite decades of our 26 experience. Extracellular potential in tissue is produced by transmembrane currents. Its 27 low-frequency part, called the Local Field Potential (LFP), is believed to mainly reflect 28 dendritic processing of synaptic inputs (Nunez and Srinivasan, 2006;Buzsáki, Anastassiou, 29 and Koch, 2012). To facilitate understanding of the processes underlying the recorded 30 signal it is useful to estimate the density of transmembrane current sources (Current Source 31 Density, CSD) (Pitts, 1952;Nicholson and Freeman, 1975;Mitzdorf, 1985;Einevoll et 32 al., 2013;Gratiy et al., 2017). CSD gives direct access to the physiologically relevant 33 information, which is often concealed in original data (Mitzdorf, 1985). The relation 34 between the CSD and the extracellular potential can be described by the Poisson equation 35 where C is the CSD, V is the extracellular potential, and σ -the conductivity tensor. 36 For isotropic and homogeneous tissue Eq. (1) reduces to which can be solved: Eq.
(3) show that knowing the potential in the whole space, we can com-39 pute the CSD, and knowing the CSD in the whole space, we can compute the potential. 40 Experimentally, we can only access the potential at discrete electrode locations, so direct 41 determination of the CSD in the whole space from Eq. (2) is impossible. To deal with this 42 problem different methods for estimation of current sources have been proposed since the 43 middle of the past century (Pitts, 1952;Nicholson and Freeman, 1975 Fig. 1 shows the different 53 experimental scenarios for which this software is applicable. 54 We also introduce some new diagnostic features which we found useful in CSD analysis 55 and illustrate the application of the package and the new diagnostic tools in typical analysis 56 workflow implemented as a Jupyter notebook (Kluyver et al., 2016) tutorial. 57 For reader convenience we first briefly restate the kCSD method (Potworowski et al.,58 2012; Chintaluri et al., 2021). Then in the Results section we introduce several new di-59 agnostic tools included in the presented package. First, we added L-curve method for 60 parameter selection. Then, to get informed on reconstruction accuracy we introduce mea-61 surement uncertainty maps, which show how measurement noise affects the estimated 62 CSD, and reliability maps, which help build intuition on reliability of estimates for classes 63 of possible sources for a given setup. We then briefly introduce the new package which 64 is extensively illustrated in the provided tutorial in the Supplementary Materials. This 65 jupyter notebook (Kluyver et al., 2016) tutorial, which is part of the toolbox, is provided 66 to facilitate understanding and usage of the kCSD method. The tutorial enables the user 67 to conduct analysis of the CSD using simulated surrogate data (electrodes positions and 68 recordings) or actual recordings. Kernel Current Source Density method uses kernel interpolation of measured potential to 71 obtain V(x) in the whole space: where x j , j = 1, . . . , N , are electrode positions, K(x, x ′ ) is a symmetric kernel function. 73 To avoid overfitting, correction for noise is made by minimizing prediction error Figure 1: Overview of experimental contexts where kCSD-python is applicable. 1D setups such as A) laminar probes and equivalent; 2D setups, such as B) multishaft silicon probes, Neuropixel or SiNAPS probes, or D) planar MEA; 3D electrode setups, such as multiple multishaft silicon probes, Utah arrays, multiple electrodes placed independently in space with controlled positions, where the sources are assumed to come C) from tissue (kCSD) or E) from single cells with known morphology (skCSD). For description of parameters see Methods.
which gives where V is the vector of measured potentials, λ is regularization parameter, and In the second step, we use function K(x, x ′ ), which we call cross-kernel, to estimate the 77 CSD: This procedure would work for arbitrary smoothing kernels K but in general it is 79 difficult to identify the relevant corresponding K, except in the simple cases. To circumvent 80 this in kCSD we introduce a large basis of CSD sources spanning the region of interest, 81 b j (x), think of many little Gaussians, and corresponding basis sources in the potential 82 space, b j (x), and construct both kernels from these basis sources (Potworowski et In these basis the CSD and the potential are given by The challenges of the method are how to construct K and K, how to select the basis, 86 the relevant parameters, and reliability of the estimation, some of which were introduced 87 before (Potworowski et al., 2012; Chintaluri et al., 2021) and some we introduce here and 88 illustrate with the provided package. For more details of kCSD method and the theory 89 behind, see (Potworowski et al., 2012;Chintaluri et al., 2021). and has a bare minimal library requirements (numpy, scipy and matplotlib). It can be 97 installed using the Python package installer (pip) or using Anaconda python package en-98 vironment (conda). Details of the installation can be found in the package documentation 99 at https://kcsd-python.readthedocs.io/en/latest/INSTALL.html.

100
The package contains a set of tools for kCSD analysis and to validate the results 101 obtained from this analysis. To facilitate uptake of this resource, the package comes 102 with extensive tutorials implemented in jupyter notebooks. These tutorials allow users 103 to test different configurations of current sources and electrodes to see the method in 104 action. The users can analyze their own data or explore the method with data generated 105 in silico. These provisions illustrate the advantages and limitations of kCSD method to 106 its users. The tutorials can also be accessed without any installation on a web browser 107 via Binder (Project Jupyter et al., 2018). The package is extensively documented (https: 108 //kcsd-python.readthedocs.io) and includes all the necessary scripts to generate the 109 figures in this manuscript.

110
An extensive tutorial overview of the kCSD-python package, is provided as an Ap-111 pendix. Its goal is to show how to use kCSD-python to perform CSD analysis, how to 112 apply the provided analysis tools, and to validate the results. We first consider a regular 113 grid of ideal (noise-free) electrodes, where we compute the potentials from a known test 114 source (the ground truth). We then use these potentials to reconstruct the sources which 115 we compare with ground truth (Basic features). Then, we explore the effects of noise on 116 the reconstruction and test the robustness of the method (Noisy electrodes). In the final 117 part of the tutorial we look at how the errors in the estimation depend on the sources and 118 1 supports also Python 2.7 version the electrode configuration by testing the effects of broken electrodes on reconstruction 119 (Broken electrodes). 120 Supplementary Fig. 11 shows error propagation maps, which we introduce below, for 121 1D regular grid of 12 electrodes. Supplementary Fig. 12 shows an example of 3D kCSD 122 source reconstruction. Supplementary Fig. 13 shows an example of skCSD reconstruction 123 (single cell kernel CSD) which corresponds to Fig. 8 from Cserpan et al., 2017. The 124 simulation, reconstruction and visualization have all been re-implemented in Python. 125 We also provide some pre-computed examples of CSD estimations using our library for 126 users to explore. We provide these as .pdf files available at http://bit.ly/kCSD-supplementary. 127 Here, the files small_srcs_3D_all.pdf and large_srcs_3D_all.pdf show 100 exam-128 ple setups of 3D kCSD reconstructions from small and large sources. Similarly, files 129 small_srcs_all.pdf and large_srcs_all.pdf show 100 example setups of 2D kCSD 130 reconstructions from small and large sources. "Smallness" and "largeness" of sources is 131 defined by the ratios of typical spatial scales of the source with respect to interelectrode 132 distances.
where the minimizing vector and where ′i means i-th column and row are removed from the given matrix (or vector). We 146 repeat this for all the electrodes i = 1, . . . , N and compare predictions from the remaining 147 electrodes against actual measurements: For final analysis, λ giving minimum prediction error is selected. It is worth checking if 149 the global minimum is also a local minimum. If the λ selected is one of the limiting values 150 this may indicate that extending the range of λ might result is more optimal result or 151 that the problem is ill-conditioned, for example too noisy, and we are either underfitting 152 or overfitting, as we discuss below for the L-curve. As a rule of thumb, the range of tested 153 λ parameters should cover eigenvalues of the (7).

154
L-curve Consider the error function, Eq. (5), which we minimize to get the regularized 155 solution, V λ = Kβ λ It is a sum of two terms we are simultaneously minimizing, prediction 156 error and the norm of the model weighted with λ. Taking λ = 0 is equivalent to assuming noise-free data. In this case we 159 are fitting the model to the data, in practice, overfitting. On the other hand, taking large 160 λ means assuming very noisy data, in practice ignoring measurements, which results in a 161 flat underfitted solution. Between these extremes there is usually a solution such that if 162 we decrease λ, the prediction error, ϱ, slightly decreases, while the norm of the model, η, 163 increases fast, and if we increase λ, the prediction error, ϱ, increases fast, while the norm 164 of the model, η, slightly decreases; see Fig. 2D. 165 This is apparent when the prediction error and the norm of the model are plotted in 166 the log-log scale. This curve follows the shape of the letter L, hence the name L-curve 167 (P. C. Hansen, 1992). Several methods have been proposed to measure the curvature of the 168 L-curve and to identify optimal parameters (P. C. Hansen, Jensen, and Rodriguez, 2007). 169 In kCSD-python, we have implemented the triangle area method proposed by Castellanos,170 Gómez, and Guerra, 2002. To distinguish between convex and concave plot, clockwise 171 directed triangle area is measured as negative Fig. 2.D shows this estimated curvature for 172 our example as a function of λ.

173
To illustrate this method in the context of CSD reconstructions, we study an example 174 of 1D dipolar current source with a split negative pole (sink; see Fig. 2.C, True CSD, 175 red dashed line). We compute the potential at 32 electrodes ( Fig. 2.A, C, black dots) 176 with additive noise at every electrode. Notice that if we want to interpret the recorded 177 potential directly (Fig. 2.A, red dots) it is difficult to discern the split sink. Fig. 2.D 178 shows the estimated curvature for our example as a function of λ. The optimal value of λ 179 is found by maximizing the curvature of the log-log plot of η versus ϱ, Fig. 2.B. The red 180 dot in Fig. 2.B, D, indicates the ideal λ parameter for this setup obtained through the 181 L-curve method. With kCSD procedure one can easily estimate optimal CSD consistent with the obtained 184 data. However, so far we have not discussed estimation of errors on the reconstruction. 185 Since the errors may be due to a number of factors -the procedure itself, measurement 186 noise, incorrect assumptions -one may approach this challenge in different ways.

187
First, to understand the effects of the selected basis sources and setup, one may consider 188 the estimation operator K(K+λ) −1 and the space of solutions it spans. This space is given 189 by the eigensources, introduced and described thoroughly in (Chintaluri et al., 2021). The 190 orthogonal complement of this space in the original estimation space, spanned byb j (x) 191 basis functions, is not accessible to the kCSD method. The study of eigensources facilitates 192 understanding which CSD features can be reconstructed and which are inaccessible. 193 Second, to consider the impact of the measurement noise on the reconstruction, for any 194 specific recording consider the following model-based procedure. Reconstruct CSD from 195 data with optimal parameters. Compute potential from estimated CSD. Add random 196 noise to each computed potential. The noise could be estimated from data, either as a 197 measure of fluctuations on a given electrode for a running signal, or from variability of 198 evoked potentials. Then, for any realization of noise, compute estimation of CSD. The 199 pool of estimated CSD gives estimation of the error at any given point where estimation 200 is made.

201
This computation can be much simplified by taking advantage of the linearity of the 202 resolvent, E = K(K+λI) −1 . Then, the i-th column (E i ) represents contribution of unitary 203 change of i-th measured potential (the i-th element of the vector V) to the estimated CSD 204 (C * ). As the contribution is proportional to the change, the column can be considered 205 an Error Propagation Map for i-th measurement (Fig. 3.A). Note that these vectors (the 206 columns of resolvent, E i ) also happen to form another basis of the solution space, an 207 alternative to the basis of eigensources.

208
If ε i is an error of i-th measurement, then its contribution to C * is ε i E i . Moreover, if 209 the measurement errors follow multivariate normal and the estimated CSD also follows multivariate normal The diagonal of EΣ V E T represents a map of CSD measurement uncertainty (uncertainty 213 attributed to the noise in the measurement, Fig. 3.B) 2 .  Figure 3: A) Error propagation maps for 3×3 regular grid of electrodes. Every panel represents the contribution of the potential measured at the corresponding electrode marked with a black circle (•) to the reconstructed CSD. Every other electrode is marked with a black cross (×). B) Map of CSD measurement uncertainty for 3 × 3 regular grid of electrodes. The CSD measurement uncertainty is represented by variance of the CSD reconstruction caused by the uncertainty in measurement of the potentials. It is assumed that measurement errors for electrodes are mutually independent and follow standard normal distribution (ε i ∼ N (0, 1)). Location of electrodes is marked with red crosses (×).
Third, one can study reconstruction accuracy for a meaningful family of test functions. 215 This could be the Fourier modes for rectangular regions or a collection of Gaussian test 216 functions, centered in different places, of single or multiple radii. For each of these test 217 functions one would compute the potential, perform reconstruction, and compare the 218 results with the original at every point. Finally, one could average this information over 219 multiple different test sources computing a single Reliability Map, which we now introduce. 220 Reliability maps Assume the standard kCSD setup, that is a region R ⊂ R n where we 221 want to estimate the sources, set of electrode positions, x i , and perhaps additional infor-222 mation, such as morphology for skCSD Cserpan et al., 2017. We now want to characterize 223 predictive power of the combination of our setup and our selected basis,b i . To do this we 224 select a family of test functions, C i (x), for example Gaussian test functions, centered in 225 different places, of multiple radii, or products of Fourier modes, etc. Then, for each C i we 226 compute V i = AC i by forward modeling, generating a surrogate dataset. Next, we apply 227 the standard kCSD reconstruction procedure obtaining estimation of the tested ground 228 truth,C i . We can then compute reconstruction error using point-wise modification of 229 Relative Difference Measure (RDM) proposed by (Meijs et al., 1988): where i = 1, 2, . . . enumerates different ground truth profiles. A simple measure of recon-231 struction accuracy is then given by the average over these profiles: Values on the map can be interpreted as follows: the closer to 0, the higher reconstruction accuracy might be achieved for a given measurement condition.
The class of functions used were the families of small and large sources mentioned 234 above. We used eight mirror symmetries of the grid in computation. 235 We can use reliability map as another source of information about the precision of 236 reconstruction, which is shown in Fig. 5. In A) we show some dipolar source which 237  Another interesting question is the effect of broken or missing electrodes on the recon-245 struction. Formally one can attempt kCSD reconstruction from a single signal but it is 246 naive to expect much insight this way. It is thus natural to ask what information can be 247 obtained from a given setup and what we lose when part of it becomes inaccessible.  In the present work we introduced a new Python package for kCSD estimation, kCSD-python. 261 This toolbox implements the basic functionality of source estimation with regularization. 262 We introduced here L-curve approach as an approach to parameter selection alternative to 263 cross-validation. This might be particularly important if regularization is attempted for 264 high-density setups, such as multi-electrode arrays with thousands of channels. We also 265 proposed several ways of diagnosing reconstruction errors, specifically Error Propagation 266 Maps and Reliability Maps, which can build intuition regarding veracity of reconstructions 267 obtained in given experimental contexts. Finally, we provided a brief tutorial to the new 268 package (see Supplementary Material) illustrating its main functions. All the figures in 269 this paper showing CSD and LFP were computed with this package and source files are 270 provided. In this section we discuss several issues related to CSD analysis in general and 271 kCSD in particular.

272
The benefits of CSD estimation over direct LFP analysis. Extracellular poten-273 tials provide valuable insight about the functioning of the brain. Thanks to recent advances 274 in multielectrode design and growing availability of sophisticated recording probes we can 275 monitor the electric fields in the nervous system at thousands of sites, often simultane-276 ously. One may wonder if this increased resolution makes CSD analysis unnecessary. In 277 our view, as we have discussed many times, it is not so. The long range nature of electric 278 potential means that even a single source contributes to every recording. Thus in princi-279 ple we should always expect strong correlation between neighboring sites. However, if the 280 separation between the electrodes becomes substantial, on the order of millimeters, the 281 level of correlation between the recordings on different electrodes will decrease. This is 282 because each electrode effectively picks up signals from a different composition of sources. 283 Even if some are shared they are overshadowed by others which may lead to small in-284 terelectrode correlations. Still, our experience shows that significant correlations can be 285 observed in the range of several millimeters (  In view of these facts we argue that it is always beneficial to attempt kCSD analysis. 292 The caveat is not to believe the reconstructed CSD blindly but always interpret it against 293 known anatomical and physiological knowledge supported by the tools such as provided 294 in the present work (eigensources, reliability maps, etc). This could be combined with 295 other decomposition methods if desired which may give more physiologically interpretable 296 results, in particular we have had very good results using independent component analysis, 297 which we recommend (Łęski, Kublik, et al., 2010;Głąbska et al., 2014).

298
Approaches to CSD estimation. Several procedures for CSD estimation were intro-299 duced over the years and are still in use. The first approach, which probably still dominates 300 today, was introduced by Walter Pitts in 1952 (Pitts, 1952) and gained popularity after 301 Nicholson and Freeman adopted it for laminar recordings (Nicholson and Freeman, 1975). 302 This was a direct numerical approximation to computation of the second derivative in the 303 Poisson equation (1). Only minor improvements were introduced over the years to stabilize 304 estimation (Rappelsberger, Pockberger, and Petsche, 1981) or handle boundaries (Vaknin,305 DiScenna, and Teyler, 1988). The first major conceptual change was introduced by Pet-306 tersen et al., 2006 who introduced model-based estimation of the sources. Their idea was 307 to assume a parametric model of sources, for example, spline interpolated CSD between 308 electrode positions, and using forward modeling to connect measured potentials to model 309 parameters. This model-based approach was generalized by Potworowski et al., 2012 who 310 proposed a non-parametric kernel Current Source Density method which is implemented 311 in the presented toolbox.

312
Apart from these main approaches several CSD variant methods were proposed. For 313 example, one may interpolate the potential first before applying traditional CSD approach, 314 or the opposite, interpolate traditionally estimated CSD. Although in some cases the 315 obtained results may look close to those obtained with kCSD, we do prefer kernel CSD 316 approach due to the underlying theory which facilitates computation of estimation errors 317 but also yields a unified framework for handling underlying assumptions, noisy data and 318 irregular electrode distributions. In our view approaches combining ad hoc interpolation 319 with numerical derivatives conceptually and computationally are less convincing to iCSD 320 and kCSD and we would not recommend them.

321
Models of tissue. Throughout this work and in the toolbox we assumed purely ohmic 322 character of the tissue. This has been debated in recent years (Bédard and Destexhe, 2011; 323 Riera et al., 2012; Gratiy et al., 2017) and it is true that more complex biophysical models 324 of the tissue, taking into account frequency dependent conductivity or diffusion currents, 325 would influence the practice of source reconstruction or its interpretation. However, the 326 available data indicate that in the range of frequencies of physiological interest these effects 327 are small. While one should keep eyes open on the new data as they become available 328 and keep in mind the different possible sources which may affect the reconstruction or 329 interpretation, we believe that the traditional view of ohmic tissue is an adequate basis for 330 typical experimental situations and going beyond that would probably require additional 331 dedicated measurement for the experiment at hand which may not always be feasible. For 332 example, as we discussed in (Ness et al., 2015), the specimen variability of the cortical 333 conductivity in the rat is much bigger than the variability between different directions 334 within a given rat (Goto et al., 2010). This means that unless we have conductivity 335 measurements for our specific rat we are probably making smaller error assuming isotropic 336 conductivity than taking different values from literature. We feel there is not enough data 337 to justify inclusion of more complex terms in the standard CSD analysis to be applied 338 throughout the brains and species.

339
In this manuscript and in the kCSD-python toolbox we also assumed constant conduc-340 tivity (with exception of MoI case below). We are convinced this is a reasonable approx-341 imation for typical depth recordings. In general, however, this approximation needs to 342 be justified or alternative models of tissue need to be considered. In principle, the kCSD 343 method can be applied for a variety of tissue models as long as the basis potentials can 344 be computed from the basis sources while incorporating the geometric and conductivity 345 changes.

346
For example, Ness et al., 2015 considered a cortical slice placed on a microelectrode 347 array (MEA) in which they included the geometry of the slice and modeled saline-slice 348 interface with changing conductivity in the forward model. They found that Method 349 of Images (MoI) gives a good approximation to the full solution obtained using finite-350 element model (FEM). This approximation was incorporated within the kCSD method as 351 MoIkCSD variant and is available in the kCSD-python package.

352
It is possible to generalize kCSD to reconstruct sources from recordings of multiple 353 electrical modalities -LFP, ECoG, EEG. In this case one needs to include the head 354 geometry and the changing tissue properties within the forward model and in the kCSD 355 method. The anisotropic (white matter tracts) and inhomogeneous (varying between skull, 356 cerebro-spinal fluid, gray matter and white matter) electrical conductivity changes can 357 be approximated using data obtained with imaging techniques such as MRI, CT or DTI. 358 Such sophisticated head models require numerical solutions such as finite element modeling 359 (FEM) to compute the basis potentials from the basis sources. We are currently working on 360 this approach to make it generic for any animal head and to eventually utilize it as a source 361 localization method for human data, for example, to localize foci of pharmacologically 362 intractable epilepsy seizures in humans. We call this extension kernel Electrical Source 363 Imaging (kESI). . Such massive high density data from thousands 371 of electrodes should greatly increase insight into the studied systems and significantly im-372 prove results of CSD reconstructions. There are two obstacles to fully benefit from kCSD 373 analysis of data from these new systems. First, kCSD involves inversion of the kernel 374 matrix which is quadratic in the number of electrodes. Combined with cross-validation 375 the necessary matrix operations quickly become overwhelming. This can be mitigated in a 376 number of ways, by subsampling the data, approximate inversions, and by switching from 377 cross-validation to L-curve method, but the challenge remains. This is the easy problem. 378 The difficult problem is physical. As we move away from a source its contribution to the 379 recorded potential goes down. In consequence, since the present version of kCSD uses all 380 recordings to estimate every source, when using remote signals to estimate local source, 381 we obtain mainly contributions from noise. In effect we get a very reliable estimation 382 of sources varying slowly in space but the sources changing fast in space are treated as 383 noise and silenced by the regularization. To take full advantage of these data a different 384 approach may give better results. One possibility could be a multiscale approach where 385 one would perform reconstructions in small windows in multiple scales to optimally recon-386 struct multiscale features of the source distribution. The challenge would be to efficiently 387 and correctly stitch them together.

388
Parameter selection. In the Results section we discussed our strategy for data-based 389 parameter selection using cross-validation or L-curve. Often, we need to tune not just λ 390 but also other parameters. For example, for Gaussian basis sources we may want to decide 391 on the width of the Gaussian used, R. To obtain the optimal set of parameters in that case 392 we compute the curvature of the L-curve or the cross-validation error for some ranges of 393 parameters considered and select parameters corresponding to the maximum curvature / 394 minimum error in the parameter space. This is a simplification of the proposition by Belge,395 Kilmer, and Miller, 2002 which in practice we found very effective.

396
As an example, in Fig. 7 we show results of such a scan for the problem shown in Figure 7: L-curve curvature (left) and CV-error (right) for the problem studied in Fig. 2. Observe that in both cases there are ranges of promising candidate parameter pairs, R, λ, which can give good reconstruction given the measured data. Red dots shows local extrema for each value of R fixed. See text for discussion of this effect. Fig. 2. The range of λ to be considered can be set by hand but by default we base it on 397 the eigenvalues of K. The smallest λ is set as the minimum eigenvalue of K which here 398 was around 1e-10. We set maximum λ at standard deviation of the eigenvalues, which 399 here was around 1e-3. The range of R values studied was from the minimum interelectrode 400 distance to half the maximum interelectrode distance. Note that for very inhomogeneous 401 distributions of electrodes this approach may be inadequate.

402
What we find is that apart from a global minimum in R, λ space there is a range of 403 R values fixing which we can find optimal λ(R) which leads to very close curvatures / 404 CV-errors / estimation results. What happens is that within some limits we may achieve 405 similar smoothing effects changing either λ or R. Bigger λ means more smoothing, but 406 bigger R means broader basis functions and effectively also smoother reconstruction space. 407 This is why the CV-error and curvature landscapes are relatively flat, or have these marked 408 valleys observed in Fig. 7. This effect supports robustness of the kCSD approach.

409
The sources of error in kCSD estimation and how to deal with them. Kernel 410 CSD method assumes a set of electrode positions and a corresponding set of recordings. 411 Additionally, single cell kCSD requires morphology of the cell which contributed to the 412 recordings and its position relative to the electrodes. Each of these may be subject to 413 errors.

414
In analysis, we assume that the electrode positions are known precisely. This is a 415 justified assumption in case of multishaft silicon probes or integrated CMOS-MEA but 416 not necessarily when multiple laminar probes are placed independently within the brain 417 or for many other scenarios. For example, for SEEG electrodes in human patients used in 418 presurgical evaluation we expect the localization errors due to the workflows used clinically 419 to be significant. We do not provide dedicated tools to study the effects of misplaced elec-420 trodes on the reconstructed CSD, however, this can be achieved easily with the provided 421 package if needed. The location of the cell with respect to the electrodes is much more 422 questionable, especially in 3D. Nevertheless, the necessary data to perform skCSD today 423 are too scarce to start addressing these issues.

424
On the other hand we do assume that the recordings are noisy and we use regulariza-425 tion to counteract the effects of noise. We have no mechanism to differentiate between 426 electrodes with varying degrees of noise to compensate this differently. However, we ob-427 served that for cases with very bad electrodes, similar results are obtained for analysis of 428 complete data and for analysis of partial data with bad electrodes removed from analysis. 429 The difference was in λ selected which was larger when broken electrodes were included 430 in the analysis. Depending on situation, if there is a big difference in the noise visible in 431 different channels, an optimal strategy may be to discard the noisy data and perform re-432 construction from the good channels only, which kCSD permits. In the end, data analysis 433 remains an art and a healthy dose of common sense is always recommended.

434
The main limitation of the method itself lies in the character of any inverse problem. 435 Here it means that there is an infinite number of possible CSD distributions each consistent 436 with the recorded potential. It is thus necessary to impose conditions which allow unique 437 reconstruction and this is what every variant of CSD method is about. In kCSD this 438 condition is minimization of regularized prediction error. In practical terms one may 439 think of the function space in which we are making the reconstruction. This space is span 440 by the eigensources we discussed before (Chintaluri et al., 2021). We feel it is useful to 441 consider both this space as well as its complement, that is the set of CSD functions whose 442 contribution to every potential is zero. This can facilitate understanding of which features 443 of the underlying sources can be recovered and which are inaccessible to the given setup. 444 While for the most common regular setups, such as rectangular or hexagonal MEA grids 445 or multishaft probes, intuitions from Fourier analysis largely carry over, in less regular 446 cases this quickly becomes non-obvious. 447 To facilitate intuition building in the provided toolbox we include tools to compute 448 the eigensources for a given setup. We also proposed here reliability maps, heuristic tools 449 to build intuition regarding which parts of the reconstructed CSD can be trusted and 450 which seem doubtful. These reliability maps are built around specific test ground truth 451 distributions and some default parameters facilitating validation for any given setup are 452 provided. Due to the open source nature of the provided toolbox more complex analysis 453 is possible if the setup or experimental context require that.

454
Appendix: kCSD-python package tutorial 618 In this section we first illustrate the use of kCSD package for CSD reconstruction in the 619 simplest case of a regular 2D square grid. This is a simplified version of a slice on a 620 microelectrode array (Ness et al., 2015), or a planar silicone probe within the brain, where 621 we assume constant conductivity in the whole space. In the following sections we show 622 how we validate our methods and what kind of diagnostics we find useful in the analysis 623 of experimental data. This tutorial is available as a jupyter notebook and can also be 624 accessed through a web-browser without installation. For more details, see 'Overview of 625 kCSD-python package' in the Discussion.

626
Basic features 627 We start with the basic CSD estimation on a regular grid. First, we define a region of 628 interest. Then, using predefined test functions for the current sources, we place a ground 629 truth current source in this region. We define the distribution of electrodes. Assuming 630 ideal electrodes, we compute the potential generated by the selected current sources as 631 measured at the electrodes. Given these potentials and the electrode locations we estimate 632 the current source density using kCSD. As a final step, we perform cross-validation to avoid 633 overfitting. Since we know the ground truth used to generate the potentials that were used 634 in the kCSD estimation, we can compare the ground truth to the estimate and see the 635 reconstruction accuracy. We define the region of interest between 0 and 1 in the xy plane with a resolution of 101 638 points in each dimension. We will assume the distance is given in mm, so we want to 639 perform a reconstruction on a square patch of 1mm 2 size.

640
Setting up the ground truth The kCSD-python library provides functions to generate 641 test sources which can be imported from the csd_profile module. Here we use the 642 gauss_2d_small function to generate two-dimensional Gaussian sources which are small 643 in the scale set by the interelectrode distance. The other implemented option for two-644 dimensional test sources is the gauss_2d_large function. To generate the exact same 645 sources in each run we must invoke this function using the same random seed which is 646 stored in the seed variable. For simplicity, these current sources are static and do not 647 change with time. We visualize the current sources as a heatmap.

648
In [2]: from kcsd import csd_profile as CSD CSD_PROFILE = CSD.gauss_2d_small true_csd = CSD_PROFILE(csd_at, seed=15) The code below displays this test source as the True CSD. For convenience we define this 649 as a function make_plot. The output for this code is shown in Fig. 8A.  figure(figsize=(7, 7)) ax = plt.subplot(111) ax.set_aspect('equal') t_max = np.max(np.abs(zz)) levels = np.linspace(-1 * t_max, t_max, 32) im = ax.contourf(xx, yy, zz, levels=levels, cmap=cmap) ax.set_xlabel('X (mm)') ax.set_ylabel('Y (mm)') ax.set_title(title) ticks = np.linspace(-1 * t_max, t_max, 3, endpoint=True) plt.colorbar(im, orientation='horizontal', format='%.2f', ticks=ticks) return ax make_plot(csd_x, csd_y, true_csd, title='True CSD', cmap=cm.bwr)  To visualize the potential, we interpolate the hundred values computed at the electrodes 663 positions with interpolate.griddata function. Note that the kCSD estimation uses only 664 the potential recorded at the electrode positions. To distinguish between the potentials 665 and CSD plots we use different colormaps. The electrodes are marked with dots in this 666 plot. The output from this step is shown in Fig. 8B. Estimated current sources are shown in Fig. 8C. Compare this to the True CSD obtained 681 before, Fig. 8A. Observe that the estimation is not very faithful. This is caused by the 682 ground truth varying significantly in the scale of a single inter-electrode distance. In the 683 next step we will use cross-validation to select better reconstruction parameters. Cross validation Leave-one-out cross-validation is performed with a single line com-685 mand. In this procedure we scan a range of R values which set the size of the Gaussian 686 basis functions and the regularization parameter λ values. At the end of this step we 687 obtain the optimal parameters that would correct for overfitting. The function outputs 688 the progress of the cross-validation step and displays the optimal candidates in the last 689 line. Alternatively, one could use the L-curve method to find these optimal parameters. 690 Fig. 8D shows the kCSD reconstruction obtained after cross-validation. We find that this 691 estimation of the current sources resembles the True CSD better. Until now we assumed noise-free data, however, experimental data are always noisy. In 694 this section we investigate how noise affects the kCSD estimation. We first show how 695 to compute the reliability map which we introduced before, Eq. (21). Then we discuss 696 reproducible generation of noisy data with varying noise amplitude. Finally, we study the 697 error in the reconstruction as a function of changing noise level.

698
Reconstruction quality measure To assess the estimation quality we measure the 699 point-wise difference between the true sources and the sources reconstructed with the kcsd. 700 We define a function point_errors which takes the true_csd and the estimated_csd 701 as the inputs, normalizes them individually, and computes the Frobenius norm of their 702 difference.
Noise definition To study resilience of the reconstruction against the noise in a con-707 trolled way we seed the random number generator in the function add_noise. We consider 708 normally distributed noise with the mean and standard deviation set by reference to the 709 recorded potentials.

710
In [13]: def add_noise(pots, noise_level=0, noise_seed=23): rstate = np.random.RandomState(noise_seed) noise = noise_level*0.01*rstate.normal(np.mean(pots), np.std(pots), size=(len(pots), 1)) return pots + noise.reshape(pots.shape) pots_noise = add_noise(pots, noise_level=15, noise_seed=23) Source reconstruction from noisy data With these tools we can study the effects 711 of noise on the reconstruction. We now generate noise for a given noise level between 0 712 and 100, add it to the simulated potential, and estimate CSD from these noisy potentials. 713 We can then use the error plots to compare the reconstruction with the True CSD. Notice 714 that the parameters giving best reconstruction obtained for noisy data in general will be 715 different from those obtained for clean potentials to compensate for noise. We can display this error with the make_error_plot plotting function which we defined 717 earlier. Changing the noise_level and the noise_seed affects the reconstruction, but the 718 error depends also on the sources, so changing the True CSD type to a gauss_2d_large or 719 changing csd_seed will lead to different results. This is illustrated in Fig. 9A-D for small 720 Gaussian sources, and Fig. 9E-H for large Gaussian sources, with varying noise levels. 721 The actual ground truth and reconstructions are shown in Fig. 8 It is often the case that due to experimental constraints, some subset of recordings are dis-724 carded in the final analysis. This can happen when some electrodes are used for stimulation 725 and cannot be used for recording, or due to bandwidth limitations requiring a compromise 726 between sampling rates and the number of simultaneously recording electrodes, or in the 727 event of an electrode break down. In this section of the tutorial we discuss how to handle 728 such cases and to estimate the errors in reconstruction despite the loss of recordings. We 729 first show how we remove recordings from selected (broken) electrodes from considered 730 data. Then we calculate the estimation error for a given source for data from such a 731 damaged setup. Finally, we compute the average error across many sources from such an 732 incomplete setup. Note that kCSD reconstruction is designed to work with arbitrary elec-733 trode setups and removing specific electrodes does not change the situation significantly. 734 We focus on broken electrodes as it is a common situation in practice and deserves consid-735 eration. This may be used to gain intuition regarding ways in which CSD reconstruction 736 may go wrong, due to slight disturbances in a familiar setup.

737
Remove broken electrodes To test the effects of removed electrodes on reconstruction 738 from a given setup we simulate this with a function remove_electrodes that takes all the 739 electrode positions for this setup and the number of electrodes that are to be removed. 740 In this example we remove the electrodes randomly. Like we did previously, to facilitate 741 repeatability we also pass a broken_seed variable, so that at each subsequent run the 742 same electrodes are discarded. By changing this seed we select a different set of electrodes 743 for removal. Error in estimation with broken electrodes After removing the broken electrodes 745 we compute the estimation error to gauge the effect of electrode removal on reconstruction. 746 Here, a fuction calculate_error takes a csd_seed as an input, which selects a specific 747 ground truth source, and all the remaining electrode positions, ele_pos. The function 748 computes the True CSD for a gauss_2d_small type source, computes the potential at 749 these electrode locations, performes kcsd estimation from these data, and computes the 750 error in the estimation of the true csd.   Figure 11: Error propagation maps for 1D regular grid of 12 electrodes Every panel represents the CSD contribution (red line) of the potential measured at the corresponding electrode, for which the potential is 1 (green line).

768
An example of 3D source reconstruction.
769 Figure 12: An example of 3D kCSD source reconstruction. Each column shows five consecutive parallel cuts through a box of size 1. A) Ground truth for the CSD seed of 16. B) Estimated potential; black dots indicate electrodes where potential is collected for further reconstruction. C) 3D kCSD reconstruction from the measured potentials, λ = 0. D) 3D kCSD reconstruction with cross-validation.