In computational neuroscience, multicompartment models are among the most biophysically realistic representations of single neurons. Constructing such models usually involves the use of the patch-clamp technique to record somatic voltage signals under different experimental conditions. The experimental data are then used to fit the many parameters of the model. While patching of the soma is currently the gold-standard approach to build multicompartment models, several studies have also evidenced a richness of dynamics in dendritic and axonal sections. Recording from the soma alone makes it hard to observe and correctly parameterize the activity of nonsomatic compartments. In order to provide a richer set of data as input to multicompartment models, we here investigate the combination of somatic patch-clamp recordings with recordings of high-density microelectrode arrays (HD-MEAs). HD-MEAs enable the observation of extracellular potentials and neural activity of neuronal compartments at subcellular resolution.

In this work, we introduce a novel framework to combine patch-clamp and HD-MEA data to construct multicompartment models. We first validate our method on a ground-truth model with known parameters and show that the use of features extracted from extracellular signals, in addition to intracellular ones, yields models enabling better fits than using intracellular features alone. We also demonstrate our procedure using experimental data by constructing cell models from in vitro cell cultures.

The proposed multimodal fitting procedure has the potential to augment the modeling efforts of the computational neuroscience community and provide the field with neuronal models that are more realistic and can be better validated.

In computational neuroscience, multicompartment models provide one of the most biophysically detailed representations of single neurons. They are built by combining a precise morphological reconstruction of neurons, obtained through imaging techniques, with electrophysiological characteristics of ion-channel dynamics and their distribution over the neuron morphology.

The use of multicompartment models has enabled researchers to explore several characteristics of neuronal dynamics, including active dendritic properties (Hay et al., 2011; Almog & Korngreen, 2014; Gidon et al., 2020) and the role of the axonal initial segment in initiating action potentials (Hallermann et al., 2012; Goethals & Brette, 2020; Fékété et al., 2021). The models also were used to generate experimentally testable hypotheses to drive research forward. In recent years, we have witnessed massive international efforts in constructing and sharing biophysically detailed multicompartment models. The neocortical microcircuit portal (Ramaswamy et al., 2015; Markram et al., 2015) of the Blue Brain Project contains thousands of publicly available cell models of the rat somatosensory cortex. A similar effort is being conducted by the Allen Institute of Brain Science, whose cell-type database includes a growing number of cell models of mice and even humans (Gouwens et al., 2018).

Historically, multicompartment models were constructed by taking ionic mechanisms from preexisting knowledge bases and by manually and iteratively tuning the many parameters of the models, until the model matched the expected or observed behavior (Almog & Korngreen, 2016). In more recent years, computer-based optimization has enabled computational neuroscientists to explore the large parameter space more thoroughly and faster—for example, by using evolutionary strategies to search for sets of parameter values that provide good fits to experimental data (Druckmann et al., 2007).

For the vast majority of available multicompartment models, experimental features used to fit the model are extracted solely from somatic patch-clamp recordings. While the soma clearly is a very important “compartment” of a nerve cell, neurons are much more than just their somata. For example, complex dynamics can arise from active dendritic properties (e.g., calcium spikes, nonlinear integration) (Larkum et al., 1999; Koch & Segev, 2000; Ujfalussy et al., 2015; Hay et al., 2011) and from the axon initial segment, which is known to be the location where spikes initiate and which has a determining role for neuron excitability (Huang & Rasband, 2018; Debanne et al., 2011). However, performing simultaneous patch-clamp recordings of multiple compartments of the same cell is extremely challenging (Larkum et al., 2001; Hay et al., 2011; Almog & Korngreen, 2014).

A strategy to capture neuronal firing dynamics over a larger spatial range includes the use of extracellular recordings. Extracellular signals are generated by transmembrane currents of all neuronal compartments (Nunez & Srinivasan, 2006; Lindén et al., 2014; Hagen et al., 2018) and provide an indirect readout of the intracellular signaling. Gold et al. (2007) first suggested that extracellular electrical recordings could provide suitable data for devising extracellular models, but their suggestion has not been followed up by the modeling community. With the advent of high-density microelectrode arrays (HD-MEAs), however, extracellular signals have been demonstrated to enable the recording, especially in vitro (and to a certain extent also in vivo; Jun et al., 2017; Jia et al., 2019), of signals of individual neurons at subcellular resolution (Müller et al., 2015). The richness of information obtained from extracellular signals enables one to observe electrical potential and signal distribution over the entire neuron. Therefore, extracellular recordings allow for parameterizing neuronal compartments, which cannot be simultaneously and directly probed by patch-clamp experiments.

Here, we present a proof-of-concept study combining patch-clamp and HD-MEA data to construct single-neuron multicompartment models. We first give a general overview of the workflow and describe the steps to include extracellular signals as an extra source of data into the models. Next, we validate our approach on a ground-truth model with known parameters to assess the benefit of including extracellular data on the final outcome. Finally, we present examples of cell models, constructed from rat cortical cell cultures, and discuss the differences and improvements of our multimodal strategy over the use of patch clamp alone.

Before diving into our multimodal fitting approach, we introduce here the general workflow for constructing multicompartment models in order to lay out the different steps and terminology that computational neuroscientists use to build biophysically detailed models of single neurons. The main objective of the construction of a multicompartment model is to build an in silico replica of a real neuron that can be used to describe the electrical activity of a neuron. The pipeline to construct multicompartment models consists of several steps (see Figure 1) that we next outline.

Figure 1:

General workflow to construct multicompartment models. In the data acquisition phase, a patched neuron is stimulated by applying different protocols, and the respective responses are recorded. Then the patched neuron is imaged, and the neuron morphology is reconstructed. Mechanisms that determine or describe the electrophysiological characteristics of the neuron need to be defined. The corresponding parameters are then used for the fitting. Features of the measured neuronal characteristics and responses are extracted to reduce the dimensionality of the recorded data. In the optimization phase, the neuron morphology and the defined mechanisms are used to emulate the applied protocols, and the simulated electrophysiological characteristics are evaluated against the experimentally obtained ones. This procedure is repeated while iteratively varying parameters until convergence is reached. Finally, the morphology, mechanisms, and optimized parameters yield the fitted model.

Figure 1:

General workflow to construct multicompartment models. In the data acquisition phase, a patched neuron is stimulated by applying different protocols, and the respective responses are recorded. Then the patched neuron is imaged, and the neuron morphology is reconstructed. Mechanisms that determine or describe the electrophysiological characteristics of the neuron need to be defined. The corresponding parameters are then used for the fitting. Features of the measured neuronal characteristics and responses are extracted to reduce the dimensionality of the recorded data. In the optimization phase, the neuron morphology and the defined mechanisms are used to emulate the applied protocols, and the simulated electrophysiological characteristics are evaluated against the experimentally obtained ones. This procedure is repeated while iteratively varying parameters until convergence is reached. Finally, the morphology, mechanisms, and optimized parameters yield the fitted model.

Close modal

2.1  Data Acquisition

The first and essential step is data collection. A target neuron is patched with a micro-pipette in current-clamp mode to record its responses to different stimulation protocols. The different protocols are designed to unveil different characteristics of neural dynamics, such as firing patterns, sub-threshold behavior, or spike adaptation. Each protocol is usually applied in several runs. The patch pipette is loaded with a dye (either biocytin or a fluorescent dye) that fills the cytosol of the target neuron during the recording session and that is used for reconstructing the neuron morphology.

2.2  Morphology Reconstruction

After the electrophysiology data acquisition, the target neuron is imaged to reconstruct its 3D morphology. The imaging consists of a z-stack of high-resolution images, which allows for a precise reconstruction of the neurites and their diameters using specialized software tools—for example, Neurolucida, Simple Neurite Tracer (Arshadi et al., 2021), and vaa3D (Peng et al., 2014).

2.3  Definition of Model Mechanisms

The neuron morphology only provides geometrical information for the model. In order to reproduce neuronal behavior, the model needs to be populated with mechanisms, which govern the equations that give rise to the neuronal dynamics. Mechanisms can be passive (leaky) or active (voltage-gated ion channels); the latter are usually modeled with a Hodgkin-Huxley formalism (Hay et al., 2011; Markram et al., 2015). Additional mechanisms can further control the dynamics of specific ion concentrations—for example, of Ca2 +, which can affect the dynamics of subpopulations of ion channels. Different mechanisms are defined in different neuronal sections (e.g., soma, dendrites, axons). Computational neuroscientists need to select the mechanisms for the different neuronal compartments, which will greatly influence the performance of the final model. The definition of the mechanisms also determines the parameters that need to be identified during the fitting procedure. Parameters include, for example, the maximum conductances of an ion channel in a certain compartment; properties of passive components, such as membrane or axial resistances; or free variables that are included in the mechanistic equations, such as the decay of calcium concentrations over time.

2.4  Feature Extraction

As the goal of the model is to optimally reproduce experimental data, one needs to define a cost function to quantify the fit of the obtained solution. While a straightforward approach could be to compute the point-to-point distance between the experimentally obtained values and the in silico counterparts, such a procedure has proven to be suboptimal in several previous studies (Druckmann et al., 2007; Weaver & Wearne, 2006) due to the intrinsic variability of neural responses even to the exact same input. A better and more viable solution includes processing the raw neuronal responses and extracting features that describe the response behavior in a compact and informative way. For example, if a protocol is designed to characterize the firing properties of a neuron, one could extract features, such as mean firing rate and spike counts, as well as bursting patterns and adaptation indices. Additional features can be used to describe the action potential (AP) waveforms (e.g., AP amplitude, AP duration) or sub-threshold behavior. For each feature, one can extract the average value and the standard deviation over different runs, namely, μexp  and σexp . After extracting all the relevant features from the raw neuronal responses, the model is ready to be optimized.

2.5  Optimization

The optimization step is aimed at finding good sets of parameters to fit the extracted features. By combining the protocols, morphology reconstruction, and mechanisms, a simulator runs the experiment for a given set of parameters. Using the simulated responses, features, analogous to those that have been defined in the feature extraction step, are computed and compared to those extracted from the experimental data. For each feature i, a score si is computed,
si=|μexpi-fi|σexpi,
(2.1)
where fi is the feature value computed from the simulated response, μexpi is the mean feature value from the experimental data, and σexpi is the experimental standard deviation of the feature. The scores of all features are then summed to compute an overall fitness for a certain solution. Several algorithms, mainly based on evolutionary strategies (Damart et al., 2020), are then used to iteratively explore the parameter space to find solutions with good fitness—solutions that fit the features extracted from the experimental data well. The optimization procedure yields, upon convergence, the solution (set of parameters) that yields the best fitness. The combination of morphology, mechanisms and the best parameters then constitutes the final multicompartment model.

3.1  Combining Patch-Clamp and HD-MEA Data for Fitting Multicompartment Models

The general framework described above was developed mainly for use of data from whole-cell patch clamp alone. In order to combine patch-clamp and HD-MEA readouts, we need to augment several steps of the model-fitting pipeline.

First, during the data acquisition, intracellular and extracellular recordings of the target neuron need to be acquired (see Figure 2A, left). Therefore, we used cultured neurons from embryonic rats plated on top of an HD-MEA with 26,400 electrodes (Müller et al., 2015; see section 5.6.2 for details). The patch-clamp and the HD-MEA systems were synchronized so that we could simultaneously acquire the patch-clamp readout and the extracellular signals from more 500 electrodes in proximity of the patch pipette. Since the HD-MEA signals featured much lower signal-to-noise ratios with respect to the patch-clamp readout, we used the patch-clamp signals to precisely detect action potential times, which we then used to average the extracellular signals (patch-triggered average; PTA). The resulting signal is referred to as an electrical “footprint” (distribution of measured electrical potentials across the array electrodes) or template (see Figure 2A, right) and is a clean signal that contains high spatiotemporal resolution information, reflecting the underlying morpho-electric properties of the patched neuron. This footprint or template is used as input data for feature extraction from the extracellular signals. Note that the template is computed using all action potentials acquired during the execution of all protocols, since the more the spikes are acquired, the cleaner the template becomes (extracellular spikes are assumed to be relatively constant across different protocols). In addition, the HD-MEA system we used features a switch-matrix architecture so that not all electrodes but only 1024 of them are recorded simultaneously (Frey et al., 2010), which may result in some holes in the template. Finally, electrodes recording low-amplitude signals (peak-to-peak amplitude below 5 μV) were excluded (see Figure 2B).

Figure 2:

Combining HD-MEA and patch-clamp data. (A) Simultaneously recorded patch-clamp (blue trace) and HD-MEA data (gray traces) were used to extract the extracellular template using patch-triggered averaging (PTA). (B) From the latter (after excluding channels with a peak-to-peak amplitude below 5 μV), several extracellular features were computed, either channel by channel (e.g., half-width in red, peak-to-trough ratio in green) or in relation to the channel with the largest signal amplitude (e.g., negative peak time difference, in purple). (C) Examples of features computed across all the channels: peak-to-trough ratio (left) and negative peak time difference (right).

Figure 2:

Combining HD-MEA and patch-clamp data. (A) Simultaneously recorded patch-clamp (blue trace) and HD-MEA data (gray traces) were used to extract the extracellular template using patch-triggered averaging (PTA). (B) From the latter (after excluding channels with a peak-to-peak amplitude below 5 μV), several extracellular features were computed, either channel by channel (e.g., half-width in red, peak-to-trough ratio in green) or in relation to the channel with the largest signal amplitude (e.g., negative peak time difference, in purple). (C) Examples of features computed across all the channels: peak-to-trough ratio (left) and negative peak time difference (right).

Close modal

The second required modification of the overall process was at the feature extraction level. While intracellular features are scalar values, extracted from the somatic patch-clamp trace (e.g., mean frequency or mean AP amplitude), the extracellular template E includes C electrodes and T time points and is therefore multidimensional in nature (ERC×T). In order to extract relevant and lower-dimensional features from such templates, we compiled a set of 11 extracellular features fext that were defined and used for each recording electrode (or channel) (fextRC). The features either describe channel-specific waveform parameters of the template (termed absolute features, such as peak-to-valley duration, peak half-width duration, and peak-to-trough ratio) or values relating the measured value in each channel to that of the channel exhibiting the largest signal amplitude (relative features, for example, negative/positive peak time difference, negative/positive peak amplitude). In total, we defined 11 extracellular features (Nfext=11) (see section 5.1.2 for details). In Figure 2B, we display the thresholded template (including channels with peak-to-peak amplitude above 5 μV) of a recorded cell, with the largest-signal-amplitude channel depicted in blue and a second channel colored in orange. The right part shows a close-up version of the template for these two channels and some examples of how channel-wise features (peak-to-trough ratio, half-width) were calculated. The negative peak difference was computed by comparing the negative-peak amplitude of the second channel (orange) with the peak of the largest-amplitude channel (blue). Figure 2C shows the feature maps across the recording electrodes for two exemplary extracellular features, peak-to-trough ratio and negative peak time difference.

