µSpikeHunter: An advanced computational tool for the analysis of neuronal communication and action potential propagation in microfluidic platforms

Understanding neuronal communication is fundamental in neuroscience, but there are few methodologies offering detailed analysis for well-controlled conditions. By interfacing microElectrode arrays with microFluidics (μEF devices), it is possible to compartmentalize neuronal cultures with a specified alignment of axons and microelectrodes. This setup allows the extracellular recording of spike propagation with a high signal-to-noise ratio over the course of several weeks. Addressing these μEF devices, we developed an advanced yet easy-to-use publically available computational tool, μSpikeHunter, which provides a detailed quantification of several communication-related properties such as propagation velocity, conduction failure, spike timings, and coding mechanisms. The combination of μEF devices and μSpikeHunter can be used in the context of standard neuronal cultures or with co-culture configurations where, for example, communication between sensory neurons and other cell types is monitored and assessed. The ability to analyze axonal signals (in a user-friendly, time-efficient, high-throughput manner) opens the door to new approaches in studies of peripheral innervation, neural coding, and neuroregeneration, among many others. We demonstrate the use of μSpikeHunter in dorsal root ganglion neurons where we analyze the presence of both anterograde and retrograde signals in μEF devices. A fully functional version of µSpikeHunter is publically available for download from https://github.com/uSpikeHunter.


Introduction
Electrical signaling is recognized as the principal modality of communication in neurons, where it is used to encode and transmit information via action potentials (APs). Electrophysiology recordings play therefore a fundamental role in understanding neuronal circuits in physiological and pathological conditions. Multielectrode arrays (MEAs) standout among the different methodologies with the appropriate spatial and temporal scales to assess neuronal circuits in well-controlled in vitro settings 1,2 . However, it is difficult to detect and characterize AP propagation (e.g., direction and velocity) using conventional MEAs in neuronal cultures: it is virtually impossible to ensure that electrodes are well placed for APs detection, source-target information is inaccessible, and the amplitude of the APs recorded from axons are typically very low and impossible to discriminate. To improve the signal-to-noise ratio (SNR) of electrophysiological recordings and the localization of neuronal processes on recording electrodes, devices combining microElectrodes and microFluidics (µEF devices) have been developed for neuroscience applications 3 . A µEF device is composed of a microfluidic device mounted on a MEA to form an enclosed culture chamber composed of two (or more) isolated compartments connected by microchannels. The reduced dimensions of the microchannels ensure that somata are excluded from these channels whereas axons are able to grow through.
Among other applications, this type of device allows for the spatially compartmentalized but functionally connected cultures of separate populations of cells. With the microfluidic device aligned such that the microchannels are positioned above multiple microelectrodes, signal propagation along the axons can be observed and analyzed. The small dimensions of the microchannels provide the added benefit of increasing the SNR of recorded axonal APs 4 .
A number of studies have been conducted using µEF devices to assess i) the directionality of communication and the origin of bursting behavior in networks of distinct populations of neurons [5][6][7][8][9] , ii) changes to the propagation velocity with culture age 10,11 , or iii) results of pharmacological, biochemical, or electrical stimulation 10,12,13 .
Compartmentalized microfluidic devices have also been used to investigate the interaction between neurons and other cells co-cultured in the separate compartments. But although studies on neuronal circuit dynamics and signal communication are intimately related with electrophysiology, the vast majority of these microfluidic studies on neuronal co-cultures with other types of cells have lacked an electrophysiological facet, which would serve to complement the fundamentally biochemical results and further elucidate the interaction between the cell types. For example, past studies have assessed various factors affecting myelination in the central nervous system [14][15][16] . Additionally, the biochemical and morphological facets of the interaction between neurons and cells from other organ systems, including osteoblasts 17 , dental pulp 18 , and myocytes 19,20 have been investigated.
In both single-and co-cultures configurations, more than the technical difficulty of combining microelectrodes and microfluidics, the challenge for most labs lies on the sheer volume and complexity of the recorded electrophysiology data. The electrophysiology data obtained from neurons cultured in µEF devices is unquestionably highly informative, allowing for example for the detection of propagating APs, calculation of their propagation direction and velocity as well as their spike times or inter-spikes interval (relevant for neuronal coding/decoding analysis). However, although the necessary experimental protocols and tools for combining MEAs and microfluidics are readily available [JoVE REF], there is a lack of user-friendly tools available to analyze recordings obtained using µEF devices (but see 21,22 , for visualization and spike sorting tools).
In this study, an advanced yet user-friendly data analysis program called µSpikeHunter was developed for the detection and characterization of APs propagating along axons confined to the microchannels of a µEF device. With this free, open-source software, the presence of traveling waves (APs) in electrophysiological recordings can be readily detected, visualized and characterized. Spike sorting can also be performed in µSpikeHunter based on the waveforms of recorded APs. In this study, µSpikeHunter was validated in data analysis of electrophysiology experiments using cortical neurons and dorsal root ganglia (DRGs), and was used to uncover the otherwise elusive presence of APs traveling towards the somata compartment in common microfluidic devices.

Results µSpikeHunter graphical user interfaces
The computational tool µSpikeHunter runs on both Microsoft Windows and Apple's macOS operating systems, and is composed of two graphical user interfaces (GUIs): the main GUI, in which the data are imported and analyzed at the single-spike level, and the spike sorting GUI, in which the user can sort spikes into source clusters associated to different neurites.
µSpikeHunter was developed to be compatible with recordings obtained using custom setups or commercial recording systems (MEA2100 from MultiChannel Systems) for 60-, 120-, and 252-electrode MEAs. Analysis can be performed for signals recorded by a series of up to 16 electrodes with a uniform inter-electrode spacing. Details on how to use µSpikeHunter can be found in the detailed user manual in the supplementary materials.
The main GUI ( Fig. 1) allows the user to select a file for analysis and, if a commercial MultiChannel Systems setup was used for the recording, the electrodes and time range for analysis. Once the desired data have been selected and the parameters set for event detection, the user is presented with a list of detected propagation sequences. These sequences are series of events detected along the electrodes which satisfy specific criteria (see Materials and Methods) to be considered an AP traveling (in an axon) through the microchannel containing the analyzed electrodes. The propagation sequences can be analyzed/characterized and the results exported to a data file.
The spike sorting GUI (Fig. 2) presents the user with a spike overlay for each of the electrodes consisting of all the detected propagation sequences, and allows the user to sort the spikes into source clusters, which are collections of propagation sequences that are presumed to have arisen from different axons inside the same microchannel. The spike overlay plots presented to the user are aligned about the peak voltage values for each detected event.

Algorithms performance assessment using synthetic data
The µSpikeHunter algorithms for detection and characterization of propagating APs were assessed/validated using synthetic data simulating electrophysiological recordings of microelectrodes along a microchannel, with different levels of SNR.

Propagation sequence detection
The propagating wave detector was evaluated based on two performance indices: the precision rate and the detection rate. The precision is defined as where is the number of true positives, i.e., the number of detected sequences that correspond to actual traveling waves, and is the number of false positives, i.e., the number of detected sequences that do not correspond to actual traveling waves. The detection rate is defined as = where is the total number of actual sequences in the recording dataset. Both of these performance indices range from 0 to 1.
A high precision indicates that the propagation sequence detector rarely yields temporally linked events that do not correspond to actual traveling waves (APs). A high detection rate indicates that the propagation sequence detector is able to recognize a high percentage of the actual traveling waves in the data. It should be noted that in spike characterization, a high precision is more important than a high detection rate. That is, it is more desirable to be certain that the analyzed spikes are actually traveling APs than to detect the majority of spikes in a recording (due to noise and detection constrains, the same spike may not be observed in all microelectrodes along a microchannel).
The precision and detection rate results are presented in Fig where it is still very close to 1 (0.96).
The detection rate drops steadily as the noise increases from a SNR of 0.7 to 0.2. Whereas the propagation sequence detector is able to detect 83% of all actual propagation sequences at a SNR of 0.7, for a SNR of 0.4 less than 10% of all propagation sequences are detected.

Propagation velocity estimation
The performance of the propagation velocity estimation of µSpikeHunter was assessed using three measures (see Materials and Methods):  cluster propagation velocity (CPV), calculated from the relative timing and shape of the voltage waveforms, measured between the two most distant electrodes;  single-sequence propagation velocity (SPV), calculated directly from the crosscorrelation between the voltage waveforms in a particular pair of electrodes;  mean SPV, calculated using the average SPV for all electrode pairs.
A comparison of the methods is shown in Fig. 3c

Detailed characterization of axonal signal propagation in controlled in vitro settings
Two sets of experiments, using microfluidic platforms with a somal and an axonal compartment separated by microchannels, were used to demonstrate the data analysis capabilities of µSpikeHunter.

Spike sorting
The spike sorting performance of µSpikeHunter was evaluated using recordings obtained from Interestingly, in addition to the distinct appearance of the waveforms, the propagation velocities of the sorted clusters can provide further confirmation that the spikes arise from different sources. For example, in Fig. 4a2 and b2 clusters 1 and 2 in channel B show significantly different CPVs with low standard deviations corresponding to errors of approximately 5% and high confidence indices owing to the highly consistent waveforms in each cluster. The same high reliability can also be observed in the CPV estimate for cluster 1 in channel M, which has an error of approximately 8% at one standard deviation. However, the standard deviation of the CPV estimates were higher for cluster 0 in both channels B and M; in these cases, the error exceeded 20%. This is likely because the SNRs of the spikes in these clusters were lower and the propagation velocity estimates were therefore more affected by noise. Alternatively, it may also be the case that these clusters actually represent more than one source, but the waveforms are indistinguishable, again because of the SNR. where DRG integrates afferent signals).

Monitoring culture activity over multiple days in vitro
As a final benchmark for the analysis capabilities µSpikeHunter a DRG explant culture was monitored in a µEF device as it aged, using 10 minutes recordings at DIVs 4, 6, and 8.

Preparation of microelectrode-microfluidic devices
The µEF devices were prepared by placing standard microfluidic devices with an appropriate microgrooves spacing against a pre-coated MEA ( Tissue fragments were then mechanically dissociated with a 5 ml plastic pipette and subsequently with 1 ml pipette tips. Viable cells were counted using the trypan blue (0.4% (w/v), Sigma-Aldrich Co.) exclusion assay, and the cell density was adjusted to 2 × 10 7 viable cells/ml. Thereafter, 5 μl of the cell suspension was seeded in the cell body compartment of a µEF device, previously treated with 0.01 mg/ml PDL as described above. Cells were cultured in Neurobasal medium supplemented with 0.5 mM glutamine, 2% B27, and 1% penicillin/streptomycin (P/S, 10,000 units/ml penicillin and 10,000 μg/ml streptomycin) and kept in a humidified incubator at 37 °C supplied with 5% CO2.
Primary embryonic mouse dorsal root ganglion (DRG) explants were isolated from wild-type C57BL/6 mice embryos (E16.5). Lumbar DRG explants were removed and placed in HBSS until use. Upon use, one DRG explant was seeded in the cell body compartment of a µEF device, previously treated with 0.01 mg/ml PDL plus laminin as described above. Cells were cultured in Neurobasal medium supplemented with 0.5 mM glutamine, 2% B27, 50 ng/ml of NGF 7S (Calbiochem®, Millipore), and 1% P/S, and kept in a humidified incubator at 37 °C supplied with 5% CO2. respectively. The electrodes are referred to using the electrode labels defined by MultiChannel Systems, with a number representing the row and a letter representing the column, and the microchannels are referred to using the letter representing the electrode column over which they are aligned (see Fig. 1, low left corner of the GUI, for a schematic example of this labeling system). Recordings were obtained at a sampling rate of 20 kHz, and all recorded activity was spontaneous activity.

Computational algorithms
µSpikeHunter was developed in MATLAB R2016b (The MathWorks Inc., USA). The graphical user interfaces were developed using MATLAB's GUI development environment (GUIDE).
With the GUIDE Layout Editor, GUI layouts were developed for compatibility with both Windows and Macintosh operating systems. For an event detected on the central th electrode to be considered part of a propagation sequence, all of the following three conditions must be met.

Propagation sequence detection
(1) Temporally linked events are detected on all other electrodes.
(2) The times of the linked events on the first and last electrodes correspond to a propagation velocity of less than 100 m/s. Here, , and , are the voltage traces of the th propagation sequence on the th and th electrodes, respectively; is the lag time between the voltage traces; and ⊗ represents the cross-correlation, which is calculated as a function of . Additionally, the condition j > i is defined so that the propagation velocity estimate described below yields positive and negative velocities for anterograde and retrograde propagation, respectively. The time window for this cross-correlation is defined as , ± , , where the time window constant is 7.5 s/m.
A change of variable is also performed for the lag timescale of each cross-correlation so that the cross-correlation is given as a function of the inverse of the velocity, as However, it should be noted that the time resolution affects the accuracy of the SPV estimates in a nonlinear way and this is not taken into consideration when the average is computed. That is, when the lag is smaller, as for more closely located electrodes, small errors in the lag produce larger errors in the SPV. Thus, the SPV is less error prone when a more distant pair of electrodes are selected. Additionally, regardless of the inter-electrode distance, an underestimation of the absolute value of the lag produces a larger error than an overestimation of the same magnitude, though the difference between the errors is larger for closer electrodes as a result of the −1 dependence. Because of this, the SPV tends to overestimate the propagation velocity, especially when a closer electrode pair or the average of all pairs is selected for the estimation.
The SPV is paired with a confidence index, which indicates the similarity of the spike shapes recorded on the two electrodes. When a single pair is selected for SPV calculation, the confidence index S is the peak value of the normalized cross-correlation: When the average of all pairs is selected, the confidence index is the peak value that is lowest among the electrode pairs.

Kymograph and audio playback
µSpikeHunter also contains two interactive elements for the detection of traveling signals: a kymograph and an audio playback function. The kymograph is an image representation of the voltage signals recorded on each electrode. The voltages are converted to pixel color map intensities and plotted in time-electrode space (see Fig. 1). This representation allows the user to readily visually assess whether there is propagation, determine the direction and speed of propagation, and observe the relative peak voltage magnitudes on each electrode. The user may draw a line on the kymograph to manually calculate the propagation velocity.
The audio playback function assigns a different tonal frequency (note) to each electrode selected for analysis and converts the voltage signals on each electrode to audio intensities. The time dilation for the playback is 500 times to allow the spikes recorded on the different electrodes to be distinguishable by the human ear. With this playback function, traveling waves are detectable as a sequence of ascending or descending frequencies for anterograde or retrograde propagation, respectively.
Together with the detailed quantitative methods, these two functions allow the users to qualitatively detect traveling waves in two different modalities: visual and auditory.

Spike sorting
Spike sorting is performed based on regions of interest (ROIs) drawn by the user for up to four source clusters (clusters 1-4) with a fifth cluster (cluster 0) comprising any spikes not sorted into clusters 1-4. ROIs are drawn on the plot for the electrode selected in the main GUI prior to the spike sorting process; this electrode is hereafter referred to as the "event electrode." Up to two ROIs can be drawn for each cluster. Spikes are sorted sequentially from cluster 1 to 4 and are removed from the sorting pool once they have been assigned to a source cluster. This means that ROIs drawn for clusters with higher cluster identification numbers (IDs) may be drawn less selectively than for clusters with lower cluster IDs.
Because the plotted events for each electrode are all part of propagation sequences, each event on the event electrode is tied to corresponding events on all other electrodes. Thus, once the spikes are sorted on the event electrode, the spikes are also sorted on all other electrodes in accordance with the propagation sequence to which they belong. The user may then visually confirm that the intra-cluster spike shapes are consistent not only on the event electrode but also on all other electrodes.

Cluster propagation velocity
The From these realignment times, the CPV of the th event is calculated as Here, the condition > is defined such that anterograde and retrograde propagation yield positive and negative CPVs, respectively. The user may select any pair of electrodes to calculate the CPV; however, errors in the CPV due to the time resolution tend to be larger and more nonlinear for closer electrodes.
As with the SPV, the CPV is also accompanied by a confidence index, which indicates the similarity of the spike shapes recorded on the two selected electrodes and with all other spikes on the same electrode in the same cluster. This is represented by where the cross-correlation is performed over the same time window as stated previously (from 1.0 ms before to 1.0 ms after the peak voltage in each event).
It should be noted that the SPV and CPV estimates of the propagation velocity are fundamentally different in their approaches. Whereas the SPV estimates the lag for a single propagation sequence by matching the spike waveforms detected on two electrodes, the CPV essentially estimates the delay between the peak voltages on two electrodes based on the timing of the peak voltage of a "master spike" representing the cluster. Therefore, when the waveforms are not consistent across all electrodes, the SPV and CPV can yield different results, as the lag that yields the highest cross-correlation for the SPV may not correspond well to the lag between the peak voltages on the two considered electrodes.

Generation of synthetic data for validation
Synthetic data were generated in MATLAB R2016b to validate the propagation sequence detection and propagation velocity estimation capabilities of µSpikeHunter. Synthetic spikes were generated as one phase of a sinusoidal wave with an amplitude of peak = 60 µV and a duration of ts = 1.

Data availability statement
The µSpikeHunter software, as well as a sample test dataset, is available for download from [SourceForge/Github links].

Acknowledgments/Funding
This work was partly financed by FEDER -Fundo Europeu de Desenvolvimento Regional was responsible for the co-culture application ideas. All authors discussed the results and contributed to the final manuscript.

Competing interests
The authors declare no competing interests.