A Multimodal Fitting Approach to Construct Single-Neuron Models With Patch Clamp and High-Density Microelectrode Arrays

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.

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.

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).

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 Ca^{2 +}, 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.

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.

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).

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 (⁠$E∈RC×T$⁠). In order to extract relevant and lower-dimensional features from such templates, we compiled a set of 11 extracellular features f_ext that were defined and used for each recording electrode (or channel) (⁠$fext∈RC$⁠). 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 (N_selected), 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 s_i,

  $si=1-μexpi×fi∥μexpi∥∥fi∥,$

  (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 N_sections 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.

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, Na_v1.2, Na_v1.6, and K_v1 (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.

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⁻⁵, soma VS sections - p < 10⁻⁵, soma VS single - p < 10⁻⁵). 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.

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—K_v7.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 K_v7.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.

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.

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.

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.

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 I_holding 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 I_rheo, the current at which a single action potential was induced. The subsequent stimulus amplitudes are reported in percent values of the estimated I_rheo current.

• IDthresh: The purpose of this protocol was to finely screen the current levels to correctly reestimate I_rheo 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.

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 (Na_v1.2 and Na_v1.6) that have been shown to be present in the AIS of L5PC (Hu et al., 2009). As shown experimentally (Hu et al., 2009), Na_v1.2 density gradually decreased, while Na_v1.6 density gradually increased along the AIS. Moreover, we also added a gradual increase of K_v1 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).

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)

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 fⁱ 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).

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

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.

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 mm² (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.

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% CO₂. All experiments were conducted between days in vitro (DIVs) 18 and 24.

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), MgCl₂·6H₂O (1 mM), NaH₂PO₄ (1.25 mM), and CaCl₂·2H₂O (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 NaHCO₃, and topping up to 2 L with nanopure water. This solution was bubbled with carbogen (5% CO₂ / 95% O₂) 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.

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).

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.

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 K_v7.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-α-K_v7.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).

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.

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.

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).

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.).