1 Introduction

Classifications based on morphology and spiking dynamics have demonstrated that cortical neurons exhibit several distinct patterns that are robust across cortical areas (Markram et al., 2004, 2015; Spruston, 2008; Gouwens et al., 2019; Kanari et al., 2019). However, the study of waveform features and dynamics of extracellularly-recorded action potentials (EAPs) has lagged. Researchers using EAP features such as trough-to-peak duration or trough half-width have shown that binary classifications (regular/fast spiking or broad/narrow spiking) provide a practical approach to studying putative excitatory and inhibitory single-unit activity (SUA) in cortex (McCormick et al., 1985; Andermann et al., 2004; Barthó et al., 2004; Mitchell et al., 2007; Marques-Smith et al., 2018; Jia et al., 2019). Novel approaches that train machine learning models on additional features (e.g., repolarization slope) have resulted in moderate improvements beyond previous binary classifications (Jia et al., 2019; Trainito et al., 2019). Despite these apparent limitations in identifying neuron-types using extracellular recordings, multiple studies suggest that particular EAP attributes may predictably vary based on the recording location relative to the specific morphology of the neuron (Gold et al., 2006, 2007; Marques-Smith et al., 2018; Teleńczuk et al., 2018; Bakkum et al., 2019).

A convention adopted for single-channel recordings is to analyze EAPs that exhibit a canonical waveform shape (McCormick et al., 1985; Andermann et al., 2004; Barthó et al., 2004; Mitchell et al., 2007; Marques-Smith et al., 2018; Jia et al., 2019). This “canonical” waveform consists of a prominent negative trough correlated with the negative of the first derivative of the intracellular action potential (IAP) at the soma (Henze et al., 2000; Anastassiou et al., 2015; Neto et al., 2016). The negative trough is typically followed by a prominent positive peak associated with the repolarization phase of the IAP. However, this canonical waveform has a non-linear relationship with the recording location along the somato-dendritic axis of the neuron, where displacement from the soma leads to increased contributions from non-somatic morphological domains (Gold et al., 2007; Pettersen and Einevoll, 2008). The increased contribution from different morphological domains influences the ratio of the negative trough to the subsequent positive peak, as well as the time between them (Gold et al., 2007). This is determined by the change in inward and outward currents at specific locations during the different electrophysiological phases of the action potential (Gold et al., 2007). In other words, a negative trough can be accounted for by the inward sodium or calcium current phase and a positive peak accounted for by an outward potassium current phase. This phase-based analysis can be extended even further. The appearance of a prominent positive peak prior to a negative trough indicates a large contribution of a local capacitive current to the EAP waveform (Gold et al., 2007; Teleńczuk et al., 2018). Moreover, the variability in amplitude of simulated EAPs at fixed distances relative to the soma is explained by variability in total cross-sectional area of neurites where they meet the soma (Pettersen and Einevoll, 2008). It is plausible that such relationships could be exploited to detect asymmetries in neurite geometry for single neurons.

Because spatiotemporal profiles of EAPs vary with electrode location and cell type, it is not possible to infer neuron-type from the snapshot obtained in a single-channel recording. However, given that subcellular differences impact phase-based features of EAP waveforms, the use of high-density recording probes has an important advantage over single-channel probes. Multichannel recordings maintain the spatial relationships among isolated EAPs, and these spatial relationships, in principle, hold additional information about neuron-types. In an experimental study using high-density recordings, Jia et al. (2019) detect backpropagation across multiple channels and demonstrate that spatiotemporal features, such as waveform amplitude spread and the velocity of the EAP trough above and below the soma, are important features for identifying putative neuron-types. Moreover, both experiments and neuron modeling studies show that spatial patterns of EAPs reflect information about the spatial relationship between neuron and recording probe (Somogyvári et al., 2005, 2012; Szymanska et al., 2013; Delgado Ruz and Schultz, 2014; Buccino et al., 2018; Marques-Smith et al., 2018), the orientation of the somato-dendritic axis (Somogyvári et al., 2012; Delgado Ruz and Schultz, 2014; Buccino et al., 2018), gross morphology (Somogyvári et al., 2012; Delgado Ruz and Schultz, 2014; Buccino et al., 2018), and nearby morphological domains (Gold et al., 2007; Radivojevic et al., 2016; Teleńczuk et al., 2018). Thus, subcellular differences in neuronal electrophysiology and morphology are observable across EAPs centered around the soma.

These results illustrate a gap in studies of neuron classification based on morphology and fast timescale electrophysiology. We address this gap by systematically studying sources of spatiotemporal EAP waveforms (EAP sources) from simulated extracellular recordings. Here, detailed biophysically-based neuron models are simulated using computational protocols that mimic high-density, linear probe configurations. These computational models provide repeatability, the ability to look at many different neuron-types, and the possibility to link EAPs to details of the neuron morphology and membrane mechanisms. Distinct spatial and multi-resolution temporal patterns are discovered using a combination of computational modeling, automated feature engineering, and machine learning.

EAP sources are typically mixed due to the superposition of EAP waveforms from subcellular domains of neuron morphologies (Figure 1A). In this study, we apply tensor components analysis (Williams et al., 2018) to features extracted through multiresolution wavelet analysis (Mallat, 1989; Quiroga et al., 2004) to demix the relative prevalence of these EAP sources in simulated extracellular recordings (Figure 1B). The result is a low-dimensional representation for recorded units that describes how EAP sources vary across single-units, relative recording site, and phases of the action potential (Figure 1B, bottom). Here, we find that four EAP sources best characterize waveform patterns shared across a diverse population of cortical neuron-type families. We demonstrate how the demixed representations can be used as features to predict population-level distributions of morphological properties and to identify reliable morpho-electrophysiological (ME) types given biophysical constraints for recordings of SUA (Figure 1C). Finally, we show how our demixing strategy performs on the problem of unsupervised discovery of excitatory and inhibitory neuron models.