Still, extracellular features were defined for each channel (C can be as high as several hundreds of electrodes), so we explored three different strategies to further reduce the dimensionality of extracellular features that were then included in the optimization:

  • Single: With this strategy, a subset of extracellular signals was manually selected (Nselected), and extracellular features were extracted separately for all selected channels. Since the template was obtained by averaging all available extracellular spikes from different runs, the standard deviation of the feature σexp , which was used to calculate the feature score (see equation 2.1), could not be computed from experimental data. In this case, σexp  was set to 5% of the feature value μexp . The number of additional extracellular features used for optimization was Nselected×Nfext. We used seven single electrodes per cell.

  • All: In order to use the entire information from all extracellular electrodes, we designed a second strategy that used the cosine distance between simulated and experimental features as the score si,
    si=1-μexpi×fiμexpifi,
    (3.1)
    where i refers to the ith feature. This strategy added Nfext objectives to the optimization (one scalar for each extracellular feature).
  • Sections: The third strategy lies in between the single and all strategies. In this case, we manually selected Nsections extracellular areas of several electrodes (10 to 20) that corresponded to different regions of the neuron (e.g., dendritic, perisomatic, axon initial segment areas). For each of these sections, the cosine distance between the simulated and experimental features was used as a score. This strategy provided Nsections×Nfext features to the optimization. We used three to four sections per cell.

The different strategies to include extracellular signals allowed us to explore how a modification of the balance between intracellular and extracellular features that were used for the optimization procedure affected the overall fitting performance. In order to provide informative features for optimization, when using the single and sections strategies, we selected electrodes/sections that covered different parts of the neuron (dendritic, perisomatic, axon initial segment). Figure 9 in the appendix shows the channel selection for the single and sections strategies for the different experimental cell models used in this article.

Figure 3:

Ground-truth model. (A) Visualization of the neuron on top of the 80-channel rectangular MEA used for the simulations. The extracellular electrical signal template is overlaid in black. The red trace depicts the channel with the largest amplitude, which is the channel closest to the AIS section. (B) Left: Morphology of the ground-truth neuron, highlighting the additional AIS and myelin sections. Right: Sample action potential signals at the soma (blue), middle AIS (orange), and distal AIS (green). The AP is initiated at the distal end of the AIS and travels back toward the soma. (C) eCode protocols used in the virtual experiment. For each protocol, the top traces show the somatic membrane potentials and the bottom ones the applied current stimuli.

Figure 3:

Ground-truth model. (A) Visualization of the neuron on top of the 80-channel rectangular MEA used for the simulations. The extracellular electrical signal template is overlaid in black. The red trace depicts the channel with the largest amplitude, which is the channel closest to the AIS section. (B) Left: Morphology of the ground-truth neuron, highlighting the additional AIS and myelin sections. Right: Sample action potential signals at the soma (blue), middle AIS (orange), and distal AIS (green). The AP is initiated at the distal end of the AIS and travels back toward the soma. (C) eCode protocols used in the virtual experiment. For each protocol, the top traces show the somatic membrane potentials and the bottom ones the applied current stimuli.

Close modal
Figure 4:

Optimization results on ground-truth model. (A) Distributions of intracellular (left) and extracellular (right) scores (i.e., sum of all feature scores) for 10 seeds for each applied strategy (soma: blue; all: orange; sections: green; single: red). (B) Distributions of intracellular (left, N = 86) and extracellular (right, N = 796) feature scores for the solution with the best intracellular total score. The *** represent a p < 10−5 with respect to soma. (C) Sample intracellular responses from validation protocols using the seed with the best intracellular score. Black traces are the ground-truth responses. (D) Normalized extracellular templates (normalized to the respective template signal-amplitude peak-to-peak value), computed from the fire pattern (120%) validation protocol, for the ground-truth model (black) and optimized models (seed with best intracellular score).

Figure 4:

Optimization results on ground-truth model. (A) Distributions of intracellular (left) and extracellular (right) scores (i.e., sum of all feature scores) for 10 seeds for each applied strategy (soma: blue; all: orange; sections: green; single: red). (B) Distributions of intracellular (left, N = 86) and extracellular (right, N = 796) feature scores for the solution with the best intracellular total score. The *** represent a p < 10−5 with respect to soma. (C) Sample intracellular responses from validation protocols using the seed with the best intracellular score. Black traces are the ground-truth responses. (D) Normalized extracellular templates (normalized to the respective template signal-amplitude peak-to-peak value), computed from the fire pattern (120%) validation protocol, for the ground-truth model (black) and optimized models (seed with best intracellular score).

Close modal
Figure 5:

Average AP distribution. Action potentials at different locations along the neuron morphology. Each inset displays the average AP at a specific location on the left (marked on the neuron morphology) and the mean relative error with respect to the ground-truth values of four commonly used features (AP amplitude, AP half-width, decay time, and AHP) for the four different optimization strategies (ground truth: black; soma: blue; all, orange: sections, green; single: red). The vertical dotted gray line marks the time of the somatic action potential peak. For the majority of locations, the extracellular strategies yielded the better AP waveforms fits than the soma strategy.

Figure 5:

Average AP distribution. Action potentials at different locations along the neuron morphology. Each inset displays the average AP at a specific location on the left (marked on the neuron morphology) and the mean relative error with respect to the ground-truth values of four commonly used features (AP amplitude, AP half-width, decay time, and AHP) for the four different optimization strategies (ground truth: black; soma: blue; all, orange: sections, green; single: red). The vertical dotted gray line marks the time of the somatic action potential peak. For the majority of locations, the extracellular strategies yielded the better AP waveforms fits than the soma strategy.

Close modal
Figure 6:

Experimental data. For each recorded neuron (cells 1 and 2), we display the maximum z-projection of the z-stack (top left), the reconstructed morphology overlaid over the extracellular template (top center: AIS in green, axon in blue), the Alexa Fluor 594 dye, and AIS-specific staining (bottom) and sample somatic traces of one experimental run (right).

Figure 6:

Experimental data. For each recorded neuron (cells 1 and 2), we display the maximum z-projection of the z-stack (top left), the reconstructed morphology overlaid over the extracellular template (top center: AIS in green, axon in blue), the Alexa Fluor 594 dye, and AIS-specific staining (bottom) and sample somatic traces of one experimental run (right).

Close modal
Figure 7:

Optimization results for the experimental cell 1 model. (A) Distributions of intracellular (left) and extracellular (right) scores (i.e., sum of all feature scores) for 10 seeds for each applied strategy (soma: blue; all: orange). (B) Sample intracellular responses using validation protocols of the seed with the best intracellular score. Black traces show patch-clamp recordings. (C) Neuron morphology (gray: soma + denrites; green: AIS; blue: axon) and normalized extracellular templates (black: experimental; blue: soma strategy; orange: all strategy), computed using the fire pattern (120%) validation protocol for experimental data (black) and the optimized model (seed with best intracellular score). The insets highlight the experimental, soma, and all extracellular waveforms on two channels.

Figure 7:

Optimization results for the experimental cell 1 model. (A) Distributions of intracellular (left) and extracellular (right) scores (i.e., sum of all feature scores) for 10 seeds for each applied strategy (soma: blue; all: orange). (B) Sample intracellular responses using validation protocols of the seed with the best intracellular score. Black traces show patch-clamp recordings. (C) Neuron morphology (gray: soma + denrites; green: AIS; blue: axon) and normalized extracellular templates (black: experimental; blue: soma strategy; orange: all strategy), computed using the fire pattern (120%) validation protocol for experimental data (black) and the optimized model (seed with best intracellular score). The insets highlight the experimental, soma, and all extracellular waveforms on two channels.

Close modal
Figure 8:

Optimization results on the experimental Cell 2 model. (A) Distributions of intracellular (left) and extracellular (right) scores (i.e., sum of all feature scores) for 10 seeds for each applied strategy (soma: blue; all: orange). (B) Sample intracellular responses using validation protocols of the seed with the best intracellular score. Black traces show patch-clamp recordings. (C) Neuron morphology (gray: soma + denrites; green AIS; blue: axon) and normalized extracellular templates (black: experimental; blue: soma strategy; orange: all strategy), computed using the firepattern (120%) validation protocol for experimental data (black) and the optimized model (seed with best intracellular score). The insets highlight the experimental, soma, and all extracellular waveforms on two channels.

Figure 8:

Optimization results on the experimental Cell 2 model. (A) Distributions of intracellular (left) and extracellular (right) scores (i.e., sum of all feature scores) for 10 seeds for each applied strategy (soma: blue; all: orange). (B) Sample intracellular responses using validation protocols of the seed with the best intracellular score. Black traces show patch-clamp recordings. (C) Neuron morphology (gray: soma + denrites; green AIS; blue: axon) and normalized extracellular templates (black: experimental; blue: soma strategy; orange: all strategy), computed using the firepattern (120%) validation protocol for experimental data (black) and the optimized model (seed with best intracellular score). The insets highlight the experimental, soma, and all extracellular waveforms on two channels.

Close modal
The final modification that we introduced to the overall pipeline concerned the simulator that was used to generate the simulated responses in the evaluation phase. The NEURON software (Carnevale & Hines, 2006) is the gold standard to simulate the intracellular dynamics of single-cell models. With the addition of a Python interface (Hines et al., 2009), the NEURON simulator has been included in the BluePyOpt (Van Geit et al., 2016) software framework, designed to implement the fitting pipeline, introduced in the previous section, with an easy-to-use and flexible API. In order to simulate extracellular signals in addition to intracellular ones, we extended the BluePyOpt framework (see section 5.1 for details). We added a simulator interface based on the LFPy (Lindén et al., 2014; Hagen et al., 2018) Python package. LFPy is a wrapper to the NEURON software that uses the transmembrane currents from each neuron compartment Ii, obtained via the cable equation, to compute extracellular potentials ϕ at the respective electrode positions,
φ(rj,t)=i14πσIi(t)drirj-ri,
(3.2)
where σ is the extracellular conductance (set to 0.3 S/m; Goto et al., 2010), and ri and rj denote the positions of compartment i and electrode j, respectively. To model the physical size of the recording electrodes, we used the so-called disk approximation (Hagen et al., 2018), which averages the extracellular potential over n points on the surface of the electrode (we used n = 10).

3.2  Validation of the Multimodal Fitting Strategy on a Ground-Truth Model

After extending the fitting framework to include extracellular signals, we sought to quantitatively assess whether the inclusion of extracellular features improved the overall model fitting performance. To do so, we replicated our planned experiment in silico, using an existing multicompartment model with known parameters as ground truth. Although model-to-model fitting is known to be easier than fitting a model to experimental data (Druckmann et al., 2008), this first validation step is required to assess the potential of the proposed multimodal fitting procedure.

The ground-truth model was built using an available morphology of a rat somatosensory cortex layer 5 thick-tufted pyramidal cell (L5PC; Hay et al., 2011). In the original model, the axon was substituted by a linear stub axon of 60 μm; we introduced instead an axon initial segment (AIS) section, expressing AIS-specific conductances (Hu et al., 2009; Shu et al., 2007; Kole et al., 2007), namely, Nav1.2, Nav1.6, and Kv1 (see section 5.3 for details). The AIS length was set to 35 μm following the reconstructed axon morphology with a subsequent 1 mm long myelinated axonal section with a diameter of 0.2 μm (see Figure 3B, left). The AIS ion channels in the model moved the initiation of action potentials in the distal part of the AIS (see Figure 3B, right), in accordance with a large body of experimental and computational evidence (Hallermann et al., 2012; Bakkum et al., 2019; Huang & Rasband, 2018). The AIS properties were also strongly reflected in the extracellular signals. The AIS was found to be the dominant contributor to extracellularly measured potentials due to the strong membrane currents generated to initiate an action potential (Teleńczuk et al., 2018; Bakkum et al., 2019). The ground-truth model was placed on a 20-row-by-4-column planar simulated MEA with 50 μm electrode-to-electrode center pitch. The extracellular template was computed with equation 3.2 and averaged over several action potentials in response to a step protocol. Figure 3A shows the ground-truth neuron sitting on top of the MEA model. The black lines are the extracellular template at the electrode locations. The red trace depicts the signal with the largest amplitude, which appears on the electrode closest to the AIS. The ground-truth model includes 33 free parameters (10 somatic, 13 AIS, 6 apical, and 4 basal parameters) that need to be optimized. See section 5.3 and Tables 2 and 3 for details on the ground-truth model definition and parameters.

Table 1:

Optimization Strategies.

strategy name# intracellular features# extracellular features# total features
Soma 76 — 76 
Single 76 77 153 
All 76 11 87 
Sections  76 33 109 
strategy name# intracellular features# extracellular features# total features
Soma 76 — 76 
Single 76 77 153 
All 76 11 87 
Sections  76 33 109 

Note: Summary of the four strategies used for fitting the models, specifying the intracellular, extracellular, and total number of features used for optimization (the numbers refer to the ground-truth model).

After checking that our model could reliably reproduce AIS-specific physiological features, which strongly influence the extracellular electrical-potential landscape, we fully replicated the data acquisition phase in silico. We utilized a set of protocols named eCode (Markram et al., 2015; Iavarone et al., 2019) to probe the behavior of a patched neuron against several different stimuli (see Figure 3C). The eCode protocols included several long (>1 s) step stimuli with different lengths and amplitudes (firepattern, IV, IDrest) to monitor firing behavior, sub-threshold responses, and bursting properties, as well as short-pulse steps (50 ms) to observe details of AP properties (APWaveform). Additional more complex stimuli included a hyperpolarizing step, followed by depolarization (HyperDepol), a three-step protocol to monitor after-hyperpolarization potentials (AHP) after several spikes (sAHP), and, finally, a series of triangular stimuli of different lengths (PosCheops). All current amplitudes were computed with respect to the neuron’s rheobase current (estimated using a 270 ms step). For our ground-truth model, the rheobase current was 150 pA. See section 5.2 for details about amplitudes and timings of the eCode protocols.

From the eCode-protocol ground-truth responses, we then extracted features before running the optimization step. We selected a subset of protocols that were used to fit the model (i.e., training), and we left the other protocols for validation to assess how general the solutions were. For both, the in silico and the experimental data, we used six protocols for training: IDrest (150%), IDrest (250%), IDrest (300%), APWaveform (290%), IV (−100%), and IV (−20%) (the amplitude is expressed in percentage of the rheobase current).

