## Abstract

Cortical information processing originates from the exchange of action potentials between many cell types. To capture the essence of these interactions, it is of critical importance to build mathematical models that reflect the characteristic features of spike generation in individual neurons. We propose a framework to automatically extract such features from current-clamp experiments, in particular the passive properties of a neuron (i.e., membrane time constant, reversal potential, and capacitance), the spike-triggered adaptation currents, as well as the dynamics of the action potential threshold. The stochastic model that results from our maximum likelihood approach accurately predicts the spike times, the subthreshold voltage, the firing patterns, and the type of frequency-current curve. Extracting the model parameters for three cortical cell types revealed that cell types show highly significant differences in the time course of the spike-triggered currents and moving threshold, that is, in their adaptation and refractory properties but not in their passive properties. In particular, GABAergic fast-spiking neurons mediate weak adaptation through spike-triggered currents only, whereas regular spiking excitatory neurons mediate adaptation with both moving threshold and spike-triggered currents. GABAergic nonfast-spiking neurons combine the two distinct adaptation mechanisms with reduced strength. Differences between cell types are large enough to enable automatic classification of neurons into three different classes. Parameter extraction is performed for individual neurons so that we find not only the mean parameter values for each neuron type but also the spread of parameters within a group of neurons, which will be useful for future large-scale computer simulations.

- dynamic threshold
- fitting method
- model characterization
- spike-frequency adaptation

cortical neurons exhibit a variety of different firing patterns in response to step currents (Connors and Gutnick 1990; Markram et al. 2004), which has led to intricate electrophysiological characterization schemes for three main neuronal cell-types [regular-spiking excitatory neurons (Exc), GABAergic fast-spiking neurons (FS), and GABAergic nonfast-spiking neurons (NFS)]. The electrophysiological characterization that reflects the biochemical composition of the cells (Toledo-Rodriguez et al. 2004) is often done manually by visual observation of the firing patterns or, more systematically, by automatic extraction of a few parameters such as the ratio of the first to the second interspike interval, the minimum voltage of the spike afterpotential (SAP), and the width of an action potential (AP). The definition of these parameters is arbitrary, and their relevance for neural information processing is questionable. The question arises whether a more principled characterization and classification of the cell types is possible based on the properties affecting the conversion of synaptic inputs into a spike.

In the community of computational neuroscience, it has been established over the last 20 yr that simplified spiking neuron models are capable of reproducing the variety of firing patterns that have been found in experimental preparations, including delayed spike onset, bursting, strong or weak adaptation, and refractoriness (Gerstner et al. 1996; Izhikevich 2004, 2007; Brette and Gerstner 2005; Naud et al. 2008; Touboul and Brette 2009; Mihalas and Niebur 2009). All of these models belong to the family of generalized integrate-and-fire (IF) models but vary in the way the standard leaky IF model (Lapicque 1907; Stein 1967) is generalized. Features to upgrade the simple IF model include spike aftercurrents (Baldissera et al. 1976; Gerstner et al. 1996; Benda and Herz 2003; Izhikevich 2007; Paninski et al. 2004; Brette and Gerstner 2005), dynamic threshold (Hill 1936; Fuortes and Mantegazzini 1962; Chacron et al. 2003; Jolivet et al. 2006; Badel et al. 2007), smooth spike initiation (Latham et al. 2000; Fourcaud-Trocme et al. 2003; Brette and Gerstner 2005), and linearized subthreshold currents (Richardson et al. 2003; Izhikevich 2007). Important questions are then: which of these features are needed for basic cortical computation? How many levels of complexity do we have to add to account for relevant features of cortical dynamics? Is the spike-frequency adaptation mediated by moving thresholds or spike-triggered currents?

To answer these questions, several groups used what one could call the “Turing test” for point-stimulated neurons (Keat et al. 2001; Paninski et al. 2005; Jolivet et al. 2006, 2008a; Kobayashi et al. 2009; Gerstner and Naud 2009): a somatic current is injected into a cortical neuron in vitro while its response is being recorded. The modeler then asks how the model response compares with the neuronal response on data that were not used for optimizing the parameters of the model. This is done both qualitatively [for instance by reproducing the firing patterns (Gerstner et al. 1996; Izhikevich 2004; Naud et al. 2008; Touboul and Brette 2009; Mihalas and Niebur 2009)] and quantitatively. The quantitative test consists of predicting the correct spike timing and subthreshold voltage of the real neuron. Obviously, a simplified neuron model is not expected to work across the whole range of stimuli that one can artificially design and apply to electrophysiological experiments. The optimal version of this test would therefore use for a neuron in vitro a stimulus that is similar (Poliakov et al. 1997; Monier et al. 2008) to the one a real neuron receives in an in vivo situation, so that the activity of the neuron is within a range that can be expected in vivo (Crochet and Petersen 2006; Poulet and Petersen 2008; Gentet et al. 2010). The optimal stimulus will appear as a noisy time-series reminiscent but not identical to the white noise previously used for characterizing the input/output relationship (Bryant and Segundo 1976; Marmarelis and Marmarelis 1978; de Ruyter van Stevenick and Bialek 1988). We say that the model is “good enough” if a neutral expert is not able to distinguish the activity of the model from that of the real neuron. We call this the Turing test of point-stimulated neuron models in analogy to the test of intelligence in computer programs suggested by Alan Turing over 60 yr ago (Turing 1950). A similar framework was also used to test the validity of neuron models in the visual pathway following light stimulation (Keat et al. 2001; Pillow et al. 2005, 2008; Carandini et al. 2007).

The aim of this study is threefold. Firstly, we present a systematic method for extracting the parameters that control the conversion of synaptic input into spike emission at the level of the soma. The method relies on the separation of the parameters affecting the subthreshold voltage and those affecting the firing threshold and its dynamics. Our method improves previously described methods for extracting spike-triggered currents (Paninski et al. 2004; Jolivet et al. 2006; Badel et al. 2007) and dynamic threshold (Azouz and Gray 2000; Chacron et al. 2003; Badel et al. 2007). We use this method to fit experimental data from GABAergic FS, GABAergic NFS, and GABAergic Exc. For each neuron type, we find the simplest model that reproduces the activity of neurons on data that were not used for parameter optimization. The models we extract reproduce the excitability type of GABAergic FS, GABAergic NFS, and Exc neurons (Connors and Gutnick 1990; Tateno et al. 2004; Markram et al. 2004). Moreover >80% of the spikes can be predicted on average while the mismatch in subthreshold voltage prediction is <2 mV.

Secondly, we ask which are the features essential for the neuron model to pass the Turing test for point-stimulated neurons. We consider subthreshold resonance and conductance- vs. current-based adaptation as well as current- vs. threshold-based adaptation. Our results reveal the importance of adaptation currents and of a moving threshold with time constants that can be as long as hundreds of milliseconds. We find that different types of neurons use moving thresholds and spike-triggered currents differently to mediate refractoriness and spike-frequency adaptation. Finally, we show that the parameters of the adaptation currents and dynamic threshold can be used for an automatic classification of the electrophysiological traces into three well-separated classes, whereas the passive parameters alone do not contain a sufficient amount of information to do so. We observe that the three neuron types have very different threshold dynamics and that efficient classification can be done using the parameters regulating the dynamics of the threshold.

We expect that applying the automatic fitting method on a larger database could allow an unsupervised classification of the extracted computational properties, which would open the possibility to detect the potential spread of parameters within a given class of neurons and therefore avoid a forced classification if in reality parameters are continuous. Finally, we expect that our method of parameter extraction once combined with similar results for synapses, dendrites, and connectivity patterns will enable us to build network models where neuronal parameters and firing patterns reproduce not only the mean “typical” firing of different neuron types but also the spread of firing characteristics within and between classes of neurons.

## MATERIALS AND METHODS

### In Vitro Two-Photon Microscopy and Whole Cell Recordings

All animal experiments were carried out under authorization from the Swiss Federal Veterinary Office. Brains of 17- to 22-day-old GAD67-green fluorescent protein (GFP) knock-mice (Tamamaki et al. 2003) were removed and quickly placed into an ice-cold modified artificial cerebrospinal fluid (aCSF; Bureau et al. 2006) containing the following (in mM): 110 choline chloride, 25 NaHCO_{3}, 25 d-glucose, 11.6 sodium ascorbate, 7 MgCl_{2}, 3.1 sodium pyruvate, 2.5 KCl, 1.25 NaH_{2}PO_{4}, and 0.5 CaCl_{2}. Three-hundred-micrometer oblique slices (parasagittal 35° away from vertical) were cut with a vibrating slicer (Leica VT1000S) and subsequently transferred into standard aCSF containing the following (in mM): 125 NaCl, 25 NaHCO_{3}, 25 d-glucose, 2.5 KCl, 1.25 NaH_{2}PO_{4}, 2 CaCl_{2}, and 1 MgCl_{2}, aerated with 95% O_{2}-5% CO_{2} at 35°C for 15 min. Afterwards slices were maintained at room temperature for at ≥30 min before use to allow for recovery from the slicing procedure. GFP-expressing neurons were visualized using a two-photon microscope (Prairie Technologies). Infrared two-photon excitation light of 880 nm was generated by a MaiTai laser (Spectra-Physics) and focused into the slice tissue through a 40 × 0.8 NA water immersion objective (Olympus). Detection of bandpass-filtered green fluorescence (525 ± 35 nm) was achieved using photomultiplier tubes (PMTs) above the objective and below the condenser. Infrared light was passed through a Dodt contrast element (Luigs & Neumann) and detected by an additional PMT to allow creation of a high contrast view of the brain tissue. Whole cell recordings were carried out at 33°C in standard aCSF. Borosilicate pipettes of 5–7 MΩ resistance were used. The patch-pipette intracellular solution contained the following (in mM): 135 K-gluconate, 4 KCl, 4 Mg-ATP, 10 Na_{2}-phosphocreatine, 0.3 Na-GTP, and 10 HEPES (pH 7.3, 280 mosmol/kgH_{2}O). A Multiclamp 700A amplifier (Molecular Devices) was used for whole cell recordings. Data were low-pass Bessel filtered at 10 kHz and digitized at 20 kHz with an ITC-18 acquisition interface (HEKA Electronics). The measured membrane potential was not corrected for the liquid junction potential.