Figure 1. High-level overview of neuron-type identification strategy. (A) Cartoon of instantaneous extracellular potential of canonical waveforms at time of spike. The shape, extent, and polarity at different locations of EAPs reflect differences in morphology. (B) Schematic of demixing using spatiotemporal EAP waveforms, X, from multichannel recordings with S channel locations, T time points, and N recorded units associated each of the probe locations for each model. The inputs to the demixing step are the W coefficients for the multiresolution wavelet representation for each channel and unit, denoted by F. Demixing reveals multiple EAP sources describing how factors (space, time, and units) contribute to EAP waveforms. Each EAP source describes the contribution of a pattern across spatial locations, S, as well as the corresponding timescales and times, W. The remaining factor N describes the prevalence of the corresponding sources across all units. The prevalences of these sources for each unit (n1, n2, n3) provide input features for neuron-type identification. (C) Construction of the latent morpho-electrophysiological (ME) space relies on N where differences in the prevalence of sources are correlated with local morphometrics (like cross-sectional area) and gross morphological properties (such as location of terminal ends of basal dendrites).

Our systematic, data-driven pipeline uses an unsupervised method to extract patterns spanning channels and timescales while also representing individual units. We simulate extracellular recordings of diverse cortical neuron-type models at many spatial locations. We find multiple unique spatial templates where combinations of these templates provide representations of single-units that are predictive of variations in the underlying neuron-types, and we link these to morphology. These results can provide higher precision analyses for in vivo extracellular electrophysiology recordings given the fundamental constraints of the field's existing techniques.

2 Methods

All analyses and simulations are implemented in the Python programming language. Standard Python modules that are critical for the pipeline include: Numpy (Harris et al., 2020), Pandas (McKinney, 2020), Scipy (Virtanen et al., 2020), and Scikit-Learn (Pedregosa et al., 2011).

2.1 Simulations of neuron-type models

Neuron-type models were obtained from the NeuroML Database (NeuroML-DB.org) using the database API (Birgiolas et al., 2023). The database hosts single neuron models in multiple model formats and provides numerical details for running models. Models were downloaded in the HOC format for simulation and the NeuroML format for morphometric analyses. We chose models that were previously developed from studies of rat primary somatosensory cortical neurons as part of the Blue Brain Project (Markram et al., 2015). NeuroML-DB model IDs for these models are provided in Supplementary Table 1. Each multicompartmental model consists of multiple, detailed morphological domains: soma, axon initial segment (AIS), basal dendrites, and possibly, apical dendrites. The AIS for each model is composed of two compartments of equal length and diameter. One compartment is the stem of the axon and is connected to the soma. The other is connected to the stem compartment and extends in the same direction as the stem. The dendritic domains are highly detailed, consisting of hundreds or even thousands of compartments. In total, we obtained 105 neuron models—five morphological variants for 21 neuron-type families. These families include 11 types of excitatory pyramidal cells and 10 types of inhibitory interneurons. These model families were chosen due to their availability, detail, and development within one research laboratory.

The neuron-type models were simulated using the Python package NetPyNE (Dura-Bernal et al., 2019) with the NEURON simulator (Carnevale and Hines, 2006). The model-based bias current amplitude and rheobase were obtained from the NeuroML Database (Birgiolas et al., 2023). Steady state membrane potentials were attained by simulating a bias current injection for 1,000 ms prior to applying a 50 ms duration current at 3x the rheobase to initiate an action potential.

For each model, the extracellular space was defined as a 3D grid centered at the soma. The Y-axis was oriented with positive and negative values reflecting more superficial and deeper locations relative to the soma, respectively (Figures 2A, C). In other words, the Y-axis reflects relative position of the soma along the cortical depth. Similarly, the XZ-plane was defined perpendicular to the Y-axis and centered at the soma. Extracellular potentials were computed within the 3D space using built-in functions from NetPyNE. NetPyNE uses a transmembrane current-based forward modeling scheme (Agudelo-Toro, 2012; Lindén et al., 2014).

Figure 2. Example spatiotemporal EAPs from simulated single-units. (A) Example model morphology from the thick-tufted pyramidal cell (TTPC) family. Inset shows EAP waveforms from 31 adjacent channels (black) with center channel (red) at a radial distance of 16.3 microns from the soma. Schematic of probe placement shows the location of the center channel (red dot) and the orientation of the channels above and below center. (B) Colormap shows average EAP waveforms from the 50 closest (top, d = 10.41–16.28 microns) and farthest (bottom, d = 65.74–154.24 microns) units across the TTPC neuron-type family. EAPs were normalized before averaging. Color scale adjusted to the maximum fluctuation of the EAP. Star indicates the center channel and the intracellular spike time for all units. Center waveform depicted in white. Amplitudes for each channel used to normalize EAPs (right) show amplitude spread for the n = 50 units. (C) Similar to (A) showing superficial pyramidal cell (PC) family. (D) Similar to (B) showing closest (top, d = 10.23–13.88 microns) and farthest (bottom, d = 41.10–74.67 microns) units across PC neuron-type family.