A total of 76 intracellular features were extracted from these selected protocols and used for optimization assuming that only the patch-clamp data were available (soma fitting strategy). In addition to these features, the all strategy (using all electrodes and the cosine distance in equation 3.1 as features) adds 11 extracellular features (total = 87). The sections strategy added 33 (three sections times 11 extracellular features), and the single strategy added 77 features (seven electrodes times 11 extra features). The number of intra- and extracellular features for each optimization strategy are summarized in Table 1. See section 5.4 for details on feature extraction according to different protocols.

For each optimization strategy, we ran 10 optimizations, each starting from a different seed using a covariance matrix adaptation (CMA) optimizer (see section 5.5 for details) (Damart et al., 2020). Figure 4A shows the fitness of the 10 optimized solutions for each strategy with respect to intracellular (left) and extracellular features (right: features were computed with the all strategy for consistency). While the soma, all, and sections strategies appeared to yield low scores (i.e., better fitness) for intracellular features, the single strategy produced solutions with worse intracellular fitness. Conversely, the single strategy yielded the best performance with respect to extracellular features, while the worst performance was achieved by using the soma strategy. These results reflect the balance between intracellular and extracellular features used for optimization (see Table 1). The single strategy used almost as many extracellular features as intracellular ones, which may have pushed the optimization to favor extracellular fitting at the expenses of intracellular fitness. In contrast, the use of all and sections strategies provided a good balance between intracellular and extracellular fitness.