The bridge balance feature of Multiclamp was not used during the recordings. Access resistance values were in the range of 10–15 MΩ. If there was an increase in access resistance or a drift in the resting membrane potential during the recording, traces were discarded. Moreover background synaptic activity was not blocked (e.g., with CNQX/AP5). Under our recording conditions, miniature and spontaneous excitatory postsynaptic potentials are present but neurons do not show spontaneous spiking.

Cortical layers 2/3 GABAergic inhibitory neurons were distinguished from Exc neurons by their expression of GFP (Gentet et al. 2010). GFP-expressing GABAergic neurons were further classified into FS and NFS neurons with respect to their AP kinetics upon somatic current pulse injection. An AP half-width <0.75 ms was used as a selection criterion for GABAergic FS neurons.

### Stimulation Protocol: Synaptic-Like Current

We construct input currents *I*_{syn}(*t*) as a weighted sum of six impinging spike trains constructed from independent inhomogeneous Poisson processes (see also Poliakov et al. 1997). Three spike-trains are convolved with a monoexponential filter with a time constant of 2 ms mimicking spike arrival at excitatory synapses and summed with weight *w*_{1,2,3}, while the three remaining ones are convolved with a 10-ms monoexponential filter for the inhibitory spike trains and summed with weight *w*_{4,5,6}. The unitary weights of the postsynaptic currents *w*_{1,2,3,4,5,6} = (20.00, 34.14, 5.86 −1.60, −3.89, −2.51) are then rescaled by a global factor *w*_{0} that was chosen for each cell individually to drive the neuron to a mean firing frequency between 2 and 15 Hz. In our data set, *w*_{0} is on average 0.89, 0.65, and 0.72 for the GABAergic FS neurons, the GABAergic NFS neurons, and the Exc neurons, respectively. Finally, all the Poisson processes shared the same time-dependent intensity (firing rate), which is a concatenation of blocks of 300- to 500-ms duration each with a constant intensity chosen randomly from a uniform distribution between 0 and 50 Hz. The duration of the blocks was drawn from a uniform distribution (Fig. 1, *A–C*). Note that this procedure produces a colored-noise current with time-dependent mean and SD. We call the time-dependent input constructed by the above procedure a synaptic-like current *I*_{syn}(*t*).

This synaptic-like current has the advantage that it produces voltage traces with interspike interval (ISI) distributions that cover a wide range of timescales, from a few milliseconds to >1 s (see Fig. 8 *G*, *H*, and *I*). This enables the extraction of threshold mechanisms that span different timescales, from milliseconds to seconds.

One minute of this synaptic-like current is injected in the cortical neurons repeatedly, interrupted by silent periods of 10 s. Each independent injection is called a “repetition.” Note that in each repetition the cell received exactly the same time course of synaptic-like current, which enables us to study peristimulus time histograms (PSTH) and reliability of neuronal spiking. To compare experimental results with our models, we used the first 30 s of all repetitions for fitting and reserved the last 30 s for testing the performance of the models. These two subsets of the data are called training and test sets, respectively.

### Other Stimulation Protocol

To assess the robustness of our fitting procedure on surrogate data, we also tested it on four different stimulation paradigms applied only to model neurons but not to real neurons. First, we use a gradually increasing ramp of current of 30 s. Second, we construct an input current as a series of 300-ms step current interleaving with 200 ms of silence, with increasing intensities. Third, we used as a stimulation current a white noise with 0 mean that lasts for 30 s. Finally, we construct a 30-s long input current made of colored noise, generated according to an Ornstein-Uhlenbeck process with correlation time constant of 4 ms. All the parameters of these currents are chosen so that the resulting input current produced an averaged firing frequency of ∼15 Hz.

### Performance Measurements

Since we are assessing the performance on the test set, the performance of the model will not increase by merely increasing the number of parameters because overfitting would occur on the training set and lower the performance on the test set. Two distinct criteria are used to evaluate the performance of our IF models: *1*) the precise spike time prediction, and *2*) the subthreshold voltage prediction.

#### Spike train similarity.

Neurons, as well as stochastic models, have some trial-to-trial variations due to intrinsic noise. In neurons, this intrinsic noise is mainly due to channel noise and spontaneous synaptic events (Faisal et al. 2008). To quantify the spike time prediction we used a method that corrects the bias due to the small number of available repetitions (Naud et al. 2011). We first use our optimal stochastic model to generate *R*_{m} = 1,000 spike-time predictions. The next step is to compute a quantity between zero and one that measures how well the set of *R*_{m} model spike trains matches the set of *R*_{r} experimental spike trains, where *R*_{r} is the number of repetitions available for the specific cell (lower index *r* for real neuron). To do so, we first count the total number of spikes of the model that fall within ±4 ms of a spike in the recorded spike train and average the count over all repetitions. We call the resulting raw measure of comparison between model (lower index *m*) and real neuron (lower index *r*): *n*_{mr}. Next, we count the average number of coincidences between distinct repetitions of the model spike trains and average the results across all available repetitions. We call this quantity *n*_{mm}. Similarly, *n*_{rr} is the analogous quantity evaluated over the available repetitions recorded in the real neuron. Finally, we combine these numbers into the measure:

#### Subthreshold voltage trace similarity.

The subthreshold voltage prediction is evaluated using the root mean square error on the voltage (RMSE) defined by
*V*_{ref}(*t|I*) is the recorded voltage given an input current *I, V _{predicted}*(

*t*|

*I*, θ, {t̂

_{ref}}) is the voltage predicted by a model with the set of parameters θ and the same input current. The index {t̂

_{ref}} indicates that for the voltage comparison we force spikes at exactly the same time as in the recorded voltage trace. In other words, we ask: how close is the voltage of the model to that of the data, given the input and the firing times in the recent past? The squared error is integrated over a subset Ψ

_{1}of the available data defined in appendix c.

*T*

_{1}is the total time for which the squared error is considered,

*T*

_{1}= ∫

_{Ψ1}

*dt*. Forcing the spikes at the observed spike times disentangles the subthreshold voltage prediction from the spike time prediction.

To estimate the RMSE solely due to intrinsic noise, we compute the RMSE between repetitions only on the subset of data Ψ_{4} (defined in appendix c), RMSE_{Ψ4}, which is restricted to voltage recordings sufficiently far away from the spikes to avoid the effect of adaptation currents due to spikes emitted at different times across repetitions (Fellous et al. 2004). RMSE_{Ψ4} is clearly an overestimation of the true subthreshold fluctuations because it still contains some long-lasting spike-triggered effects. Nevertheless, it can give some insights on how well our models have to predict the subthreshold voltage. To quantify the ability of the models to replicate subthreshold voltage fluctuations, we use the normalized “root mean square error ratio”:

## RESULTS

We stimulated a neuron with a time-dependent input current that mimics stochastic spike arrival (see Fig. 1 and materials and methods) and drives the neuron through episodes of high and low firing rates. We adapted parameters of generalized IF models so as to explain the data and predict spike times and subthreshold voltage for novel sitmuli not used during parameter optimization. results are organized as follows. First, we describe the model class and parameter extraction strategy. Second, we quantify the performance of the extracted neuron model and discuss differences between different cell types. Third, we perform automatic classification of the different neuron types based on the extracted parameters.

### Model Dynamics

To describe subthreshold voltage and spike generation, we studied a family of IF models augmented with a spike-triggered current η(*t*), a moving threshold γ(*t*), and stochastic spike emission via escape noise.

#### Subthreshold voltage dynamics.

The somatic voltage is deterministic and follows the differential equation:
*C*, *g*_{l}, and *E*_{l} are the passive membrane parameters of the neuron: the membrane capacitance, leak conductance, and the resting potential, respectively. *I*(*t*) is a time-dependent input current. Adaptation of the subthreshold membrane potential is mediated by an adaptation current η that is triggered at the firing time *t̂ _{j}*; contributions from previous spikes in the spiking history accumulate by summation of the contributions of η(

*t*−

*t̂*) over all spike times {

_{j}*t̂*} = {

_{j}*t̂*

_{1},

*t̂*

_{2},…}. By convention, the current η in

*Eq. 3*is depolarizing when its amplitude is positive and hyperpolarizing otherwise. A short hyperpolarizing current mediates refractoriness, while a current with a long time constant leads to spike frequency adaptation because the effect of multiple spikes can accumulate. Even though the reset of the voltage after a spike is equivalent to a short hyperpolarizing pulse, we use an explicit reset that is discussed below. We approximate the adaptation current η(

*t*−

*t̂*) sometimes by a single exponential (IF

_{j}_{1}) or by two exponentials (IF

_{2}), or we keep the time course arbitrary (IF

_{η}).

#### Stochastic spike emission.

The deterministic voltage dynamics of *Eq. 3* is integrated so as to yield a voltage *V*(*t*) = *V*(*t*|*I*, θ, {*t̂ _{j}*}), which depends on the input current

*I*(

*t′*) for

*t′*≤

*t*; on the past firing times

*t̂*

_{1},

*t̂*

_{2},

*t̂*

_{3},… <

*t*, as well as the set of parameters (denoted as θ) chosen for the model. This deterministic voltage is the central variable for the stochastic spike generation. Using the escape-rate picture (Plesser and Gerstner 2000; Gerstner and Kistler 2002; Paninski et al. 2005), the conditional probability intensity of emitting a spike is a nonlinear function of the instantaneous distance from the voltage threshold. More precisely the stochastic spiking process follows an inhomogeneous point process with conditional firing intensity λ(

*t|V*,

*V*

_{T}) given by:

_{0}has units of s

^{−1}, so that λ(

*t*) is in