2.2 Constraining extracellular space

In typical extracellular recording studies, single units are often excluded if they do not meet an amplitude threshold sufficient to overcome background noise (Gray et al., 1995; Henze et al., 2000; Segev et al., 2004). In this study we first wanted to identify the region of modeled extracellular space where simulated EAPs are consistent with experimental conditions. To that end, we characterized the detectability of neuron-types and identified those with EAP amplitudes <20 μV at distances of ~10 μm from the soma. This detectability threshold was conservative relative to the accepted convention of using signal-to-noise levels that are above a value that is 3–4 times the noise floor or 50–60 μV (Abeles, 1991; Henze et al., 2000; Segev et al., 2004). Thus, simulation studies of EAPs were performed in two stages: (1) EAPs were simulated to uncover the detectability limit for each neuron-type model and (2) EAPs for the studies described here were simulated within this detectability limit.

2.2.1 Detectability limit

For the first stage, using a forward modeling scheme for each of the models, we computed the extracellular potential, Ve, at 1,050 locations. The XZ-coordinates of these locations were circularly distributed at 10 angles along the XZ-plane that were chosen uniformly at random. For each angle, we used 21 radial distances from the origin along the XZ-plane, ranging between 10–120 μm at 5 μm increments. This process was repeated, resulting in radial grids at 20, 10, 0, −10, and −20 μm from the origin along the Y-axis. Each neuron-type model was subjected to the current injection protocol described above. We then computed a detectability limit for each separate model based on amplitudes computed from individual EAPs at these 1,050 locations. Here, the amplitude was taken to be the absolute difference between the trough and peak of an EAP, where the trough and peak were defined as the minimum and maximum values of an EAP, respectively.

The detectability limit was taken to be the average radial distance resulting in an amplitude of at least 20 μV. The average detectability range across all neuron-type models was ~35 μm (μ = 35.48, σ = 14.71). By fitting the average maximum amplitudes to exponential curves, <2% (2/105) of models were unable to generate EAPs with amplitudes >20 μV at locations one micron from the soma. These two models were from the Layer 2/3 small basket cell and Layer 2/3 double bouquet cell neuron-type families and had peak amplitudes of 8.4 and 1.95 μV, respectively. We excluded these two models from subsequent analyses.

2.2.2 Simulating EAPs for further study

The model-specific detectability limits for the remaining 103 models from 21 neuron-type families were used to construct the simulation dataset used for all further analyses. For each model, we first chose 100 different locations where we place the linear probe for recording simulated EAPs. For these locations, the radial distances from the origin in the XZ-plane were uniformly distributed between 10 microns and the associated detectability limit. The angle for each location was also selected from a uniform distribution. Each of these 100 locations were taken as the center of a 64-channel linear probe with channel locations spaced every 10 μm along the Y-axis, where center channels were jittered using a uniform distribution with bounds of ±5 μm along the Y-axis centered at the origin. The current injection protocols described above were used to generate EAPs at channel locations. Additionally, we recorded the somatic membrane potential for each simulation.

2.3 Waveform pre-processing

Action potentials are classically divided into intervals, or phases, of interest. These intervals are associated with the depolarization, repolarization, and recovery phases, reflecting the activation and inactivation of inward and outward currents. A non-classical phase of interest is the capacitive phase associated with spike initiation in the axon initial segment and propagation into the soma (Gold et al., 2006; Radivojevic et al., 2016; Teleńczuk et al., 2018). The effects of these phases also can be seen in EAPs (Figure 5A). Using the average EAP (white curve, bottom panel Figure 5A), we determined the average spike width (~1.4 ms), defined as the time between the trough and following peak. The duration of snippet windows for further analysis of EAPs was taken to be a multiple of the average spike width. To include all phases of interest, the extracellular action potential was extracted using a total snippet window of 5.6 ms with a 1.4 ms pre-spike interval, where the spike time is defined as the time the somatic membrane potential achieved its peak value during an action potential. The extracted spike snippets were then used to construct the mean spike-aligned waveform.

Extracellularly recorded signals around a neuron appear to exponentially decay with recording distance (Rall, 1962; Gray et al., 1995; Segev et al., 2004; Pettersen and Einevoll, 2008). Thus, recorded extracellular action potentials have signal-to-noise ratios that decrease as background activity and electrical noise become more prominent with increasing distance. This is the case for both probes situated at large distances from the soma, as well as, more distally located channels within probes situated closer to the soma. In order to capture these noise relationships, we modeled the noise component of extracellular action potentials as Gaussian distributed noise with a 1/f power law spectrum, where f is the frequency. A total of 200 independent noise signals were generated for each channel, then scaled and averaged to match noise levels from Neuropixel recordings of rat somatosensory cortex (Marques-Smith et al., 2018). This process was independently repeated across all channels resulting in uncorrelated background noise. Finally, each channel's waveform was shifted by subtracting the median of the waveform.

2.4 Selecting canonical multichannel waveforms

Typically, studies using multichannel recordings of EAPs with linear probes designate a single channel as nearest to the soma then orient additional channels symmetrically relative to this channel, i.e., centering the waveform along the probe. Waveforms are often excluded from analysis if they do not exhibit a canonical shape, i.e., if the waveform does not exhibit a prominent negative trough preceding a prominent positive peak. Our goal is to select centered, canonical waveforms for further analysis.

