Functional imaging of conduction dynamics in cortical and spinal axons
Abstract
Mammalian axons are specialized for transmitting action potentials to targets within the central and peripheral nervous system. A growing body of evidence suggests that, besides signal conduction, axons play essential roles in neural information processing, and their malfunctions are common hallmarks of neurodegenerative diseases. The technologies available to study axonal function and structure integrally limit the comprehension of axon neurobiology. High-density microelectrode arrays (HD-MEAs) allow for accessing axonal action potentials at high spatiotemporal resolution, but provide no insights on axonal morphology. Here, we demonstrate a method for electrical visualization of axonal morphologies based on extracellular action potentials recorded from cortical and motor neurons using HD-MEAs. The method enabled us to reconstruct up to 5-cm-long axonal arbors and directly monitor axonal conduction across thousands of recording sites. We reconstructed 1.86 m of cortical and spinal axons in total and found specific features in their structure and function.
Editor's evaluation
This work provides fundamental new insight into fine axonal morphologies based solely on extracellular action potential recordings. They provide compelling evidence of fine resolution in mapping functional connections between neurons. The work may have broad use in neurobiology, bioengineering, stem cell biology, as well as tissue engineering in functional characterization.
https://doi.org/10.7554/eLife.86512.sa0Introduction
Axons are neuronal processes specialized for the conduction of action potentials (APs). Cortical axons serve as communication cables between various types of neurons that are synaptically connected and, accordingly, arranged in multiple layers to receive, process, and convey neural information between different regions in the brain. Motor neurons located in the ventral horn of the spinal cord project their axons outside the central nervous system (CNS) and are responsible for the contraction of effector muscles in the periphery. Spinal axons are specialized to innervate and precisely control different types of muscle fibers, thus ensuring refined coordination of complex body movements (Stifani, 2014).
Mainly due to difficulties to experimentally access axonal conduction, axonal information processing has been neglected, and axons are classically seen as conductive cables that do nothing more than faithfully transmit APs in an all-or-none, ‘binary’ fashion (Hodgkin and Huxley, 1952). Later studies have challenged this view and suggested that axons have much more complex roles than traditionally thought (Alcami and El Hady, 2019). Contrary to classical concepts, reports show that, besides conducting binary APs, hippocampal and cortical axons can also transmit analog currents in a passive fashion (Shu et al., 2006; Alle and Geiger, 2008). It has been shown that analog currents integrate with ongoing axonal APs, changing their waveforms and, consequently, affecting synaptic transmission in a graded, ‘analog’ manner (Shu et al., 2006; Alle and Geiger, 2006; Zbili and Debanne, 2019). Additionally, it was found that axons modulate the waveforms of APs contingent on neuronal activity and, accordingly, adjust the synaptic release (Geiger and Jonas, 2000). Moreover, during high-frequency regimes of neuronal activity, axons are able to reduce their conduction velocities by more than 20% and, thereby, tune the timing of AP arrival at synapses (Radivojevic et al., 2017). Taken together, these findings suggest that axons passively and actively process APs to tune the amount of information transmitted by synapses, but also imply that axons play a crucial role in the temporal coding of the neural information.
Morphological complexity of mammalian axons generally depends on their downstream targets’ spatial disposition. Thus, for example, cortical neurons form extensively branched axons to convey APs to numerous postsynaptic neurons located in different cortical layers and regions in the brain. Due to their small diameter of 0.08–0.4 μm (Debanne, 2004) and absence of surrounding myelin sheaths, cortical axons provide relatively low conduction velocities of 0.1–2 m/s (Radivojevic et al., 2017; Bakkum et al., 2013), which are, nevertheless, sufficient to ensure rapid communication between closely spaced neurons. As opposed to cortical neurons, motor neurons transmit signals to distant targets outside the CNS and, for that purpose, develop considerably longer (up to 1 m), thicker (0.5–10 μm), and myelinated axons that provide considerably higher conduction velocities of up to 100 m/s (Debanne et al., 2011; Saliani et al., 2017).
Despite being acknowledged as reliable conduction cables, axons can in certain cases fail to propagate APs faithfully. Such cases are referred to as conduction failures and are attributed to particular axonal morphology. For instance, experimental and theoretical studies suggest that conduction failures are likely to occur at axonal branching points and local swellings due to abruptly increased axon diameter (Debanne, 2004). In addition, the AP propagation fidelity and temporal precision depend on axon diameter, biophysical properties of various types of ion channels, and thermodynamic noise inherent to their gating dynamics (Faisal and Laughlin, 2007; Faisal et al., 2005; White et al., 2000). According to theoretical studies, thin distal axons (diameter <1 μm) are prone to ‘channel noise’, which can introduce variability in axonal AP waveforms (Neishabouri and Faisal, 2014), increase the jitter of AP propagation (Radivojevic et al., 2017; Faisal and Laughlin, 2007), and compromise the reliability of the conduction itself (Faisal and Laughlin, 2007; Skaugen and Walløe, 1979; Schneidman et al., 1998). Besides geometrical factors, neuronal activity can also affect AP conduction, which is consistent with our previous finding that high-frequency neuronal activity increases the jitter of AP propagation in cortical axons (Radivojevic et al., 2017).
Ultrastructural aberrations and malfunctioning conduction in human axons are often early hallmarks of neurodegenerative diseases. For instance, in Alzheimer’s disease (AD), abnormal protein aggregates cause local swellings in axons which have been shown to reduce conduction velocity and even block AP propagation (Blazquez-Llorca et al., 2017). Human post-mortem studies suggest that axonal degeneration may be the earliest feature of Parkinson’s disease and, therefore, an appropriate target for early therapeutic interventions (Burke and O’Malley, 2013). Motor neuron degeneration is the hallmark of amyotrophic lateral sclerosis (ALS), where axonal dysfunction begins long before symptom onset and motor neuron death (Suzuki et al., 2020).
Comprehension of axon neurobiology in health and disease is generally limited by the technologies available to study axonal function and structure integrally. The whole-cell patch-clamp technique, complemented with fluorescence microscopy, is commonly used to correlate electrophysiological data with the morphological properties of the neuron. Namely, a patch-clamp pipette can inject fluorescent dyes directly into the neuron and, therefore, allow visualization of the axonal arbor during the electrophysiological experiment. However, patch-clamp is typically limited to recording APs from a single axonal site (Shu et al., 2006; Alle and Geiger, 2006; Forsythe, 1994) and does not allow for tracking AP propagation across axons. Moreover, the technique is invasive and destructive, which constrains the duration of a recording session to about an hour. Alternatively, fluorescent indicators sensitive to voltage (Peterka et al., 2011) or calcium (Grienberger and Konnerth, 2012) can be used to observe axonal AP propagation and, at the same time, to visualize axonal morphology. However, fluorescent indicators exhibit photobleaching and phototoxicity and may perturb the physiology of the cell to the point of affecting AP conduction (Peterka et al., 2011). Complementary-metal-oxide-semiconductor (CMOS)-based high-density microelectrode arrays (HD-MEAs) have been designed to record extracellular APs from neuronal cultures (Müller et al., 2015) and allow tracking axonal signals across hundreds of microelectrodes (Radivojevic et al., 2017; Bakkum et al., 2013). Thanks to a low-noise CMOS design, HD-MEAs enable detection of APs across entire arbors of cortical axons, including tiny axon terminals (Radivojevic et al., 2017). HD-MEAs provide noninvasive access to axonal APs and impose no constraints on the duration of the recording sessions. Nevertheless, HD-MEA technology does not provide direct insights into axonal morphology. It has to be complemented with auxiliary optical techniques to allow correlation between axonal function and structure. Live-cell imaging techniques can be used to optically visualize axonal morphologies directly on HD-MEA surfaces (Radivojevic et al., 2017; Bakkum et al., 2013; Bullmann et al., 2019). These techniques, however, entail introduction of fluorescent reporters into the cell, which induces chemotoxicity, phototoxicity, and cell death (Icha et al., 2017).
Our objective was to enable the reconstruction of axonal morphologies based solely on extracellular APs recorded from in vitro mammalian neurons using the CMOS-based HD-MEA system. The present study reports a method for tracking AP propagation across tens-of-millimeter-long nonmyelinated axonal arbors of primary cortical and motor neurons cultured on HD-MEAs. The method allows for label-free electrical visualization of axonal conduction trajectories and, at the same time, provides noninvasive access to axonal APs waveforms recorded across hundreds of microelectrodes. Using the developed method, we investigated (I) morphological features of cortical and spinal axons, (II) fluctuations in AP amplitudes across axonal arbors, and (III) temporal dynamics of the axonal conduction.
Results
We developed a method for the automatic reconstruction of axonal morphology based on extracellular APs recorded from cortical and motor neurons using the HD-MEA system (Figure 1A). The reconstructed morphology is here referred to as ‘functional morphology’, and it reveals the location of the axon initial segment (AIS), axonal trunk, and higher-order branches (Figure 1B). The method was used to trace extracellular APs propagating across axonal arbors of cortical (Figure 1—video 1) and motor neurons (Figure 1—video 2). The experimental design and overview of the method are outlined in Figure 2 and are described below.
Experimental design and overview of the method
We cultured rat primary cortical (Figure 2—figure supplement 1) and motor neurons (Figure 2—figure supplement 2) on a CMOS-based HD-MEA system that comprises 26,400 densely packed microelectrodes (Figure 2—figure supplement 3). The cultures were grown in monolayers with estimated thicknesses of ~5–40 and ~5–50 μm for cortical and motor neuron cultures, respectively (Figure 2—figure supplement 4). Both cortical and motor neuron cultures yielded cell densities of ~500–2000 neurons per mm2 (Figure 2—figure supplement 4). The cultures matured in controlled conditions (see Methods), and their extracellular electrical activity was recorded between 12 and 24 days in vitro (DIV). Mature neurons grew their axonal arbors efficiently across sensing areas of the HD-MEA chips and provided a tight interface between axons and microelectrodes (Figure 2—figure supplement 5). The relative proximity of axons to the sensing area varied across axonal arbors within ranges of ~1–13 and ~1–18 μm for cortical and motor neurons, respectively (Figure 2—figure supplement 5).
The CMOS-based HD-MEA system enabled us to map extracellular APs across entire cultures and to electrically identify individual neurons (Figure 2—figure supplement 6). Spontaneous neuronal activities were recorded across all microelectrodes for 2 min, and average amplitudes of the recorded voltage traces were used to produce network-wide activity maps (Figure 2—figure supplement 6A). Because extracellular APs with the largest amplitudes arise from the AIS (Bakkum et al., 2019), local maxima found within these maps indicated the location of individual neurons in the network (Radivojevic et al., 2016).
We used spike-sorting algorithms to discern signals among individual neurons in the culture (see Methods). The dense arrangement of the microelectrodes allowed us to access the electrical activity of a single neuron at high spatial resolution (Figure 2—figure supplement 6A, Figure 2—videos 1 and 2); however, overlapping signals recorded from multiple neighboring neurons were observed in most of the cases (Figure 2—figure supplement 6B). Spike-sorting procedures enabled us to discern APs among individual neurons reliably and extract relative times of their activities (‘spike times’) (Radivojevic et al., 2017; Radivojevic et al., 2016; Jäckel et al., 2012). Sorted APs were then averaged over adjacent electrodes to reveal the spatiotemporal distribution of a single neuron’s activity. The spatiotemporal distribution of APs recorded from proximal neuronal compartments (near the AIS) is called an ‘extracellular AP footprint’. Electrical footprints reconstructed for eight neighboring neurons are presented in Figure 2—figure supplement 6C. Z-stack image series of the corresponding culture is shown in Figure 2—videos 1 and 2.
We used spike-triggered averaging for electrical imaging of axonal arbors (Figure 2—figure supplement 7). High-amplitude APs detected near the AIS were used to trigger the averaging of single voltage traces (‘single trials’) recorded across all electrodes in the array (Figure 2—figure supplement 7A). Because the AIS signal represents the first occurring (initial) trace of the neuron’s activity, the averaging reveals spatial and temporal shifts in propagating axonal signals (Figure 2—figure supplement 7B). The spatiotemporal distribution of axonal signals reconstructed using spike-trigger averaging is called an ‘axonal electrical image’. Axonal electrical images of three neighboring neurons are presented in Figure 2—figure supplement 7B and in Figure 2—video 3.
We developed a method for reconstructing axonal functional morphologies based purely on features extracted from axonal electrical images. Adaptive thresholding was used to map signal peaks across axonal arbors (Figure 3), and the moving object tracking technique was used to reconstruct axonal conduction trajectories (Figure 4). The method for reconstruction of axonal functional morphologies is described in the following two sections.
Adaptive thresholding enables mapping of axonal signal peaks in space and time
Adaptive thresholding applied to axonal electrical images enabled us to map AP peaks propagating across axonal arbors (Figure 3). Axonal electrical images were obtained by averaging 200 voltage traces per electrode. Time derivatives of averaged traces (μV/μs) were computed for each of the electrodes, and the resulting data were divided into 400 consecutive timeframes (with 50 μs inter-frame intervals). Adaptive thresholds were next applied at each timeframe to detect AP peaks at different time points during axonal propagation.
Our thresholding technique generally utilized a ‘greedy algorithm’ principle (see Discussion). The algorithm allowed for recurrent adaptations of the threshold based on parameters updated along the detection process. High-amplitude APs were mapped first, and parameters obtained from the mapped locations were used to tailor the detection of low-amplitude signals in the later steps. The thresholds were initially determined based on electrical noise observed from the HD-MEA chip. Electrical noise was estimated from voltage traces sampled across an entire array during periods when the observed neuron was inactive, and the noise was estimated for each neuron separately.
Adaptive thresholding involves three steps (Figure 3, Figure 3—videos 1–3). In the first step, a simple planar threshold, set to 9 STD of the estimated noise, was used to detect high-amplitude signal peaks (Figure 3A, Figure 3—video 1). A high threshold level allowed us to detect AP peaks far above the background noise. However, axonal APs with lower amplitudes remained undetected, leaving gaps along axonal conduction trajectories. In the second step, confined thresholds, set to 2 STD of the estimated noise, were applied locally to detect low-amplitude AP peaks (Figure 3B, Figure 3—video 2). The confined thresholds were positioned on spatial and temporal coordinates of previously mapped peaks. The thresholds were confined spatially to 50 μm radii and temporally to periods encompassing three consecutive timeframes (tprevious, tcurrent, tnext, see Figure 3B). The spatiotemporal confinement enabled us to detect low-amplitude signals in close proximity to previously mapped peaks and to fill local gaps along the conduction trajectories. The third step utilized the same detection strategy but with differently tuned parameters (Figure 3C, Figure 3—video 3). Namely, the threshold level was further lowered to 1 STD of the estimated noise, and the spatial confinement was broadened to a 100 μm radius. These parameters enabled the detection of low-amplitude APs near axon terminals.
A moving object tracking technique enables reconstruction of axonal conduction trajectories
We developed an algorithm for tracking mapped axonal signals across consecutive timeframes (Figure 4). The algorithm enabled us to reveal signal conduction trajectories and reconstruct functional morphologies of cortical (Figure 1—video 1) and spinal axons (Figure 1—video 2).
Our tracking algorithm was designed to predict axonal conduction trajectories based on three factors: (I) spatiotemporal proximities of the mapped peaks, (II) topology of skeletonized axonal signal, and (III) conduction velocities estimated from previously reconstructed trajectories.
The tracking procedure involved three iterative steps (Figure 4, Figure 4—videos 1–3). (I) Direct interconnection – the closest signal peaks found within a 100 μm Euclidean distance and mapped in two consecutive timeframes were interconnected via direct links (Figure 4A, Figure 4—video 1). Direct interconnection revealed fragments of axonal conduction trajectories but failed to reconstruct axonal branching forks in most cases. Reconstructed fragments were used to calculate axonal conduction velocities. Obtained velocities served as criteria to tailor the tracking procedure in the next two steps. (II) Skeletonization-assisted interconnection – consecutive signal peaks that could not be interconnected directly, but were found within a 200 μm Euclidean range, were interconnected via skeletonized remnants of the axonal signal (Figure 4B, Figure 4—video 2). Signal amplitudes extracted from the corresponding electrical images were averaged over two consecutive timeframes (Δt=100 μs) and mapped over interconnection areas. Obtained maps were next skeletonized to infer propagation trajectory between the consecutive signal peaks. Conduction velocities estimated in the previous step were used as criteria for the selection of optimal propagation trajectories. Trajectories whose velocities deviated from previously estimated values by more than 50% were discarded. Skeletonization-assisted interconnection enabled reconstruction of some, but not all, axonal branching forks. It was designed to predict axonal trajectories over consecutive (continuous) signal peaks but could not predict trajectories between discontinuous signal peaks. Axonal conduction velocities estimated from the reconstructed trajectories served as criteria to tailor the tracking procedure in the next step. (III) Indirect interconnection – discontinuous signal peaks mapped in every other timeframe and found within a 400 μm Euclidean range were indirectly interconnected (Figure 4C, Figure 4—video 3). The propagation trajectory between the discontinuous signal peaks was reconstructed using remnants of skeletonized signals. Signal amplitudes extracted from the corresponding electrical images were averaged over three consecutive timeframes (Δt=150 μs) and mapped over interconnection areas. Obtained maps were next skeletonized to infer propagation trajectory between signal peaks found in the first and third timeframe. Conduction velocities estimated in the previous steps were used as criteria for selecting optimal propagation trajectories and predicting spatial coordinates of data for the second timeframe. Reconstruction of axonal trajectories over discontinuous signal peaks revealed axonal branching forks that could not be reconstructed in the previous steps.
Performance of the method for reconstructing axonal functional morphology
The Bayes optimal template-matching technique was used to estimate the algorithm’s performance for detecting axonal AP peaks (Figure 5). The technique can reliably discriminate axonal APs from the background noise (Radivojevic et al., 2017) and therefore provide a ground truth for the validation of the algorithm. Herein, template-matching was used to discriminate between ‘false peaks’ caused by the noise and ‘true peaks’ that resulted from axonal electrical activity (Figure 5A). Discrimination criteria were based on similarities (‘matching’) between waveforms of averaged signals (‘templates’) and corresponding single trials. High similarities were expected in cases of accurate axonal signals (APs) and low similarities in cases where signals were derived from the background noise. We constructed templates for each signal peak in an electrical image. Constructed templates were next compared with waveforms of corresponding single trials, and a percentage of ‘matching’ trials was computed for each of the peaks. Signal peaks that yielded a match of >70% were classified as true peaks. Classified (true and false) peaks were used as the ground truth to estimate performances of the algorithm for detecting axonal AP peaks (Figure 5B and C). Data obtained from 20 cells (10 cortical and 10 motor neurons) were used for this analysis.
The performance of the algorithm for the detection of axonal AP peaks is shown in Figure 5B. We estimated the performance of the three detection steps (see Figure 3) with varying threshold levels. We were able to detect 45%, 74%, and 98% of the actual peaks after the first, second, and third steps, respectively. We observed no false peak detections for thresholds set to 9, 2, and 1 STD of the noise in first, second, and third steps, respectively. The detectability of AP peaks across axonal arbors of cortical and motor neurons is shown in Figure 5C.
Axonal conduction trajectories, obtained from stimulation-triggered neuronal activity, were applied to estimate the performance of the algorithm for tracking axonal AP peaks (Figure 6). The key concepts of the stimulation protocol used in this study are presented in Figure 6—figure supplement 1 and described in the Methods section. We used targeted stimulation to reveal and verify the same axon’s conduction trajectory by altering its conduction’s direction (Figure 6A, Figure 3—video 3). Stimulation-triggered APs were mapped spatially across all microelectrodes and temporally over discrete timeframes (with 50 μs inter-frame interval) to produce the movie. Observing the spatial movements of AP peaks in consecutive movie frames enabled us to track axonal conduction in different directions visually and to reconstruct axonal conduction trajectories (as shown in Figure 6A and Figure 6—videos 1 and 2).
The algorithm’s performance for tracking axonal conduction is shown in Figure 6B. We estimated the performance of the three tracking steps (see Figure 4) with varying diameters of corresponding interconnection areas. We could interconnect 70%, 85%, and 91% of the mapped AP peaks after the first, second, and third steps, respectively. We observed no false links when using the interconnection areas with diameters of 100, 200, and 400 μm in first, second, and third steps, respectively. The efficiency of the signal tracking across axonal arbors of cortical and motor neurons is shown in Figure 6C.
Functional morphologies of cortical and spinal axons
The reconstructed functional morphologies revealed axonal arbors of cortical and motor neurons and indicated positions of their AISes and somas (Figure 7). They further provided insights into total axonal lengths, the spatial distribution of axonal branching points and terminals, lengths of inter-branching segments and their branching orders. We used electrical images to extract areas occupied by axonal electrical activity (‘active areas’).
Representative examples of functional morphologies reconstructed for cortical and motor neurons are shown in Figure 7A, and displayed in Figure 7—video 1. The reconstructed arbor of the cortical axon yielded a total length of 27.12 mm, comprising 101 inter-branching segments and 53 axon terminals. Axonal electrical activity was detected on 7295 electrodes, occupying an active area of 2.23 mm2. The axonal arbor of the motor neuron yielded a total length of 15.26 mm, comprising 43 inter-branching segments and 23 axon terminals. Axonal activity was detected on 6663 electrodes, occupying an active area of 2.04 mm2. We reconstructed functional morphologies for 50 cortical and 50 motor neurons and analyzed morphological features of their axonal arbors, a total length of 1.04 m and 0.81 m of cortical and spinal axons, respectively.
We found that cortical and motor neurons grew axons comparable in their total lengths. However, ratios between axonal lengths and their corresponding active areas were significantly higher for cortical than spinal axons (Figure 7B). The average axonal lengths were 16.73±1.20 and 14.63±0.88 mm for cortical and motor neurons, respectively (p=0.25). The average sizes of active areas were 1.38±0.08 and 1.61±0.08 mm2 for cortical and motor neurons, respectively (p=0.23). The average length-to-area quotients were 11.72±0.32 and 9.49±0.24 mm–1 for cortical and motor neurons, respectively (p<10–6).
We found cortical axons to develop branching points in more proximal parts of their arbors and to form shorter inter-branching segments as compared to spinal axons (Figure 7C). Average axial distances of axonal branching points were 0.27±0.01 and 0.43±0.01 mm for cortical and motor neurons, respectively (p<10–6). Average lengths of axonal inter-branching segments were 0.34±0.01 and 0.48±0.01 mm for cortical and motor neurons, respectively (p<10–6).
We found significantly more axon terminals in cortical than in spinal axons; however, spinal axons projected their terminals at significantly greater distances as compared to cortical axons (Figure 7D). Average numbers of axon terminals were 29.92±1.64 and 17.82±1.03 for cortical and motor neurons, respectively (p<10–6). Average axial distances of axon terminals were 1.17±0.02 and 2.41±0.04 mm for cortical and motor neurons, respectively (p<10–6).
Signal amplitude fluctuations during axonal conduction
In addition to exposing the axonal structure, functional morphologies also reveal APs as propagating across axonal arbors (Figure 8). They allow observation of AP waveforms at any axonal location and mapping of AP amplitudes across entire arbors. Representative AP waveforms extracted from functional morphologies of cortical and spinal axons are shown in Figure 8A, and displayed in Figure 8—video 1. We investigated fluctuations in AP amplitudes across functional morphologies of 50 cortical and 50 spinal axons (same neurons as in Figure 7). The analysis included 45,232 data points obtained from 1.044 m of cortical axons and 36,286 data points obtained from 0.812 m of spinal axons.
We found that APs recorded from the most proximal parts of axons had much higher amplitudes than APs recorded from more distal axonal locations (Figure 8A–C). To further compare AP amplitudes across different cortical and spinal axon domains, we segregated the reconstructed morphologies into proximal and distal axons, primary trunks, and axonal trees. Proximal axons pertain to locations found within the first 0.2 mm of axial length; all other locations (axial distance >0.2 mm) were considered as distal (see Figure 8B, right). Primary trunk entails locations found in the axonal domain between the soma and the first branching fork, all other locations were assigned to the axonal tree (see Figure 8C, right).
APs recorded from proximal axons had significantly higher amplitude in cortical (44.58±4.81 µV) than in motor neurons (19.07±2.53 µV; p<10–6). On the contrary, APs recorded from distal axons had significantly smaller amplitude in cortical (1.50±0.10 µV) than in motor neurons (3.30±0.38 µV; p<10–6; Figure 8B).
Similarly, APs obtained from primary axonal trunks had significantly higher amplitude in cortical than motor neurons, whereas, APs obtained from axonal trees had significantly lower amplitude in cortical than motor neurons (Figure 8C). Average amplitudes of APs recorded from primary axonal trunks were 33.29±4.78 and 15.06±2.78 µV for cortical and motor neurons, respectively (p=0.002). Average amplitudes of APs recorded from axonal trees were 1.25±0.13 and 2.17±0.30 µV for cortical and motor neurons, respectively (p=0.024).
Temporal dynamics of axonal conduction
Functional morphologies carry information about times at which APs arrived at any axonal site and, as such, provide direct insights into temporal dynamics of axonal conduction (Figure 9). Axonal functional morphology reveals time at which axonal conduction is initiated (‘initial time’) and allows to obtain times at which axonal APs arrive at axon terminals (‘arrival time’). It enables extraction of the interval between the earliest and the latest arrival time (‘interval of arrivals’) and inspection of the total duration of the axonal conduction (‘active timespan’) (see Figure 9A, left). Functional morphologies allow for the inspection of AP conduction dynamics across entire axonal arbors and over individual axonal paths (see Figure 9A, right). Finally, structural and temporal parameters from functional morphologies enable computation of conduction velocities for any axonal segment.
We investigated temporal aspects of axonal conduction across functional morphologies of 50 cortical and 50 spinal axons (same neurons as in Figure 7 and Figure 8). The analysis included 45,232 data points obtained from 1.044 m of cortical axons and 36,286 data points obtained from 0.812 m of spinal axons.
Cortical neurons had significantly shorter active timespans and higher synchronization of the arrival times as compared to motor neurons (Figure 9B). Active timespans were 4.60±0.33 and 7.40±0.32 ms for cortical and motor neurons, respectively (p<10–6). Intervals of signal arrivals at axon terminals were 3.44±0.31 and 5.40±0.34 ms for cortical and motor neurons, respectively (p=0.0024). Variance of signal arrival times were 0.95±0.23 and 2.50±0.35 ms2 for cortical and motor neurons, respectively (p<10–4).
We next computed axonal conduction velocities using structural and temporal parameters obtained from the functional morphologies. Conduction velocities were calculated from axial distances versus AP propagation times (Figure 9A, right), and were computed over 100 µm long axonal chunks stepped by 17.5 µm.
We found cortical axons to exhibit significantly slower and more uniform conduction velocities compared to spinal axons (Figure 9C). Average conduction velocities were 0.46±0.01 and 0.56±0.02 m/s for cortical and motor neurons, respectively (p<10–4). Variances of conduction velocities were 0.03±0.002 and 0.04±0.003 m/s for cortical and motor neurons, respectively (p<10–4).
Discussion
We developed a method for noninvasive functional imaging of unmyelinated mammalian axons in vitro using the CMOS-based HD-MEA system (Figure 1, Figure 1—videos 1–2). The method yields an axonal ‘functional morphology’, comprising multidimensional data derived from extracellular APs recorded during axonal conduction. Functional morphology contains information about axonal conduction trajectories mapped at high spatial and temporal resolution (Figure 3—videos 1–3 and Figure 4—videos 1–3). It allows the reconstruction of axonal ‘electrical morphology’ at different timepoints during the conduction and, at the same time, to expose waveforms of the propagating APs (Figure 8, Figure 8—video 1).
The presented method has been developed and validated on primary rodent cortical and motor neurons cultured directly over HD-MEA surfaces (Figure 2—figure supplements 1–3). The cultures matured within 2 weeks and exhibited spontaneous electrical activities (Figure 2—figure supplement 6). They efficiently adhered to HD-MEA surfaces and formed monolayers that varied in local thicknesses and cell densities (Figure 2—figure supplement 4). Neurons developed complex axonal arbors and provided a tight interface between axons and microelectrodes. Because axons grew in a 3D manner, their distances from the HD-MEA surfaces varied locally (Figure 2—figure supplement 5), but remained within the sensing range of the microelectrodes (Obien et al., 2019).
Spontaneous electrical activity enabled us to map individual neurons in the cultures (Figure 2—figure supplement 6) and expose the spatiotemporal distribution of APs recorded across their axonal arbors (Figure 2—figure supplement 7, Figure 2—video 3). Mapping individual neurons was possible, thanks to high-amplitude APs recorded at their AISes (Bakkum et al., 2019; Radivojevic et al., 2016). The AIS signals broadly surpass the background noise and hence could be easily detected (Radivojevic et al., 2017). Because they colocalized with the AIS, local maxima found within activity maps (Figure 2—figure supplement 6A) indicate proximal regions of individual neurons in the network (Bakkum et al., 2019; Radivojevic et al., 2016). Owing to the high-density arrangement of recording microelectrodes (Figure 2—figure supplement 1; Ballini et al., 2014), spike-sorting enabled discerning mixed signals recorded from neighboring neurons (Figure 2—figure supplement 6B) and thereby extracting APs arising from individual neurons (Figure 2—figure supplement 6C). Averaging of extracted signals over all electrodes revealed the spatiotemporal distribution of a single neuron’s activity. Such signal representation was referred to as an ‘axonal electrical image’ (Figure 2—figure supplement 7, Figure 2—video 3) and it carries information about neuronal subcellular elements (Bakkum et al., 2019; Radivojevic et al., 2016). Thus, the largest and the earliest signals in the footprint arise from the AIS, much smaller and fast propagating signals arise from axons, and slow back-propagating signals arise from the soma and proximal dendrites (Bakkum et al., 2019; Radivojevic et al., 2016). Thanks to the low-noise CMOS design, the HD-MEA chip allows detection of AP propagation across large portions of entire axonal arbors (Radivojevic et al., 2017). The soma and proximal dendrites provide minor contributions to the electrical image, typically masked by larger AIS signals (Bakkum et al., 2019). Owing to their low-amplitude extracellular signals, distal dendrites do not seem to be detectable with the HD-MEA system used here.
Adaptive thresholding applied to axonal electrical images enabled us to map extracellular APs propagating across axonal arbors (Figure 3, Figure 3—videos 1–3). Signal peaks were traced at high spatial and temporal resolution, thanks to the dense arrangement of the electrodes and high-frequency sampling rate (Ballini et al., 2014). The thresholding scheme was constructed to detect signals with various amplitudes while minimizing detection errors. Simple planar thresholds, set to 9 STD of the background noise, were used to detect high-amplitude signal peaks (Figure 3A, Figure 3—video 1). Since these thresholds provided detection high above the level of the background noise, they could be applied across an entire array while yielding no detection errors. They, however, failed to detect low-amplitude signal peaks and only provided fragmentary detection. Instead, spatiotemporal confinement of low-level thresholds, contingent on coordinates of high-amplitude peaks, enabled the mapping of low-amplitude signals (Figure 3B, Figure 3—video 2). This was possible, thanks to the very nature of axonal conduction, which can be represented using a simplified Markov chain model (Gagniuc, 2017). Namely, because signal propagation in unmyelinated axons is continuous, it represents a cascade of successive events in which the probability of each event depends solely on the state attained in the preceding event. Therefore, there was a high probability of detecting signal peaks in the immediate spatiotemporal proximity of previously mapped peaks. While confining the reach of the thresholds to local areas minimized the influence of the background noise, temporal confinement to preceding and following timeframes excluded signals that occurred much earlier and later than the reference peak.
Further relaxation of the threshold confinements enhanced the detection of low-amplitude signals but inevitably increased the risk of detection errors (Figure 3C, Figure 3—video 3). In principle, our adaptive thresholding scheme followed a greedy algorithm paradigm (Cormen et al., 2022), and as such, it is vulnerable to detection errors. Generally, a greedy algorithm builds up a global solution stepwise by making the locally optimal choice at each step. However, erroneous solutions made at earlier steps can proliferate and thereby compromise the global solution (Gutin et al., 2002). Therefore, carefully tuned thresholding parameters are key to the optimal detection of signal peaks.
The approach of multiple-moving object tracking enabled us to reconstruct axonal conduction trajectories through iterative associations of the mapped signal peaks (Figure 4, Figure 4—videos 1–3). The tracking algorithm was designed to estimate inter-frame correspondences between the mapped peaks in a data-driven fashion. It operates within predefined spatiotemporal boundaries, but also directly learns data association criteria from sequences reconstructed along the tracking process. Thanks to strict boundaries, the algorithm discarded trajectories yielding conduction velocities that are unlikely for mammalian axons. While such boundaries ensure the exclusion of slow-propagating dendrites (Bakkum et al., 2019), they could also omit super-fast saltatory conduction inherent to myelinated axons (Debanne, 2004), which were most likely not present in our preparations. Owing to the data-driven refinement of the association criteria, the algorithm allowed customization of the tracking process for a specific neuron. Such customization provides applicability of the algorithm to axons with different conduction velocities, allowing observation of variations between neurons of the same or different types.
Considering the Markov model of axonal conduction, AP peaks mapped in immediate spatiotemporal proximities can be associated directly using a greedy nearest neighbor matching (Figure 4A, Figure 4—video 1). Due to the limited reach of the matching caliper, direct interconnection enabled only reconstruction of short fragments of the axonal conduction trajectory. Reconstructed fragments, in return, enabled estimation of the conduction velocity and thereby refining the association criteria. Because conduction velocities vary across a single axon (Radivojevic et al., 2017; Bakkum et al., 2013), greedy matching was limited to peaks mapped closely over sections of the trajectory with relatively constant conduction velocities. Relaxation of the association criteria could potentially extend the reach of the matching. However, it would inevitably introduce erroneous interconnections between spatially segregated peaks stretched over notably faster sections of the trajectory.
Skeletonization of an axonal electrical image enabled reproduction of conduction paths between distant peaks (Figure 4B, Figure 4—video 2). Because skeletal remnants reflect the spatiotemporal landscape of the axonal signals, they can optimize the association of mapped peaks graphically, thus sidestepping complex heuristic techniques needed for maintaining computational tractability. In general, topology and structural complexity of skeletal remnants depend on spatiotemporal boundaries within which the axonal signals were skeletonized. Narrow boundaries encompassing two consecutive timeframes captured brief sequences of the axonal conduction and yielded skeletal remnants with fairly simple structure. Simple remnants enabled interconnection of the peaks across consecutive timeframes, however, discontinuous peaks remained beyond their reach.
Skeletonization within wider spatiotemporal boundaries allowed interconnecting peaks across discontinuous timeframes (Figure 4C, Figure 4—video 3), but yielded more complex remnants, especially in regions where axons branched or intertwined. Complex remnants provided multiple possible solutions for the peak associations, thus aggravating the selection of an optimal conduction path between corresponding peaks. By restricting the skeletonization to three consecutive timeframes and using the refined association criteria, we could select the optimal path in most cases. However, further widening of the skeletonization boundaries increased the number of possible solutions exponentially and inevitably led to a deterioration of tracking performance.
The Bayes optimal template-matching technique was used to validate the algorithm for signal peak detection (Figure 5). We chose this technique because it has previously been demonstrated to provide reliable detection of extracellular APs recorded using the HD-MEA system (Radivojevic et al., 2017). It enables the detection of axonal APs within single recording trials and microsecond differences in axonal propagation (Radivojevic et al., 2017; Franke et al., 2015). Bayes optimal template-matching is sensitive enough to discriminate low-amplitude APs from the background noise (Radivojevic et al., 2017), and thus can be used to validate our algorithm.
Stimulation-triggered axonal activities allowed us to trace axonal conduction trajectories accurately and thus provided the ground truth for validation of our tracking algorithm (Figure 6, Figure 6—videos 1 and 2). We have previously shown that both spontaneous and stimulation-triggered neuronal activities recorded by HD-MEAs can be used to outline axonal conduction trajectories (Bullmann et al., 2019; Radivojevic et al., 2016). Moreover, we demonstrated a congruence between the reconstructed trajectories and actual axonal morphology revealed optically (Radivojevic et al., 2017; Bakkum et al., 2013; Bullmann et al., 2019; Radivojevic et al., 2016).
Previous research on HD-MEA systems provided methods for the electrical imaging of axons (Radivojevic et al., 2017; Bakkum et al., 2013; Müller et al., 2015; Bullmann et al., 2019; Radivojevic et al., 2016) but did not provide algorithms for the reconstruction of axonal morphology. Therefore, optical imaging of axonal morphology or manual tracking of AP propagation was necessary to assign the recorded signals to their respective axonal locations (Radivojevic et al., 2017; Bakkum et al., 2013; Bullmann et al., 2019; Radivojevic et al., 2016). Recent studies have proposed algorithms for deducing axonal arbors and propagation velocity based on extracellular APs (Yuan et al., 2020; Buccino et al., 2022). However, these algorithms have solely been validated on simulated morphologies, and no experimental validation has been done to assess their performances (Buccino et al., 2022). The present algorithm enables precise reconstruction of axonal electrical morphology and has been validated on experimental data obtained from cortical and motor neurons (see Figures 5 and 6). The algorithm’s performances were investigated, and optimal parameters for accurate reconstruction were identified. The algorithm’s adaptive functionality enables it to adjust and optimize its parameters in a data-driven fashion, which allows for customizing the tracking process for a specific neuron.
We investigated functional morphologies of cortical and spinal axons and found significant differences between their structures, AP amplitudes, and conduction dynamics (Figures 7—9), which potentially reflect biological hallmarks of different neuronal subtypes. Although comparable in overall length, cortical and spinal axons exhibited substantially different functional morphologies. Cortical axons had a more complex branching pattern than spinal axons, but spinal axons projected their terminals to much longer distances than cortical axons. The amplitude of APs recorded from spinal axons was higher than those recorded from cortical axons. However, when comparing signals recorded from the AIS, cortical neurons displayed significantly higher amplitude than motor neurons. Cortical neurons had shorter active timespans and more synchronized arrival times than motor neurons, while cortical axons had slower and more uniform conduction velocities than spinal axons.
The structural complexity of cortical axons was reflected in their extensive branching, followed by the abundance of the axon terminals that were projected locally and distally (Figure 7, Figure 7—video 1). Such morphological features of cortical axons can be attributed to their role in providing synaptic connectivity within local neuronal assemblies as well as across large neuronal networks, thus enabling higher-order operations performed in the brain (Buzsáki, 2010). Spinal motor neurons, however, have evolved much simpler morphology to provide fundamentally different biological roles. Motor neurons innervate and precisely control muscles in the periphery outside the CNS and are, for that purpose, equipped with the longest known axons in the body (Stifani, 2014). Our results obtained in in vitro conditions are consistent with such a notion. They indicate that spinal axons indeed have much simpler morphologies when compared to cortical axons, yet still retain relatively long total lengths (Figure 7, Figure 7—video 1). We found that spinal axons had moderately branched arbors, but projected their axon terminals at great distances, thus resembling morphological adaptations observed in vivo. Thanks to their extensively branched axons, cortical neurons enabled the delivery of APs to numerous axon terminals, surprisingly synchronously and within short timespans. On the contrary, spinal axons with less branches at their disposal required considerably more time to forward their APs to axon terminals in a less synchronized fashion (Figure 9B). In this study, the assessment of the structural complexity of cortical and spinal axons relies entirely on their functional morphologies. It would be valuable for future studies to compare the morphological features obtained through electrical and optical imaging.
We found that extracellular APs recorded from proximal axons near the putative AIS exposed disproportionately larger amplitudes than signals obtained from other axonal parts (Figure 8, Figure 8—video 1). It has been previously shown that the AIS is the dominant contributor to the neuron’s extracellular electrical landscape (Bakkum et al., 2019), which can be attributed to the high densities of voltage-gated ion channels expressed in the AIS (Kole et al., 2008; Hu et al., 2009; Lorincz and Nusser, 2010). Interestingly, we found that APs obtained from proximal axons had significantly higher amplitudes in cortical than in motor neurons (Figure 8B and C). The relatively low AP amplitudes at the AIS in motor neurons could result from the potentially disrupted organization of the AIS in unmyelinated motor neurons. Such disruption pertains to altered distribution and combination of voltage-gated ion channels as well as incomplete segregation of the AIS, para-AIS, and juxtapara-AIS compartments inherent to myelinated spinal axons (Duflocq et al., 2011). The diameter, length, and positioning of the AIS in relation to recording electrodes may also substantially impact shaping the extracellular AP profile (Bakkum et al., 2019; Kumar et al., 2022), which in turn can influence the AP amplitude. We found that signals obtained from distal axons and axonal trees had significantly lower amplitudes in cortical than motor neurons (Figure 8B and C). This difference could be explained by the fact that spinal axons have a considerably larger diameter than cortical axons (Saliani et al., 2017). As such, they engage a greater number of voltage-gated ion channels and create stronger transmembrane currents during an AP electrogenesis (Stein and Pearson, 1971). However, a comprehensive modeling study that incorporates biophysical, structural, and spatial parameters is necessary to explain the observed differences in AP amplitudes between cortical and spinal axons. Large diameters of spinal axons could also explain their significantly faster conduction velocities compared to cortical axons (Figure 9C). In general, axonal diameter and the presence of myelin sheaths are crucial factors that control conduction velocity in mammalian axons (Debanne et al., 2011; Waxman, 1980). On the one hand, conduction velocity of unmyelinated axons is proportional to the square root of the axon diameter (Hodgkin and Huxley, 1952). On the other hand, myelination provides a fundamentally different mechanism of AP propagation known as ‘saltatory’ conduction (Huxley and Stämpeli, 1949). In this case, myelin sheaths spatially restrict the distribution of voltage-gated ion channels to nodes of Ranvier and thus impose discontinuous AP propagation up to 100-fold faster than propagation in unmyelinated axons (Debanne et al., 2011; Saliani et al., 2017). The culture protocols used in this study do not support the formation of myelin sheaths around axons. Therefore, it is unlikely that complete myelination occurred in our cultures. Namely, the myelination of primary neurons in vitro requires either supplement that supports well-balanced development of both neurons and myelin-producing cells (Pang et al., 2012) or, in the case of motor neurons, co-culture with a feeder layer of Schwann cells (Hyung et al., 2015). In both cases, the timely addition of specific neurotrophic factors is required to initiate the myelination program. The culture media used in this study were not specifically designed to support the growth of myelin-producing cells, and no myelination-promoting factors were included in the media. The protocols did not involve co-culturing motor neurons with Schwann cells that would provide myelination. Moreover, motor neurons were grown in a culture medium containing NT3, a neurotrophic factor known to significantly inhibit myelination in vitro (Chan et al., 2001). Finally, we did not detect any saltatory high-speed signal conduction in any recorded neurons, suggesting that myelin sheaths were not formed around the axons. We, however, cannot exclude potential cases where axons were partially myelinated. Moreover, incomplete myelinization or discontinuous distribution of voltage-gated ion channels could explain relatively large variances of conduction velocities in spinal axons (Figure 9C).
Generally, standard deviations of the data presented in Figures 7—9 could result from morphological and biophysical differences between various neuronal subtypes that may have existed in our cultures. Thus, for example, developed cortical cultures comprise excitatory glutamatergic and inhibitory GABAergic neurons which are known to differ in size and morphological complexities (Björklund et al., 2010). Even greater diversity of neuronal subtypes can be found among motor neurons, including alpha, beta, and gamma motor neurons that are structurally specialized to innervate different muscle fiber types in the body (Stifani, 2014). However, another question is to what extent such specializations can be recapitulated in vitro and whether electrophysiological parameters observed in cultures allow for discriminating between neuronal subtypes.
The presented method for functional imaging of cortical and spinal axons provides direct insights into the biophysical properties of axonal conduction. It allows investigation of interdependences between axonal function and structure, noninvasively and over extended periods. Functional morphologies contain data obtained from large portions of axonal arbors and, as such, allow studying the fidelity of signal conduction in different axonal elements, including conduction failures at branching points (Debanne, 2004) and signal attenuation in tiny axon terminals (Faisal and Laughlin, 2007). They also allow the inspection of activity-dependent modulation of extracellular AP waveforms (Radivojevic et al., 2017; Lewandowska et al., 2016) and analysis of consequent changes in times at which modulated signals arrive at axon terminals (Radivojevic et al., 2017). Owing to its noninvasive nature, the method enables the study of structural and functional changes in axons over long periods. Thus, for example, comparative analysis of functional morphologies obtained at different timepoints allows investigation of biophysical changes during axonal outgrowth and pathfinding – developmental processes that guide axons toward specific targets and enable wiring within neuronal networks (Myers et al., 2011; Uesaka et al., 2005). Complemented with techniques for electrical microstimulation (see Figure 6—figure supplement 1), the method can be used to study activity-dependent axonal plasticity, including structural alteration such as the positional shift of the AIS (Grubb and Burrone, 2010; Kumar et al., 2022) as well as functional adaptations that involve changes in axonal conduction velocities (Bakkum et al., 2013). Because the HD-MEA system allows observing electrical activities across entire neuronal networks (see Figure 2—figure supplement 6), the method can serve as a complementary tool for studying the functional interplay between network dynamics and plastic adaptations in axons (Grubb and Burrone, 2010; Kuba et al., 2010).
In summary, the presented method is designed for specific tissue-on-a-chip platforms and may serve as an attractive tool for pharmacological and preclinical studies of neurodegenerative diseases (Ronchi et al., 2021). The key potential of the method lies in its ability to gain multilevel knowledge of neuronal functionality, which may serve as a valuable asset in drug discovery and safety pharmacology. The method holds great potential for developing disease-specific bioassays, especially considering the advent of human induced pluripotent stem cell technology that allows for utilizing patient-derived neurons.
Methods
Animal use
All experimental protocols were approved by the Uppsala Animal Ethical Committee under animal license C97/15 and follow the guidelines of the Swedish Legislation on Animal Experimentation (Animal Welfare Act SFS 2009:303) and the European Communities Council Directive (2010/63/EU).
HD-MEA and signal processing
A CMOS-based HD-MEA system (MaxWell Biosystems AG) was used for extracellular neuronal recording and stimulation (Figure 2—figure supplement 1). The array comprises 26,400 platinum microelectrodes (9.3×5.45 μm2) packed within an area of 3.85×2.10 mm2, providing a density of 3150 electrodes per mm2 (17.5 μm center-to-center pitch). Owing to a flexible switch matrix technology, up to 1024 readout and/or stimulation channels could be routed to the desired electrodes and reconfigured within a few milliseconds. On-chip circuitry was used to amplify (0–80 dB programmable gain), filter (high pass: 0.3–100 Hz, low pass: 3.5–14 kHz), and digitize (8-bit, 20 kHz) the neuronal signals. Digitized signals were sent to a field-programmable gate array board and further streamed to a host PC for real-time visualization and data storage. Recorded signals were up-sampled to 200 kHz following the Whitaker-Shannon interpolation formula. Python 3.7 and Matlab R2020a were used for data analysis and to design experimental protocols.
Cortical cultures
We used a culturing protocol for the long-term maintenance of neural cultures (Hales et al., 2010: Figure 2—figure supplement 2). Minor adaptations were introduced to the protocol to constrain culture growth to the sensing area of the array and maintain optimal growth media conditions during long-term experimentation (Radivojevic et al., 2016). Cortices from embryonic day 18 Sprague Dawley rat were dissociated enzymatically in trypsin supplemented with 0.25% EDTA (Thermo Fisher Scientific) and physically by trituration. The remaining cell aggregates and debris were filtered using 40 µm cell strainer (Corning). For cell adhesion, a layer of 0.05% polyethyleneimine (Sigma-Aldrich) in borate buffer (Chemie Brunschwig), followed by a layer of 0.02 mg/ml human recombinant laminin (BioLamina) in DPBS containing Ca2+ and Mg2+ (Thermo Fisher Scientific) was deposited on the electrode array. To constrain culture growth to the electrode array, a cell drop containing ~50,000 cells and covering ~4 mm2 was seeded in the center of the array. The plating media were changed to growth media after 24 hr and regularly changed every 6 days. Plating media consisted of Neurobasal, supplemented with 10% horse serum (HyClone), 0.5 mM GlutaMAX, and 2% B27 (Thermo Fisher Scientific). The growth media consisted of DMEM, supplemented with 10% horse serum (HyClone), 0.5 mM GlutaMAX, and 1 mM sodium pyruvate (Thermo Fisher Scientific). Cultures were maintained inside an incubator under controlled environmental conditions (36°C and 5% CO2). The culturing chambers were sealed with an ~1 mm layer of light mineral oil (Sigma-Aldrich) floating above the growth medium. The sealing provided selective permeability to gases, such as O2 and CO2, and prevented evaporation and consequent changes in the growth media’s osmolarity during long-term experiments.
Motor neuron cultures
Protocol for growing motor neuron cultures on HD-MEA surface is illustrated in Figure 2—figure supplement 3. Spinal cords isolated from embryonic day 14 Sprague Dawley rat were dissociated enzymatically in trypsin supplemented with 0.25% EDTA (Thermo Fisher Scientific) and physically by trituration. The remaining cell aggregates and debris were filtered out using 40 µm cell strainer (Corning). A density gradient medium composed of 15% OptiPrep (Sigma-Aldrich) in Leibovitz’s L-15 medium (Thermo Fisher Scientific) was used to fractionate motor neurons. For cell adhesion, a layer of 0.05% polyethyleneimine (Sigma-Aldrich) in borate buffer (Chemie Brunschwig), followed by a layer of 0.02 mg/ml human recombinant laminin (BioLamina) in DPBS containing Ca2+ and Mg2+ (Thermo Fisher Scientific) was deposited on the electrode array. To constrain culture growth to the electrode array, a cell drop containing ~50,000 cells and covering ~4 mm2 was seeded in the center of the array. The growth media consisted of Neurobasal supplemented with 20% MyoTonic differentiation medium (Cook MyoSite), 1% fetal bovine serum, 2% horse serum, 2% B27, 1% antibiotic-antimycotic, 0.2 mM gentamicin, 0.7 mM L-glutamine, 1.5 mM sodium pyruvate (Thermo Fisher Scientific), 0.8 nM brain-derived neurotrophic factor, 0.2 nM glial cell line-derived growth factor, 0.2 nM ciliary neurotrophic factor, 1.5 nM neurotrophin 3 (NT3), 1.4 nM neurotrophin 4 (NT4, R&D Systems), and 1.3 nM human insulin-like growth factor 1 (hIGF-1, Thermo Fisher Scientific). The growth media were changed after 24 hr and further regularly changed every 4 days. Cultures were maintained inside an incubator under controlled environmental conditions (36°C and 5% CO2). The culturing chambers were sealed with an ~1 mm layer of light mineral oil (Sigma-Aldrich) floating above the growth medium.
Immunocytochemistry
Neuronal cultures were fixed in 4% paraformaldehyde (Thermo Fisher Scientific) in PBS (Sigma-Aldrich) at pH 7.4 for 15 min at room temperature, washed twice with ice-cold PBS, permeabilized with 0.25% Triton X-100 (Sigma-Aldrich) in PBS for 10 min and washed three times in PBS. Fixed cultures were exposed to phosphate-buffered saline with Tween 20 (1% bovine serum albumin and 0.1% Tween 20 in PBS; Sigma-Aldrich) for 30 min to prevent unspecific binding of antibodies. The primary antibodies Anti-Map2 (Abcam, ab5392), Anti-β-III tubulin (Abcam, ab18207), and Anti-mCherry (Abcam, ab125096), diluted in phosphate buffered saline with Tween 20 to ratios of 1:1000, 1:500, and 1:500 respectively, were added and left overnight at 4°C on a shaker. Cultures were washed three times in PBS for 5 min each time on the shaker. The secondary antibodies Alexa Fluor 488 (Abcam, ab150173), Alexa Fluor 488 (Abcam, ab150073), and Alexa Fluor 594 (Abcam, ab150116) diluted to ratio of 1:200 in PBS with 1% BSA were added and left for 60 min in the dark at room temperature. Samples were washed three times in PBS for 5 min in the dark and then stored at 4°C. Cultures were immunostained with anti-β-III tubulin (Figure 2—figure supplements 4 and 5, Figure 1—videos 1 and 2), anti-Map2 (Figure 2—figure supplements 6 and 7, Figure 2—videos 1–3), and a combination of anti-mCherry and anti-β-III tubulin antibodies (Figure 2—figure supplement 5).
Live imaging
Live-cell visualization of neurons was performed by transfection using pAAV-hSyn-mCherry plasmid from Karl Deisseroth (Addgene plasmid # 114472) and Lipofectamine 3000 (Thermo Fisher Scientific), jetPRIME (Polyplus) or TurboFect (Thermo Fisher Scientific) transfection reagent following the manufacturer’s protocols. Micrographs of cortical and motor neurons transfected using Lipofectamine 3000 are shown in Figure 2—figure supplement 5.
Microscopy and 3D image reconstruction
A Nikon Eclipse LVDIA-N microscope, Nikon DS-Fi2 camera, and the Nikon NIS-Elements imaging software were used to produce micrographs. Epifluorescence microscopy was used to collect Z-stack image series from immunostained cultures shown in Figure 2—figure supplements 4 and 5 and Figure 2—videos 1 and 2. ImageJ and Matlab custom designed codes were used to create 3D surface plots (Figure 2—figure supplement 4) and to reconstruct 3D neuronal morphologies (Figure 2—figure supplement 5) based on the intensities of immunofluorescent signals.
Network-wide activity mapping
Spontaneous extracellular APs were sampled across the entire microelectrode array by sequential scanning over 28 recording configurations. Up to 1024 randomly selected electrodes recorded neuronal activity in each configuration for 2 min. The average voltage traces, recorded by each electrode, were used to reconstruct the network-wide activity map (Figure 2—figure supplement 6A). Since the largest extracellular APs occur near the AIS, and APs arising from axonal arbors have much smaller amplitude, regions with high-amplitude APs in the activity maps indicated the locations of the AISes (Radivojevic et al., 2017; Bakkum et al., 2019; Radivojevic et al., 2016).
Electrical identification of individual neurons
Simultaneous access to signals arising from the AIS region enabled us to reconstruct the spatiotemporal distribution of extracellular APs, referred to as the ‘electrical footprint’ (Figure 2—figure supplement 6C). To obtain these signals, we used high-density recording configurations covering the AIS locations revealed in the activity maps (Figure 2—figure supplement 6A). In each configuration, blocks of 13×13 electrodes were connected to readout channels to sample neuronal activity for 2 min. Recorded signals were sorted by using the Spyking-Circus algorithm (Yger et al., 2018), and electrical footprints were reconstructed by using custom-designed Matlab code. Because the first-occurring (initial) trace found in the electrical footprint colocalizes with the neuron’s AIS, it was used to trigger the averaging of voltage traces recorded across an entire axonal arbor (see below).
Electrical imaging of axonal arbors
Array-wide averaging of voltage traces, synchronized with the initial trace (recorded from the AIS), reveals the spatiotemporal distribution of extracellular APs across an entire axonal arbor (Figure 2—figure supplement 7 and Figure 2—video 3). The first step in obtaining these data was selecting 9 electrodes that were closest to the putative AIS. We next designed multiple recording configurations covering the entire array – in each configuration, 9 of the 1024 readout channels were routed to the 9 preselected electrodes, and remaining available channels were routed to randomly selected electrodes. Each configuration was used to sample neuronal activity during 2 min. Signals recorded by the 9 preselected electrodes in each configuration were sorted using the Spyking-Circus algorithm (Yger et al., 2018). Timestamps of the sorted signals were used to trigger the averaging of voltage traces recorded across all other electrodes in the array. The spatiotemporal distribution of averaged signals was reconstructed using a custom-designed Matlab code.
Electrical stimulation of individual neurons in the network
Stimulation performances of the HD-MEA system allow for stimulating any neuron in the culture with subcellular spatial precision (Radivojevic et al., 2017; Bakkum et al., 2013; Radivojevic et al., 2016; Ronchi et al., 2019). Key concepts of the stimulation protocols that were used in this study are presented in Figure 6—figure supplement 1. Extracellular electrical stimulations targeted at predefined axonal locations allow to experimentally control the direction of the axonal conduction (Radivojevic et al., 2016; Figure 6—figure supplement 1A). Electrical stimulation directed at the AIS was used to elicit orthodromic neuronal activation (Figure 6 – figure supplement 1A, left). Stimulation directed at distal axons was used to elicit antidromic neuronal activation (Figure 6 – figure supplement 1A, right). All stimulation protocols utilized balanced positive-first biphasic voltage pulses, with phase durations of 200 μs, because of their proven effectiveness in electrical stimulation (Wagenaar et al., 2004). Stimuli were applied to one electrode at a time at a frequency of 4 Hz. Neuronal responses to orthodromic stimulations were estimated by observing APs recorded from distal axons (Figure 6 – figure supplement 1A, left). Responses to antidromic stimulation were estimated by observing APs recorded from proximal axons and the AIS (Figure 6 – figure supplement 1A, right). To find optimal parameters for stimulation of individual neurons in the network, we used approaches thoroughly described in our previous study (Radivojevic et al., 2016). In brief, we applied neuron-wide stimulation over a range of voltages to reveal sites with the lowest activation threshold (Figure 6—figure supplement 1B). Stimulation was applied at 4 Hz for voltages from ±10 to ±300 mV, with steps of ±10 mV. Each stimulation voltage was applied 60 times per site, and activation thresholds were defined as the minimum voltage to trigger an AP in 100% of the trials. To get more detailed excitability profiles, the most sensitive sites were then stimulated with voltages stepped by ±1 mV (Figure 6—figure supplement 1C).
Statistical analysis
All quantitative data presented in Figures 7—9 are expressed as mean ± SEM. Numbers of biological replicates (N) used for analyses presented in Figures 5—9 are denoted in the figures and are also stated in their legends. We used non-parametric tests for comparing distributions of parameters presented in Figures 7—9, since normal distribution of the underlying data could not be determined unequivocally. The two-sided Mann-Whitney U-test was applied and a p-value <0.05 was considered significant. Logarithmic scaling of graphs shown in Figure 8B and C was used to display a wide range numerical data in a compact way.
Data availability
Datasets generated for this study, entitled 'Electrical imaging of cortical and spinal axons using high-density microelectrode arrays', are stored in Dryad.
-
Dryad Digital RepositoryElectrical imaging of cortical and spinal axons using high-density microelectrode arrays.https://doi.org/10.5061/dryad.gxd2547r1
References
-
Axonal ComputationsFrontiers in Cellular Neuroscience 13:413.https://doi.org/10.3389/fncel.2019.00413
-
Analog signalling in mammalian cortical axonsCurrent Opinion in Neurobiology 18:314–320.https://doi.org/10.1016/j.conb.2008.08.006
-
A 1024-Channel CMOS microelectrode array with 26,400 electrodes for recording and stimulation of electrogenic cells in vitroIEEE Journal of Solid-State Circuits 49:2705–2719.https://doi.org/10.1109/JSSC.2014.2359219
-
High plasticity of axonal pathology in Alzheimer’s disease mouse modelsActa Neuropathologica Communications 5:14.https://doi.org/10.1186/s40478-017-0415-y
-
An automated method for precise axon reconstruction from recordings of high-density micro-electrode arraysJournal of Neural Engineering 19:026026.https://doi.org/10.1088/1741-2552/ac59a2
-
Large-scale mapping of axonal arbors using high-density microelectrode arraysFrontiers in Cellular Neuroscience 13:404.https://doi.org/10.3389/fncel.2019.00404
-
Axon degeneration in Parkinson’s diseaseExperimental Neurology 246:72–83.https://doi.org/10.1016/j.expneurol.2012.01.011
-
Information processing in the axonNature Reviews. Neuroscience 5:304–316.https://doi.org/10.1038/nrn1397
-
Bayes optimal template matching for spike sorting – combining fisher discriminant analysis with optimal filteringJournal of Computational Neuroscience 38:439–459.https://doi.org/10.1007/s10827-015-0547-7
-
Traveling salesman should not be greedy: domination analysis of greedy-type heuristics for the TSPDiscrete Applied Mathematics 117:81–86.https://doi.org/10.1016/S0166-218X(01)00195-0
-
How to culture, record and stimulate neuronal networks on micro-electrode arrays (MEAs)Journal of Visualized Experiments e2056.https://doi.org/10.3791/2056-v
-
A quantitative description of membrane current and its application to conduction and excitation in nerveThe Journal of Physiology 117:500–544.https://doi.org/10.1113/jphysiol.1952.sp004764
-
Evidence for saltatory conduction in peripheral myelinated nerve fibresThe Journal of Physiology 108:315–339.https://doi.org/10.1113/jphysiol.1949.sp004335
-
Applicability of independent component analysis on high-density microelectrode array recordingsJournal of Neurophysiology 108:334–348.https://doi.org/10.1152/jn.01106.2011
-
Tracking axon initial segment plasticity using high-density microelectrode arrays: A computational studyFrontiers in Neuroinformatics 16:957255.https://doi.org/10.3389/fninf.2022.957255
-
Regulation of axonal outgrowth and pathfinding by integrin-ECM interactionsDevelopmental Neurobiology 71:901–923.https://doi.org/10.1002/dneu.20931
-
Axonal noise as a source of synaptic variabilityPLOS Computational Biology 10:e1003615.https://doi.org/10.1371/journal.pcbi.1003615
-
Single-cell electrical stimulation using CMOS-Based high-density microelectrode arraysFrontiers in Neuroscience 13:208.https://doi.org/10.3389/fnins.2019.00208
-
Axon and myelin morphology in animal and human spinal cordFrontiers in Neuroanatomy 11:129.https://doi.org/10.3389/fnana.2017.00129
-
Firing behaviour in a stochastic nerve membrane model based upon the Hodgkin-Huxley equationsActa Physiologica Scandinavica 107:343–363.https://doi.org/10.1111/j.1748-1716.1979.tb06486.x
-
Predicted amplitude and form of action potentials recorded from unmyelinated nerve fibresJournal of Theoretical Biology 32:539–558.https://doi.org/10.1016/0022-5193(71)90155-x
-
Motor neurons and the generation of spinal motor neuron diversityFrontiers in Cellular Neuroscience 8:293.https://doi.org/10.3389/fncel.2014.00293
-
Effective parameters for stimulation of dissociated cultures using multi-electrode arraysJournal of Neuroscience Methods 138:27–37.https://doi.org/10.1016/j.jneumeth.2004.03.005
-
Channel noise in neuronsTrends in Neurosciences 23:131–137.https://doi.org/10.1016/s0166-2236(99)01521-0
-
Past and future of analog-digital modulation of synaptic transmissionFrontiers in Cellular Neuroscience 13:160.https://doi.org/10.3389/fncel.2019.00160
Article and author information
Author details
Funding
Vetenskapsrådet (2016-02184)
- Anna Rostedt Punga
Vetenskapsrådet (2014-02048)
- Anna Rostedt Punga
Vetenskapsrådet (2014-07603)
- Anna Rostedt Punga
Vetenskapsrådet (2020-02040)
- Anna Rostedt Punga
Göran Gustafssons Stiftelse för Naturvetenskaplig och Medicinsk Forskning
- Anna Rostedt Punga
Hjärnfonden
- Anna Rostedt Punga
Olle Engkvists Stiftelse (197-0235)
- Anna Rostedt Punga
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Olga Netsyk, Bryndis Birnir, Tanel Punga, and Evgenii Bogatikov for help with experiments; Markus de Ruijter and Stefan Engblom for help with software developments; Marta K Lewandowska and Yu-Fang Huang for critical discussions. This study was supported by Swedish Research Council grant 2016-02184, Swedish Research Council grant 2014-02048, Swedish Research Council grant 2014-07603, Swedish Research Council grant 2020-02040, Göran Gustafsson foundation for medical research, Hjärnfonden (the Swedish Brain Foundation), Olle Engkvist Byggmästare grant 197-0235, Margarethahemmets foundation and Bissen Brainwalk Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Ethics
All experimental protocols were approved by the Uppsala Animal Ethical Committee under animal license C97/15 and follow the guidelines of the Swedish Legislation on Animal Experimentation (Animal Welfare Act SFS 2009:303) and the European Communities Council Directive (2010/63/EU).
Copyright
© 2023, Radivojevic and Rostedt Punga
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,904
- views
-
- 147
- downloads
-
- 2
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Computational and Systems Biology
- Neuroscience
Diffusional kurtosis imaging (DKI) is a methodology for measuring the extent of non-Gaussian diffusion in biological tissue, which has shown great promise in clinical diagnosis, treatment planning, and monitoring of many neurological diseases and disorders. However, robust, fast, and accurate estimation of kurtosis from clinically feasible data acquisitions remains a challenge. In this study, we first outline a new accurate approach of estimating mean kurtosis via the sub-diffusion mathematical framework. Crucially, this extension of the conventional DKI overcomes the limitation on the maximum b-value of the latter. Kurtosis and diffusivity can now be simply computed as functions of the sub-diffusion model parameters. Second, we propose a new fast and robust fitting procedure to estimate the sub-diffusion model parameters using two diffusion times without increasing acquisition time as for the conventional DKI. Third, our sub-diffusion-based kurtosis mapping method is evaluated using both simulations and the Connectome 1.0 human brain data. Exquisite tissue contrast is achieved even when the diffusion encoded data is collected in only minutes. In summary, our findings suggest robust, fast, and accurate estimation of mean kurtosis can be realised within a clinically feasible diffusion-weighted magnetic resonance imaging data acquisition time.
-
- Neuroscience
The retina transforms patterns of light into visual feature representations supporting behaviour. These representations are distributed across various types of retinal ganglion cells (RGCs), whose spatial and temporal tuning properties have been studied extensively in many model organisms, including the mouse. However, it has been difficult to link the potentially nonlinear retinal transformations of natural visual inputs to specific ethological purposes. Here, we discover a nonlinear selectivity to chromatic contrast in an RGC type that allows the detection of changes in visual context. We trained a convolutional neural network (CNN) model on large-scale functional recordings of RGC responses to natural mouse movies, and then used this model to search in silico for stimuli that maximally excite distinct types of RGCs. This procedure predicted centre colour opponency in transient suppressed-by-contrast (tSbC) RGCs, a cell type whose function is being debated. We confirmed experimentally that these cells indeed responded very selectively to Green-OFF, UV-ON contrasts. This type of chromatic contrast was characteristic of transitions from ground to sky in the visual scene, as might be elicited by head or eye movements across the horizon. Because tSbC cells performed best among all RGC types at reliably detecting these transitions, we suggest a role for this RGC type in providing contextual information (i.e. sky or ground) necessary for the selection of appropriate behavioural responses to other stimuli, such as looming objects. Our work showcases how a combination of experiments with natural stimuli and computational modelling allows discovering novel types of stimulus selectivity and identifying their potential ethological relevance.