*Hz*,

*V*(

*t*) is the membrane potential,

*V*

_{T}(

*t*) is the voltage threshold, and Δ

*V*describes the sharpness of the exponential function. In principle, any function of

*V*−

*V*

_{T}would be possible, but it was shown previously that the exponential function explains experimental data well (Jolivet et al. 2006). Note that the value of λ

_{0}can be chosen arbitrarily since a scaling λ

_{0}=

*a*λ

_{1}= λ

_{1}exp(log

*a*) can be compensated in

*Eq. 4*by a shift of the threshold

*V*

_{T}by Δ

*V*log(

*a*).

When a spike is emitted the numerical integration is stopped and after a short absolute refractory period *T*_{refr}, *V* is reset to *E*_{reset}. This voltage reset is typical for a large class of IF models.

#### Moving threshold.

The voltage threshold *V*_{T}(*t*) can be either a constant threshold, with *V*_{T}(*t*) = *V*_{0} or a time-dependent sliding threshold that implements an additional source of refractoriness (Gerstner and Kistler 2002; Chacron et al. 2003; Badel et al. 2007). In our model, the dynamic threshold is a cumulative function of previous spike times:
*V*_{0} is the threshold baseline, {*t̂ _{j}*} is the set of spike times that were emitted before

*t*, and γ(

*t*) is the spike-triggered shift in voltage threshold (see

*Shape of the Dynamic Threshold*for details). When a dynamic threshold

*V*

_{T}(

*t*), as specified in

*Eq. 5*, is used instead of the constant threshold, we denote it with “+dyn”. For example, IF

_{2}+ dyn is an IF model with two exponential spike-triggered currents and the dynamic threshold defined by γ(

*t*).

### Fitting Procedure

To extract all the parameters of the neuron model a four-step procedure is applied. The main steps of the method are illustrated in Fig. 2. For the details of the fitting procedure and the definition of appropriate subsets of the available data, see appendix d.

#### Step 1.

We measure the average spike shape by computing the spike-triggered average voltage (Jolivet et al. 2004). The refractory period *T*_{refr} and reset potential *E*_{reset} are defined by the value of the minimum of the afterhyperpolarization (AHP; see Fig. 2*B*, *i*). For Exc neurons, the AHP trace does not have any local minimum, and we arbitrarily chose *T*_{refr} = 4 ms and fixed *E*_{reset} to the corresponding voltage at this moment. As long as it remains short, the exact value of *T*_{refr} has no impact on the quality of the fit.

#### Step 2.

We extract from the experimental voltage traces the first-order estimate of the time derivative using the finite difference: *V̇*_{t} = (*V*_{t+1} − *V*_{t})*/dt*. Here, *V*_{t} denotes the binned voltage time-series as obtained from the recordings, using a bin size *dt* of 0.05 ms. We use the data set Ψ_{2} where spikes are removed (see appendix c) so as to optimize the parameters by minimizing the sum square error between the observed voltage time-derivative *V̇*_{t} and that of the model. The IF model given in *Eq. 3* with η(*t*) expressed as a linear combination of rectangular basis functions is linear in the parameters, and it is straightforward to obtain the optimal set of parameters with a multilinear regression on the derivative of the voltage (Weisberg 2005; Paninski et al. 2005; Huys et al. 2006; see *Shape of the Spike-Triggered Current* and *E**stimation of the Subthreshold Voltage Paramters* for details).

#### Step 3.

To extract the cumulative dynamic threshold γ(*t*) from the data, we maximize the likelihood of generating the experimental spike train by our model. The log-likelihood for a spike train can be written in terms of the probability *p*_{t} of observing no spike in an experimental time bin by using Bayes theorem recursively in time (Paninski 2004):
*t̂ _{j}*} is taken to be 0.5 ms before the peak of the spike and Ψ̃

_{2}contains segments of the voltage trace with the spike times removed (see appendix c). The optimal set of parameters θ̂

_{3}is obtained by convex maximization of the log-likelihood function with respect to the parameters (Paninski et al. 2005). The probability to spike (1 −

*p*

_{t}) or not to spike

*p*

_{t}is calculated from the parameters of the model via

*Eq. 4*(see

*Shape of the Dynamic Threshold*and

*Estimation of the Voltage Threshold*for details).

#### Step 4.

At this point, all the parameters have been extracted, but to obtain an optimal spike time prediction in terms of the spike train similarity measure *M*_{d}* (see *Method Spike Train Similarity*), we recompute the baseline threshold *V*_{0} so that it maximizes *M*_{d}* (*V*_{T}). To do so, we find the parameter *V*_{0} that maximizes *M*_{d}* through an exhaustive search over a large range of parameters *V*_{0}.

### Efficiency and Accuracy of the Fitting Method on Surrogate Data

Before turning to experimental data obtained from cortical neurons, we checked the consistency of our parameter extraction method on artificially generated data. We first generated 60 s of surrogate data from an IF model augmented with a spike-triggered current decaying exponentially with two time constants and a dynamic threshold (IF_{2} + dyn) and used our fitting method to retrieve its parameters. Thus we used an IF_{2} + dyn model to fit data from another IF_{2} + dyn model. The error in the estimated parameters and thus the prediction performance of the model depends on the amount of data used for fitting.

Figure 3 shows how the fitting quality evolves as a function of the amount of data used for parameter extraction. We observe that the voltage prediction of the fitted model becomes better when a larger amount of data is used (see Fig. 3, *A* and *C*). However, even with a small amount of data (1 s), the predicted voltage is relatively accurate, producing a RMSE of 0.43 mV (the RMSE obtained with 15 s of data is 0.26 mV). The spike time prediction also depends on the size of the training set, leading to a spike train similarity *M*_{d}* of 0.79 with 1 s of data while *M*_{d}* = 0.99 when 15 s are used. Furthermore the relative error ε in the parameter estimate (computed according to *M*_{d}*, and ε at each training set size. We find that if <14 s of recordings are used (Fig. 3*C*, shaded area), parameters are significantly different from their steady-state value as quantified by ε (two-sample *t*-test, α = 0.05). Therefore, when our reference model is fitted to itself, 14 s of recordings with a 10-Hz firing frequency are sufficient to obtain a good fit in terms of RMSE, *M*_{d}*, and ε.

Since fitting a model to data (model extraction) is more difficult than fitting to itself (model identification), for the fits to real neurons (as done in the next subsection) we always use a training set of ≥20 s and a separated test set that contains ≥100 spikes with firing frequencies from 0.1 to 40 Hz.

We checked that the fitting procedure works for other stimulation protocols. To study this, we used the stimulated IF_{2} + dyn model with various kinds of widely used currents: series of steps, ramp, white noise, and colored noise. We tested the fitted models on a test set made with the synaptic-like current. The results are shown in Table 1. We conclude that step currents or synaptic-like currents are more informative and enable a more reliable parameter extraction than ramp currents and (white or colored) noise injection. This can be explained by two facts. First, constant mean noise yields models with almost constant firing rate. The produced ISI distributions do not span a sufficient range of timescales to allow a fair extraction of the threshold parameters that potentially act on long timescales. Second, since we use a training set of fixed length, the distribution of the spike times (i.e., the spiking pattern) plays an important role in the total amount of available data for the parameter extraction. Indeed, for the estimation of the model parameters we have to remove a period of time around each spike; thus there are more available data when the model fires irregularly with long period of silence interleaved with short period of high activity (which is the case with synaptic-like current and series of steps) than when the model fires regularly (white and colored noise). Both facts together explain why parameter extractions fails with white noise stimuli with constant mean and SD. The problem occurring with the ramp is of another nature. Indeed the ramp stimuli do not carry enough information to allow a fair parameter extraction: the time derivative is almost constant, which is insufficient to evaluate the membrane time constant τ and the passive conductance *g*_{l}.

### Quantitative and Qualitative Accuracy of Fitted Models on FS, NFS, and Exc Neurons

The accuracy of the fitting method is summarized in Fig. 4 for exemplars of FS, NFS, and Exc neurons. With voltage prediction consistently below the intrinsic RMSE of the data (estimated across several repetitions, see materials and methods) and spike time prediction above a spike train similarity of *M*_{d}* = 0.78, we conclude that the simple neuron models, in combination with our fitting method, accurately model the three different cell types.

Our neuron models extracted by the above procedure are able to reproduce the typical behavior of the different cell types in terms of their firing patterns in response to step currents. To show this, we reproduce an experiment done by Tateno et al. (2004), where the authors stimulate GABAergic FS and Exc neurons with step currents of 1 s at various amplitudes and classify cells as a function of their frequency-intensity curves (*f-I* curve). Using this method, the authors conclude that GABAergic FS neurons have a step in the *f-I* curve (type-II excitability) whereas Exc neurons exhibit a smooth *f-I* curve (type-I excitability).

Figure 5, *A*, *B*, and *C*, which have been generated using our optimal models for GABAergic FS, GABAergic NFS, and Exc neurons, are analogous to Fig. 4 of the study of Tateno et al. (2004). Our model of FS neurons is of type-II exhibiting a minimal frequency at a critical input amplitude, with *f*_{c} = 15.5 Hz (when *f*_{c} is computed as the inverse of the first interspike interval) and a steady-state critical frequency of 5.45 Hz. When constant currents are used, smaller frequencies are possible but the *f-I* curve always exhibits a finite jump. The critical frequency and emergence of type-II behavior can be traced back to the facilitating component of the effective moving threshold (Fig. 5*A*, *inset*, and discussion for details). If the amplitude of the stimulation current is sufficient to evoke a spike, then the facilitating part of η(*t*) causes repetitive firing at nonzero frequency that is maintained as long as the stimulating current is maintained. For NFS and Exc models, we obtained a smooth transition of the *f-I* curve between silence and repetitive firing, which is the behavior of the type-I neurons consistent with experiments of Tateno et al. (2004). Links between the type of *f-I* curve and the bifurcation theory have been previously established (Izhikevich 2007; Naud et al. 2008). These relations are discussed in the *Links with Bifurcation Theory*. Finally, Fig. 5, *D–F*, provides examples of firing patterns from our optimal models for each cell type, analogous to the examples provided in Fig. 5 of Tateno et al. (2004).

We conclude that the fitted models predict quantitatively the subthreshold voltage and the spike times as well as qualitatively features such as the firing patterns and the type of excitability.

### Essential Features for Subthreshold Voltage Prediction

The membrane filter κ(*t*) (described in appendix b) for each neuron type is well approximated by a single exponential function (Fig. 6, *A–C*). Since the extraction method of the filter (see appendix b) is flexible enough to extract resonances, or multiple exponentials, the absence of resonances and the presence of a single timescale shows that voltage-dependent subthreshold currents (Richardson et al. 2003; Izhikevich 2007) are small and can be neglected. Moreover, it follows from this finding that, for the neuron types studied here, subthreshold resonance is not the most important factor for accurate prediction of the voltage traces. Specific spike-triggered currents, however, are necessary to reproduce the subthreshold voltage.

The spike-triggered current η(*t*) corresponds to the stereotypical current that flows into the neurons after a spike. After the onset of the spike, the dominant features of the spike-triggered currents consist of *1*) the current that produces the spike (i.e., the spike shape), *2*) a short refractory current that follows the end of the AP (just a few milliseconds after the spike onset), and *3*) a long cumulative current that can adapt the spike frequency of the neurons. Here, we consider only *parts 2* and *3* since these currents are the most important part for the processing of information done by the neuron (Koch 1999; Hille 1992; Jolivet et al. 2008a). To investigate the shape of the spike-triggered currents, we measured the cumulative adaptation current η(*t*) for each neuron and group these by cell type. Figure 6, *D–F*, shows the mean adaptation current η(*t*) (black lines), one example for each neuron class.