Usually, the channel designated as the center channel is the channel with the largest amplitude waveform (Buzsáki and Kandel, 1998; Somogyvári et al., 2005, 2012; Delgado Ruz and Schultz, 2014; Jia et al., 2019). However, high-resolution empirical findings and modeling of EAPs near the axon initial segment (AIS) indicate that the current dipole created by AP propagation from the AIS to the soma contributes more to the amplitude of the EAP than the soma alone (Teleńczuk et al., 2018; Bakkum et al., 2019). Thus, the largest amplitude waveforms occur between the distal end of the AIS and the soma.

Because we know the location of the channel closest to the soma in these modeling studies, we also assessed the reliability of selecting the channel with the largest amplitude waveform as the center. This process is complicated by the fact that there may be multiple channels with comparable large amplitudes. For example, this can occur when the center channel is sufficiently far from the soma but a different channel is close to the axon initial segment (AIS). In this case, the waveforms from the channel located close to the axon initial segment have spike troughs that occur earlier than the intracellular voltage peak used for spike detection. With this in mind, we developed a set of heuristics to center multichannel waveforms using the amplitude profile from across a linear probe (Figures 2B, D). We collected each peak in the amplitude profile that exceeded 60% of the maximum amplitude of the probe. If there was only a single peak, we took this to be the center channel. If there were high amplitudes associated with waveforms from multiple channels, we considered the two locations with the largest peak. The two waveforms were characterized as inverted or not, based on Ve at the intracellular spike time. If the potential at the spike time was greater than the median value of the waveform, then the waveform was taken to be inverted and rejected as the center channel. After the center channel was selected, we retained 31 channels for further analysis—the center channel along with 15 channels above and 15 channels below.

Spike detection often involves defining a negative threshold for recorded extracellular potentials and taking the nearby trough as the spike time. Since we have access to the intracellular spike time, the selected center channel EAP waveforms were classified as canonical or non-canonical based on potential values around the intracellular spike time (t = 0 ms) so that non-canonical waveforms could be excluded from further analysis. In a typical experimental setting, such waveforms would not be excluded from analyses if detected; however, we made this choice to ensure that the center channel waveform was not dominated by the axonal spike nor inverted due to the soma-AIS current dipole (Teleńczuk et al., 2018). The known spike time was compared to the time when the extracellular potential achieved its minimum value. We computed upper and lower tolerances based on the mean and standard deviation of the waveform (μ±σ). We defined waveforms as non-canonical if they: (1) had a prominent positive deflection greater than the upper tolerance during the pre-spike interval (−0.21 to 0 ms), (2) had a minimum value that occurred earlier than the pre-spike interval, or (3) did not have a prominent negative deflection less than the lower tolerance during the post-spike interval (0–0.42 ms). Non-canonical waveforms comprised <3% of all waveforms based on these heuristic criteria.

Finally, to ensure the same number of simulated probe locations for each neuron-type family in our further analysis, the final dataset consisted of 400 probe locations for each neuron-type family. These locations were chosen randomly from all simulated data associated with a given neuron-type family with center channel waveforms that have the canonical waveform shape.

2.5 Demixing steps

Next we use a machine learning approach to demix the underlying sources of these spatiotemporal EAP waveforms. Our demixing approach relies on a multiresolution wavelet representation of the waveforms followed by tensor decomposition. Here we describe these steps in more detail.

2.5.1 Organization of the data

Simulated EAPs were organized into a single tensor with three axes (Figure 1B). The first axis corresponds to the relative spatial locations and has dimension S = 31 for the number of channels associated with the simulated probe. The second axis corresponds to the times associated with each EAP and has dimension T = 128 for the number of time points for each recorded waveform. As described above, we generate data of dimension S×T for each probe location for each model, which we now refer to as a single-unit or simply unit. The third axis corresponds to the final set of units across the multiple neuron-type models and has dimensions N for the total number of units (N = 400 × 21). We use X to denote the data tensor with dimensions S×T×N.

2.5.2 Wavelet representation

Most feature analyses of single-channel EAP waveforms rely on readily apparent waveform features, such as trough-to-peak duration or slope during an interval of interest, but this approach limits the analysis to the coarsest of timescales. In order to discover patterns that are observable across multiple spatial locations but possibly hidden within different timescales, we use a multiresolution wavelet representation for each waveform. Such representations have been shown to result in optimal spike sorting when used as features for clustering (Quiroga et al., 2004). Automated feature extraction is applied to each of the S rows of X producing a wavelet-based feature representation of dimension W for each single-channel waveform. This feature representation is denoted by the tensor F with dimensions S×W×N.

Our multiresolution wavelet analysis is implemented using the PyWavelets Python package (Mallat, 1989; Lee et al., 2019) and requires a choice of mother wavelet and number of levels. We use the Haar wavelet because it is real-valued with both orthogonal and symmetric properties. This symmetry enables similar representation of EAPs generated by current source dipoles (e.g., axon-soma dipoles, Teleńczuk et al., 2018) that contribute to waveform inversion (Figure 3A, left). Further, this symmetry is conserved by taking the absolute value of all multiresolution coefficients (Figure 3A, right). At each level, the multiresolution wavelet analysis reduces the length of a waveform by half in an iterative process. The total number of levels was determined based on physiological intervals of interest relative to the spike time, including the prespike, depolarization, repolarization, and recovery associated phases of interest. In our final analysis, we include additional levels to capture less apparent, faster timescale subdivisions. F is the input for the tensor decomposition at the next stage of analysis.

