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

Interpretation of extracellular recordings can be challenging due to the long range of electric field. This challenge can be mitigated by estimating the current source density (CSD). Here we introduce kCSD-python, an open Python package implementing Kernel Current Source Density (kCSD) method and related tools to facilitate CSD analysis of experimental data and the interpretation of results. We show how to counter the limitations imposed by noise and assumptions in the method itself. kCSD-python allows CSD estimation for an 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 which illustrates a typical analytical workflow and main functionalities useful in validating analysis results.

Setting up the ground truth The kCSD-python library provides functions to generate test sources which can be imported from the csd_profile module.Here we use the gauss_2d_small function to generate two-dimensional Gaussian sources which are small in the scale set by the interelectrode distance.The other implemented option for two-dimensional test sources is the gauss_2d_large function.To generate the exact same sources in each run we must invoke this function using the same random seed which is stored in the seed variable.For simplicity, these current sources are static and do not change with time.We visualize the current sources as a heatmap.
In [2]: from kcsd import csd_profile as CSD CSD_PROFILE = CSD.gauss_2d_smalltrue_csd = CSD_PROFILE(csd_at, seed=15) The code below displays this test source as the True CSD.For convenience we define this as a function make_plot.The output for this code is shown in Fig 1A.
In Place electrodes We now define the virtual electrodes within the region of interest.We place them between 0.05 mm and 0.95 mm of the region of interest, with a resolution of 10 (as indicated by 10j in mgrid) in each dimension, totalling 100 electrodes.Notice that the electrodes do not span the entire region of interest.Although in this example the electrodes are distributed on a regular grid, this is not required by the kCSD method as it can handle arbitrary distributions of electrodes.
In [4]: ele_x, ele_y = np.mgrid[0.05:0.95: 10j, 0.05: 0.95: 10j] ele_pos = np.vstack((ele_x.flatten(),ele_y.flatten())).T Compute potential To obtain the potential, pots, at the given electrode positions due to the current sources that were placed in the previous steps we use the function forward_method.We assume the sources are localized within a slab of tissue of thickness 2h on top the MEA (See [1,2] and Methods).We also assume an infinite homogeneous medium of conductivity sigma equal to 1 S/m.Finally, we assume that the electrodes are ideal, point-size, and noise-free.kCSD method Here we illustrate the most basic estimation of CSD with the kcsd library.Since our example is two dimensional the relevant method is KCSD2D.For convenience, we encapsulate the actual method call with parameters being set inside a function do_kcsd.We first set the h and sigma parameters of the forward model.Then we restrict the potentials to the first time point of the recording.For typical experimental data the shape of this matrix would be N ele × N time , where N ele is the number of electrodes and N time is the total number of recorded time points.Next, we call the KCSD2D class with the relevant parameters.The only required parameters are the electrode positions, ele_pos, and the potentials they see, pots.We can also provide here the parameters for the forward model, h and sigma.We define a rectangular region of estimation by setting the values xmin, xmax and ymin, ymax.The number of basis functions, n_src_init is set to 1000, basis functions are of the type gauss, and the width of the Gaussian basis source R_init is set to be 1.Finally, the estimated CSD is stored as est_csd.Observe that the estimation is not very faithful.This is caused by the ground truth varying significantly in the scale of a single inter-electrode distance.In the next step we will use cross-validation to select better reconstruction parameters.
In [8]: make_plot(k.estm_x,k.estm_y, est_csd[:, :, 0], # First time point title='Estimated CSD without CV', cmap=cm.bwr) Cross-validation Leave-one-out cross-validation is performed with a single-line command.In this procedure we scan a range of R values which set the size of the Gaussian basis functions and the regularization parameter λ values.At the end of this step we obtain the optimal parameters that would correct for overfitting.The function outputs the progress of the cross-validation step and displays the optimal candidates in the last line.Alternatively, one could use the L-curve method to find these optimal parameters.Fig 1D shows the kCSD reconstruction obtained after cross-validation.We find that this estimation of the current sources resembles the True CSD better. In

Noisy electrodes
Until now we assumed noise-free data, however, experimental data are always noisy.In this section we investigate how noise affects the kCSD estimation.We first show how to compute the reliability map which we introduced before, Eq (26).Then we discuss reproducible generation of noisy data with varying noise amplitude.Finally, we study the error in the reconstruction as a function of changing noise levels.
Reconstruction quality measure To assess the estimation quality we measure the point-wise difference between the true sources and the sources reconstructed with the kcsd.We define a function point_errors which takes the true_csd and the estimated_csd as the inputs, normalizes them individually, and computes the Frobenius norm of their difference.
In We visualize this difference as before, except we use a greyscale colormap to display the intensity of the reconstruction error.For convenience, we define the plotting in a function called make_error_plot.The output from this step is shown in Fig 2A.A) The error between the True CSD and the estimation obtained with kCSD for a 2 dimensional small Gaussian current source, using the csd seed of 15.The electrodes in this case are assumed to be noise-free.B, C, D) Same as A, however, noise is added to the recorded potentials, whose magnitude is 5%, 10%, or 30%, respectively.E-H) Analogous to A-D, except in this case large Gaussian sources with seed 6 were used.
Noise definition To study resilience of the reconstruction against the noise in a controlled way we seed the random number generator in the function add_noise.We consider normally distributed noise with the mean and standard deviation set by reference to the recorded potentials.
In Source reconstruction from noisy data With these tools we can study the effects of noise on the reconstruction.We now generate noise for a given noise level between 0 and 100, add it to the simulated potential, and estimate CSD from these noisy potentials.We can then use the error plots to compare the reconstruction with the True CSD.Notice that the parameters giving best reconstruction obtained for noisy data in general will be different from those obtained for clean potentials to compensate for noise.In [15]: make_error_plot(k_noise.estm_x, k_noise.estm_y,error_noise, title='Error CSD, with noise')

Broken electrodes
It is often the case that due to experimental constraints, some subset of recordings are discarded in the final analysis.This can happen when some electrodes are used for stimulation and cannot be used for recording, or due to bandwidth limitations requiring a compromise between sampling rates and the number of simultaneously recording electrodes, or in the event of an electrode break down.In this section of the tutorial we discuss how to handle such cases and to estimate the errors in reconstruction despite the loss of recordings.We first show how we remove recordings from selected (broken) electrodes from considered data.Then we calculate the estimation error for a given source for data from such a damaged setup.Finally, we compute the average error across many sources from such an incomplete setup.Note that kCSD reconstruction is designed to work with arbitrary electrode setups and removing specific electrodes does not change the situation significantly.We focus on broken electrodes as it is a common situation in practice and deserves consideration.This may be used to gain intuition February 26, 2024 7/10 regarding ways in which CSD reconstruction may go wrong, due to slight disturbances in a familiar setup.
Remove broken electrodes To test the effects of removed electrodes on reconstruction from a given setup we simulate this with a function remove_electrodes that takes all the electrode positions for this setup and the number of electrodes that are to be removed.In this example we remove the electrodes randomly.As we did previously, to facilitate repeatability we also pass a broken_seed variable, so that at each subsequent run the same electrodes are discarded.By changing this seed we select a different set of electrodes for removal.
In Error in estimation with broken electrodes After removing the broken electrodes we compute the estimation error to gauge the effect of electrode removal on reconstruction.Here, a fuction calculate_error takes a csd_seed as an input, which selects a specific ground truth source, and all the remaining electrode positions, ele_pos.The function computes the True CSD for a gauss_2d_small type source, computes the potential at these electrode locations, performs kcsd estimation from these data, and computes the error in the estimation of the true csd.

Fig 1 .
Fig 1. Basic features tutorial A) Shows the ground truth (True CSD), here, two-dimensional small Gaussian current sources, for the CSD seed of 15.B) The interpolated potentials generated by this current source are shown, the electrodes are displayed as black dots.C) CSD estimated with kCSD using the potentials from the electrode positions, without cross-validation.D) Same as C but cross-validation was used.E-H) Analogous to A-D, except large Gaussian current sources for seed 6 were used.

Fig 2 .
Fig 2. Noisy electrodes.A) The error between the True CSD and the estimation obtained with kCSD for a 2 dimensional small Gaussian current source, using the csd seed of 15.The electrodes in this case are assumed to be noise-free.B, C, D) Same as A, however, noise is added to the recorded potentials, whose magnitude is 5%, 10%, or 30%, respectively.E-H) Analogous to A-D, except in this case large Gaussian sources with seed 6 were used.
In [14]: k_noise = do_kcsd(ele_pos, pots_noise) k_noise.cross_validate(Rs=np.linspace(0.01,0.15, 15)) estm_csd_noise = k_noise.values('CSD')error_noise = point_errors(true_csd, estm_csd_noise) No lambda given, using defaults Cross validating R (all lambda) : 0.01 Cross validating R (all lambda) : 0.02 ... Cross validating R (all lambda) : 0.15 R, lambda : 0.01 0.00110069417125 We can display this error with the make_error_plot plotting function which we defined earlier.Changing the noise_level and the noise_seed affects the reconstruction, but the error depends also on the sources, so changing the True CSD type to a gauss_2d_large or changing csd_seed will lead to different results.This is illustrated in Fig 2A-D for small Gaussian sources, and Fig 2E-H for large Gaussian sources, with varying noise levels.The actual ground truth and reconstructions are shown in Fig 1.
Fig 3. Broken electrodes.A) Shows the average error between the True CSD and the CSD estimated with kcsd for 100 random small Gaussian current sources.B) The same average error as in A, except in this case 5 electrodes were discarded in the estimation.Likewise for C and D, where 10 and 20 electrodes out of the 100 were considered broken.E-H) analogous to A-D, except in this case we show the averages for 100 large Gaussian current sources.