Once we have η(*t*) in terms of rectangular basis functions, we want to test whether a simpler model made of a combination of a small number of exponential processes also works. To this end, we fit *N* exponentially decaying functions with different amplitude *b*_{i} and time constant τ_{i} on the extracted η(*t*). An IF model made with *N* exponential adaptation currents will be called IF_{N}, so that IF_{2} means an IF model with two adaptation currents. A model where the shape η comes from the raw fit of η is called IF_{η}. We explore the voltage prediction in terms of RMSER mentioned in materials and methods, where the lower the RMSER, the better the prediction. A RMSER < 1 is possible, because we systematically overestimate the intrinsic RM̃SE. Here we investigate models with static thresholds, since the dynamic threshold has no impact on the predicted voltage when the spikes of the model are enforced at the correct spike times.

Figure 6*D* shows the average spike-triggered current of the nine GABAergic FS neurons. We observe that this spike-triggered current has two main parts, a strong and fast hyperpolarization that prevents repetitive firing during the first 40 ms, followed by a weaker but longer depolarization that lasts for 350 ms. This resonance is distinct from strictly subthreshold (resonating) membrane currents (Hutcheon and Yarom 2000) since, as discussed above, no resonance was observed in the membrane filter κ(*t*). Note that the shape of the spike-triggered current is similar to the feedback kernel of the FS neurons observed in (Tateno and Robinson 2009). The time course η(*t*) of the spike-triggered current can be well approximated by a double exponential decay (IF_{2} with *b*_{1} = −111.61 pA, τ_{1} = 36.86 ms, and *b*_{2} = 72.64 pA, τ_{2} = 61.76 ms). The membrane potential prediction also shows that *N* = 2 time constants are necessary and sufficient for optimal RMSER (Fig. 6*G*).

The NFS GABAergic neurons show a simpler spike-triggered current that only mediates hyperpolarization and that can be fitted with a single exponential function (IF_{1} with amplitude and time constant *b*_{1} = −29.02 pA, τ_{1} =34.58 ms; Fig. 6*E*, blue trace). This current produces the relatively weak adapting behavior of the GABAergic NFS neurons, characteristic of their firing patterns (Kawaguchi and Hama 1987). The membrane potential prediction also shows that a single time constant (*N* = 1) is necessary and sufficient for optimal RMSER (Fig. 6*H*).

Exc neurons have a stronger and longer adaptation current η(*t*) (Fig. 6*F*) than the GABAergic NFS and the GABAergic FS cells, which mediates the regular spiking (accommodating) behavior of these cells. Again this current is well approximated with a monoexponential function, with *b*_{1} = −48.35 pA, τ_{1} = 44.89 ms. Moreover, we observed more variability in η(*t*) across individual cells for Exc neuron than for the two other groups. The membrane potential prediction also shows that a single time constant (*N* = 1) is necessary and sufficient for optimal RMSER (Fig. 6 *I*).

From these results, we conclude that the shape and dimensions of the adaptation current η(*t*) are cell-type specific. Moreover, we observe that η(*t*) in GABAergic NFS and Exc differ only by their time scale and amplitude whereas η(*t*) in GABAergic NFS and Exc have a shape distinct from η(*t*) in GABAergic FS. We also notice that independent of the neuron type, spike-triggered currents extend for a few hundred milliseconds, so that due to the cumulative effect, the spike-triggered current influence the spike-frequency adaptation of the neurons on long timescales. This suggests that adaptation currents might mediate some aspect of cell-type specific behavior (i.e., firing patterns).

In summary, we found that *1*) the adaptation currents of different neuron classes reflect their typical firing behavior (see Fig. 5), *2*) adaptation currents act on multiple timescales (FS: 37 and 62 ms, NFS: 35 ms, and Exc: 45 ms), and *3*) NFS and Exc have strictly hyperpolarizing currents (leading to spike-frequency adaptation), but FS have both hyperpolarizing and depolarizing currents. From these results, it is clear that spike-triggered currents can cause different types of spike-frequency adaptation. However, some other mechanisms can also contribute to adaptation such as a fatigue of AP initiation mechanisms (Kobayashi et al. 2009; Benda et al. 2010). This is discussed in the next section.

In the previous paragraphs, we discussed the importance of spike-triggered currents. However, we know that these mechanisms are mediated by a spike-triggered change in conductance rather than current (Schwindt et al. 1988a,b). To address this issue, we fit spike-triggered conductances instead of spike-triggered currents and look at the magnitude of the improvements that follow (see *Conductance-Based Adaptation* and *Estimation of the Subthreshold Voltage Parameters* for a description of the conductance-based spike-triggered adaptation). Figure 7 *A* shows the movement of the cumulative change in conductance η_{cond} following each spike for Exc neurons. The RM̃SE depends on the reversal potential *E*_{rev} and shows a minimum at *Ê*_{rev} = −51.89mV (Fig. 7 *B*). Note that *E*_{rev} could be attributed to the reversal potential of potassium since most of the spike-triggered currents are mediated by potassium ion channels. However, the high value we found indicates that *E*_{rev} reflects a mixed ion reversal, which, in addition to potassium, also includes sodium and calcium.

The shape of the spike-triggered conductances η_{cond} is more difficult to interpret than the standard spike-triggered current η, because of its dependency on the reversal potential *E*_{rev}. In fact, the effect of η_{cond} on voltage depends on the instantaneous difference between the actual membrane potential and the reversal potential. Furthermore one can observe that given the mean voltage reset *E*_{reset} = −35.56 mV of the Exc neurons, that is above *E*_{rev}, the effect of the spike-triggered conductance after a spike is to hyperpolarize the membrane potential and thus η_{cond} also mediate spike-frequency adaptation. Figure 7, *C* and *D*, shows the averaged RMSER of voltage prediction as well as the spike train similarity measure *M*_{d}* for the three cell types and with the two model variants, IF_{η} and IF_{ηcond}. One can observe that conductance-based adaptation does not lead to any significant improvements in terms of spike train similarity measure *M*_{d}* and the subthreshold voltage prediction RMSER (two-sample *t*-test, *P* > 0.2 for *M*_{d}* and *P* > 0.07 for RMSER). Thus we conclude that conductance- or current-based spike-triggered events make equally valid models.

### Essential Features for Spike Time Prediction

To explore the spike time prediction of our models, we compute the spike train similarity measure *M*_{d}* (see materials and methods) averaged across all neurons of a given type. For example, as mentioned in the materials and methods, a value of *M*_{d}* = 0.8 indicates that 80% of the PSTH is correctly predicted by the model. To compare the effects of adaptation current η and sliding threshold γ, we show the spike time prediction for models with and without spike-triggered currents as well as with or without sliding threshold. Figure 8, *D*, *E*, and *F*, shows our results for GABAergic FS, GABAergic NFS, and Exc neurons, respectively. The dynamic thresholds were in each case described as a double exponential, consistent with the results of Fig. 8, *A*, *B*, and *C* (green traces).

For the GABAergic FS cells (Fig. 8 *A*), we find that there is almost no movement of the threshold. The small fluctuations in the extracted moving thresholds are presumably due to noise in the estimation. This implies that all the adaptive behavior of GABAergic FS cells is mediated by the spike-triggered current and not by any changes in the AP threshold. We did not find that adding a moving threshold yields any significant improvement in terms of spike-time prediction (*P* > 0.07 for all pairs, with two-sample *t*-test with α = 0.05). Nevertheless when IF_{0} is augmented with bin-based spike-triggered currents η (IF_{η}), we obtained a minor gain in spike time predictions (Δ*M*_{d}* = 0.05). There is negligible increase in *M*_{d}* when augmenting the model further with sliding threshold (IF_{η} + dyn; Δ*M*_{d}* = 0.007). The optimal model for GABAergic FS neurons is IF_{η} + dyn producing a *M*_{d}* of 0.87 ± 0.06, but the gain compared with other model variants is marginal.

The moving threshold γ(*t*) of GABAergic NFS neurons follows a double exponential decay, with parameters *b*_{1} = 3.64 mV, τ_{1} = 21.88 ms for the fast component and *b*_{2} = 1.24 mV, τ_{2} = 336.50 ms for the late component (Fig. 8*B*, green trace). After half a second to a second, γ(*t*) is weak but by cumulating over multiple spikes the late component can contribute to spike-frequency adaptation. When comparing the value of spike time prediction *M*_{d}*, a model IF_{0} performs always worse than models augmented with spike-triggered currents η or dynamic AP threshold γ. More precisely, adding a sliding threshold produces a highly significant gain in the spike prediction measure of Δ*M*_{d}* = 0.10. Adding a spike-triggered current produces a net gain of Δ*M*_{d}* = 0.07 compared with IF_{0}. This leads us to the optimal model for GABAergic NFS cells being IF_{η} + dyn with a value of spike time prediction *M*_{d}*(IF_{η} + dyn) that indicates that >90% of the spikes in the PSTH are indeed predicted by the model.

The Exc neurons also exhibit a dynamic threshold with a double exponential decay with *b*_{1} =12.45 mV, τ_{1} = 37.22 ms and *b*_{2} = 1.98 mV, τ_{2} = 499.80 ms (Fig. 8*C*, green trace) with an amplitude that is at least twice as strong as for GABAergic NFS cells. The effect of the spikes on the threshold lasts for >1 s. The case of the Exc neurons is special because they show a voltage reset to a value above the threshold baseline (*E*_{reset} > *V*_{T}). Therefore, all models with static threshold produce repetitive firing at very high and unrealistic frequencies, which leads to a very low *M*_{d}*, as one can observe on Fig. 8*F*. Thus models upgraded with a sliding threshold always generate a significantly higher *M*_{d}* than similar models with static threshold (Δ*M*_{d}* = 0.48). We also observe a small increase of Δ*M*_{d}* = 0.07 for models upgraded with spike-triggered currents. Thus the adaptation process is dominated by the effect of the sliding threshold (also see Fig. 11*C*). The optimal model for Exc neurons is an IF_{η} + dyn that produces *M*_{d}* = 0.81 ± 0.04. The ISI distributions of the data (Fig. 8, *G-I*) agree with the ones coming from the optimal models but not with those from the simple model without spike-triggered adaptation. We summarize results from Figs. 6 and 8 by observing that the minimal optimal model for GABAergic FS cells must have two spike-triggered currents and a static threshold (IF_{2}), whereas for GABAergic NFS and Exc neurons the minimal model must have only one spike-triggered current but a dynamic threshold (IF_{1} + dyn). These minimal models reproduce the subthreshold voltage traces of the experiments (RMSER ≤ 1) and ≥80% of the spikes (scaled to experimental reliability).

We observe that GABAergic NFS and Exc neurons have AP threshold dynamics that extend from milliseconds to >500 ms. Moreover, due to its cumulative property, the moving threshold can tune the neuron's firing frequencies, and thus the PSTH, on timescales beyond 1 s. Finally, AP threshold dynamics are only present in some cell types and, when present, act on very long timescales. We also note that it is the effect of the dynamic threshold γ(*t*) combined with the spike-triggered current η(*t*) that produces the effective adaptation behavior of a given neuron.

### Dependency of the Extracted Parameters on the State of the Neurons

Given the intricate voltage dependence of ion channels, it is clear that spike-triggered currents as well as passive parameters of neurons change as a function of the instantaneous state of the neurons. It is conceivable that the intensity and the time course of the spike-triggered currents would be affected by the instantaneous firing frequency of the neuron (Tateno and Robinson 2009). To investigate this point, we applied our method to subsets of the data that differ in mean firing rate and in mean membrane potential. In other words, we studied the dependence of the extracted parameters on the state of the neuron. To do so, we split the 60-s long experiments in 6 classes of 10 s and computed the mean firing frequencies and the mean voltage of these classes. Note that we used only six classes because as mentioned in *Efficiency and Accuracy of the Fitting Method on Surrogate Data* we need ≥10 s of recordings to obtain a fair estimate of the parameters. Moreover, as mentioned in *Fitting Procedure* our method required data that covers a large range of ISIs to obtain a good estimate of the threshold parameters γ. With six classes both requirements are fulfilled. In Fig. 9*A*, one can observe that a single 10-s class contains high and low rate episodes interleaved with periods of silence. This produces an ISI distribution covering a wide range of instantaneous rates (from 0 to 25 Hz; Fig. 9*B*, *bottom*) estimated from episodes of 0.5 s. If averaged over segments of 10 s, the apparent range of firing frequency is smaller (from 1 to 6 Hz; Fig. 9 *B*, *top*) but contains the exact same episodes. The six classes of 10 s enable us to investigate the dependency of extracted parameters on the neuronal regime.

Using these six classes, we could not find important dependencies of the spike-triggered current η (Fig. 9, *C* and *D*) or the passive parameters *C*, τ, and *E*_{l} (Fig. 9, *E-J*) on the firing rate or on the mean voltage. Moreover, it is impossible to distinguish parameters extracted on a single class from parameters extracted on the whole 60-s long experiment.

These results suggest that over the range of firing frequencies that we have studied, the parameters of the neurons do not change significantly. This tends to justify our assumptions that a unique set of parameters per neuron is sufficient to capture its dynamics. This does not exclude that a neuron behaves differently when driven to some exotic regime (in fact, it definitely will), but over the range studied here, our method did not indicate important parameter changes. Other experiments, where rates vary more drastically would be needed to investigate these points in a more systematic way.

### Cell-Type Classification

In the last sections, we showed that the membrane filter, the time course of adaptation, and the AP dynamics strongly depend on the neuron type. Here we ask whether we can characterize cell types solely on their extracted parameters. We classify cell types based on *1*) their passive parameters *C*, *g*_{l}, and *E*_{l}, *2*) the parameters describing the adaptation current including the value of voltage reset *E*_{reset} (note that here we use a combination between η and *E*_{reset}, because the extracted spike-triggered current depends on the voltage reset, and in this view, *E*_{reset} is a part of the adaptation process), *3*) the parameters describing the shape of the dynamic threshold γ, or *4*) all the parameters. To do so we perform standard principal component analysis (PCA) using *1*) only the passive parameters for each cells, *2*) only the parameters of adaptation current, *3*) only the parameters of the dynamic threshold, and *4*) using all the parameters. Using any subset of parameters is sufficient for classification with the first two principal components (Fig. 10, *B–D*) except when passive parameters only are used (Fig. 10*A*). Classification based on passive parameters fails because GABAergic FS and NFS neurons do not differ in a significant way in their capacitance, their leak conductance, or their reverse potential. However, if PCA is applied on the parameters that characterize the adaptation current and/or the dynamic threshold, we can successfully classify neurons. Moreover, we also observe that the variance between the different cells is mainly explained by the reset value *E*_{reset} and the dynamic threshold (Fig. 10*D*, *right*).

## DISCUSSION

### Automatic Fitting Method

There is a rich history of fitting neuron models to intracellular recordings of real neurons (Vanier and Bower 1999; Rauch et al. 2003; Keren et al. 2005; Achard and De Schutter 2006; Huys et al. 2006; Jolivet et al. 2006; Kobayashi and Shinomoto 2007; Druckmann et al. 2007; Badel et al. 2008; Kobayashi et al. 2009). The variety of approaches arises from the choice of neuron model and the fitting method. Still, not all methods yield models that can predict the spike times and membrane potential with high accuracy (Jolivet et al. 2008a; Gerstner and Naud 2009). To predict the membrane potential with Hodgkin-and-Huxley compartmental models, one needs prior knowledge on the dynamics of the ion channels present in the recorded cell (Huys et al. 2006; Druckmann et al. 2008). Without knowledge of the ion channels, fitting Hodgkin-and-Huxley compartmental models becomes plagued with local minima, and there are often multiple parameters that share the same fitting quality (Achard and De Schutter 2006). The only hope for a fitting method that can easily be applied to multiple systems for which we have insufficient knowledge of the ion channel dynamics is to use convex fitting methods in combination with IF models (Paninski 2004). Earlier work (Jolivet et al. 2006) had an efficient fitting method for the subthreshold voltage, but the black-box fitting of the adaptive threshold was not convex. The method of Badel et al. (2008) may have been convex, but it applied only to models without spike-frequency adaptation. The method of Paninski (2004) was convex but did not use the information contained in the voltage trace while the method of Paninski et al. (2005) was convex and used the voltage trace but lacked the moving threshold required for efficient spike prediction. In this study, we have used a method that improves on the earlier multilinear regression method (Paninski et al. 2005) by adding a second fitting step for the moving threshold taken from the literature on generalized linear models (McCullagh and Nelder 1998). The method is sure to find only one set of optimal parameters because it is made of two convex fitting methods [multiple linear regression and generalized linear model with Poisson or Bernoulli probability distribution, but see constraints for the convexity in Paninski et al. (2004)]. We expect the method to generalize well to many cell types because the total time of the spike-triggered current, number, and size of the basis functions are not expected to depend on the cell types.

The notion of a dynamic threshold affecting neuronal computation also has a long history (Hill 1936; Azouz and Gray 1999, 2000; Gerstner and Kistler 2002; Chacron et al. 2003). So far, the methods for fitting the dynamics of the firing threshold have relied on the measurement of the effective threshold for each spike (Azouz and Gray 2000; Chacron et al. 2003). Instead, the method presented in this article uses the whole voltage trace, providing information about the firing threshold each time a transient increase in the membrane potential is not followed by a spike. We expect this method to be more precise since the number of data points used to constrain the moving threshold is not proportional to the number of spikes but to the number of data points constituting the subthreshold voltage trace in the regime close to threshold.

Since the choice of model will affect the prediction performance the question arises, why we omit the nonlinearity responsible for spike initiation in IF-type models (Fourcaud-Trocme et al. 2003; Badel et al. 2008)? An exponential nonlinearity in IF models was shown to be crucial for accurate processing of the inputs at high frequency (Fourcaud-Trocme et al. 2003). Such a nonlinear term can interact with subthreshold currents to produce a variety of firing patterns (Izhikevich 2004, 2007; Naud et al. 2008). The IF models used here assume strictly linear voltage dynamics and are fitted away from the spikes so that the increased nonlinearity close to a spike does not bias the parameter estimation. It is not trivial to generalize the present convex method to include an exponential nonlinearity for spike initiation. However, in the present model the spike initiation is stochastic with an exponential nonlinearity for the probability of spiking as a function of voltage as extracted from experimental data (Jolivet et al. 2006). This exponential nonlinearity should not be confused with the exponential nonlinearity in the adaptive exponential IF (AdEx) model (Brette and Gerstner 2005) or the model of (Badel et al. 2008). Nevertheless, there are some links between the exponential spike initiation of the AdEx model and the exponential transfer function of the generalized IF models discussed here, so that the AdEx model can be approximately mapped onto a generalized IF model (Mensi et al. 2011). The escape-noise IF models discussed here can reproduce all the main firing patterns, except delayed spike initiation upon pulse current input (see below).

Another feature of some neurons that is not present in our model is a subthreshold resonance. A strong subthreshold resonance has been observed in the dendrites of large Exc neurons (Cook et al. 2007) and in “mes V” neurons (Izhikevich 2007). Subthreshold resonance is thought to be mediated by an additional current linearly coupled with the membrane potential. Similarly, delayed spiking upon step current injection has also been attributed to an additional current linearly coupled with the membrane potential (Naud et al. 2008). One can check the presence of such an additional effect by looking at the shape of the membrane filter κ(*t*) extracted with a method that does not force an exponential shape [i.e., the Wiener-Hopf optimal filter method described in appendix b and by Jolivet et al. (2004)]. Indeed, a resonance corresponds to a filter with a negative undershoot, while a delayed onset corresponds to a filter with a double exponential decay. Both cases can be described by an additional current having linear coupling with the soma to take into account resonance or delayed spiking. We have tested a method which involves adding the term *ae*^{−t/τw} * *V*(*t*) to *Eq. 3*. Multilinear regression can still be applied to determine the strength of the coupling *a* when a time constant τ_{w} is assumed. Iterating through a large range of possible τ_{w} by repeating the multilinear regression usually yields a convex function of the mean square error as a function of τ_{w}. This method enables us to fit the parameters mediating subthreshold resonance or delayed spiking, but this was not necessary here since the neuron types studied have no resonance or do not display delayed spiking.

The choice of an appropriate input stimulus is of crucial importance for all model identification methods. Here we have successfully tested our method on different kinds of input currents injected in current-clamp mode. However, in vivo neurons receive excitatory and inhibitory synaptic inputs acting in conductance. It is possible to mimic this complex input scenario in vitro by patching neurons in dynamic-clamp mode, where excitatory and inhibitory conductances are dynamically injected in a patched neuron. We do not directly investigate the robustness of our method in this particular input scenario in this study. However, we have shown that our fitting procedure can handle dynamic-clamp input by participating in an international competition in 2008, where the goal was to predict the spike times of a single neuron stimulated in dynamic clamp [see Jolivet et al. (2008b) and Jolivet et al. (2008a) for details and results]. For this competition, we used a related model and fitting procedure and obtained similar results.

### Effective Moving Threshold

We have seen that the effects of spike-triggered currents and dynamic threshold merge to produce spike-frequency adaptation. The effect of the spike-triggered adaptation current η(*t*) on the voltage is simply given by the convolution of the spike-triggered current η(*t*) with the membrane filter κ (η_{v}(*t*) = −κ * η). The effective moving threshold ζ(*t*) arises from the combination of the threshold dynamics γ(*t*) and the effect in voltage of the spike-triggered current. Since spikes are triggered when the membrane potential hits the threshold, the relevant variable is the difference ζ(*t*) between the change in dynamic threshold (increasing after a spike) and the contribution of the adaptation current to the voltage trace.

To answer the question of whether the spike-frequency adaptation is dominantly mediated by spike-triggered currents or moving threshold, we can look at the respective contribution of γ(*t*) and η_{v}(*t*) towards the effective adaptation ζ(*t*). Figure 11*A* shows that the effective moving threshold of GABAergic FS cells are clearly dominated by spike-triggered currents. For the Exc cells, ζ(*t*) is dominated by the moving threshold (see Fig. 11*C*), whereas for the GABAergic NFS cells the effective adaptation process ζ(*t*) (Fig. 11*B*) is a combination of η and γ, where each spike-triggered mechanism mediates approximatively half of the effective moving threshold. Therefore, we conclude that the adaptation is mediated by different processes in different cell types: adaptation is dominated by the moving threshold for the Exc neurons; caused entirely by spike-triggered currents for the FS neurons; and consists of an equal mix of threshold and current for the NFS neurons.

The effective moving threshold for GABAergic FS cells has the particularity of crossing zero after 30–80 ms and remaining negative (i.e., facilitating) thereafter. The zero-crossing then determines the type of excitability: under constant current injection, after the neuron model fires its first spike, the effective threshold is first high, decreases, and eventually crosses the membrane potential, hence forcing a spike after a period approximately equal to the first zero-crossing in the function ζ(*t*). One can see that under constant current injection this type of neuron will not fire with a period longer than the time of the minimum of the effective threshold ζ(*t*) : *f*_{c} ≈ 12–33 Hz, in agreement with Fig. 5. The facilitating tail of the effective adaptation must lead to spike-frequency facilitation as observed in the firing patterns of Fig. 5. Moreover, the peak of the facilitating part in η(*t*) will indicate a preferred frequency around 10 Hz. On the other hand, strictly decaying functions ζ(*t*) as in NFS and Exc will produce adapting firing patterns and type-I excitability.

### Links with Bifurcation Theory

We observe that our FS models have a type-II exctitability, whereas our Exc and NFS models are of type-I. In the literature these type-I and type-II behaviors are known to occur via different types of bifurcations (Koch 1999; Gerstner and Kistler 2002; Izhikevich 2007; Naud et al. 2008; Touboul and Brette 2009). For instance, type-I may occur via a saddle-node bifurcation onto an invariant circle, whereas type-II typically occurs through an Hopf bifurcation or a saddle-node bifurcation off invariant circle. The presence of the hard reset and multiple timescales in the spike-triggered adaptation means that standard theorems of bifurcation theory in continuous two dimensions do not apply. In spite of this, we can draw analogies with two-dimensional bifurcation theory.

Hopf bifurcations are associated with subthreshold resonances: when a neuron model is stimulated with a slowly increasing current ramp, subthreshold oscillations would be observed in response to a short current pulse before stability is lost through the Hopf bifurcation. These oscillations happen even before the neuron has fired a first spike and are therefore strictly subthreshold. Our fitted FS models do not generate such subthreshold oscillations as can be deduced from the absence of any resonance in the membrane filter κ(*t*). Nevertheless, the FS neuron model exhibits resonance in the spike-triggered currents. These spike-triggered currents, which are summarized in the adaptation current η(*t*) are responsible for the minimal firing frequency that is characteristic of type-II behavior. We note that spike-triggered currents result in membrane potential oscillations after an AP, i.e., they give rise to a nonmonotonic SAP. Such an oscillatory SAP should not be confused with pure subthreshold oscillations. Let us describe in more detail how the type-II behavior arises in a simplified version of our fitted models. We simplify our models such that adaptation is mediated by a single exponential spike-triggered current of amplitude *b* that mediates spike-frequency adaptation (*b* > 0) for Exc and NFS neurons and facilitation (*b* < 0) for FS neurons. This can be done by coupling a second differential equation τ_{w}*ẇ* = −*w* to the main voltage *Eq. 3* and to add a constant value *b* to *w* each time a spike is emitted (Izhikevich 2007). Note that this makes a *w*-nullcline that is independent of the membrane potential *V*, as expected from the fact that the membrane filter is a single exponential (as shown in Fig. 6, *A–C*). Under these assumptions, the only difference between type-I and type-II models is the sign of the spike-triggered current. A schematic of the phase plane of these simplified type-I and type-II models is shown in Fig. 12. Before application of a depolarazing step current, the membrane potential is at rest (Fig. 12*A*). The stable fixed point looses its stability through a saddle-node bifurcation. The simple model can still create type-I or type-II excitability depending on whether the saddle node bifurcation is on or off the limit cycle. Facilitating spike-triggered adaptation (*b* < 0) implies that the limit cycle stays away from the location of the saddle-node bifurcation in phase space, because the reset of the *w* variable restarts at a lower value (type II; Fig. 12*B*). True adaptation (*b* > 0), however, leads to the classic saddle-node bifurcation onto invariant circle (type-I; Fig. 12*C*).

We have argued that our fitted FS models loose stability via a saddle-node bifurcation, which causes type-II excitability because it is off the invariant circle. We have discarded the Hopf bifurcation because of the absence of subthreshold resonance. It is possible, however, that the fitting procedure did not capture a subthreshold resonance that appears only at voltages close to the threshold, and there are exactly the type of resonances that would be expected when stability is lost via a Hopf bifurcation. Therefore, the present analysis cannot draw conclusions on the type of bifurcation responsible for firing in real FS neurons.

### Interpretation of Model Parameters

We extract a set of parameters for a typical model of each neuron class. The parameters are related to the underlying biophysics (density, distribution and dynamics of ion channels, membrane capacitance, and resistance), but the exact relationship is unclear. For instance, the extracted membrane time constants (τ_{m} ranged from 5 ms to 20 ms) are slightly shorter than the typical membrane time constant measured in experiments with voltage recording of the response to subthreshold step current injection at rest. This discrepancy can be explained by the fact that we measure passive parameters of neurons when they are stimulated instead of at rest. The different working regime of membrane potential will activate differently the subthreshold conductances (spike-triggered or not) and thus modify the effective membrane time constant (Richardson et al. 2003; Jolivet et al. 2004; Kobayashi et al. 2009). The presence of the electrode may also bias our estimate of the membrane capacitance as discussed by Badel et al. (2008). Similarly, the reversal potential, voltage threshold, and voltage reset may depend on the bath solution and the intrapipette solution. The amplitude and time scale of the spike-triggered adaptation should not be affected by the electrode or the bath solution, but the temperature at which the experiment was performed can affect the dynamics of the underlying ion channels.

We also investigated the dependency of the extracted parameters (i.e., the membrane filter and the spike-triggered current) on the regime of the neuron. Somewhat surprisingly, we did not observe significant dependencies, in spite of what has been already observed (Tateno and Robinson 2009). Our results could be explained by the limited range of firing rates that we studied and the nature of our stimulation protocol. Indeed, we have restricted our analysis to current-clamp experiments and it is known that neurons exhibit different behaviors under in vivo conditions (Prescott et al. 2006, 2008). Those in vivo conditions can be approximated in vitro by the dynamic clamp. The dynamic-clamp allows the study of the neuron behavior as a function of the relative strength and nature of the input conductances (i.e., balance between excitatory and inhibitory conductances, and shunting). For instance, the integration properties of Exc neurons change drastically depending on the conductance state (Jolivet et al. 2004; Prescott et al. 2008). Moreover, these changes also affect the adaptation mechanisms of the neurons. The method presented in this article can be applied to dynamic-clamp recordings to address these issues.

The effect of the spike-triggered current on the voltage η_{v}(*t*) is closely related to the SAP (Sah 1996). The SAP may differ from η_{v}(*t*) since it is measured around the resting potential while η(*t*) is an average of the spike-triggered current under synaptic-like current injection. Furthermore, the SAP is measured after a spike that was artificially triggered by a large and short current injection. The amount of charge that was injected to produce the spike will leak out of the membrane on a time scale given by the membrane time constant. The spike after current extracted by standard experimental protocols (Sah 1996) is thus biased by the current used for stimulating the spike. Yet, the close relationship between the SAP and η_{v}(*t*) indicates that η(*t*) should be mediated by the same ion channels mediating the SAP, namely: the M-type current *I*_{M} (Adams et al. 1982), the afterhyperpolarization current *I*_{AHP} (Madison and Nicoll 1984), or any other calcium-dependent ion channels (Hille 1992; Sah 1996; Koch 1999; Wang et al. 2003). Moreover, spike-triggered events in the dendrites can also shape the spike-triggered current (Doiron et al. 2007).

The movement of the threshold after a spike has been proposed to depend on sodium channel de-inactivation (Fleidervish et al. 1996; Fleidervish and Gutnick 1996; Kobayashi et al. 2009). Following a spike, a portion of the Na^{+} channels responsible for the spike initiation stays inactivated, which leads to a higher effective threshold. The sodium channels then de-inactivate, which results in a gradual decay of the spiking threshold (Henze and Buzsáki 2001; Platkiewicz and Brette 2011). It has been proposed that only a subtype of sodium channels are inactivating (Martina and Jonas 1997). Thus, in our framework, the observed dynamic threshold must be related to the proportion and the type (inactivating or noninactivating) of the sodium channel. Our results corroborate this hypothesis since only the GABAergic NFS and Exc types have a moving threshold, which suggest that the GABAergic FS neurons do not express the inactivating sodium channels.

Finally, it would be possible to investigate the biological mechanisms underlying spike-triggered adaptation by the use of specific pharmacological experiments. For instance, by blocking some specific ion channels one can study how the shape of the spike-triggered current is affected by calcium channels or high-voltage activated potassium channels. Similarly, if one can block specific sodium channels, like the inactivating Na^{+} channels, it is possible to investigate the dependency of the dynamic threshold on the type of sodium channels expressed in a given neuron type.

### Classification

Classification of neuron types can be done on multiple features: firing pattern, spike shape, morphology, and expression of molecular markers (Markram et al. 2004). Here we classify based on the computational properties of the neurons, that is, the parameters regulating how the neuron encodes the incoming current into spike trains. These computational properties are determined by the expression of ion channels and lead to firing patterns that depend on the neuron type. We have shown that the classification of neuron types relates to a classification of the computational properties beyond the classification of firing patterns (GABAergic NFS and Exc cell types are both regular spiking and accommodating neurons). In other terms, classification is possible even if the shape of the spike used traditionally for the distinction between GABAergic NFS and Exc neurons is not taken into account, since the different neuron types encode the incoming current differently.

We found that the passive properties of the neurons (capacitance, input resistance, and membrane time constant) are not sufficient to efficiently distinguish between the neurons types. It is the adaptation properties (voltage reset, spike-triggered current, and moving threshold) that distinguish the different cell types. These results would indicate that the strength and time scale of adaptation are important parameters of cortical network computation, indicating a direction for investigating the importance of microcircuitry in cortical networks.

We have studied only three types of neurons, but more types of Exc neurons (Connors and Gutnick 1990) and GABAergic neurons (Markram et al. 2004) have been described. Further work is needed to check that such extensive classification can be done on the properties affecting the conversion of synaptic inputs to a spike. Classification on a greater pool of neurons would also enable one to study how distinct the computational properties of different neuron types are and if a finer classification can be inferred. We expect that one would need on the order of 100 recorded neurons under synaptic-like current injection to study in detail the classification of cortical neurons and perform unsupervised clustering on the computational properties.

Neurons recorded in vitro can show very different properties from their alter ego in the intact, awake and behaving animal. In the awake animal, it is not yet possible to know the input a neuron receives from its synaptic connection, but the somatic voltage can be recorded (Crochet and Petersen 2006). In such an experiment, it is not possible to apply this part of the fitting method for the spike-triggered current and passive properties because the method requires the knowledge of the stimulating current arriving at the soma. The moving threshold, on the other hand, can be extracted since our fitting method only requires the voltage trace and the time of the spike. It remains to be tested if moving thresholds can effectively be extracted from in vivo recordings. It is an interesting avenue for further research since this would allow to study the correspondence of in vitro and in vivo threshold dynamics and its classification across cell types.

### How Good is Good?

Depending on the neuron type, the optimal IF model is able to predict between 81 and 91% of the spikes and reproduce the subthreshold voltage fluctuations with a precision in the range of millivolts. One can then ask: Why can we not achieve a perfect prediction (i.e., *M*_{d}* = 100% and a RMSE close to 0 mV)? What is missing? We can think of three possible explanations.

First, the experimental data suffer from some unavoidable drifts that are not due to the neuron itself and that we do not model. These drifts are presumably due to some additional currents that flow out of the neuron near the patch junction and so affect the recorded membrane potential in a nonsystematic manner. The experimental drifts can greatly limit the maximum prediction performance.

Second, it is known that injecting current through the same electrode used for recording the voltage corrupts the recorded voltage (Brette et al. 2008). This artifact manifests itself as a high-frequency component of the recorded voltage that is correlated with the current being injected. Since we are fitting on the voltage trace, the artifactual component of the voltage will bias the estimated parameters. Mainly, the electrode artifact will affect the estimation of the membrane time-constant (Brette et al. 2008). We can also speculate that the average current triggering a spike (the so-called spike-triggered averaged current) will contribute erroneously to the measured spike-triggered current. However, we expect these effects to be small and to affect only very small time scales since the time constant of the electrode contribution was measured to be below the millisecond range (Badel et al. 2008). Similarly, the erroneous high-frequency component of the voltage can bias our estimation of the threshold and its dynamics. The extent to which the bias in membrane time-constant and spike-triggered adaptation affects the prediction performance would have to be studied in detail, but we have preliminary results showing that the effect is negligible. In any case, the artifact due to simultaneous current injection can be greatly reduced by the use of active electrode compensation methods (Brette et al. 2008; Badel et al. 2008).

Third, by modeling neurons with simple IF models, we neglect some nonlinearities that are present in the neurons. For instance, the voltage dynamics close to a spike become strongly nonlinear due to the activation of sodium channels (Naundorf et al. 2006; Badel et al. 2008). Saturation in the open/close fraction of the ion channels can cause the spike-triggered current of spikes in a burst to differ from the spike-triggered current of isolated spikes, leading to higher order dependencies on the spiking history. Furthermore, we know that ion channels mediating the spike-triggered currents have time constants which depend nonlinearly on the voltage, while our formalism imposes voltage-independent time constants. Similarly, the escape-noise formalism is only an approximation to the full dynamics entailed by stochastic activation of a limited number of voltage-dependent ion channels. For injection of synaptic-like current into the soma of a cortical neuron, all these approximations prove to be very good since the spike-time prediction is high and would be even higher if experimental drifts and electrode artifacts could be completely removed.

## GRANTS

This work was supported by Swiss National Science Foundation Project no. 200020-132871 and the BrainScales Project.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## APPENDIX A: VARIANTS OF THE FAMILY OF IF MODELS

We studied different variants of the main model described by *Eqs. 3*, *4*, and *5* (see *Model Dynamics*). Here we describe more precisely how these models are built and some variants of the standard models considered in this study.

### Shape of the Spike-Triggered Current

As mentioned in the results *Subthreshold Voltage Dynamics*, the spike-triggered current η(*t*) could be any function of time, and for fitting purposes we implemented this spike-triggered current as a linear combination of rectangular basis function, so that η(*t*) is given by:
*f* ^{(k)}(*t*) are rectangular basis functions centered at *T*_{k} and of width δ*T*_{k} and where Θ denotes the Heaviside step function. Since we expect the spike-triggered current to change faster close to the previous spike and slower far from it, we spaced the rectangular basis function logarithmically. The index *k* runs from 1 to *K* where *K* is the total number of rectangular functions used for fitting. The parameters *T*_{k}, δ*T*_{k}, and *K* are metaparameters fixed a priori so that the shape of the spike aftercurrents is controlled by the coefficient *a*_{k} [see also Paninski et al. (2005)]. Parameters *T*_{k}, δ*T*_{k}, and *K* are chosen to give sufficient freedom to span all plausible shapes of η while avoiding overfitting. *Equation 8* will be used to extract the precise shape of the spike aftercurrent η(*t*) without assumptions with respect to its shape (single exponential, double exponential, etc.).