Figure 3. Feature representation and selection of number of EAP sources. (A) (left) Normalized EAP waveforms at two example channel locations from the same simulated unit illustrate a canonical waveform (top) and an associated inverted waveform (bottom). Such waveforms indicate a common source, e.g., a soma-axon current dipole. As indicated in the text, each temporal waveform represents 5.6 ms. (right) Multiresolution analysis using Haar wavelets results in multiresolution wavelet representations (black curves) shown grouped by resolution levels (vertical gray bars). Resolution at each level between gray bars is given by dw where the exponent k is indicated by the dot color legend. Each level consists of multiple coefficients ordered in time and corresponding to the convolution of the Haar wavelet associated with that level. The absolute values of the multiresolution wavelet representations are used as the input to TCA (red curves). (B) When performing TCA, we must select the parameter value for the number of sources, R, based on similarity and performance. The result of TCA is a set of matrices with dimensions S×R, W×R, and N×R. For each source count parameter value from 1 to 19, four independent demixings are trained to the data. The model with the lowest reconstruction error is taken as the target model for computing similarity scores among the remaining three models (similarity scores gray circles). Source counts are selected for further testing if the average source similarity (red curve) is >99%. (C) Each source count with sufficient similarity is used to classify morphological classes using a random forest classifier, and the performance accuracy for each is plotted with 95% bootstrapped confidence intervals.

2.5.3 Tensor decomposition

Next we applied tensor components analysis (TCA) to the multiresolution wavelet representation, F. TCA is a higher-dimensional generalization of principal components analysis with special properties to optimally handle tensor data and demix data generated from mixed contributions originating from multiple, non-disjoint sources (Harshman, 1970; Williams et al., 2018). TCA has the additional benefit of not requiring orthogonality, and it allows for interpretations of results in terms of the original axes of the feature representation.

Our goal is to use TCA to demix EAP sources across the population of units. We implemented TCA as a non-negative, canonical polyadic (CP) tensor decomposition with the block coordinate descent (BCD) method using the Python package TensorTools (tensortools.ncp_bcd, Williams et al., 2018). This decomposition results in R non-orthogonal components (sr, nr, and wr for r∈1, …, R) that approximate the feature representation as

fi,nk,j≈∑r=1Rsirnnkrwjr,

where i∈, nk∈, and j∈ (see Figure 1B). Here, the column vectors sr and wr for a given r represent the population-level spatial distribution of an EAP source across single-units and the relative contribution of the EAP source across multiple timescales, respectively (Figure 1B). The column vector nr represents how prevalent a given EAP source is for each of the N single-units. Finding the neuron-type source average using the nnkr values for all single-units belonging to the same neuron-type family yields a neuron-type representation in terms of the EAP sources and thus represents the probability that the EAP sources would be observed for that neuron-type. The reconstruction error and decomposition similarity scores were also computed using TensorTools and tested for multiple values of R (Supplementary Figure 1).

In what follows, the column vectors sr and wr are termed spatial sources and multiresolution temporal sources, respectively. When referring to the spatial and multiresolution temporal sources corresponding to the same r, we use the term EAP source. The numerical values of these two components are referred to as their contributions, e.g., the spatial source s1 has a contribution of si1 at the ith channel. Thus, in our results, we show contributions for both spatial courses (Figure 4, left) and multiresolution temporal sources (Figure 5B). In contrast, the column vectors nr are termed the single-unit source prevalences and their numerical values are referred to as the prevalence for the corresponding unit.

Figure 4. Summary of final demixed EAP sources. (A–D) The normalized contribution of each EAP source at the 64 channels along the probe. The center channel is taken to be closest to the putative perisomatic region of all neuron models. The multiresolution temporal sources (middle) describe the contribution of each EAP source across multiple timescales (black curves, gray lines show grouped timescales), as in Figure 3A. In other words, the values of each spatial and multiresolution temporal source reflect the contribution of that EAP source in space and time to the EAP waveforms. Source prevalences (right) were averaged within neuron-type families. Error bar depicts 1 standard deviation. Shown prevalences were normalized by their maximum values.

Figure 5. Multiresolution patterns reveal underlying electrophysiological phases. (A) Normalized canonical IAP (top) and EAP (bottom) waveforms show variability within electrophysiological phases of interest (I–IV). The average (black) of the action potentials represent the canonical IAP shape. The corresponding EAPs (gray) and canonical EAP (white) at the center channel across the probes (bottom) are labeled based on the phases of the canonical IAP. (B) The contributions of the multiresolution temporal sources (right) were reconstructed, yielding their temporal dynamics, to illustrate the contribution relative to the spike time (gray dashed line). Intervals with larger contributions of multiresolution temporal sources are depicted in gray. Insets show the corresponding spatial sources.

2.6 Morphometric analysis

All morphometrics were computed from the NeuroML format description for each model (Gleeson et al., 2010), where each reproduces the morphological description from the original Blue Brain Project model (Markram et al., 2015). This XML-based language enables flexible, automated extraction of compartment details using the ElementTree XML API in Python, as well as computation of model morphological features. NeuroML specifies the structure of morphologies within a nested hierarchy of parent-child relationships between compartments (Crook et al., 2007). Each compartment (segment element) is specified by the xyz-coordinates and radius for the proximal and distal ends of the compartment. All neuron-type model segments were specified as belong to somatic, axonal, basal dendrites, or apical dendrites segment groups, using the segment labels provided in the original Blue Brain Project models.