Next, we report the results of the solution with the best intracellular fit for each strategy. Figure 4C displays intracellular responses for some validation protocols. All of these solutions yielded a good fitting of the intracellular traces used for validation (firepattern, HyperDepol, sAHP, and PosCheops). Interestingly, for the HyperDepol stimuli, the soma, all, and sections strategies failed to reproduce the absence of spikes observed in the depolarization phase (second row: HyperDepol (−160%). Only the single strategy correctly did not generate a spike for HyperDepol (−160%) and did generate one for HyperDepol (−40%). The overall distribution of scores for all intracellular features, computed on validation protocols (total = 86 features), is shown in Figure 4B, left. No significant differences were detected across strategies (post hoc Mann-Whitney tests with Holm p-value correction). To quantify extracellular performance, we computed all features using the single strategy on all 80 electrodes (total of 796 features after excluding features that resulted in invalid values—for example, when a clear signal peak was not available). The right panel of Figure 4B shows the distribution of these extracellular features across the different strategies. All strategies making use of extracellular data significantly outperformed the soma strategy (post hoc Mann-Whitney tests with Holm p-value correction: soma VS all - p < 10−5, soma VS sections - p < 10−5, soma VS single - p < 10−5). Figure 4D shows the normalized extracellular templates (each channel value has been normalized to the respective template signal-amplitude peak-to-peak value) for the solution yielding the best extracellular fit. It is evident at first glance that the three extracellular strategies achieve a better fit for almost all electrodes, while the solution from the soma strategy, despite providing a very good fit for intracellular features and traces, failed to capture the complex dynamics of extracellular signals in many recording channels.

Extracellular signals indirectly reflect membrane potentials across the entire neuron morphology, as the latter drive transmembrane currents, which are responsible for the generation of extracellular potentials. Having a ground-truth model at our disposal, we could also analyze the intracellular AP distribution, that is, the AP waveforms across different neurites. We first took the fire pattern (120%) response of the ground-truth model and the optimized models, detected action potentials in the somatic trace, and finally used the detected AP peaks to average the action potentials of different neurites (the first and last spikes were excluded as they could contain artifacts originating from stimulus onset and offset). Once the average AP was calculated, we computed the mean relative error of four commonly used features with respect to the ground-truth AP: AP amplitude, AP half-width, decay time, and AHP (afterhyperpolarization). Figure 5 shows the morphology of the ground-truth neuron and the recording positions on the neurites. Each inset displays the AP at the respective location obtained with the ground-truth (black), soma (blue), all (orange), sections (green), and single (red) models (left) alongside the mean relative errors with respect to the ground-truth AP (right). The selected locations included 4 points along the apical dendrite (apical left, right, middle, proximal), the distal and proximal AIS locations, and two locations on the basal dendrites (basal left, right). For all tested locations, the all strategy produces lower errors than the strategy. For all dendritic locations (both apical and basal), all extracellular strategies (all, single, sections) outperform the soma strategy. This suggests that combining intracellular and extracellular features could provide a better representation of the AP waveforms across the entire neuronal arbor.

In this section, we performed an in silico experiment using a ground-truth model with known parameters to investigate whether the inclusion of extracellular features (using different strategies) could be beneficial to construct multicompartment models. We have shown that adding extracellular features, in particular by using the all and sections strategies, improved the optimization performance in terms of fitting intracellular signals, extracellular templates, and action potential distribution across the entire neuronal morphology. Note that this scenario is ideal because we imposed that the ground-truth solution lies within the parameter space (i.e., all the mechanisms and respective parameter boundaries are “correct”) and the forward model used in the optimization to generate extracellular signals was the same as the one used to generate ground-truth signals.

3.3  Cell Models from Cultured Neurons

After assessing the benefits of combining intra- and extracellular signals with in silico experiments, we looked into experimental data. We acquired paired patch-clamp and extracellular signals from two cells of cortical cultures in embryonic rats. After the electrophysiological data acquisition, cells were imaged with a confocal microscope and fixated. We then stained the cultures with AIS-specific antibodies to locate the axon initial segment in the morphology. Figure 6 shows the two cells that we used to construct multicompartment models, including a maximum intensity z-projection of the z-stack (top left), the reconstructed morphology overlaid on the extracellular template (top middle: AIS is green, axon is blue, other neurites in gray), sample somatic patch-clamp traces from one of the runs (right), and AIS stainings (bottom: left panel shows Alexa Fluor, the center panel displays the AIS-specific channel yielding the best signal—Kv7.3 for cell 1 and ankyrin-G for cell 2—and the right panel shows the merged image with an indication of the AIS location). Note that we stained with both Kv7.3 and ankyrin-G, and we report here the brightest indicator. Given the early developmental stage of the cultures (experiments performed between DIVs 18 and 21), there was no clear and unique apical dendrite. For this reason, when reconstructing the morphology, we did not differentiate between apical and basal dendrite, and we applied the same mechanisms to all dendritic sections. Moreover, for both cells, the AIS did not originate directly from the soma but from one of the dendrites, which is termed an axon-bearing dendrite (ABD; Thome et al., 2014; Goaillard et al., 2020). Despite these differences, we used the same dendritic mechanisms of the ground-truth model detailed in the previous sections and considered the ABD as a “normal” dendrite section. Another important difference with respect to the ground-truth L5PC model is that cultured cells do not exhibit myelinated axons. We therefore used active axonal mechanisms without myelination (see section 5.6 for details).

Figures 7 and 8 show the results of the optimizations performed on cells 1 and 2 in Figure 6. Panel A displays the intracellular (left) and extracellular (right) scores for the optimized models starting from 10 different seeds for the soma (blue) and all (orange) strategies (we report results only for the all strategy for the extracellular options, as it showed the best performance in the in silico tests). For both cell models, the application of the soma strategy yielded better intracellular scores, while the use of the all strategy resulted in better extracellular fitting. The scores obtained were almost twice the ones we got for the fitting of the ground-truth model in Figure 4 (intracellular score: GT fitting: 62.8 ± 13.6 for soma and 59.3 ± 17.0 for all; cell 1 fitting: 123.9 ± 11.2 for soma and 149.2 ± 47.2 (excluding one seed that did not converge) for all; cell 2 fitting: 152.8 ± 9.22 for soma and 177.6 ± 20.4 for all). These discrepancies can be partly explained by the fact that for the in silico validation, the optimization procedure was performed with a correct cell model using the same mechanisms and distributions for optimization that have been used for the ground-truth model. On the contrary, for the experimental cells, we did not know the exact ion channels expressed by the neurons and used the mechanisms present in the ground-truth model. In fact, the experimental cells (putative excitatory neurons recorded in a DIVs 18-24 cell culture) and the ground-truth model (designed to fit a L5PC recorded in a P14-P20 slice) are likely to express different ion channels. Despite these apparent differences, the models fitted on the experimental data could reproduce the main features of the intracellular spiking patterns of the validation protocols (panels B, with the exception of the PosCheops (300%) response for the soma strategy in cell 2; see Figure 8B).

Panels C in Figures 7 and 8 show the experimental extracellular templates (black) and the resulting simulated templates using the soma (blue) and all (orange) fitting strategies on top of the reconstructed cell morphologies. The insets show a close-up of two representative channels for each cell. The shapes of the fitted extracellular templates of both strategies are qualitatively similar to the experimental templates. The shape characteristics of the waveforms appear to be reproduced by the optimized cell models, but mismatches are still apparent, in particular for cell 2 (e.g., in the area below the soma in Figure 8). The worse fit of cell 2 with respect to cell 1 is also reflected by the overall higher extracellular scores (right box plots of panels B in Figures 7 and 8).

To investigate this further, we took a closer look at the nature of the mismatch for cell 2. Specifically, focusing on the left panel of Figure 8C, one can see that neither of the strategies, soma and all, can reproduce the initial negative phase followed by the strong, positive upstroke observed experimentally. We hypothesize that this is due to a wrong weight of different compartment contributions to the extracellular signal. We therefore virtually patched the all cell model at the soma, AIS, and basal dendrite locations and showed that indeed the mismatch could be explained by an underestimation of the somatic contribution and an overestimation of the AIS one to the EAP (extracellular action potential; for further details, see Figure 10). From a modeling perspective, the mismatch could be caused by several factors: (1) an underestimation of the size of the soma, resulting in lower transmembrane currents (the overall current is proportional to the compartment’s area); (2) the optimization boundaries chosen for the AIS and the soma are not “correct,” resulting in an AIS that is too excitable; and (3) the distance between the dendrite and the electrode plane is overestimated, yielding lower dendritic contributions to the extracellular potential, which in turn reduces the positive upstroke.

Considering the confounding factor of probable mismatches between experimentally recorded signal traces and the model traces based on the mechanisms included in the model fits, it is virtually impossible to state whether any of the strategy (soma or all) is better than the other in this case. Nevertheless, the inclusion of extracellular features in the fitting procedure does not impede the optimization to find acceptable models.

We have presented a multimodal strategy to construct biophysical multicompartment models using a combination of patch clamp and high-density microelectrode arrays (HD-MEAs). After introducing the general concepts of fitting multicompartment models, we presented the extensions that were required to include extracellular signals in the procedure. We then validated the approach using an existing ground-truth model with known parameters and showed that the inclusion of extracellular data improves the overall fitting performance, specifically in terms of fitting extracellular signals and the membrane potentials over the entire neuron morphology. Finally, we applied the validated method to experimental data from two cultured cortical neurons and demonstrated how to build cell models either fitted on the base of patch-clamp data alone or by applying a combined patch-MEA fitting strategy.

While the suggestion to use extracellular data as a source to construct multicompartment models dates back more than 10 years ago (Gold et al., 2007), the presented work is the first concise attempt to include high-quality extracellular data into the fitting of single-neuron models. Our in silico validation provides convincing evidence that extracellular data can be used to improve the optimization performance and may provide better-validated models. Nevertheless, including extracellular measurements as an additional data source is challenging and may entail various issues that merit discussion. First, a setup to perform patch clamp and simultaneous extracellular recordings needs to be available. It requires the use of an upright microscope due to the opaqueness of the MEA substrate. Moreover, the presence of the MEA within the patch-clamp setup may provide additional noise sources for both systems, potentially yielding lower-quality data, although several approaches to reduce noise are available. Second, in order to record reliable and large-amplitude extracellular signals, the patched cell needs to be located in close proximity to the extracellular electrodes (Buzsáki, 2004). The closer the cell, the larger the signal and, therefore, the signal-to-noise ratio. While proximity is not an issue when working with 2D cell cultures where (ideally) cells form a monolayer on top of the MEA, the use of thicker preparations, such as brain slices, may add complications: the patch experimenter will need to blindly target deep cells that lie close to the surface. However, techniques to use the HD-MEA recordings to estimate the precise patch pipette location help to render the cell targeting easier (Allen et al., 2018; Obien et al., 2019). Another attractive possibility to use slices is to use organotypic slice cultures (Egert et al., 1998; Gong et al., 2016), where slices are cultured for several weeks on an MEA device and thin down over time. Finally, for imaging the patched neurons, a precise registration between the 3D electrode locations and the neuron morphologies is required in order to correctly assign extracellular signals. Due to the relatively low resolution in the z-dimension of the microscope (in our case 0.4 μm), the electrode registration may include errors that may cause mismatches between experimentally recorded and reconstructed extracellular signals. In order to minimize this mismatch, we imaged the cell right after the patch experiment without a staining step, which causes shrinkage. However, there could still be a movement of the soma with respect to the MEA due to the retraction of the pipette. We tried to minimize this effect by carefully and slowly removing the pipette from the cell and targeting neurons that were not too isolated, so that the neighbor cells could help maintain the neuron during the retraction.

The forward-modeling framework used to simulate extracellular signals from transmembrane currents is well established and widely accepted (Einevoll et al., 2013). However, this framework makes several assumptions that may not be fully satisfied under experimental conditions. The extracellular milieu is assumed to be infinite, isotropic, homogeneous, and ohmic. Clearly, the sample is not infinite, and it is bounded by the well and the substrate on which the cells are cultured. However, considering that the reference electrodes are at the borders of the sensing area and we targeted cells lying in the center of the HD-MEA (Frey et al., 2009), the infinity conditions can be considered to be partially met. The isotropy and homogeneity assumption are harder to relax: For in vitro cultures, the cells mainly are arranged in a monolayer on top of the substrate. Therefore, the arrangement of cells and medium is not homogeneous and the conductivity of the tissue is not isotropic. Cell density is high and conductivity is low along lateral directions within the cell layer. Cell density becomes lower, and conductivity increases on moving perpendicular to the cell layer into the medium. Moreover, the MEA substrate, on which the cell layer is arranged, is an insulating material, while the medium is conductive. However, different analytical solutions can be used to correct for the anisotropy of the tissue (Hagen et al., 2018) and discontinuities, for example, using the method of images (Ness et al., 2015), which corrects for the amplitudes of the extracellular potentials. We did not explicitly include these modeling extensions in the forward model, but we did account for the presence of the MEA substrate in computing extracellular features by using relative amplitudes instead of absolute ones. The adoption of the method of images would only affect the overall amplitude of the extracellular signals; therefore, the computation of features based on relative amplitudes would not be influenced by different modeling schemes. The assumption that the medium shows ohmic behavior seems to be justified in the frequency range of extracellular signals (Goto et al., 2010). Besides assumptions concerning the extracellular medium, the cable equation that has been used to calculate transmembrane currents assumes that extracellular potentials outside the cell membrane are constant (i.e., extracellular axial currents are zero), and self-ephaptic phenomena (the effects of the self-generated extracellular potentials on the neuron itself) are therefore ignored (Buccino et al., 2019). This assumption is widely accepted due to the very low amplitude of extracellular signals in comparison to intracellular ones. Some of the simplifications and assumptions mentioned above could partially account for the mismatch between simulated and experimentally recorded signals. The use of more advanced and sophisticated modeling frameworks, such as the finite element method (Agudelo-Toro & Neef, 2013; Tveito et al., 2017; Buccino et al., 2019, 2021), could help to overcome some of the simplifications in that MEA substrate, and intracellular and extracellular space can be specifically modeled. However, the use of such methods would be computationally very expensive for performing model optimizations with several thousands of simulations until convergence to a stable solution. A probably major contributor to the discrepancies observed between experimental and fitted data (see Figures 7 and 8) may originate from the model definition that we used to reconstruct the data of the experimental cells. We used the model structure of a layer 5 pyramidal cell recorded in a P20 slice to fit data recorded of putative pyramidal cells recorded in culture at DIV 18-24. Clearly, the developmental path of neurons in a living animal (from which the slice was obtained) differs from that of cultured neurons in a dish starting from embryonic state. Therefore, the use of an L5PC model to fit neurons of an embryonic cell culture is not ideal. In addition, we did not explicitly model the axon-bearing dendrite (ABD), but rather treated it as a normal dendrite. It is likely that the ABD is more excitable than the rest of the dendritic tree, and therefore modeling it explicitly could produce better-fitting models. Nevertheless, to the best of our knowledge, there are no models for cell-cultured neurons available in the literature. Another possible line of improvement would be to extend on the optimization itself. The extended feature space proposed here could benefit from recent optimization approaches that, instead of finding a single “optimal” solution, aim at finding solution subspaces (Gonçalves et al., 2019; Schneider et al., 2023) and/or from including additional constraints to reduce parameter degeneracy, such as energy efficiency (Jedlicka et al., 2022).

As the main goal of this work was not to build models from cultured neurons but to investigate and develop a methodology and software tools to combine multimodal data sources in a fitting procedure, the application to experimental data serves only as a proof of concept.

The main intention of introducing the multimodal method here was to better constrain the many parameters of multicompartment models by using other information-rich data sources. While most of the available models are constructed by using electrophysiological data of a single somatic patch-clamp experiment, a simultaneous patch recording from several neurites, especially dendrites, can shed light on activities and nonlinear properties of dendritic regions (Larkum et al., 2001; Hay et al., 2011; Almog & Korngreen, 2014). Despite such multi-pipette patch techniques being experimentally very challenging, one of their main advantages is that one can stimulate at both the soma and the dendrite location to induce calcium spikes and fit the model to reproduce the associated phenomena. While the combined use of an HD-MEA with a somatic patch clamp does not provide intracellular access to dendritic neurites, one can exploit the extracellular stimulation capability of HD-MEAs to selectively stimulate dendrites to induce calcium spikes, the effects of which then can be observed in the soma (Hay et al., 2011). Here, we used high-density extracellular microelectrode arrays as an additional recording modality to probe the spatiotemporal distribution of electrophysiological signals across the entire neuron morphology. However, other multimodal procedures could potentially constitute attractive alternatives. For example, one could combine the patch-clamp technique with optical electrophysiology readouts to obtain a proxy of the intracellular signals from many different regions of a neuron simultaneously. Genetically encoded voltage indicators (GEVIs; Abdelfattah et al., 2019) provide indirect membrane potential readout of neuronal compartments that can be used to extract region-specific features for optimization. One drawback of this approach, however, is that neurons would need to be genetically modified to express the fluorescent indicators, which possibly affects electrophysiological characteristics and spontaneous firing of the cell. In considering the low yield of patch-clamp experiments, the question arises whether it is possible to build multicompartment neuronal models without patch clamp. The two main advantages of the patch-clamp technique are the capability of applying precise stimuli to the neuron through multiple protocols (both supra- and subthreshold) and filling the neuron with fluorescent markers for subsequent imaging to reconstruct its morphology (an essential component of the model). The combination of GEVI imaging and HD-MEA could represent another attractive multimodal approach that can provide similar features. Electrical stimulation with extracellular electrodes can be used to apply different stimuli to the cell, with the advantage that subcellular compartments of the neuron (Ronchi et al., 2019) can be targeted; in addition, extracellular stimulation could also be used to generate patch-clamp-like sustained stimuli (Abbott et al., 2020), but probably with lower control and precision than patch-clamp stimulation. The GEVI readout can provide indirect membrane potential measurements across multiple neural compartments, as outlined above, and can also be used to image the target neuron at high resolution for morphology reconstruction owing to the brightness and stability of recently developed indicators (Abdelfattah et al., 2019).

In summary, we developed and investigated a novel approach to combine multiple data acquisition modalities, specifically patch clamp and HD-MEAs, to construct more accurate multicompartment models that can be better validated. As novel neurotechnologies are developed, it can be foreseen that new multimodal strategies will become available to acquire even richer data sets for building and constraining single-cell multicompartment models.

5.1  Extensions to the Fitting Framework

As part of this project, we extended the open-source BluePyOpt Python framework (Van Geit et al., 2016), designed to fit neuroscience models, to enable us to use extracellular signals in the fitting procedure.

5.1.1  Simulator Back End

BluePyOpt originally implemented a NEURON-based cell model and simulator (CellModel and NrnSimulator classes). We implemented two additional classes, the LFPyCellModel and LFPySimulator, that use instead LFPy (Lindén et al., 2014; Hagen et al., 2018) as a back end to construct and simulate neuron models. The LFPySimulator can be instantiated with a recording probe object from the MEAutility package (Buccino & Einevoll, 2021) and computes extracellular signals generated by the neuron’s transmembrane currents (see equation 3.2). In addition, a stimulus (LFPySquarePulse), response (TimeLFPResponse), and recording (LFPRecording) classes were added to the framework to support extracellular simulations over the entire pipeline.

5.1.2  Extracellular Feature Extraction

We additionally extended the feature-extraction capabilities of BluePyOpt to include extracellular features. These features are based on the extracellular template, the average extracellular action potential. We introduced a new class, extraFELFeature, that preprocesses the responses to obtain the extracellular template and computes extracellular features from a specified protocol (in our case, we used IDrest (300%); see the eCode protocols section below). Note that in the ideal and noise-free simulation, it is sufficient to use one protocol to compute the extracellular template, while for experimental data, we combine several protocols to obtain a cleaner template. The extracellular signals at the respective electrode locations are first interpolated to match the sampling frequency of experimentally recorded extracellular traces (in our case, 20 kHz). Optionally, the interpolated traces can also be filtered with a zero-phase filter to mimic the filter applied to experimental traces (in our case, we used a Butterworth bandpass filter with cutoff frequencies at 300 and 6000 Hz). The somatic response was used to reliably identify peak times using the eFEL “peak_time” feature (the first and last spike can be discarded because they may contain artifacts of stimulus onsets and offsets). Peak times were used to cut out snippets of the extracellular traces, which were averaged and optionally upsampled to provide more precision in feature values. We used 2 ms before and 5 ms after the peak time as cutouts and upsampled the template 10 times.

When the template was extracted, extracellular features could be calculated. We compiled a list of 11 features that have been included in the BluePyOpt package (Anastassiou et al., 2015). Extracellular features are either absolute, computed for each channel separately, or relative, computed with respect to the channel with the largest extracellular signal amplitude:

  • Peak-to-valley duration (absolute): Time in seconds between the negative and positive peaks. If the negative peaks precede the positive one, the value of the feature is positive. Conversely, when the positive peak precedes the negative one, the value is negative.

  • Halfwidth (absolute): Width of waveform at half of its amplitude in seconds. If the positive peak precedes the negative one, the value is negative. This procedure helps to maximize the shape information carried by the feature value.

  • Peak-to-trough ratio (absolute): The ratio between positive and negative peaks.

  • Recovery slope (absolute): After depolarization, the neuron repolarizes until the signal peaks. The recovery slope is the slope of the action potential after the peak, returning to the baseline in dV/dT. The slope is computed within a user-defined window after the peak (default = 0.7 ms).

  • Repolarization slope (absolute): After reaching its maximum depolarization, the neuronal potential will recover. The repolarization slope is defined as the dV/dT of the action potential between the negative peak and the baseline.

  • Negative peak relative amplitude (relative): The relative amplitude of the negative peak with respect to the negative signal peak of the channel with the largest amplitude. For the largest-amplitude channel, this feature has a value of 1.

  • Positive peak relative amplitude (relative): The relative amplitude of the positive signal peak with respect to the positive signal peak of the channel with the largest amplitude. For the largest-amplitude channel, this feature has a value of 1.

  • Negative peak time difference (relative): The time difference between the negative signal peak with respect to the negative signal peak of the channel with the largest amplitude. For the largest-amplitude channel, this feature has a value of 0. Note that values can also be negative if the respective negative signal peak occurs earlier than the negative signal peak on the largest-amplitude channel.

  • Positive peak time difference (relative): The time difference between the positive peak with respect to the occurrence of the positive peak of the channel with the largest amplitude. For the largest-amplitude channel, this feature has a value of 0. Note that values can also be negative if the respective positive signal peak occurs earlier than the positive signal peak on the largest-amplitude channel.

  • Negative image (relative): Voltage values at the time of the negative signal peak on the largest-amplitude channel. The values are normalized by the negative signal-amplitude value on the largest-amplitude channel.

  • Positive image (relative): The voltage values at the time of the positive signal peak on the largest-amplitude channel. The values are normalized by the positive signal-amplitude value on the largest-amplitude channel.

5.2  eCode Protocols

The eCode protocols used for both the in silico validation and experimental data acquisition included a series of stimuli designed to probe the neuron’s behavior under a wide range of conditions (Markram et al., 2015; Iavarone et al., 2019). After a neuron was patched, the holding current Iholding was manually adjusted so that the read-out membrane potential was −70 mV (which corresponded to approximately −84 mV after liquid junction potential correction (see section 5.6.4). A series of step stimuli of 270 ms duration with increasing amplitudes were then applied to roughly estimate the neuron’s rheobase current Irheo, the current at which a single action potential was induced. The subsequent stimulus amplitudes are reported in percent values of the estimated Irheo current.

  • IDthresh: The purpose of this protocol was to finely screen the current levels to correctly reestimate Irheo before processing. It consisted of 21 sweeps with a 270 ms square pulse of increasing amplitudes from 50% to 130% with step increments of 4%.

  • Firepattern: This protocol was aimed at characterizing the firing properties of the neuron. It consisted of two step pulses of 3.6 s at 120% and 200%.

  • IV: The purpose of this protocol was used to check on sub-threshold properties of the neuron (e.g., the sag). The protocol consisted of 11 step stimuli of 3 s ranging from −140 to 20% with increments of 20%.

  • IDrest: This protocol was used to characterize the input/output curve of the neuron. Eleven step stimuli of 1.35 s were delivered with amplitudes ranging from 50% to 300% in increments of 25%.

  • APWaveform: This protocol was defined to precisely record one to two action potentials at high sampling rates (typically 200 kHz). It consisted of short 50 ms pulses with amplitudes from 200% to 350% and step increments of 30%.

  • HyperDepol: This protocol was aimed at looking at cell excitability at a hyperpolarization potential. It consisted of four sweeps with hyperpolarizing steps of 450 ms from −40% to −160% (increment step 40%), followed by a 270 ms 100% step.

  • sAHP: This protocol was aimed at characterizing after-hyperpolariza-tion potentials (AHPs) after several spikes. It consisted of four sweeps starting with a 250 ms depolarizing step at 40%, immediately followed by a larger-amplitude phase of 225 ms, ranging from 150% to 300% with 50% step increments, and a third 450 ms phase back at 40% amplitude.

  • PosCheops: This final protocol was used to characterize the neuron’s response to ramp stimuli. It contained three ramps of up to 300% amplitude of different duration: the first ramp featured ramp-up and -down phases of 4 s, the second ramp 2 s, and the third ramp 1.33 s.

All protocols (excluding IDthres, which was not used for training or validation) are shown in Figure 3D. For experimental data, the eCode protocols were run four to six times.

Table 2:

Ground-Truth Model Parameters (Nondistribution Type) and Optimization Bounds.

ParameterGT valueSection listBoundsUnits
v_init −72  — mV 
celsius 34  — °C 
e_pas −74.64792 all — mV 
cm somatic, ais, basal — μF/cm2 
cm apical — μF/cm2 
Ra 100 all — Ω-cm 
g_pas 0.00000226 myelinated — S/cm2 
cm 0.09091 myelinated — μF/cm2 
g_pas 0.0000248 somatic, ais, apical, basal [1e−5, 6e−5] S/cm2 
ek −90 somatic, ais, apical — mV 
ena 50 somatic, ais, apical — mV 
decay_CaDynamics_DC0 52.7290615 somatic [20, 300] ms 
gamma_CaDynamics_DC0 0.01682752 somatic [5e−3, 5e−2] — 
gCa_LVAstbar_Ca_LVAst 0.00412676 somatic [0, 0.01] S/cm2 
gCa_HVAbar_Ca_HVA2 0.00098005 somatic [0, 0.001] S/cm2 
gSKv3_1bar_SKv3_1 0.25953206 somatic [0, 1] S/cm2 
gSK_E2bar_SK_E2 0.08579915 somatic [0, 0.1] S/cm2 
gK_Tstbar_K_Tst 0.0108473 somatic [0, 0.1] S/cm2 
gK_Pstbar_K_Pst 0.14280937 somatic [0, 0.2] S/cm2 
vshiftm_NaTg 13 somatic — mV 
vshifth_NaTg 15 somatic — mV 
slopem_NaTg somatic — mV 
gNaTgbar_NaTg 0.27576985 somatic [0, 0.3] S/cm2 
g_pas 0.0000248 basal [1e−5, 6e−5] S/cm2 
gamma_CaDynamics_DC0 0.03485579 basal [0.005, 0.05 — 
gCa_LVAstbar_Ca_LVAst 0.00055044 basal [0, 0.001 S/cm2 
gCa_HVAbar_Ca_HVA2 0.0000959 basal [0, 0.0001 S/cm2 
g_pas 0.0000248 apical [1e−5, 6e−5] S/cm2 
cm apical — μF/cm2 
gamma_CaDynamics_DC0 0.00625451 apical [5e−3, 5e−2] — 
gSKv3_1bar_SKv3_1 0.0000133 apical [0, 3e−3] S/cm2 
vshiftm_NaTg apical — mV 
vshifth_NaTg apical — mV 
gCa_LVAstbar_Ca_LVAst 0.00016523 apical [0, 1e−3] S/cm2 
gCa_HVAbar_Ca_HVA2 0.0000361 apical [0, 5e−4] S/cm2 
gCa_LVAstbar_Ca_LVAst 0.00014743 ais [0, 0.01] S/cm2 
gCa_HVAbar_Ca_HVA2 0.00029333 ais [0, 0.001] S/cm2 
gSKv3_1bar_SKv3_1 0.26081743 ais [0, 1] S/cm2 
gSK_E2bar_SK_E2 0.07743135 ais [0, 0.1] S/cm2 
gK_Tstbar_K_Tst 0.16512831 ais [0, 0.2] S/cm2 
gK_Pstbar_K_Pst 1.03159443 ais [0, 2] S/cm2 
gNap_Et2bar_Nap_Et2 0.00231049 ais [0, 0.02] S/cm2 
decay_CaDynamics_DC0 288.16009 ais [20, 300] ms 
gamma_CaDynamics_DC0 0.02491751 ais [5e−3, 5e−2] — 
ParameterGT valueSection listBoundsUnits
v_init −72  — mV 
celsius 34  — °C 
e_pas −74.64792 all — mV 
cm somatic, ais, basal — μF/cm2 
cm apical — μF/cm2 
Ra 100 all — Ω-cm 
g_pas 0.00000226 myelinated — S/cm2 
cm 0.09091 myelinated — μF/cm2 
g_pas 0.0000248 somatic, ais, apical, basal [1e−5, 6e−5] S/cm2 
ek −90 somatic, ais, apical — mV 
ena 50 somatic, ais, apical — mV 
decay_CaDynamics_DC0 52.7290615 somatic [20, 300] ms 
gamma_CaDynamics_DC0 0.01682752 somatic [5e−3, 5e−2] — 
gCa_LVAstbar_Ca_LVAst 0.00412676 somatic [0, 0.01] S/cm2 
gCa_HVAbar_Ca_HVA2 0.00098005 somatic [0, 0.001] S/cm2 
gSKv3_1bar_SKv3_1 0.25953206 somatic [0, 1] S/cm2 
gSK_E2bar_SK_E2 0.08579915 somatic [0, 0.1] S/cm2 
gK_Tstbar_K_Tst 0.0108473 somatic [0, 0.1] S/cm2 
gK_Pstbar_K_Pst 0.14280937 somatic [0, 0.2] S/cm2 
vshiftm_NaTg 13 somatic — mV 
vshifth_NaTg 15 somatic — mV 
slopem_NaTg somatic — mV 
gNaTgbar_NaTg 0.27576985 somatic [0, 0.3] S/cm2 
g_pas 0.0000248 basal [1e−5, 6e−5] S/cm2 
gamma_CaDynamics_DC0 0.03485579 basal [0.005, 0.05 — 
gCa_LVAstbar_Ca_LVAst 0.00055044 basal [0, 0.001 S/cm2 
gCa_HVAbar_Ca_HVA2 0.0000959 basal [0, 0.0001 S/cm2 
g_pas 0.0000248 apical [1e−5, 6e−5] S/cm2 
cm apical — μF/cm2 
gamma_CaDynamics_DC0 0.00625451 apical [5e−3, 5e−2] — 
gSKv3_1bar_SKv3_1 0.0000133 apical [0, 3e−3] S/cm2 
vshiftm_NaTg apical — mV 
vshifth_NaTg apical — mV 
gCa_LVAstbar_Ca_LVAst 0.00016523 apical [0, 1e−3] S/cm2 
gCa_HVAbar_Ca_HVA2 0.0000361 apical [0, 5e−4] S/cm2 
gCa_LVAstbar_Ca_LVAst 0.00014743 ais [0, 0.01] S/cm2 
gCa_HVAbar_Ca_HVA2 0.00029333 ais [0, 0.001] S/cm2 
gSKv3_1bar_SKv3_1 0.26081743 ais [0, 1] S/cm2 
gSK_E2bar_SK_E2 0.07743135 ais [0, 0.1] S/cm2 
gK_Tstbar_K_Tst 0.16512831 ais [0, 0.2] S/cm2 
gK_Pstbar_K_Pst 1.03159443 ais [0, 2] S/cm2 
gNap_Et2bar_Nap_Et2 0.00231049 ais [0, 0.02] S/cm2 
decay_CaDynamics_DC0 288.16009 ais [20, 300] ms 
gamma_CaDynamics_DC0 0.02491751 ais [5e−3, 5e−2] — 

Note: Only parameters with optimization bounds are fitted.

Table 3:

Ground-Truth Model Parameters with Distributions and Optimization Bounds.

ParameterDistributionGT valueSection listRef.BoundsUnits
gNa12bar_Na12 11+expdistance-253×Value 6.15808 ais ais[0](0) [0, 8] S/cm2 
gNa16bar_Na16 1-11+expdistance-202.5×Value 4.39817 ais ais[0](0) [0, 8] S/cm2 
gkbar_Kd 1-11+expdistance-202.5×Value 0.04356 ais ais[0](0) [0, 2] S/cm2 
gIhbar_Ih (− 0.8696 + 2.087 × 0.00001 somatic, soma(0.5) — S/cm2 
 exp ((distance × 0.0031)) × Value  apical, basal    
gNaTgbar_NaTg exp (distance × (− 0.00352)) × Value 0.06642 apical soma(0.5) [0, 0.1] S/cm2 
ParameterDistributionGT valueSection listRef.BoundsUnits
gNa12bar_Na12 11+expdistance-253×Value 6.15808 ais ais[0](0) [0, 8] S/cm2 
gNa16bar_Na16 1-11+expdistance-202.5×Value 4.39817 ais ais[0](0) [0, 8] S/cm2 
gkbar_Kd 1-11+expdistance-202.5×Value 0.04356 ais ais[0](0) [0, 2] S/cm2 
gIhbar_Ih (− 0.8696 + 2.087 × 0.00001 somatic, soma(0.5) — S/cm2 
 exp ((distance × 0.0031)) × Value  apical, basal    
gNaTgbar_NaTg exp (distance × (− 0.00352)) × Value 0.06642 apical soma(0.5) [0, 0.1] S/cm2 

Note: Only parameters with optimization bounds are fitted.

5.3  Ground-Truth Model

To build the ground-truth model, we used the morphology of a previously reconstructed layer-5 pyramidal neuron (L5PC; Hay et al., 2011). The axon was removed except for the first 35 micrometers (axon initial segment or AIS). A 1-mm-long myelinated cylindrical axon with a 0.2 μm diameter was attached to the AIS. The model was composed of four sections, which contained specific voltage-gated conductances (apical dendrites, basal dendrites, soma, AIS) and one section (myelinated axon) that contained only passive channels. The voltage-gated mechanisms were derived from the L5PC model of Markram et al. (2015). In order to improve the simulation of action potential waveforms in the AIS, we replaced the axonal sodium channels present in Markram et al. (2015) by two different sodium channels (Nav1.2 and Nav1.6) that have been shown to be present in the AIS of L5PC (Hu et al., 2009). As shown experimentally (Hu et al., 2009), Nav1.2 density gradually decreased, while Nav1.6 density gradually increased along the AIS. Moreover, we also added a gradual increase of Kv1 channels along the AIS (Kole et al., 2007; Shu et al., 2007). The parameters of the model were optimized in the different sections using BluePyOpt (Van Geit et al., 2016) in order to reproduce 58 electrical features of L5PC. The somatic features (e.g., action potential waveforms, firing rates) were taken from previously recorded L5PC (Markram et al., 2015). The dendritic features (action potential backpropagation and calcium entry) were adapted from Stuart and Sakmann (1994), Helmchen et al. (1996), and Schiller et al. (1995). The AIS features (action potential waveform) were taken from Ritzau-Jost et al. (2021).

5.4  Feature Selection and Extraction

Relevant electric features (e-features) were extracted from the somatic patch-clamp recordings of the BluePyEfe Python package (Blue Brain Project, 2022a), which itself relies on the eFEL Python package (Blue Brain Project, 2022b) for e-feature computation.

To perform e-feature extraction, the e-features were first computed on all voltage recordings independently. Then the rheobase of each cell was computed as the lowest current amplitude, triggering at least one spike in the majority of recordings of the cell. Each recording was then associated with an amplitude, expressed as a percentage of the rheobase of the cell to which it belongs. Finally, the mean and standard deviation of the e-features were computed across the recordings of the same stimulus and amplitude. To perform this operation, a tolerance of 10% was added to the amplitudes. This operation can be better understood as a running average of width 20% of the e-feature value along the amplitude. The end result was a set of e-features where combos e-feature name + stimulus + amplitude were associated with a mean and standard deviation. Note that for the ground-truth model described above, due to its deterministic nature and therefore lack of variance, the e-feature standard deviation was set to 5% of the e-feature mean.

Following is a list of the e-features used as targets for optimizing and validating the models. These e-features were chosen to fully describe the firing pattern and shape of the APs seen in the experimental data. Amplitudes were expressed as percentages of the rheobase of the cell. Refer to the documentation of eFEL for a description of each e-feature.

The e-features used for optimization for different protocols are (total of 79 features):

  • IDrest (amplitudes: 150%, 250%, 300%): mean_frequency, burst_number, adaptation_index2, ISI_CV, ISI_log_slope, inv_time_to_first_spike, inv_first_ISI, inv_second_ISI, inv_third_ISI, inv_fourth_ISI, inv_fifth_ISI, AP_amplitude, AHP_depth, AHP_time_from_peak, voltage_base, Spikecount (before step), Spikecount (after step)

  • IV (amplitudes: −100%, −20%): Spikecount, voltage_base, voltage_deflection, voltage_deflection_begin, steady_state_voltage_stimend, ohmic_input_resistance_vb_ssse, decay_time_constant_after_stim, sag_amplitude (only for −100%), sag_ratio1 (only for −100%), sag_ratio2 (only for −100%)

  • APWaveform (amplitude: 290%): AP_amplitude, AP1_amp, AP2_amp, AP_duration_half_width, AP_begin_width, AP_begin_voltage, AHP_depth, AHP_time_from_peak

The e-features used for validations for different protocols are (total of 110 features):

  • firepattern (amplitudes: 120%, 200%): mean_frequency, burst_number, adaptation_index2, ISI_CV, ISI_log_slope, inv_time_to_first_spike, inv_first_ISI, inv_second_ISI, inv_third_ISI, inv_fourth_ISI, inv_fifth_ISI, AP_amplitude, AHP_depth, AHP_time_from_peak

  • HyperDepol (amplitudes: −160%, −120%, −80%, −40%): Spikecount, burst_number, AP_amplitude, ISI_values for the depolarized phase of the stimulus and sag_amplitude, sag_ratio1, sag_ratio2 for the hyperpolarized phase

  • sAHP (amplitudes: 150%, 200%, 250%, 300%): Spikecount, AP_amplitude, ISI_values, AHP_depth, AHP_depth_abs, AHP_time_from_peak, steady_state_voltage_stimend

  • Poscheops (amplitude: 300%): Spikecount, mean_frequency, burst_number, adaptation_index2, ISI_CV, ISI_log_slope, inv_time_to_first_spike, inv_first_ISI, inv_second_ISI, inv_third_ISI, inv_fourth_ISI, inv_fifth_ISI (computed during the three pyramids independently)

5.5  Optimization Algorithm

The optimization was performed using the BluePyOpt python package (Van Geit et al., 2016) whose optimization module relies on the DEAP Python package (Fortin et al., 2012). We extended the package and implemented a hybrid CMA optimization strategy (Hansen, 2016; Damart et al., 2020) tasked to both:

  • Minimize the sum of the scores defined as si=|μexpi-fi|/σexpi, where μexpi and σexpi were the experimental mean and standard deviation for e-feature i, respectively, and fi was the feature value computed on the model to evaluate.

  • Maximize the hypervolume of the Pareto front formed by the current population of models (Bader & Zitzler, 2011).

At each generation, all models in the population of size λ = 20 were ranked for both criteria, and a mixed rank was obtained following the formula:
rankmixedj=(whv×rankhvj)+((1-whv)×rankscoresj)
The weight assigned to the hypervolume ranking, whv, was set to 0.4 in our study. Following this ranking, the μ = λ/2 first individuals were selected to update the CMA kernel for the next generation. Note that during the computation of the scores, the scores computed for the extracellular e-features by cosine similarity (all and sections strategies) were weighed by an empirically obtained factor of 2.5, compared to the intracellular ones, in order to balance the values of intracellular and extracellular scores.

For each model, we ran 10 optimizations of 600 generations, each for a different starting seed of the random number generator.

5.6  Experimental Data Acquisition

5.6.1  Ethics Statement

All animal experimental protocols were approved by the Basel-Stadt veterinary office (cantonal no. 2358) according to Swiss federal laws (national no. 30692) on animal welfare and were carried out in accordance with the approved guidelines.

5.6.2  HD-MEA System

We used an in-house-developed high-density microelectrode array (HD-MEA) system with 26,400 electrodes covering a sensing area of around 2 × 4 mm2 (Frey et al., 2010; Müller et al., 2015). The center-to-center electrode distance was 17.5  μm, and the electrodes were coated with Pt-black to lower impedance and improve the signal-to-noise ratio. The HD-MEA system can be configured to record from up to 1024 electrodes simultaneously at 20 kHz. Electrode configurations were chosen according to the position of the target neuron.

5.6.3  Embryonic Rat Cortical Cultures

All experimental protocols involving animals were approved by the Basel-Stadt veterinary office according to Swiss federal laws on animal welfare and were carried out in accordance with the approved guidelines. We used rat primary neurons, obtained from dissociated cortices of Wistar rats at embryonic day 18, following the protocol described in Ronchi et al. (2019).

Before cell plating, HD-MEA chips were sterilized with 70% ethanol for 30 minutes. After the ethanol was removed, the chips were rinsed three times with sterile tissue-culture-grade water. The HD-MEA chips were coated with a layer of 0.05% polyethylenimine (Sigma-Aldrich, Buchs, Switzerland) in borate buffer to make the surface more hydrophilic. Next, we added a layer of laminin (Sigma-Aldrich, Buchs, Switzerland, 0.02 mg/mL) in neurobasal medium (Thermo Fisher Scientific, Waltham, MA) on the HD-MEA and incubated for 30 minutes at 37°C to facilitate cell adhesion. We used trypsin with 0.25% EDTA (Gibco, Thermo Fisher Scientific), followed by trituration, to dissociate embryonic cortices of E-18 Wistar rats enzymatically and then seeded 15,000 to 20,000 cells in 7 μL on top of the electrode arrays. The plated chips were incubated for 30 minutes at 37°C before adding 2 mL of plating medium. The plating medium consisted of neurobasal, supplemented with 10% horse serum (HyClone, Thermo Fisher Scientific), 0.5 mM Glutamax (Invitrogen, Thermo Fisher Scientific), and 2% B-27 (Invitrogen, Thermo Fisher Scientific). After three days, 50% of the plating medium was replaced by the BrainPhys neuronal medium (Stem Cell Technologies) growth medium. The procedure was repeated twice a week. The chips were kept inside an incubator at 37°C and 5% CO2. All experiments were conducted between days in vitro (DIVs) 18 and 24.

5.6.4  Patch-Clamp Solutions

The intracellular patch-clamp solution consisted of potassium-gluconate (110 mM), phosphocreatine (10 mM), KCl (10 mM), hepes (10 mM), GTP (0.3 mM), ATP-Mg (4 mM) dissolved in nanopure water (chemicals purchased from Sigma-Aldrich, Buchs, Switzerland). The pH was adjusted to 7.2–7.3 by adding potassium hydroxide at 5 M. On each experimental day, an aliquot was thawed, Alexa Fluor 594 (50 μM; Sigma-Aldrich) was added, and the final solution was filtered using a Millex GV 0.22 μm filter.

The extracellular artificial cerebrospinal fluid (aCSF) solution contained NaCl (125 mM), KCl (2.5 mM), MgCl2·6H2O (1 mM), NaH2PO4 (1.25 mM), and CaCl2·2H2O (2 mM). On the day of the experiment, 2 L of 1× aCSF was prepared by combining 200 mL of a 10× stock solution with 1 L nanopure water, dissolving 9 g of glucose and 4.2 g of NaHCO3, and topping up to 2 L with nanopure water. This solution was bubbled with carbogen (5% CO2 / 95% O2) and heated to a temperature of around 34°C throughout the experiment.

The combination of intra- and extracellular solutions yielded a liquid junction potential (LJP) of approximately 14 mV that we used to correct the values of the recorded membrane potentials.

5.6.5  Simultaneous HD-MEA and Patch-Clamp Recording

The experiments were conducted with a custom patch-clamp rig with an upright microscope equipped with the HD-MEA recording unit. The patch-clamp system included a MultiClamp 700B amplifier (Axon Instruments, Victoria, Australia), Axon Digidata 1440A (Axon Instruments), PatchStar Micromanipulator (Scientifica, Uckfield, U.K.), and an Olympus BX61WI microscope (Olympus, Tokyo, Japan). The WinWCP software was used to control the patch-clamp system and acquire the intracellular recordings. Glass micropipettes (4–8 MΩ) were pulled using a P-1000 micropipette puller (Sutter Instruments, Novato, CA). Whole-cell current-clamp recordings were low-pass-filtered at 10 kHz and digitized at 20 kHz. After whole-cell current-clamp mode was established with a target cell, the holding current was manually adjusted to obtain a resting membrane potential of around −70 mV (note that with LJP correction, this corresponds to −84 mV). Brief current pulses to induce action potentials were then delivered to the neuron in order to assess that the extracellular signals reliably showed spikes. We discarded cells with extracellular spikes lower than 50 μV. The patch-clamp and HD-MEA systems were synchronized by sending TTL pulses from the patch-clamp output channel to the field-programmable-gate-array (FPGA) board controlling the HD-MEA device at the beginning of each sweep. Next, we roughly estimated the patched neuron’s rheobase current using increasing step pulses of 270 ms duration. The estimated rheobase current was used to automatically generate the eCode protocol files that were then injected into the neuron. The eCode protocols were repeated during four to six runs (each run lasted approximately 200 s).

5.6.6  Confocal Imaging

After the electrophysiology data acquisition, the chip was mounted under a Nikon NiE upright microscope, equipped with a Yokogawa W1 spinning disk scan head, an ORCA-Flash4.0 V2 digital CMOS camera (Hamamatsu Photonics, Tokyo, Japan), and a 60×/1.00 NA water-objectives (Nikon, Tokyo, Japan). Prior to imaging, we co-registered the location of the electrodes to the microscope stage position using the stage positions of three electrodes at the vertices of the sensing area. The confocal channel for imaging had an excitation laser of 561 nm and an emission filter of 609/54 nm. The patched neurons were then imaged using multiple tiles covering their entire morphology. For each tile, a z-stack with 0.4  μm z-step was acquired. The x-y resolution of the images was 112.5 nm. We acquired 6 × 4 tiles and 66 z images for cell 1, and 7 × 4 tiles and 61 z images for cell 2.

5.6.7  Fixation and Staining

After imaging, the culture was fixated by removing the remaining medium and adding 1 mL of a 4% paraformaldehyde (PFA, Thermo Fisher Scientific, MA) solution for 10 minutes. Subsequently, the chip was washed three times with PBS for 10 minutes with gentle agitation. Permeabilization of the membrane and blocking of unspecific binding sites were performed in one step by adding 0.5% Triton X-100 (Sigma-Aldrich, Buchs, Switzerland) to the blocking solution (PBS; Thermo Fisher Scientific) with 10% normal donkey serum (Jackson ImmunoResearch, West Grove, PA), 0.02% Na-Az (Sigma-Aldrich), 1% bovine serum albumin (Sigma-Aldrich) for 30 minutes. The cultures were then stained for AIS-specific antibodies (Ankyrin-G and Kv7.3) to identify the AIS tract for the morphology reconstruction. The primary antibodies were diluted in the reaction solution (PBS with 3% normal donkey serum, 0.02% Na-Az, 1% bovine serum albumin, 0.05% Triton X-100) and added to the culture. After 2 hours at room temperature on a shaker, the primary antibody solution was aspirated, and the culture was washed thrice with PBS. The secondary antibodies were diluted in fresh reaction solution and added to the culture. After 4 hours at room temperature on a shaker, the secondary antibody solution was removed, and the culture was washed thrice with PBS. Until imaging, the culture was kept in PBS at 4°C. The antibodies used in this study included rabbit-α-Kv7.3 (Alomone Labs, Jerusalem, Israel), guinea pig-α-AnkG (Synaptic Systems, catalog number 386 004), donkey-α-rabbit Alexa Fluor Plus 488 (Thermo Fisher Scientific), donkey-α-guinea pig Alexa Fluor 647 (Jackson ImmunoResearch).

5.7  Experimental Data Analysis

5.7.1  Electrophysiology

The combined patch-clamp and HD-MEA signals were analyzed using SpikeInterface (Buccino et al., 2020) and NEO (Garcia et al., 2014) to interface with extracellular and intracellular recordings, respectively.

The patch-clamp and extracellular recordings for each run were processed separately. The TTL pulses, sent by the patch-clamp system to the HD-MEA FPGA, were detected in both data streams and used to synchronize the intra- and extracellular traces (all patch-clamp protocols for each run were concatenated at this time). We removed extracellular channels that displayed saturation of the amplifiers (denoising) and filtered the signals with a third-order Butterworth filter with cutoff frequencies at 300 and 6000 Hz. Next, the patch-triggered average was used to compute extracellular templates: intracellular action potentials were detected from the patch signals (the PosCheops protocol was discarded because of the strong spike waveform modulation due to the high-frequency firing rate) and used to extract extracellular waveforms, which were further centered around the extracellular peaks on the channel with the largest signal amplitude to correct for possible mismatches between the intra- and extracellular peak times. To obtain a cleaner template, we discarded waveforms, whose amplitudes were over ±2 standard deviations away from the median amplitude on the channel featuring the largest signal amplitude. The template was then computed as the sample-wise median of the waveforms (we used median because it is less sensitive to noise outliers). Templates were extracted for each run, and the template, which yielded the largest number of channels after denoising, was used for downstream analyses.

5.7.2  Imaging and Morphology Reconstruction

We employed Huygens Professional (version 21.10, Scientific Volume Imaging, the Netherlands, http://svi.nl) to perform the image deconvolution and stitching. Specifically, we first deconvolved images using the CMLE algorithm with SNR:12 and 40 iterations. Subsequently, the deconvolved images were stitched together with an overlap of 10%, using the circular vignetting correction model.

We used the SNT plug-in in Fiji (Arshadi et al., 2021) to reconstruct the 3D morphology of the neurons. The soma position was reconstructed as a single point by merging all the branches originating from it to the same root. Neurites were tagged as “soma,” “dend,” and “axon.” The AIS tract was identified using the staining images and it was labeled as “apic” since the SWC format that we used to export the morphology did not formally support the axon initial segment tag. Subsequently we relabeled the AIS correctly at instantiation in BluePyOpt. The radii of the reconstructed paths were fitted using the Refine/Fit tool and exported to SWC format. The raw morphology was preprocessed with a custom Python function (https://github.com/alejoe91/multimodalfitting/blob/master/multimodalfitting/imaging_tools/correct_swc.py) that interpolated missing radii (below or equal to 0.1 μm) and smoothed the radii of each path with a 15-sample moving average. The imaging data were also used to estimate the distance between the center of the soma and the underlying HD-MEA plane.

5.7.3  Experimental Model Definition and Feature Selection

The models based on experimental data were built using a similar set of mechanisms and boundaries used for the ground-truth model. However, due to the early-development stage of the cultured neurons and the absence of clearly defined basal and apical dendrites, all dendrites were considered the same and expressed all the basal mechanisms (no apical dendrites were defined). The axon-bearing dendrite (ABD) was also treated as a basal dendrite. Furthermore, an axonal section was added to the model to represent the unmyelinated axon. In total, the experimental models had 42 free parameters to optimize. Tables 4 and 5 list all the parameters for the experimental cell models, including their values (for fixed parameters), boundaries (for free parameters), and distributions (see Table 5).

Table 4:

Experimental Model Parameters (Nondistribution Type) and Optimization Bounds.

ParameterValueSection listBoundsUnits
v_init −84  — mV 
celsius 34  — °C 
cm somatic, axonal, ais — μF/cm2 
cm basal — μF/cm2 
Ra 100 all — Ω-cm 
ek -90 all — mV 
ena 50 all — mV 
e_pas — all [−95, −60] mV 
g_pas — all [1e−05, 6e−05] S/cm2 
decay_CaDynamics_DC0 — somatic, axonal, ais [20, 300] ms 
gamma_CaDynamics_DC0 — all [5e−3, 5e−2] — 
gCa_HVAbar_Ca_HVA2 — somatic, axonal, ais [0, 0.001] S/cm2 
gCa_HVAbar_Ca_HVA2 — basal [0, 0.0001] S/cm2 
gCa_LVAstbar_Ca_LVAst — somatic, axonal, ais [0, 0.01] S/cm2 
gCa_LVAstbar_Ca_LVAst — basal [0, 0.001] S/cm2 
vshiftm_NaTg basal — mV 
vshifth_NaTg basal — mV 
gSKv3_1bar_SKv3_1 — basal [0, 3e−3] S/cm2 
vshiftm_NaTg 13 somatic — mV 
vshifth_NaTg 15 somatic — mV 
slopem_NaTg somatic — mV 
gNaTgbar_NaTg — somatic [0, 0.3] S/cm2 
gK_Pstbar_K_Pst — somatic [0, 0.2] S/cm2 
gK_Tstbar_K_Tst — somatic [0, 0.1] S/cm2 
gSKv3_1bar_SKv3_1 — somatic, axonal, ais [0, 1] S/cm2 
gSK_E2bar_SK_E2 — somatic, axonal, ais [0, 0.1] S/cm2 
gNap_Et2bar_Nap_Et2 — axonal, ais [0, 0.02] S/cm2 
gK_Pstbar_K_Pst — axonal, ais [0, 2] S/cm2 
gK_Tstbar_K_Tst — axonal, ais [0, 0.2] S/cm2 
gNa16bar_Na16 — axonal [0, 8] S/cm2 
gkbar_Kd — axonal [0, 2] S/cm2 
constant_decay — meta [−0.1, 0] — 
ParameterValueSection listBoundsUnits
v_init −84  — mV 
celsius 34  — °C 
cm somatic, axonal, ais — μF/cm2 
cm basal — μF/cm2 
Ra 100 all — Ω-cm 
ek -90 all — mV 
ena 50 all — mV 
e_pas — all [−95, −60] mV 
g_pas — all [1e−05, 6e−05] S/cm2 
decay_CaDynamics_DC0 — somatic, axonal, ais [20, 300] ms 
gamma_CaDynamics_DC0 — all [5e−3, 5e−2] — 
gCa_HVAbar_Ca_HVA2 — somatic, axonal, ais [0, 0.001] S/cm2 
gCa_HVAbar_Ca_HVA2 — basal [0, 0.0001] S/cm2 
gCa_LVAstbar_Ca_LVAst — somatic, axonal, ais [0, 0.01] S/cm2 
gCa_LVAstbar_Ca_LVAst — basal [0, 0.001] S/cm2 
vshiftm_NaTg basal — mV 
vshifth_NaTg basal — mV 
gSKv3_1bar_SKv3_1 — basal [0, 3e−3] S/cm2 
vshiftm_NaTg 13 somatic — mV 
vshifth_NaTg 15 somatic — mV 
slopem_NaTg somatic — mV 
gNaTgbar_NaTg — somatic [0, 0.3] S/cm2 
gK_Pstbar_K_Pst — somatic [0, 0.2] S/cm2 
gK_Tstbar_K_Tst — somatic [0, 0.1] S/cm2 
gSKv3_1bar_SKv3_1 — somatic, axonal, ais [0, 1] S/cm2 
gSK_E2bar_SK_E2 — somatic, axonal, ais [0, 0.1] S/cm2 
gNap_Et2bar_Nap_Et2 — axonal, ais [0, 0.02] S/cm2 
gK_Pstbar_K_Pst — axonal, ais [0, 2] S/cm2 
gK_Tstbar_K_Tst — axonal, ais [0, 0.2] S/cm2 
gNa16bar_Na16 — axonal [0, 8] S/cm2 
gkbar_Kd — axonal [0, 2] S/cm2 
constant_decay — meta [−0.1, 0] — 

Note: Only parameters with optimization bounds are fitted. The constant_decay parameter is a meta-parameter that controls the gNaTgbar_NaTg distribution in Table 5.

Table 5:

Experimental Model Parameters with Distributions and Optimization Bounds.

ParameterDistributionSection listRef.BoundsUnits
gNa12bar_Na12 11+expdistance-253×Value ais ais[0](0) [0, 8] S/cm2 
gNa16bar_Na16 1-11+expdistance-202.5×Value ais ais[0](0) [0, 8] S/cm2 
gkbar_Kd 1-11+expdistance-202.5×Value ais ais[0](0) [0, 2] S/cm2 
gIhbar_Ih (− 0.8696 + 2.087 × somatic, soma(0.5) [0, 2e−4] S/cm2 
 exp ((distance × 0.0031)) × Value basal    
gNaTgbar_NaTg exp(distance×constant_decay)×Value basal soma(0.5) [0, 0.1] S/cm2 
ParameterDistributionSection listRef.BoundsUnits
gNa12bar_Na12 11+expdistance-253×Value ais ais[0](0) [0, 8] S/cm2 
gNa16bar_Na16 1-11+expdistance-202.5×Value ais ais[0](0) [0, 8] S/cm2 
gkbar_Kd 1-11+expdistance-202.5×Value ais ais[0](0) [0, 2] S/cm2 
gIhbar_Ih (− 0.8696 + 2.087 × somatic, soma(0.5) [0, 2e−4] S/cm2 
 exp ((distance × 0.0031)) × Value basal    
gNaTgbar_NaTg exp(distance×constant_decay)×Value basal soma(0.5) [0, 0.1] S/cm2 

Note: Ref. is the reference point from which NEURON distance is measured for the respective formula.

Figure 9:

Visualization of channels selected for the single strategy (top row) and sections strategy (bottom row) for the ground truth (A–D), cell 1 (B–E), and cell 2 (C–F) models.

Figure 9:

Visualization of channels selected for the single strategy (top row) and sections strategy (bottom row) for the ground truth (A–D), cell 1 (B–E), and cell 2 (C–F) models.

Close modal
Figure 10:

Analysis of cell 2 EAP mismatch. To further characterize the EAP mismatch, we virtually patched the cell 2 all model at the soma, AIS, and at basal dendrite. The left panel shows the location of the virtual pipettes, the cell morphology, and the overlaid EAP signals. The associated transmembrane currents Im, at the top right (normalized), show different characteristics: soma (red), initial small positive phase, followed by a strong negative signal and a repolarization positive phase; basal dendrite (brown), strong initial positive upstroke and a long negative phase afterward; AIS (purple), early strong negative phase going back to zero without a positive phase. Looking at the zoomed-in EAP signals, we believe that the mismatch can be explained by an underestimation of the somatic contribution and an overestimation of the AIS contribution. The upstroke of the experimental EAP close to the basal dendrite (black trace in the green box on the right) seems to coincide with the somatic negative peak, rather than the AIS one. The initial negative phase instead is time-locked to the AIS negative peak. A smaller AIS contribution paired with a stronger somatic current could therefore explain the experimental trace: the positive upstroke in the green box would be the result of the larger return current associated with the somatic inward current. To further motivate this hypothesis, the blue box shows the EAP on an electrode close to the AIS. The experimental EAP, different from the modeled ones, shows a stronger contribution of the soma (late negative peak), despite being closer to the AIS compartment.

Figure 10:

Analysis of cell 2 EAP mismatch. To further characterize the EAP mismatch, we virtually patched the cell 2 all model at the soma, AIS, and at basal dendrite. The left panel shows the location of the virtual pipettes, the cell morphology, and the overlaid EAP signals. The associated transmembrane currents Im, at the top right (normalized), show different characteristics: soma (red), initial small positive phase, followed by a strong negative signal and a repolarization positive phase; basal dendrite (brown), strong initial positive upstroke and a long negative phase afterward; AIS (purple), early strong negative phase going back to zero without a positive phase. Looking at the zoomed-in EAP signals, we believe that the mismatch can be explained by an underestimation of the somatic contribution and an overestimation of the AIS contribution. The upstroke of the experimental EAP close to the basal dendrite (black trace in the green box on the right) seems to coincide with the somatic negative peak, rather than the AIS one. The initial negative phase instead is time-locked to the AIS negative peak. A smaller AIS contribution paired with a stronger somatic current could therefore explain the experimental trace: the positive upstroke in the green box would be the result of the larger return current associated with the somatic inward current. To further motivate this hypothesis, the blue box shows the EAP on an electrode close to the AIS. The experimental EAP, different from the modeled ones, shows a stronger contribution of the soma (late negative peak), despite being closer to the AIS compartment.

Close modal

All the data and code used in this article are openly available. The contributions to the BluePyOpt (https://github.com/BlueBrain/BluePyOpt) package are summarized in the pull request (https://github.com/BlueBrain/BluePyOpt/pull/385) and included in the 1.13.88 version release. The script, notebooks, and preprocessed data used for optimizations (including cell morphologies, feature and protocol files, extracellular templates and probe files) are available on this GitHub repo: https://github.com/alejoe91/multimodalfitting. The reconstructed cell morphologies are also available through the NeuroMorpho.org (Ascoli et al., 2007) portal: https://neuromorpho.org/neuron_info.jsp?neuron_name=Buccino_Hierlemann_Cell1 and https://neuromorpho.org/neuron_info.jsp?neuron_name=Buccino_Hierlemann_Cell2. Finally, the raw electrophysiology data (simultaneous patch-clamp and extracellular HD-MEA recordings) are available in Neurodata Without Borders (NWB) format (Teeters et al., 2015; Rübel et al., 2021) on the DANDI archive: https://dandiarchive.org/dandiset/000294.

This study was supported by the ETH Zurich Postdoctoral Fellowship 19-2 FEL-17 (A.P.B.), the ERC Advanced Grant 694829 “neuroXscales” (J.B., X.X., T.G., V.E., A.H.), the China Scholarship Council, by funding to the Blue Brain Project, a research center of the École polytechnique fédérale de Lausanne, from the Swiss government’s ETH Board of the Swiss Federal Institutes of Technology (T.D., D.M., M.Z., A.J., H.M., W.V.G.), and by the European Union’s Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreement 945539 (Human Brain Project SGA3) (T.D., A.J.).

Abbott
,
J.
,
Ye
,
T.
,
Krenek
,
K.
,
Gertner
,
R. S.
,
Ban
,
S.
,
Kim
,
Y.
, . . .
Ham
,
D
. (
2020
).
A nanoelectrode array for obtaining intracellular recordings from thousands of connected neurons
.
Nature Biomedical Engineering
,
4
(
2
),
232
241
.
Abdelfattah
,
A. S.
,
Kawashima
,
T.
,
Singh
,
A.
,
Novak
,
O.
,
Liu
,
H.
,
Shuai
,
Y.
, . . .
Schreiter
,
E. R
. (
2019
).
Bright and photostable chemigenetic indicators for extended in vivo voltage imaging
.
Science
,
365
(
6454
),
699
704
.
Agudelo-Toro
,
A.
, &
Neef
,
A
. (
2013
).
Computationally efficient simulation of electrical activity at cell membranes interacting with self-generated and externally imposed electric fields
.
Journal of Neural Engineering
,
10
(
2
),
026019
.
Allen
,
B. D.
,
Moore-Kochlacs
,
C.
,
Bernstein
,
J. G.
,
Kinney
,
J.
,
Scholvin
,
J.
,
Seoane
,
L.
, . . .
Boyden
,
E. S.
(
2018
).
Automated in vivo patch clamp evaluation of extracellular multielectrode array spike recording capability
.
Journal of Neurophysiology
,
120
(
5
).
Almog
,
M.
, &
Korngreen
,
A
. (
2014
).
A quantitative description of dendritic conductances and its application to dendritic excitation in layer 5 pyramidal neurons
.
Journal of Neuroscience
,
34
(
1
),
182
196
.
Almog
,
M.
, &
Korngreen
,
A.
(
2016
).
Is realistic neuronal modeling realistic?
Journal of Neurophysiology
,
116
(
5
),
2180
2209
.
Anastassiou
,
C. A.
,
Perin
,
R.
,
Buzsáki
,
G.
,
Markram
,
H.
, &
Koch
,
C
. (
2015
).
Cell type and activity-dependent extracellular correlates of intracellular spiking
.
Journal of Neurophysiology
,
114
(
1
),
608
623
.
Arshadi
,
C.
,
Günther
,
U.
,
Eddison
,
M.
,
Harrington
,
K. I.
, &
Ferreira
,
T. A
. (
2021
).
SNT: A unifying toolbox for quantification of neuronal anatomy
.
Nature Methods
,
18
(
4
),
374
377
.
Ascoli
,
G. A.
,
Donohue
,
D. E.
, &
Halavi
,
M
. (
2007
).
Neuromorpho.org: A central resource for neuronal morphologies
.
Journal of Neuroscience
,
27
(
35
),
9247
9251
.
Bader
,
J.
, &
Zitzler
,
E
. (
2011
).
HypE: An algorithm for fast hypervolume-based many-objective optimization
.
Evolutionary Computation
,
19
(
1
),
45
76
.
Bakkum
,
D. J.
,
Obien
,
M. E. J.
,
Radivojevic
,
M.
,
Jäckel
,
D.
,
Frey
,
U.
,
Takahashi
,
H.
, &
Hierlemann
,
A
. (
2019
).
The axon initial segment is the dominant contributor to the neuron’s extracellular electrical potential landscape
.
Advanced Biosystems
,
3
(
2
),
1800308
.
Blue Brain Project
(
2022a
).
Blue brain python e-feature extractionbluepyefe
. https://github.com/BlueBrain/BluePyEfe/tree/BPE2.
Blue Brain Project
(
2022b
).
Electrophys feature extraction library (EFEL)
. https://github.com/BlueBrain/eFEL.
Buccino
,
A. P.
, &
Einevoll
,
G. T
. (
2021
).
MEArec: A fast and customizable testbench simulator for ground-truth extracellular spiking activity
.
Neuroinformatics
,
19
(
1
),
185
204
.
Buccino
,
A. P.
,
Hurwitz
,
C. L.
,
Garcia
,
S.
,
Magland
,
J.
,
Siegle
,
J. H.
,
Hurwitz
,
R.
, &
Hennig
,
M. H
. (
2020
).
Spikeinterface, a unified framework for spike sorting
.
eLife
,
9
,
e61834
.
Buccino
,
A. P.
,
Kuchta
,
M.
,
Jæger
,
K. H.
,
Ness
,
T. V.
,
Berthet
,
P.
,
Mardal
,
K. A.
, . . .
Tveito
,
A.
(
2019
).
How does the presence of neural probes affect extracellular potentials?
Journal of Neural Engineering
,
16
(
2
).
Buccino
,
A. P.
,
Kuchta
,
M.
,
Schreiner
,
J.
, &
Mardal
,
K.-A
. (
2021
).
Improving neural simulations with the EMI model
. In
A.
Tveito
,
K.-A.
Mardal
, &
M. E.
Rognes
(Eds.),
Modeling excitable tissue
(pp.
87
98
).
Springer
.
Buzsáki
,
G
. (
2004
).
Large-scale recording of neuronal ensembles
.
Nature Neuroscience
,
7
(
5
),
446
451
.
Carnevale
,
N. T.
, &
Hines
,
M. L.
(
2006
).
The NEURON book
.
Cambridge University Press
.
Damart
,
T.
,
Van Geit
,
W.
, &
Markram
,
H
. (
2020
).
Data driven building of realistic neuron model using IBEA and CMA evolution strategies
. In
Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion
(pp.
35
36
).
Debanne
,
D.
,
Campanac
,
E.
,
Bialowas
,
A.
,
Carlier
,
E.
, &
Alcaraz
,
G
. (
2011
).
Axon physiology
.
Physiological Reviews
,
91
(
2
),
555
602
.
Druckmann
,
S.
,
Banitt
,
Y.
,
Gidon
,
A. A.
,
Schürmann
,
F.
,
Markram
,
H.
, &
Segev
,
I.
(
2007
).
A novel multiple objective optimization framework for constraining conductance-based neuron models by experimental data
.
Frontiers in Neuroscience
,
1
(
1
).
Druckmann
,
S.
,
Berger
,
T. K.
,
Hill
,
S.
,
Schürmann
,
F.
,
Markram
,
H.
, &
Segev
,
I
. (
2008
).
Evaluating automated parameter constraining procedures of neuron models by experimental and surrogate data
.
Biological Cybernetics
,
99
,
371
379
.
Egert
,
U.
,
Schlosshauer
,
B.
,
Fennrich
,
S.
,
Nisch
,
W.
,
Fejtl
,
M.
,
Knott
,
T.
, . . .
Hämmerle
,
H
. (
1998
).
A novel organotypic long-term culture of the rat hippocampus on substrate-integrated multielectrode arrays
.
Brain Research Protocols
,
2
(
4
),
229
242
.
Einevoll
,
G. T.
,
Kayser
,
C.
,
Logothetis
,
N. K.
, &
Panzeri
,
S
. (
2013
).
Modelling and analysis of local field potentials for studying the function of cortical circuits
.
Nature Reviews Neuroscience
,
14
(
11
),
770
785
.
Fékété
,
A.
,
Ankri
,
N.
,
Brette
,
R.
, &
Debanne
,
D.
(
2021
).
Neural excitability increases with axonal resistance between soma and axon initial segment
.
Proceedings of the National Academy of Sciences
,
118
(
33
).
Fortin
,
F.-A.
,
De Rainville
,
F.-M.
,
Gardner
,
M.-A.
,
Parizeau
,
M.
, &
Gagné
,
C
. (
2012
).
DEAP: Evolutionary algorithms made easy
.
Journal of Machine Learning Research
,
13
,
2171
2175
.
Frey
,
U.
,
Egert
,
U.
,
Heer
,
F.
,
Hafizovic
,
S.
, &
Hierlemann
,
A.
(
2009
).
Microelectronic system for high-resolution mapping of extracellular electric fields applied to brain slices
.
Biosensors and Bioelectronics
,
24
(
7
),
2191
2198
.
Frey
,
U.
,
Sedivy
,
J.
,
Heer
,
F.
,
Pedron
,
R.
,
Ballini
,
M.
,
Mueller
,
J.
, . . .
Hafizovic
,
S
. (
2010
).
Switch-matrix-based high-density microelectrode array in CMOS technology
.
IEEE Journal of Solid-State Circuits
,
45
(
2
),
467
482
.
Garcia
,
S.
,
Guarino
,
D.
,
Jaillet
,
F.
,
Jennings
,
T. R.
,
Pröpper
,
R.
,
Rautenberg
,
P. L.
, . . .
Davison
,
Y. P
. (
2014
).
Neo: An object model for handling electrophysiology data in multiple formats
.
Frontiers in Neuroinformatics
,
8
,
10
.
Gidon
,
A.
,
Zolnik
,
T. A.
,
Fidzinski
,
P.
,
Bolduan
,
F.
,
Papoutsi
,
A.
,
Poirazi
,
P.
, . . .
Larkum
,
M. E
. (
2020
).
Dendritic action potentials and computation in human layer 2/3 cortical neurons
.
Science
,
367
(
6473
),
83
87
.
Goaillard
,
J.-M.
,
Moubarak
,
E.
,
Tapia
,
M.
, &
Tell
,
F.
(
2020
).
Diversity of axonal and dendritic contributions to neuronal output
.
Frontiers in Cellular Neuroscience
(
January
), 570.
Goethals
,
S.
, &
Brette
,
R
. (
2020
).
Theoretical relation between axon initial segment geometry and excitability
.
eLife
,
9
,
e53432
.
Gold
,
C.
,
Henze
,
D. A.
, &
Koch
,
C
. (
2007
).
Using extracellular action potential recordings to constrain compartmental models
.
Journal of Computational Neuroscience
,
23
(
1
),
39
58
.
Gonçalves
,
P. J.
,
Lueckmann
,
J.-M.
,
Deistler
,
M.
,
Nonnenmacher
,
M.
,
Öcal
,
K.
,
Bassetto
,
G.
, . . .
Macke
,
Jakob
. (
2019
).
Training deep neural density estimators to identify mechanistic models of neural dynamics
.
bioRxiv
.
Gong
,
W.
,
Senčar
,
J.
,
Bakkum
,
D. J.
,
Jäckel
,
D.
,
Obien
,
M. E. J.
,
Radivojevic
,
M.
, &
Hierlemann
,
A. R
. (
2016
).
Multiple single-unit long-term tracking on organotypic hippocampal slices using high-density microelectrode arrays
.
Frontiers in Neuroscience
,
10
,
537
.
Goto
,
T.
,
Hatanaka
,
R.
,
Ogawa
,
T.
,
Sumiyoshi
,
A.
,
Riera
,
J.
, &
Kawashima
,
R
. (
2010
).
An evaluation of the conductivity profile in the somatosensory barrel cortex of Wistar rats
.
Journal of Neurophysiology
,
104
(
6
),
3388
3412
.
Gouwens
,
N. W.
,
Berg
,
J.
,
Feng
,
D.
,
Sorenson
,
S. A.
,
Zeng
,
H.
,
Hawrylycz
, . . .
Arkhipov
,
A.
(
2018
).
Systematic generation of biophysically detailed models for diverse cortical neuron types
.
Nature Communications
,
9
(
1
),
710
.
Hagen
,
E.
,
Næss
,
S.
,
Ness
,
T. V.
, &
Einevoll
,
G. T.
(
2018
).
Multimodal modeling of neural network activity: Computing LFP, ECOG, EEG, and MEG signals with LFPy2.0
Frontiers in Neuroinformatics
,
12
.
Hallermann
,
S.
,
De Kock
,
C. P.
,
Stuart
,
G. J.
, &
Kole
,
M. H
. (
2012
).
State and location dependence of action potential metabolic cost in cortical pyramidal neurons
.
Nature Neuroscience
,
15
(
7
),
1007
1014
.
Hansen
,
N.
(
2016
).
The CMA evolution strategy: A tutorial
.
arXiv
.
Hay
,
E.
,
Hill
,
S.
,
Schürmann
,
F.
,
Markram
,
H.
, &
Segev
,
I
. (
2011
).
Models of neocortical layer 5b pyramidal cells capturing a wide range of dendritic and perisomatic active properties
.
PLOS Computational Biology
,
7
(
7
),
e1002107
.
Helmchen
,
F.
,
Imoto
,
K.
, &
Sakmann
,
B
. (
1996
).
Ca2+ buffering and action potential-evoked Ca2+ signaling in dendrites of pyramidal neurons
.
Biophysical Journal
,
70
(
2
),
1069
1081
.
Hines
,
M.
,
Davison
,
A. P.
, &
Muller
,
E
. (
2009
).
Neuron and Python
.
Frontiers in Neuroinformatics
,
3
,
1
.
Hu
,
W.
,
Tian
,
C.
,
Li
,
T.
,
Yang
,
M.
,
Hou
,
H.
, &
Shu
,
Y
. (
2009
).
Distinct contributions of Na(v)1.6 and Na(1.2) in action potential initiation and backpropagation
.
Nature Neuroscience
,
12
(
8
),
996
1002
.
Huang
,
C. Y.-M.
, &
Rasband
,
M. N
. (
2018
).
Axon initial segments: structure, function, and disease
.
Annals of the New York Academy of Sciences
,
1420
(
1
),
46
61
.
Iavarone
,
E.
,
Yi
,
J.
,
Shi
,
Y.
,
Zandt
,
B.-J.
, & O’Reilly C.,
Van Geit
,
W.
, . . .
Hill
,
S. L
. (
2019
).
Experimentally-constrained biophysical models of tonic and burst firing modes in thalamocortical neurons
.
PLOS Computational Biology
,
15
(
5
),
e1006753
.
Jedlicka
,
P.
,
Bird
,
A. D.
, &
Cuntz
,
H
. (
2022
).
Pareto optimality, economy–effectiveness trade-offs and ion channel degeneracy: improving population modelling for single neurons
.
Open Biology
,
12
(
7
),
220073
.
Jia
,
X.
,
Siegle
,
J. H.
,
Bennett
,
C.
,
Gale
,
S. D.
,
Denman
,
D. J.
,
Koch
,
C.
, &
Olsen
,
S. R
. (
2019
).
High-density extracellular probes reveal dendritic backpropagation and facilitate neuron classification
.
Journal of Neurophysiology
,
121
(
5
),
1831
1847
.
Jun
,
J. J.
,
Steinmetz
,
N. A.
,
Siegle
,
J. H.
,
Denman
,
D. J.
,
Bauza
,
M.
,
Barbarits
,
B.
, . . .
Harris
,
T. D
. (
2017
).
Fully integrated silicon probes for high-density recording of neural activity
.
Nature
,
551
(
7679
),
232
.
Koch
,
C.
, &
Segev
,
I
. (
2000
).
The role of single neurons in information processing
.
Nature Neuroscience
,
3
(
11
),
1171
1177
.
Kole
,
M. H.
,
Letzkus
,
J. J.
, &
Stuart
,
G. J
. (
2007
).
Axon initial segment Kv1 channels control axonal action potential waveform and synaptic efficacy
.
Neuron
,
55
(
4
),
633
647
.
Larkum
,
M. E.
,
Zhu
,
J. J.
, &
Sakmann
,
B
. (
1999
).
A new cellular mechanism for coupling inputs arriving at different cortical layers
.
Nature
,
398
(
6725
),
338
341
.
Larkum
,
M. E.
,
Zhu
,
J. J.
, &
Sakmann
,
B
. (
2001
).
Dendritic mechanisms underlying the coupling of the dendritic with the axonal action potential initiation zone of adult rat layer 5 pyramidal neurons
.
Journal of Physiology
,
533
(
2
),
447
466
.
Lindén
,
H.
,
Hagen
,
E.
,
Leski
,
S.
,
Norheim
,
E. S.
,
Pettersen
,
K. H.
, &
Einevoll
,
G. T
. (
2014
).
LFPy: A tool for biophysical simulation of extracellular potentials generated by detailed model neurons
.
Frontiers in Neuroinformatics
,
7
,
41
.
Markram
,
H.
,
Muller
,
E.
,
Ramaswamy
,
S.
,
Reinmann
,
M. W.
,
Abdellah
,
M.
,
Sanchez
,
C. A.
,
Ailamaki
,
A.
, . . .
Schürmann
,
F
. (
2015
).
Reconstruction and simulation of neocortical microcircuitry
.
Cell
,
163
(
2
),
456
492
.
Müller
,
J.
,
Ballini
,
M.
,
Livi
,
P.
,
Chen
,
Y.
,
Radivojevic
,
M.
,
Shadmani
,
A.
, . . .
Hierlemann
,
A
. (
2015
).
High-resolution CMOS MEA platform to study neurons at subcellular, cellular, and network levels
.
Lab on a Chip
,
15
(
13
),
2767
2780
.
Ness
,
T. V.
,
Chintaluri
,
C.
,
Potworowski
,
J.
,
Łȩski
,
S.
,
Gła̧bska
,
H.
,
Wójcik
,
D. K.
, &
Einevoll
,
G. T
. (
2015
).
Modelling and analysis of electrical potentials recorded in microelectrode arrays (MEAs)
.
Neuroinformatics
,
13
(
4
),
403
426
.
Nunez
,
P. L.
, &
Srinivasan
,
R.
(
2006
).
Electric fields of the brain: The neurophysics of EEG
.
Oxford University Press
.
Obien
,
M. E. J.
,
Hierlemann
,
A.
, &
Frey
,
U.
(
2019
).
Accurate signal-source localization in brain slices by means of high-density microelectrode arrays
.
Scientific Reports
,
9
(
1
),
1
19
.
Peng
,
H.
,
Bria
,
A.
,
Zhou
,
Z.
,
Iannello
,
G.
, &
Long
,
F
. (
2014
).
Extensible visualization and analysis for multidimensional images using Vaa3D
.
Nature Protocols
,
9
(
1
),
193
208
.
Ramaswamy
,
S.
,
Courcol
,
J.
,
Abdellah
,
M.
,
Adaszewski
,
S. R.
,
Antille
,
N.
,
Arasever
,
S.
, . . .
Markram,
H.
(
2015
).
The neocortical microcircuit collaboration portal: A resource for rat somatosensory cortex
.
Frontiers in Neural Circuits
,
9
.
Ritzau-Jost
,
A.
,
Tsintsadze
,
T.
,
Krueger
,
M.
,
Ader
,
J.
,
Bechmann
,
I.
,
Eilers
,
J.
, . . .
Hallermann
,
S
. (
2021
).
Large, stable spikes exhibit differential broadening in excitatory and inhibitory neocortical boutons
.
Cell Reports
,
34
(
2
),
108612
.
Ronchi
,
S.
,
Fiscella
,
M.
,
Marchetti
,
C.
,
Viswam
,
V.
,
Müller
,
J.
,
Frey
,
U.
, &
Hierlemann
,
A
. (
2019
).
Single-cell electrical stimulation using CMOS-based high-density microelectrode arrays
.
Frontiers in Neuroscience
,
13
,
208
.
Rübel
,
O.
,
Tritt
,
A.
,
Ly
,
R.
,
Dichter
,
B. K.
,
Ghosh
,
S.
,
Niu
,
L.
, . . .
Bouchard
,
K. E.
(
2021
).
The neurodata without borders ecosystem for neurophysiological data science
.
bioRxiv
.
Schiller
,
J.
,
Helmchen
,
F.
, &
Sakmann
,
B
. (
1995
).
Spatial profile of dendritic calcium transients evoked by action potentials in rat neocortical pyramidal neurones
.
Journal of Physiology
,
487
(
3
),
583
600
.
Schneider
,
M.
,
Bird
,
A. D.
,
Gidon
,
A.
,
Triesch
,
J.
,
Jedlicka
,
P.
, &
Cuntz
,
H
. (
2023
).
Biological complexity facilitates tuning of the neuronal parameter space
.
PLOS Computational Biology
,
19
(
7
),
e1011212
.
Shu
,
Y.
,
Yu
,
Y.
,
Yang
,
J.
, &
McCormick
,
D. A
. (
2007
).
Selective control of cortical axonal spikes by a slowly inactivating K+ current
.
Proceedings of the National Academy of Sciences
,
104
(
27
),
11453
11458
.
Stuart
,
G. J.
, &
Sakmann
,
B
. (
1994
).
Active propagation of somatic action potentials into neocortical pyramidal cell dendrites
.
Nature
,
367
(
6458
),
69
72
.
Teeters
,
J. L.
,
Godfrey
,
K.
,
Young
,
R.
,
Dang
,
C.
,
Friedsam
,
C.
,
Wark
,
B.
, . . .
Sommer
,
F. T
. (
2015
).
Neurodata without borders: Creating a common data format for neurophysiology
.
Neuron
,
88
(
4
),
629
634
.
Teleńczuk
,
M.
,
Brette
,
R.
,
Destexhe
,
A.
, &
Teleńczuk
B.
(
2018
).
Contribution of the axon initial segment to action potentials recorded extracellularly
.
Eneuro
,
5
(
3
).
Thome
,
C.
,
Kelly
,
T.
,
Yanez
,
A.
,
Schultz
,
C.
,
Engelhardt
,
M.
,
Cambridge
,
S. B.
, . . .
Egorov
,
A. V
. (
2014
).
Axon-carrying dendrites convey privileged synaptic input in hippocampal neurons
.
Neuron
,
83
(
6
),
1418
1430
.
Tveito
,
A.
,
Jæger
,
K. H.
,
Lines
,
G. T.
,
Paszkowski
,
L.
,
Sundnes
,
J.
,
Edwards
,
A. G.
, . . .
Einevoll
,
G. T
. (
2017
).
An evaluation of the accuracy of classical models for computing the membrane potential and extracellular potential for neurons
.
Frontiers in Computational Neuroscience
,
11
,
27
.
Ujfalussy
,
B. B.
,
Makara
,
J. K.
,
Branco
,
T.
, &
Lengyel
,
M
. (
2015
).
Dendritic nonlinearities are tuned for efficient spike-based computations in cortical circuits
.
eLife
,
4
,
e10056
.
Van Geit
,
W.
,
Gevaert
,
M.
,
Chindemi
,
G.
,
Rössert
,
C.
,
Courcol
,
J.-D.
,
Muller
,
E. B.
, . . .
Markram
,
H
. (
2016
).
BluePyOpt: Leveraging open source software and cloud infrastructure to optimise model parameters in neuroscience
.
Frontiers in Neuroinformatics
,
10
,
17
.
Weaver
,
C. M.
, &
Wearne
,
S. L
. (
2006
).
The role of action potential shape and parameter constraints in optimization of compartment models
.
Neurocomputing
,
69
(
10–12
),
1053
1057
.

Author notes

Andreas Hierlemann and Werner Van Geit share senior authorship.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode