ORIGINAL RESEARCH article

Front. Cell. Neurosci., 06 September 2019

Sec. Cellular Neurophysiology

Volume 13 - 2019 | https://doi.org/10.3389/fncel.2019.00404

Large-Scale Mapping of Axonal Arbors Using High-Density Microelectrode Arrays

  • 1. RIKEN Quantitative Biology Center, RIKEN, Kobe, Japan

  • 2. Graduate School of Informatics, Kyoto University, Kyoto, Japan

  • 3. Carl Ludwig Institute for Physiology, University of Leipzig, Leipzig, Germany

  • 4. Department of Biosystems Science and Engineering, ETH Zurich, Basel, Switzerland

  • 5. Graduate School of Frontier Biosciences, Osaka University, Osaka, Japan

  • 6. MaxWell Biosystems AG, Basel, Switzerland

Abstract

Understanding the role of axons in neuronal information processing is a fundamental task in neuroscience. Over the last years, sophisticated patch-clamp investigations have provided unexpected and exciting data on axonal phenomena and functioning, but there is still a need for methods to investigate full axonal arbors at sufficient throughput. Here, we present a new method for the simultaneous mapping of the axonal arbors of a large number of individual neurons, which relies on their extracellular signals that have been recorded with high-density microelectrode arrays (HD-MEAs). The segmentation of axons was performed based on the local correlation of extracellular signals. Comparison of the results with both, ground truth and receiver operator characteristics, shows that the new segmentation method outperforms previously used methods. Using a standard HD-MEA, we mapped the axonal arbors of 68 neurons in <6 h. The fully automated method can be extended to new generations of HD-MEAs with larger data output and is estimated to provide data of axonal arbors of thousands of neurons within recording sessions of a few hours.

Introduction

The classical view of axons is that of mere transmission cables (Hodgkin and Huxley, ), while dendrites integrate distinct synaptic inputs (Spruston, ), and learning and memory are perceived to be a consequence of synaptic plasticity (Redondo and Morris, ). Recent data, however, indicate that the functional capacity of axons may be much more complex (Debanne, ; Ohura and Kamiya, ; Rama et al., ).

According to the classical theory describing giant squid axons (Hodgkin and Huxley, ), voltage-gated sodium and potassium channels support self-sustained action potentials that propagate along the axon. That is also true for mammalian axons, but at least two additional types of cationic channels have been described. These channels are activated by G-protein dependent receptors or hyperpolarization and modify the shape and propagation of the action potentials (Elgueta et al., ; Ko et al., ).

Axons can be surrounded by a myelin sheet that leads to saltatory conduction of action potentials (Tasaki, ), instead of the continuous conduction that has been observed in non-myelinated axons, such as the giant squid axons (Hodgkin and Huxley, ). Mammalian axons also show considerable variations in their diameter (Gaussian distribution around a peak at 200 nm diameter) and their number of branch points and varicosities, which affect the conduction velocity of action potentials. Differences in axon lengths have been found to produce defined temporal delays for spike trains for coincidence detection in the auditory system (Seidl et al., ; Stange-Marten et al., ). Furthermore, axonal delays in the neocortex show sub-millisecond precision (Swadlow, ), which is compatible with mechanisms of spike-timing-dependent synaptic plasticity (Dan and Poo, ). These characteristics lead, according to computational studies (Izhikevich, ), to the emergence of precisely timed firing patterns.

Early computational studies have suggested that some axonal arbors can act as filters for spike patterns by the selective failure of action potential propagation (Lüscher and Shiner, ). Whereas cerebellar mossy fibers in the cerebellum can reliably transmit spike trains up to 1.6 kHz (Ritzau-Jost et al., ), occasional failures have been observed in the axons of CA3 pyramidal neuron at firing frequencies of 30–40 Hz (Meeks and Mennerick, ). Failure depends on the diameter and the branching morphology as well as on the activation of presynaptic A-type potassium channels. The morphology and molecular composition of the individual axonal arbors has been found to influence the timing and shape of the action potential (Bischofberger et al., ; Alle and Geiger, ; Cho et al., ) arriving at the pre-synapse, which affects synaptic transmission and plasticity. For example, broadening of action potentials due to slow inactivation of voltage-gated potassium channels during high-frequency spike trains was found to facilitate synaptic release (Geiger and Jonas, ). Furthermore, axonal signaling and synaptic connectivity have been found to constitute important parameters in induced-pluripotent-stem-cell (iPSC) models of Parkinson's (Kouroupi et al., ) and amyotrophic lateral sclerosis (Wainger et al., ).

The understanding of information processing in axons is a fundamental question in neuroscience; however, the availability of experimental data is severely limited due to the small axon diameter. Most data sets have been gathered by performing patch clamp recordings of axonal membranes, at axon terminals and boutons, at the intact axon shaft or at the axonal bleb that is formed upon cutting axons (Ohura and Kamiya, ). However, these patch-based methods are not very well-suited to track the propagation of action potentials at more than two sites or even across the full axonal arbor. Moreover, the patch recording time is limited to a few hours while axonal recording is a serial process and has to be done axon after axon. Finally, patch-based methods constitute endpoint measurements and following the development of axonal arbors over extended periods is not possible. Another possibility to study axons includes the use of imaging tools and optogenetics (Chen et al., ), and it has been shown that imaging of single axon terminals is possible (Hoppa et al., ). However these methods are still limited in signal-to-noise ratio as well as in temporal resolution (Emmenegger et al., ).

A viable alternative to the methods described above includes the use of high-density microelectrode arrays (HD-MEAs), based on complementary-metal-oxide-semiconductor (CMOS) technology (Eversmann et al., ; Berdondini et al., ; Frey et al., ; Ballini et al., ; Bertotti et al., ; Viswam et al., ; Tsai et al., ). HD-MEAs can be used to capture neuronal activity at high temporal resolution across spatial scales, including networks, dendrites and most importantly, axons (Obien et al., ). CMOS-based HD-MEAs have been used to study the action potential propagation along axonal arbors (Bakkum et al., ) and the initiation of action potentials at the axon initial segment (Bakkum et al., ), as well as to track single action potentials along axons (Radivojevic et al., ) and stimulate single axon initial segments (Ronchi et al., ). Furthermore, HD-MEA recordings can be combined with classical patch clamp recording to study postsynaptic currents in response to pre-synaptic stimulation (Jäckel et al., ).

The objective of this work was to provide a method for high-throughput scanning of axonal arbors and mapping their axonal delays with minimal or no need to adjust parameters for the detection of axonal signals. Using a standard HD-MEA (Frey et al., ), we mapped the axonal arbors of more than 68 neurons in <6 h.

Materials and Methods

Animal Use

Timed pregnant rats (Wistar) were obtained from a commercial vendor (Nihon SLC, Japan). Animals were sacrificed on the day of arrival to obtain embryos for primary neuron cultures. All experimental procedures on animals were carried out in accordance with the European Council Directive of 22 of September 2010 (2010/63/EU) and have been approved by the local authorities in Japan (Animal Care and Use Committee of RIKEN; QAH24-01).

High-Density Microelectrode Array (HD-MEA)

Sub-cellular resolution extracellular recordings were obtained using a CMOS-based HD-MEA (Frey et al., ) with 11,011 electrodes, arranged in a hexagonal pattern and featuring an electrode density of 3,150 electrodes/mm2. The culture chamber of the HD-MEA was prepared as described before (Heer et al., ) with minor modifications: after attaching the chamber ring (polycarbonate, 19 mm inner diameter, 8 mm high) using epoxy resin (EPO-TEK 301-2, Epoxy Technology Inc.), GlobTop (G8345D-37, Namics Inc.) was used to cover the bond wires while keeping the electrode area clean, and the remaining area was covered by a thin film of PDMS (Sylgard 184, Dow Corning). Platinum black was electrochemically deposited (Marrese, ) [with modifications of the original procedure (Heer et al., )] on the electrodes to decrease their impedance in order to improve signal-to-noise characteristics (Viswam et al., ). Before plating, the surface of the HD-MEAs was rendered hydrophilic by oxygen-plasma treatment (40 s, 20 W), incubated for 4 h with of 50 μg/ml Poly-D-Lysine (Sigma-Aldrich, P7280) in PBS, washed twice with aqua dest and air-dried for 1 h.

Primary Neuron Cultures

Adult rats were anesthetized with isofluorane and killed using a guillotine. The embryos were removed from the uterus and decapitated. Their neocortex was dissected in ice-cold dissection medium (HBSS without Ca2+ and Mg2+; Gibco, NO.14175) and incubated for 20 min at 37°C in Trypsin/EDTA (Sigma-Aldrich). After washing twice with plating medium (Neurobasal A, supplemented with 10% Fetal bovine serum, 2% B27 Supplement, 1:100 GlutaMax, all from Gibco, Japan, and 10 μg/ml Gentamicin, Sigma-Aldrich, Japan), the tissue was mechanically dissociated, passed through a 40 μm nylon mesh, and centrifuged 6 min at 200 g. The supernatant was removed; the cells were suspended and counted. A 20 μl drop containing 10,000 cells was placed on the electrode area of the HD-MEA in the middle of the culture chamber. The cultures were covered with a membrane, permeable to gas but not to water vapor, and placed in a standard incubator (37°C, 5% CO2, 80% relative humidity). The neurons were allowed to settle and attach to the surface during 30 min. Thereafter, the culture chamber was filled with 600 μl serum-free, astrocyte-conditioned DMEM/Hams's F12 medium (Nerve Culture Medium, Sumitomo, Japan, #MB-X9501). Medium was completely exchanged with 600 μl conditioned medium after 4 days and then every 7 days until day 17 in vitro.

Recordings and Spike Event Detection

For recording, HD-MEAs were placed in a bench-top incubator (TOKAI HIT, Japan, INU-OTOR-RE) with temperature control, and 5% CO2 was supplied by a gas-mixer and humidified by a water bath. In order to avoid the evaporation of medium during prolonged recording intervals, the water bath and the lid temperature set point was set 1K and 3K above the sample temperature set point, which was 35°C. The HD-MEA recordings were performed using custom scripts, written in LabView (National Instruments, US), Matlab (Mathworks, US), C++ and Python running on a standard PC with a Linux operating system. Data underwent lossless data compression and were directly stored on a server on the local LAN.

Offline analysis of the recordings included filtering, event detection and averaging. First, a band-pass filter (2nd order Butterworth filter, 100–3,500 Hz) was used to remove slowly changing field potentials as well as high frequency noise. The remaining (background) noise was characterized by the median absolute deviation (MAD), which is resilient to outliers in the data but represents a consistent estimator of the standard deviation, sV = 1.4826 MAD(Vsig). Using a voltage-threshold method for event detection (Lewicki, ), negative signal peaks below a threshold of Vthr = 5sV (Vthr > 50 μV in all cases) were identified (Figure 1C). To avoid multiple detection of the same spike, successive events within <0.5 ms were discarded.

Figure 1

).

Mapping of Axonal Initial Segments

To initially identify the location of axonal initial segments (AISs), the whole array was scanned using configurations in which non-overlapping blocks of 6 × 17 electrodes (Figure 1D) were connected to the amplifiers through the switch matrix (Frey et al., ) (schematic in Figure 1B). After recording for 30 s from every electrode in the array, spike detection was performed for each electrode and taking the median value summarized the amplitude of the negative peaks. The median of the negative peak for each electrode was plotted as a map showing some areas with large negative peaks (Figure 2A). Such local minima were assumed to indicate the putative (proximal) AIS locations (Bakkum et al., ). We did not distinguish between inhibitory and excitatory neurons according to the AIS waveform (Mita et al., ).

Figure 2

). The asterisc (*) indicates an artifact probably from a large spike at the AIS of a presynaptic neuron.

Mapping of Electrical “Footprints”

To map the electrical “footprint” (spatial distribution of extracellularly measured electrical potentials obtained with the densely packed electrodes) of neuronal units and their axonal arbors over the entire array, we used a series of configurations in which so called “fixed electrodes” were always connected to the amplifiers through the switch matrix (Frey et al., ), whereas the remaining “variable electrodes” were connected in sequential configurations (schematic in Figure 1E). After recording the spontaneous activity for 115 s for each of these configurations, spike sorting was performed and the electrical “footprint” was calculated using the spike-triggered average as previously described Bakkum et al., ; Müller et al., ; Deligkaris et al., ; Radivojevic et al., , ;Jäckel et al., .

To perform spike-triggered averaging we need to reliably record spiking activity of the neurons. Therefore, we selected the fixed electrodes at the putative locations of the AISs but imposed a spatial restriction in order to not record the same neuron twice. For selection of fixed electrodes, all electrodes were ranked according to their median negative peak amplitude (see previous section). The electrode with the highest rank was selected, afterwards all electrodes in its proximity (within 100 μm distance) were discarded from the list, and the procedure was repeated.

Spike Sorting

Even in sparse cultures some electrodes will pick up spikes from multiple neurons in their vicinity. Therefore, spike sorting was performed for each fixed electrode after recording. Waveforms were extracted for a period of 1.5 ms before to 1.5 ms after the negative peak of each event comprising 61 samples for each event. In order to extract those features that best separate the different clusters of spikes, we performed a principal component analysis. We choose the first 10 principal components as spike features (Abeles and Goldstein, ), containing more than 85% of the energy of the signal. KlustaKwik (Kadir et al., ) was used to fit a mixture of Gaussians with unconstrained covariance matrices and automatically selected the number of mixture components. Clusters containing more than 2,500 spiking events in total were assigned to distinct neurons. This number corresponds to an average of 14 spikes in each of the 179 recording configurations. With 14 samples the spike triggered averaging reduces the noise by approximately a factor . This procedure omits neurons with a firing rate below 14 spikes/115 s ≈ 7 spikes/min for which we could not compute a reliable footprint. The uniqueness of the detected neurons was confirmed by assessing their different footprints (Figure 6). Clusters leading to identical footprints were manually merged to single neurons, and the footprints of the merged neurons were re-calculated.

Optimal Recording Configurations for High-Throughput Scanning of Electrical Footprints

The extracellular signals originating from axons and dendrites are very small with respect to the background electrical activity and noise, so that spike-triggered averaging was applied. We developed a set of recording configurations to map the electrical footprint of several neurons in parallel by utilizing the switch matrix of our HD-MEA. The switch matrix can be dynamically configured to connect a large number, e, of electrodes to a smaller number, a, of amplifiers. In a first-order approach, one could use one electrode as trigger and the remaining electrodes to scan the neuronal footprint, which would result in a large number of configurations, cw, needed to scan the whole array, cw = e/a. For recording axonal arbors of n neurons, electrodes near the AISs of these n neurons have to be always connected to amplifiers (n fixed electrodes). The remaining amplifiers can then be connected to the remaining electrodes in several successive configurations (variable electrodes). Following this procedure, the whole array can be scanned with

configurations. For n neurons we need c(n) configurations, which means on average C(n) = c(n)/n configurations for one neuron. An optimal strategy means to choose n in such a way that C(n) → min for 0 < n < a. With n < a ≪ e, the number of neurons being much smaller than the number of electrodes, we can approximate:

The right-hand side has a minimum for n = a/2. Therefore, approximately half of the amplifiers should be connected to fixed electrodes. The other half of the amplifiers can then be used to scan the whole array in 2cw configurations, or on average configurations per neuron. The axonal arbors of a single neuron may extend over the whole array, but by scanning the axonal arbors of many neurons in parallel, the average number of configurations, Ca, per neuron is much less than the number of configurations, cw, required for scanning the whole array for large a:

This is due to the fact that increasing the number of amplifiers quadratically decreases the average time to scan a single neuron, which allows for high-throughput acquisition of axonal delay maps.

Live Imaging

Live-cell visualization of whole neurons was performed by transfection (Bakkum et al., ; Radivojevic et al., ). Transfection was performed using a pLV-hSyn-RFP plasmid from Edward Callaway (The Salk Institute, US; Addgene, #22909) and Lipofectamine 2000 (Life Technologies) in accordance with the manufacturer's protocol. A Leica DM6000 FS microscope, a Leica DFC 345 FX camera, and the Leica Application Suite software were used to produce the micrographs.

Results

To initially identify the location of axonal initial segments (AISs), the whole array was scanned by using configurations in which non-overlapping blocks (Figure 1D). The median of the negative peak for each electrode was plotted as a map showing some areas with large negative peaks (Figure 2A). Such local minima were assumed to indicate the putative (proximal) AIS locations (Bakkum et al., ). In low-density cultures, several AISs could be distinguished (Figure 1A). So called “fixed electrodes” were then selected as trigger electrodes at the putative locations of the AISs, while the remaining “variable electrodes” were selected in sequential configurations to map the electrical footprint of neuronal units and their axonal arbors (Bakkum et al., ; Müller et al., ; Deligkaris et al., ; Radivojevic et al., , ; Jäckel et al., ) over the entire array (Figure 1E). After spike-triggered averaging, the amplitude and delay of axonal signals could be mapped across the whole recording area for each individual neuron (Figure 2B).

High-Throughput Scanning of Electrical “Footprints”

We used HD-MEA with 11,011 electrodes and 126 amplifiers (Frey et al., ). We connected 62 amplifiers to fixed electrodes and used the remaining amplifiers to scan the whole array selecting 64 electrodes at random positions, which would yield a total of configurations (see section Spike Sorting). However, due to the design of the switch matrix in our HD-MEA, some electrodes cannot be selected in the same configuration, which increased the number of recording configurations to 175 (see Table 1).

Table 1

This workExtrapolated performance using other HD-MEA
HD-MEA(Frey et al., )(Ballini et al., )(Dragas et al., )(Yuan et al., )
ModeSMSMSMSMAPS
Electrodes11,01126,40059,7608,6408,640
Amplifiers1261,0242,0481129,216
Configurations17951571531d
Neuronsa625121,024561,000e
Configurations/neuron2.80.10.12.70.001
Circuit noise (AP band)μVRMS2,42,42,42,312,4
Total noise (AP band)bμVRMS5555,013,2
Factorc11116,9
Time/configurations115115115113—d
Total timemin3439711028813
Average time/neurons3321163091

Comparison of high-throughput mapping of axonal delays using different types of CMOS-based HD-MEA.

AP band = 300 Hz−10 kHz. SM, switch matrix; APS, active pixel sensor. The expected number of neurons which can be mapped simultaneously in a single recording session are highlighted in boldface.

a

Assuming one neuron per fixed electrode. Note that spike sorting should be used in dense cultures to identify single-neuron activity, which would increase this number.

b

Total (RMS) noise consists of circuit noise, noise from dendritic Pt black electrodes (1.8 μV) (Viswam et al., ) as well as from background neuronal activity (4 μV). These noise sources are uncorrelated and they add up to 3.1 μV in saline (Frey et al., ) and to about 5 μV in the neuronal cultures cultured on our HD-MEA. From these values, we can estimate the total noise while recording from such cultures on other HD-MEA with similar electrodes.

c

Factor by which the recording time per configuration increases to obtain a similar noise reduction after spike-triggered averaging. This factor is the square of the ratio between the total noises for each HD-MEA.

d

Full-frame readout does not require selection configuration with different electrodes.

e

Number of neurons was fixed at 1,000 for comparison with HD-MEAs with larger number of electrodes and switch-matrix design.

In our culture, 53 of the 62 fixed electrodes recorded single-unit activity. However, only 28 axonal arbors could be identified, because 25 electrodes did not record enough spikes to reliably determine arbor structures (see Methods). The remaining 9 fixed electrodes recorded multi-unit activity, so that spike sorting was used to identify single-neuron activity and to obtain additional 40 axonal arbors. In total, we mapped the axonal arbors of 68 neurons. For further analysis we only used the n = 46 neurons with axonal arbors extending over more than 50 electrodes (Figure 6).

Identification of Axonal Arbors

Previously (see Figure 6 in Bakkum et al., ), axonal arbors were traced according to the occurrence of negative peaks in signal amplitudes in the spike-triggered averages of extracellularly recorded electrical signals that exceeded 5 times the background noise (sV). This method (hereafter termed “method I”) evaluates the spike-triggered averages for each electrode separately, without considering their spatial arrangement and signal correlations between neighboring electrodes. We used another method (hereafter termed “method II”) to benefit from the fact that the delay of the negative signal amplitude peak in the extracellularly measured potential originating from a single axon is very similar at neighboring electrodes it passes by. When we mapped the delay of these negative signal peaks in the spike-triggered averages (Figure 3A), the resulting map showed a distinct region of similar delays against a background of random delays (Figure 3B). This finding is due to the fact that signals originating from a common source, e.g., from an axon of the same neuron, are very similar across neighboring electrodes, whereas in the case of random signals the negative peak could occur anywhere in the interval [Δtpre, Δtpost] for which the spike triggered averaging was performed (see histogram in Figure 3C). Δtpre and Δtpost represent the boundaries relative to the spike triggering, so that the total spike-triggered average has the length T = Δtpost − Δtpre. To quantify the smoothness of the delay map, the delay was sampled in a neighborhood (compare Figure 4H) of hexagonal shape around each electrode, and the sample standard deviation Sτn was calculated. Note, that in the case of a uniform distribution over the interval [0, 1], the standard deviation of a sample is bounded by 0 ≤ s ≤ 0.5 and shows a characteristic distribution, depending on the number of observations, N, in each sample.

Figure 3

Figure 4

This distribution shows a sharp peak around the mean standard deviation:

There is no analytical expression for this distribution, but after a coordinate transformation of the interval, it can be approximated by a beta distribution .

In case that axonal signals are present, the delays in the neighborhood of an electrode with a negative peak at t are distributed in the interval , depending on the velocity, c, of the action potential and the distance, r, between the electrodes. Therefore, the mean of the sample standard deviation for axonal delays is:

For an HD-MEA with r = 18 μm and a typical conduction velocity for short-range-projecting axons in the rat neocortex of 0.3−0.44 m/s (Lohmann and Rörig, ; Telfeian and Connors, ), a standard deviation of around 30 μs can be expected (compare with Figure 3E). At the boundary, more and more neighboring electrodes will no more pick up the axonal signal, so that the sample standard deviation shifts toward:

In our case, with T = 8 ms, a standard deviation of 2.3 ms for background signals would be expected (compare with Figure 3D). Empirically, the distribution of saxon can be approximated by an (truncated) exponential distribution (see below). Therefore, axons have a distribution of sτn with a peak close to zero, which is clearly distinguishable from the distribution for the background. A threshold smin, placed at the local minimum in the sτn distribution (Figure 3G), can be used to separate both populations. The electrodes with a sτn below this threshold represent negative peaks that are consistent across neighboring electrodes (Figure 3I). If these peaks appear after the negative peak at the AIS (Figure 3H), they are assumed to originate from the axonal arbor of the neuron (Figure 3J).

For a limited number of neurons, the ground truth in the form of fluorescence images was available, so that we could compare the performance of the new method (method II; Figure 4B) with the old method (method I; Figure 4A) employing a fixed threshold at 5sV (Bakkum et al., ). An example of a transfected neuron is shown in Figure 4C. More electrodes were selected by method II than by method I, which could be compensated for by lowering the threshold, e.g., to 3sV (Figure 4D). In order to compare the electrode selection (a set of electrode coordinates, E) with the ground-truth axon information (a set of pixels in the image, A), we used the Haussdorff distance , which is commonly used in computer vision to measure how strongly two shapes differ from each other. After registering the fluorescence image to the electrode coordinate system, we calculated

using the Euclidean distance d(a, e) between the coordinate of an electrode recording an axonal signal and the coordinate of the pixel representing an axon in the fluorescence image. For smaller threshold for method I lead to a larger deviation from the ground truth than method II. This was mainly due to the fact that, more electrodes far away from the axonal arbors were selected. The new method inherently relies on adjacency and rejected these “outliers” and produced more compact maps that more closely followed the ground truth. In other words, the distributions of the feature used to classify the electrical activity as either axonal signals or as background, showed a larger overlap for method I than method II. We tested the robustness of method II against the spatial distance of the electrodes by increasing the spatial extension of the neighborhood while keeping the number of electrodes in each neighborhood constant (N = 7). When hexagonal patterns with r = 2 × 18 μm (Figure 4E) or r = 3 × 18 μm distance between electrodes were selected, the distance to the ground truth only slightly increased (Figure 4G). However, it seemed that for a carefully chosen threshold (e.g., around 4.5sV, H≈100 μm, Figure 4F), the original method performed as well as the new method.

Due to the limits of the Haussdorff distance in estimating the quality of the extracted mappings, we also calculated the receiver-operator characteristics (ROC) (Fawcett, ), which is possible, as both segmentation methods are binary classifiers. We fitted the respective empirical distributions of their scores with a mixture of two partially overlapping distributions representing axonal signals (“positive” class, P) and background activity (“negative” class, N) with:

  • two normal distributions (Figure 5A) for the segmentation with method I (Figure 5D)

    with , μN < μP

  • a beta distribution and a truncated exponential distribution (Figure 5B) for the segmentation with method II (Figure 5E):

    with , for 0 ≦ x ≦ 1

Figure 5

For both methods, we used the fitted distributions to estimate the probability observing axonal signals (P+) or background activity (P−).

  • and for method I,

  • and for method II.

We then numerically calculated the cumulative distributions in order to calculate the true positive rate (TPR) and false positive rate (FPR) for each threshold, x, and plotted them as an ROC curve (Figure 5C). The area under the curve (AUC) showed that the new method consistently had a better performance than the original method for a total of n = 46 neurons (Figure 5F). Furthermore, the automatic threshold procedure yielded a much better TPR (Figure 5G) at the expense of a slightly increased false-positive rate FPR in comparison to the original method (Figure 5H).

Discussion

HD-MEAs with more than 3,000 electrodes per square millimeter and dedicated low-noise on-chip amplifiers are suitable tools to record the electrical activity of individual axonal arbors. We first optimized a recording scheme for the switch matrix HD-MEA that relied on combinations of fixed and variable recording sites for high-throughput parallel mapping of as many neurons as possible per total recording time. The method described here shows a promising way to obtain axonal arbors at large scale from potentially all neurons during a single recording session. As an example, 68 neurons were mapped in parallel, 48 neurons of which featured large axonal arbors (Figure 6). This yield can be further improved using an HD-MEA design with an increased number of simultaneously active recording channels, decreasing the number of necessary measurement configurations and, hence, the on average the required measurement time per neuron to a few seconds (Table 1).

Figure 6

We then developed a method to distinguish axonal signals from the background noise. Axonal arbors reveal themselves by typical waveforms of the extracellular electric field potentials (Bakkum et al., , ; Petersen et al., ; Deligkaris et al., ). Previously (Bakkum et al., ), axonal arbors were traced according to the occurrence of negative peaks in signal amplitudes in the spike-triggered averages of extracellularly recorded electrical signals that exceeded 5 times the background noise. With this threshold more than 50% of the smaller axonal signals are lost (Figure 5G). Lowering the threshold can alleviate this problem, however, the likelihood that background signals are falsely assigned to axonal arbors also increases. Some electrodes measuring spurious signals can easily be detected, if they are located at comparably large distance from electrodes that record genuine axonal signals. Typically, these false positives are removed manually after visual inspection of the recorded waveforms and their isolated location in the recording area. However, this procedure is not possible for more than a few recorded axonal arbors. Therefore, we developed a new method that compares the axonal delays on neighboring electrodes. If these delays are similar, the corresponding signals most likely originated from the same axon. If there are large differences in the delay, the corresponding signals may represent background noise. Both methods, amplitude thresholding (Bakkum et al., ) and the new method, worked well when we validated them by comparing their results with the morphologies of axonal arbors obtained by sparse transfection (Bakkum et al., ; Radivojevic et al., ). Upon comparing both methods according to their ROCs on a larger data set of 46 axonal arbors, the new method based on local correlations outperformed the classic amplitude thresholding. Furthermore, the new method does not require manual intervention, e.g., for setting a threshold.

The example workflow demonstrated here enables high-throughput scanning of axonal arbors and mapping of their axonal delays without the need to adjust parameters for the detection of axonal signals. This method can be extended to new generations of HD-MEAs (Ballini et al., ; Viswam et al., ; Yuan et al., ) and can be used to obtain data from axonal arbors of thousands of neurons within a recording session of a few hours (Table 1). Our example workflow yielded axonal arbors of 68 neurons, 46 of which featured axonal arbors extending over more than 50 electrodes (Figure 6C). These axonal arbors showed considerable variation in total length, but also in local branching patterns and axonal delays (Supplemental Figure 1). Interestingly, some neurons featured footprints in the form of disconnected patches, which could be the signature of saltatory conduction of action potentials. This observation is consistent with the reported formation of nodal components in neuronal cultures (Freeman et al., ), but requires further experimental validation.

Our method can be used for automated selection of neurons with suitable axonal arbors for stimulation experiments (Jäckel et al., ), single action potential tracking (Radivojevic et al., ), automated patching (Obien et al., ), for mapping of ion receptors by local drug application (Sasaki et al., ) and even for automated single-cell phenotyping of axonal conduction in human iPSC-derived neuronal cultures.

Statements

Data availability statement

The Hana (high density microelectrode array recording analysis) analysis pipeline is open source. All source code as well as example data to replicate the figures are available at: http://github.com/tbullmann/hdmea_axon. The example data consists of spike triggered-averages that were extracted from the raw recordings.

Ethics statement

All experimental procedures on animals were carried out in accordance with the European Council Directive of 22 September 2010 (2010/63/EU) and had been approved by the local authorities (Animal Care and Use Committee of RIKEN; QAH24-01).

Author contributions

TB designed the study, performed experiments, wrote the software, analyzed data, assembled figures, interpreted the results, prepared, and revised the manuscript. MR performed recording and imaging experiments. SH implemented and tested analysis algorithms. KD performed cell culture experiments. AH interpreted the results, prepared, and revised the manuscript. UF planned the study, supported the experiments, interpreted the results, and revised the manuscript.

Funding

The HD-MEA work at ETH Zurich was financially supported by the European Community through the European Research Council Advanced Grant 694829 neuroXscales (Horizon 2020) and the Swiss National Science Foundation Grant 205321_157092/1 (Axons). Financial support through the Swiss Commission for Technology and Innovation project 25933.2 PFLS-LS is also acknowledged.

Acknowledgments

We thank Alexander Stettler and Peter Rimpf for post-processing CMOS chips as well as Manuel Schröter, Roland Diggelmann and Felix Franke for helpful discussions about spike sorting.

Conflict of interest

UF is a co-founder of MaxWell Biosystems AG, Mattenstrasse 26, Basel, Switzerland. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2019.00404/full#supplementary-material

Supplementary Figure 1

Sholl analysis for number of electrodes (A) and axonal delay (B) for n = 46 neurons shown in Figure 6.

References

Summary

Keywords

axons, high-density microelectrode array, extracellular electrical field, action potential, axonal arborizations, action potential propagation, high-throughput screening

Citation

Bullmann T, Radivojevic M, Huber ST, Deligkaris K, Hierlemann A and Frey U (2019) Large-Scale Mapping of Axonal Arbors Using High-Density Microelectrode Arrays. Front. Cell. Neurosci. 13:404. doi: 10.3389/fncel.2019.00404

Received

23 May 2019

Accepted

20 August 2019

Published

06 September 2019

Volume

13 - 2019

Edited by

Dominique Debanne, INSERM U1072 Neurobiologie des canaux Ioniques et de la Synapse, France

Reviewed by

Adelaide Fernandes, University of Lisbon, Portugal; Valentina Carabelli, University of Turin, Italy

Updates

Copyright

© 2019 Bullmann, Radivojevic, Huber, Deligkaris, Hierlemann and Frey.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Andreas Hierlemann

This article was submitted to Cellular Neurophysiology, a section of the journal Frontiers in Cellular Neuroscience

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Share on WeChat

Scan with WeChat to share this article

Article metrics