### Conductance-Based Adaptation

As a variant of the model defined in *Eq. 3*, we can turn the spike-triggered current η(*t*) into a spike-triggered conductance. This is done by replacing η in *Eq. 3* by η_{cond}(*V* − *E*_{rev}), where η_{cond} has now units of conductance and *E*_{rev} is the reversal potential. If the model's adaptation is based on conductance as opposed to an adaptation current, we call the model IF_{ηcond}.

### Shape of the Dynamic Threshold

The dynamic threshold γ(*t*) described in *Eq. 5* (see *Model Dynamics* for details) is expressed as a linear combination of rectangular basis functions:
*t*) is defined by the coefficients *c*_{p}, which control the amplitude of a set of the rectangular basis *f* ^{(p)}(*t*).

## APPENDIX B: EXTRACTION OF MEMBRANE FILTER

The passive properties of a neuron, such as the reversal potential *E*_{l}, the leak conductance *g*_{l}, and the capacitance *C*, characterize the linear response of a neuron to current injection and define the membrane filter κ of the neurons. The standard exponential membrane filter arises from the solution of *Eq. 3* using the initial condition *V*(0) = *E*_{l}:
*t*) acts as a low-pass filter on the current *I* and on the spike-triggered currents:
*Equation 10* is known as the spike-response model (Gerstner et al. 1996). It summarizes the dynamics of the membrane potential as made of two parts: the effect of the input current (first term on the right-hand side of *Eq. 10*) and the effect of the afterpotential following each spike (second term on the right-hand side *of Eq. 10*). The functions κ(*t*) and η(*t*) are response kernels that depict the membrane filter and the shape of the spike aftercurrents, respectively. *Equation 10* taken without restrictions on the shape on η(*t*) and κ(*t*) is very general and can take into account the effect of multiple subthreshold or spike-triggered ion channels.

In practice, the linear filter κ(*t*) may or may not consist of a single exponential decay. It is known that subthreshold currents, like *I*_{m} or *I*_{h} (Sabah and Leibovic 1969; Mauro et al. 1970; Koch 1984) generate subthreshold resonances or delayed spiking responses to steps (Naud et al. 2008; Izhikevich 2004).

These subthreshold currents give rise to a current-to-voltage filter κ exhibiting a negative bump (for resonances) or an exponential decay with two time constants (for delayed firing onset). Thus, to make sure that our assumption of an exponential filter is not too restrictive, we extract the shape-free membrane filter κ defined in *Eq. 10*. To do so, we compute this filter by extracting the Wiener-Hopf optimal filter (Jolivet et al. 2004) on the spike-free subset of data Ψ_{3} (see appendix c for details), such that we are left with only the last term on the right-hand side of *Eq. 10*. Note that we take only subset of the data that are far from the spike to ensure that the resulting κ filter does not take into account some spike-triggered effects, or nonlinearities of the spike onset, that will corrupt the estimation of the membrane filter.

## APPENDIX C: DATA PREPROCESSING

To extract all the parameters of the models it is convenient to define appropriate subsets of the available data. Since we do not want to model the exact shape of the AP we cut out, around each spike time *t̂*, a small segment of the data which we ignore. In the experimental data, the spike time *t̂* is defined to be the time when the membrane potential crosses a given voltage *V*_{detect} from below; here we set the detection threshold *V*_{detect} to 0 mV.

Let us, at each moment *t*, refer to the last spike time as *t̂*^{(last)} and to the forthcoming spike time as *t̂*^{(next)}. The subsets of the data that we use for fitting are:
*T*_{refr} the absolute refractory period. Thus Ψ_{1} represents a set of recording times where the voltage is subthreshold and outside the absolute refractory period. The Ψ_{2} further removes 2 ms of data before each spike. We will also use two other subsets of the data:
_{3} removes a period of time *T*_{adapt} after the spikes where *T*_{adapt} > *T*_{refr}. We use for *T*_{adapt} a period of 200 ms when the recording has a high average firing frequency (>5 Hz), but we use *T*_{adapt} = 500 ms otherwise. For our recordings made of multiple repetitions, we only consider the subset of times Ψ_{4} that are separated by a period of at least *T*_{adapt} from any previous spike, from any repetition. This subset is therefore the intersection between the subsets Ψ_{3}^{(k)} of all repetitions 1 ≤ *k* ≤ *R*. Each subset Ψ_{1}, Ψ_{2}, Ψ_{3}, and Ψ_{4} will be used in different steps of the fitting procedure.

## APPENDIX D: ESTIMATION OF THE MODELS PARAMETERS

The fitting procedure to extract all the parameters of the model from a single voltage trace and the input current is a four-step method presented in *Fitting Procedure*. Here we describe in detail two critical steps of this method, *1*) the linear regression method that allows the estimation of the optimal parameters governing the dynamics of subthreshold voltage, and *2*) the maximum likelihood method used to estimate the optimal parameters governing the firing activity of the model. Along with these descriptions some possible variants of the standard fitting protocol are briefly explained.

### Estimation of the Subthreshold Voltage Parameters

As discussed in the *step 2* of the Section *Fitting Protocol*, the model parameters *V̇*_{t}^{(mod)} so that:
*V*, the time-derivative of the voltage *V̇*, the input current *I* and the basis function are known the only unknown in *Eq. 16* are the parameters. For notational convenience, we will write the above equation as a matrix equation, defining *V→̇*^{(mod)} as the vector of the binned voltage time-derivatives and **X** as a matrix with the vector of voltage *V→ _{t}* (binned as a function of time) in the first column, a vector l→ of ones in the second column, the vector of input current

*I→*(binned as a function of time) in the third column, and the value ∑

_{t}*t̂*for 1 ≤

_{j}f_{t−t̂}^{(k)}*k*≤

*K*evaluated at the known spike times

*t̂*in the remaining

*K*columns, such that the differential equation 16 becomes:

*SSE*(θ→

_{2}) = ‖

*V→̇*

_{Ψ2}−

*X*

_{Ψ2}θ→

_{2}‖

^{2}between the voltage derivative of the experimental trace

*V*

_{Ψ2}and that of the model,

*X*

_{Ψ2}θ

_{2}, summed over all points in the data set Ψ

_{2}that comprise the voltage trace in the subthreshold regime. According to multilinear regression theory (Weisberg 2005), the optimal set of parameters is then given by

*V̇*time series. This step gives the passive parameters of the neurons and the adaptation current η (Fig. 2

*B*,

*ii*,

*inset*). Note that here we voluntarily discard the spikes from the data (because we only consider the subset Ψ

_{2}, see appendix c), but it is straightforward to apply the same linear regression method on the whole recording (including the spikes). To do this, we put

*T*

_{refr}= 0 so that the first bins in η(

*t*) model the shape of the spike and an explicit reset of the voltage is no longer be necessary.

We presented a method to extract the coefficients *a*_{k} that govern the shape of the spike-triggered current η from the data. However, when a model IF_{N} is considered, the time constants τ_{i} and the amplitude *b*_{i} of the adaptation currents *w*_{i} are extracted from η(*t*) = Σ_{k}*a*_{k} *f* ^{(k)}(*t*) by fitting *N* exponential functions to the time course η(*t*). When a model IF_{ηcond} is considered, we observe that if *E*_{rev} is known a priori, then exactly the same protocol can be applied. So we perform the linear regression defined by *Eq. 18* iteratively for a set of {*E*_{rev}} and chose the optimal *Ê*_{rev} to be the one that minimizes the SSE of the regression (see Fig. 7, *A* and *B*).

### Estimation of the Voltage Threshold Parameters

As discussed in *step 3* of the *Fitting Protocol*, it is possible to extract the cumulative dynamic threshold γ from the data by maximizing the likelihood of generating the experimental spike train by our model. The log-likelihood for a spike train can be written in terms of the probability *p*_{t} of observing no spike in an experimental time bin by using Bayes theorem recursively in time (Paninski 2004):
*t̂ _{j}*} is taken to be 0.5 ms before the peak of the spike and

*Eqs. 3*,

*4*, and

*5*, the probability of having no spike in a time bin [

*t, t*+ Δ

*T*] is (Gerstner and Kistler 2002):

*T*(here Δ

*T*= 0.05 ms). The λ

_{t}is the discretized version of λ(

*t*). Using

*Eq. 4*, we have:

*X→*is made of

_{t}*V*

_{t}, one, and ∑{

*t̂*}

_{j}*f*

_{t}

^{(p)}, and the parameters are

*V*

_{t}and not the modeled voltage, as it would be done with a purely generative model.

Now, using λ_{t} defined in *Eq. 21* with λ_{0} = 1/Δ*T* and using the fact that λ_{t}Δ*T* is small, the optimal set of parameter θ→̂_{3} is simply given by:
*Eq. 21*, we are sure that the log-likelihood is a convex function of the parameters θ→_{3} (Paninski et al. 2005).

In *Fitting Procedure*, we mentioned that to extract the optimal parameters *V*_{0}, Δ*V*, and *c*_{p} governing the threshold dynamics γ (*t*), one has to maximize the log-likelihood of a spike train with a gradient ascent. To perform the gradient-ascent of the log-likelihood function, the simplest method is perhaps to use a preprogrammed script (for instance fminunc.m in Matlab), but we used the iteratively reweighted least-square method, also called Fischer's scoring method (McCullagh and Nelder 1998). All numerical computations have been done in Matlab (The Mathworks, Natwick, MA) on a desktop computer. In practice our fitting procedure is straightforward and fast; it takes only a few minutes on a desktop computer to extract all the parameters of a model from 10 s of voltage recordings and current injection producing a firing frequency of 10 Hz. The Matlab code used to extract the model parameters along with a subset of our data will be available on ModelDB: http://senselab.med.yale.edu/modeldb/.

- Copyright © 2012 the American Physiological Society