We computed morphological features of the most proximal compartments of the neuron model for use in further analysis. This approach is based on modeling results demonstrating that overall differences in the cross-sectional area and distance from the soma of local neurites are reflected by the amplitude variability of EAPs (Pettersen and Einevoll, 2008; Teleńczuk et al., 2018). For the somatic group, we computed the overall soma height in the Y-dimension (i.e., depth axis of cortex), the number of segments with the soma as a parent (known as stems), and the total cross-sectional area (CSA) of all stems. The CSA was computed as π(d/2)2 where d was the diameter of the proximal end of the stem (Figure 6A). Many neuron-types have multiple basal dendritic stems that project in several directions. We divided the basal stems into upper and lower divisions by bisecting the somatic segments into upper and lower groups based on the soma mid-point in the Y-dimension (Figure 6A, soma Y-mid). For the basal dendrite group, we computed the total number of basal stems, the total CSA of all stem segments as well as the CSA for the upper and lower stem divisions. We also computed the average Y-locations of the basal terminal ends (Figure 6A, yterm) for the upper vs. lower stem groups. Those neuron-type models that have apical domains have at most two, projecting upward, downward or in both directions. We also computed the CSA of the upper and/or lower apical stems. For the axonal group, we computed the stem CSA of the most proximal segment in addition to noting whether it belong to the upper or lower division of soma. Finally, we determine the proximal and distal Y-locations of each neurite stem (Figure 6A, defined as compartments starting at the soma and extending up to the first branch, yprox and ydist) as well as the distal Y-locations of neurite terminals (Figure 6A, yterm).

Figure 6. Relationships between model morphologies and demixed EAP sources. (A) Schematic of multi-compartment models showing different morphological domains. (B) Bar plots depict spatial profiles for different stem geometries relative to discovered spatial sources (as depicted in Figure 4, left column). The spatial sources are scaled to illustrate shared patterns. (C) Bootstrapped canonical correlation analysis shows the relationship between morphometrics and SU source prevalences. (CSA, cross sectional area; U, upper; L, lower). (D) Scatter plots illustrate the differences in SU source prevalences across recording distances, radial distance r along XZ-plane, for each neuron-type family and relative to the nearest recording distance (~10 microns for each). Linear regression (red) and binned averages (black, bin size = 10 microns) indicate whether the source tended to increase or decrease with distance from the putative perisomatic region.

2.7 Bootstrapped canonical-correlation analysis

Canonical correlation analysis (CCA) is a cross decomposition method for finding the linear combination of two random variables that maximizes correlations between the variables. We applied CCA to determine how predictive the single-unit source prevalences were of local morphometrics. To implement CCA, we use the SciKit-Learn cross_decomposition library in Python. We fit R independent CCA models to the matrix of morphometrics for each neuron-type model, M, and the matrix of randomly selected single-unit source prevalences for each neuron-type model and corresponding to source r, Pr for r∈1, …, R. For all CCA models, we took n_components to be 1. We adopted a bootstrapped approach to CCA where the matrices M and Pr were independent constructed on each bootstrap iteration. We describe the matrices, their construction and bootstrap algorithm below.

Let M be the X×M matrix of morphometrics for each neuron-type model where the total number of neuron-type families is given by X and the total number of morphometrics to be predicted is given by M. In other words, mij of matrix M is the jth morphometric of the ith neuron-type model. Let Pr be the X×P matrix of single-unit source prevalences where P is a fixed number of source prevalences that are randomly selected from the total number of source prevalences corresponding to the same neuron-type model and r corresponds to the rth source, as above. In this case, pik of Pr corresponds to the kth source prevalence of source r and the ith neuron-type model.

At every bootstrap iteration, one neuron-type model is randomly selected from each of the X neuron-type families. The matrix of morphometrics M is then populated by the M morphometrics of the X selected neuron-type models. For the source prevalence matrix, we populated the matrix Pr with 50 randomly sampled (with replacement) source prevalences for each of the 21 selected neuron-type models and demixed source r. Matrices M and Pr were then standardized using the Python StandardScaler before fitting the CCA models. We computed a total of 1,000 canonical-correlations using the above sampling approach and repeated the bootstrapped CCA for each demixed source.

2.8 Random forest classification

We estimate the limits of using single-unit source prevalences for each unit by training random forest classifiers. Random forest classification uses several decision trees as estimators and averages the classification accuracy by averaging over the estimators. We implemented the random forest classifier from the SciKit-Learn ensemble library. Random forest classifiers were trained using the single-unit source prevalences as features. The training labels were one of three groupings: (1) neuron-type family, (2) discovered morpho-electrophysiological groups, or (3) excitatory-inhibitory types. For all random forest classifiers, we use five-fold cross-validation. We used the built-in grid search method hyperparameters max_depth and n_estimators which control the depth of the decision tree and the number of decision tress to average over, respectively (Supplementary Figure 2). From the grid search, we selected the parameter combination resulting in the greatest classifier performance (using the out-of-bag score). The final classification accuracy for any given label grouping was scored based on withholding 20% of the units for the test set prior to performing the grid search. All confusion matrices were based on the test set of the final classifiers. The importance of each feature was obtained from the final random forest classifier models using built-in methods. We computed the impurity-based importance (Gini importance) using the default feature_importances_ method for RandomForestClassifier in Scikit-learn.

In order to assess the relative performance of our demixed source prevalences in classification tasks, we compare to a more traditional approach using multiple, single-channel features extracted from the center waveform and two multichannel features. The single-channel features include peak-to-trough duration, depolarization half-width, peak-to-trough ratio, repolarization slope, and recovery slope. The trough is defined as the minimum value of the waveform, and the waveform peak is the maximum value following the trough. Depolarization half-width is defined as the time difference between the waveform values corresponding to 50% of the waveform trough voltage. The repolarization and recovery slopes (μV/ms) were calculated using the value of the waveform at the trough and peak, respectively, along with the values 30 μs after. We included two multichannel features: multichannel spread and total propagation velocity. The spread is defined as the range in distance in microns across all channels with amplitudes >12% of the maximum amplitude. For the total propagation velocity, first we extract the time of the trough for all channels. We then compute the absolute value of the median slope (μm/ms) between the center channel and all channels above (propagation velocity above) and also below the center channel (propagation velocity below). Slope is defined as the distance between channels divided by the time difference of the troughs. The total propagation velocity is then taken to be the sum of the two slopes from above and below. We applied a random forest classifier as defined previously with the added constraint that max_features is taken to be four, the number of sources we are comparing against. This has the effect that only four randomly chosen (without replacement) features are used to construct a given random tree within the ensemble.

2.9 Clustering methods

To benchmark the performance of our methods, we performed standard unsupervised learning to group single-unit source prevalences into putative excitatory and inhibitory types to compare to their true class labels. We applied k-means and agglomerative clustering from the SciKit-Learn cluster library and Gaussian mixture models from the mixture library. In the case of agglomerative clustering, we applied both Ward's method and the average linkage method for merging the hierarchical clusters. In all cases, Euclidean distance was used as the distance metric. We also took the number of groups (either n_clusters or n_components) to be two. To assess how well each method performed on this benchmark, we computed the percent correct clustered using the the true labels.

3 Results 3.1 Unsupervised discovery of EAP spatial and multi-timescale dynamic patterns of EAP sources

Here we describe the overall demixing strategy for discovering reliable morpho-electrophysiological (ME) types across multichannel recordings of EAPs. The approach is validated using biophysical simulations from morphologically-realistic cortical neuron models. Overall, the strategy finds a reduced set of population-level spatial patterns (spatial sources) and multiresolution dynamics (temporal sources) that describe variability across simulated or recorded EAP waveforms. These EAP sources expose the contribution of specific morphological domains and electrophysiological phases to EAPs. Each unit can be represented by different relative contributions (prevalences) of these spatial and temporal sources. A major strength of this approach is its ability to uniquely represent each unit within a given dataset as a combination of EAP sources which are predictive of ground-truth morphological properties. We can also find the prevalence of these sources at different spatial locations, times, and timescales across neuron-types.

Our strategy consists of four stages. First, automated feature extraction is used to represent multiresolution dynamics of EAPs recorded on individual channels. Second, a tensor decomposition is obtained that represents each of the isolated units in terms of shared multiresolution dynamic patterns across channel locations (i.e., demixing of EAP spatial and temporal sources). Third, a supervised method demonstrates the limits of this approach for multiple parameterizations of the tensor decomposition. Finally, a latent morpho-electrophysiological space is constructed to represent neuron-types using single-unit (SU) source prevalences. We show that neuron-type representations in the morpho-electrophysiological space are correlated to local model morphometrics. We also use these representations for unsupervised clustering of neuron-types.

Our strategy emphasizes the analysis of demixed EAP sources that summarize the trial-by-trial variability across the single-unit EAP spatiotemporal waveforms (Figure 1B) without trial averaging. These demixed EAP sources describe how each single-unit differs based on shared temporal dynamics across multiple timescales and where those patterns are more localized relative to the soma. These EAP sources can be examined to generate hypotheses about how populations of latent neuron-types can be discovered from multichannel recordings. Interpretations of EAP sources are supported by correlating the single-unit source prevalences with local morphometrics to predict the morphological properties discovered through demixing EAP sources.

3.2 Four source model approximates multichannel EAP waveforms and predicts neuron-type families

When applying tensor components analysis (TCA), we must choose an appropriate number of components to describe the data (see Section 2.5.3). This approach discovers a low-dimensional representation of the dataset. Additionally, random initiations can lead to quantitatively different solutions. We found that there were qualitatively similar patterns across random initiations despite differences in the exact solutions. In order to constrain the demixed source model, we studied the applicability of TCA by comparing independent trainings of the model, determining a similarity score for each number of sources (or components).

We first varied the number of sources and tested four independent trainings. We found that the average similarity scores were larger than 99% for up to four EAP sources before rapidly decaying (Figure 3B). While the standard experimental setting would not have access to ground-truth information on morphological types, the use of computational models from defined neuron-type families facilitated further constraining the number of desirable EAP sources. To that end, we determined the predictive power of the demixed source models for up to four sources using a random forest classifier. Each classifier was trained to predict the neuron-type family of a simulated unit based on the single-unit source prevalences. We found that the classification accuracy steadily increased with the addition of each new source (Figure 3C).

These results identify four unique sources that most contribute to spatiotemporal patterns shared across multiple neuron-types recorded on linear probes (Figure 4). Note that we ordered the EAP sources by decreasing importance of each source to classifying the neuron-type families using the random forest classifier (Figure 7D, top). By visualizing the spatial and multiresolution patterns associated with each source, we found that patterns discovered using fewer sources were qualitatively reproduced in the final four source model. This suggests that the discovered EAP sources may reflect fundamental differences across neuron-types that could be reproduced across many conditions.

Figure 7. Source prevalences for units predict neuron-type model families. (A) Each point in the scatter plot corresponds to a unit (a model and probe location), where the color indicates the neuron-type family (color legend on the right). The location of the point corresponds to the source prevalence for the sources shown on the axes. This 4-dimensional space illustrates the latent morpho-electrophysiological (ME) space and shows the contribution of each source to recorded multichannel EAPs. The expected demixed representation of each neuron-type family (neuron-type representation, colored squares) are depicted by the bootstrapped centroids using 1,000 random subsamples. (B) Average prediction accuracy for random forest classifier (out-of-bag score) is shown for each neuron-type family. Dashed line depicts 50% threshold. (C) Confusion matrix from random forest classifier across all neuron-types. Neuron-type labels show excitatory (red) and inhibitory (blue) families. Ground truth is comprised of model neuron-types, and predictions are predicted model neuron-types. Color indicates probability of predicted neuron-type family normalized by row. (D) Feature importances estimated by random forest classifier using 100 random subsamples of simulated units (top). Dashed line depicts the threshold for equally weighted importances. 95% confidence intervals are too small to show. Coefficient of variation for change in prevalence with distance (bottom) across all units (gray bars) and average values for each neuron-type family (black lines, error bars are one standard deviation).

3.3 Multiresolution temporal sources reflect different action potential phases

We proposed that the discovered four source model reflects fundamental differences across neuron-types. In order to evaluate that proposition, we investigated how the multiresolution temporal and spatial patterns of the EAP sources could be interpreted in terms of the electrophysiology and morphology of the neuron-type models, respectively. We first reconstructed the four EAP sources in terms of the multiresolution temporal patterns and interpreted the approximated dynamics of each source.

We divided the dynamics of the canonical waveform of the somatic action potential into electrophysiological phases (Figure 5A, top). We took the canonical waveform to be the average of all somatic IAPs (normalized and aligned by their peaks). The elecrophysiological phases were defined as: capacitive (I), depolarization (II), repolarization (III), and recovery (IV) phases. We then mapped these phases onto the canonical waveform of the EAPs closest to the perisomatic region (Figure 5D, bottom). We next reconstructed the dynamics by taking the product of the contributions of the discovered multiresolution patterns and Haar wavelets at the corresponding times and timescales. Since the EAP sources are all expressed as non-negative contributions, the resulting source reconstructions primarily reflect the magnitude of the contribution in time for a reconstructed temporal source. We found that the reconstructed source dynamics exposed electrophysiological phases with the most variability across neuron-type models, and that the dynamics predominately differed in the repolarization and recovery phases across all models (Figure 5B, gray regions).

More specifically, Source #1 is associated with changes in the recovery phase (centered at 2.1 ms), as well as fast changes during the capacitive phase (Figure 5B). Recall that this source was the most important source for classifying neuron-type families. The variability for all IAPs peaked during the initial part of the recovery phase at ~1.75 ms. Due to the delay between the time of peak variability of IAPs and the time of peak contribution after the spike time of Source #1, we conclude that this source is likely associated with propagation into dendrites. Further, the high contribution during the capacitive phase also points to plausible propagation of strong axial currents into large dendritic compartments (Gold et al., 2007). Together these results suggest that Source #1 may reflect active dendritic processes activated by backpropagation, similar to results found by Jia et al. (2019). This is the subject of ongoing further research.

The next most important sources were Sources #2 and #3 in that order. Source #2 is associated with increasingly rapid fluctuations starting at the beginning of the recovery phase. This source monotonically increases in the magnitude of its contribution with time. As it is least present during the active phases of the action potential and it contributes at channel locations most distal from the soma, we concluded that Source #2 reflects the constant background noise level present in all extracellular recordings. In contrast to Source #1, Source #3 is associated with the middle of the recovery phase (centered at 2.8 ms) and its contribution is observed on a slower timescale. Additionally, Source #3 peaks at greater distances (± 60 μm) from the soma while also contributing the least nearby the soma. This, together with the smaller contribution from the capacitive phase, points to propagation of the action potential into more distal regions of dendrites. As these distal dendrites have weaker axial currents due to their smaller diameters, their contribution can likely only be observed when the more dominant contributions of larger dendritic compartments and the soma have dissipated.

Finally, Source #4 is associated with fast timescale differences during the repolarization phase that partially carry over into the recovery phase. This source mostly contributes at channels nearest to the soma. Unlike other EAP sources, this source retains nearly 50% of the peak spatial contribution across all channels (Figure 5B inset, Figure 4D). Because of this wide reaching contribution, the peak at the center channel, and its association with the repolarization phase, it's likely that Source #4 is largely due to differences in the somatic and AIS currents which dominate EAPs (Bakkum et al., 2019). The two additional peaks aren't interpretable at this stage of analysis. Interestingly, the variability in this phase has historically been exploited for classifying units from single-channel recordings as putative excitatory and inhibitory types by computing spike half-widths or spike durations. This variability is apparent by visual examination of the EAP waveforms at the center channel during the repolarization phase (Figure 5A, bottom). However, this was the least important source for classifying neuron-type families as previously indicated.

3.4 Spatial sources reflect the population density of different morphological domains

We now visit the problem of interpreting the spatial sources in terms of the morphological properties of the neuron-type models. Cortical neurons can exhibit both symmetry in morphological domains and certain deep-superficial layer biases. For example, deep layer cortical pyramidal cells often have large apical dendrites that project superficially. However, some cortical pyramidal cells exhibit apical dendrites that project toward deeper layers or even bidirectionally. Difference