## Abstract

Neuronal response properties are typically probed by intracellular measurements of current-voltage (*I-V*) relationships during application of current or voltage steps. Here we demonstrate the measurement of a novel *I-V* curve measured while the neuron exhibits a fluctuating voltage and emits spikes. This dynamic *I-V* curve requires only a few tens of seconds of experimental time and so lends itself readily to the rapid classification of cell type, quantification of heterogeneities in cell populations, and generation of reduced analytical models. We apply this technique to layer-5 pyramidal cells and show that their dynamic *I-V* curve comprises linear and exponential components, providing experimental evidence for a recently proposed theoretical model. The approach also allows us to determine the change of neuronal response properties after a spike, millisecond by millisecond, so that postspike refractoriness of pyramidal cells can be quantified. Observations of *I-V* curves during and in absence of refractoriness are cast into a model that is used to predict both the subthreshold response and spiking activity of the neuron to novel stimuli. The predictions of the resulting model are in excellent agreement with experimental data and close to the intrinsic neuronal reproducibility to repeated stimuli.

## INTRODUCTION

Accurate models of electrically active cells and their interactions are central requirements for the understanding of the computational processes taking place in nervous tissue. The construction of network models, even at the level of cortical columns, requires the identification of cell classes and the quantification of both their typical behavior and the heterogeneities within a population. The volume of data that is required for this tissue-level modeling demands a high-throughput approach in which response properties can be routinely measured.

Electrophysiology provides an array of techniques for the extraction of neuronal response properties. Standard methods involve probing the response to step-change stimuli leading to current-voltage (*I-V*) curves for the steady-state or instantaneous response. Used systematically with pharmacology, they can yield a full conductance-based description (Hodgkin and Huxley 1952; Koch 1999) although the time required is prohibitive for routine neuron-by-neuron classification. More recently, an elegant optimization method (Huys et al. 2006) has been proposed that promises to significantly facilitate the construction of biophysically detailed models, given some prior knowledge of the kinetics of the channels present.

Detailed models, comprising hundreds of compartments, are important for understanding the biophysical properties probed during electrophysiological and pharmacological manipulations. Such models can be used for network simulations on high-performance computers, but the complexity associated with a high level of detail means that it can be difficult to get a deep understanding of the network states that emerge in these simulations. At the other end of the spectrum of neuron models are single-variable integrate-and-fire type models. The mathematical analysis of this class of models has given a great deal of insight into the emergence of dynamic network states (Abbott and van Vreeswijk 1993; Brunel and Hakim 1999; Gerstner 2000; Gerstner and van Hemmen 1993), and they are readily extended to include further biological details such as subthreshold currents (Brunel et al. 2003; Richardson et al. 2003) and adaptation (Gigante et al. 2007).

An interesting link can be made between one-variable models and *I-V* curves; both involve the reduction of multi-variable dynamics, comprising voltage and channel activation states, to a relation between the net membrane current and voltage. Standard stimuli used for measuring *I-V* curves comprise step-change current and voltage pulses (Hodgkin et al. 1952) or slow voltage ramps (see, for example, Swensen and Marder 2000). However, if a more naturalistic stimulus were used to generate an *I-V* curve, then a reduced model derived from that *I-V* curve could provide a more efficient description of the neuronal dynamics. Ultimately the success or not of such a method should be judged by the extent to which it can predict experimental data.

Here we introduce a method that measures the *I-V* curve during ongoing activity. An attractive feature is that it can be used to measure response properties in time slices triggered to events in the voltage trace—an aspect that will be used to quantify the refractory properties of pyramidal cells. A model derived from the *I-V* curve measurements will then be tested against novel stimuli, which were not used for parameter extraction, and its predictive power evaluated.

## METHODS

### Slice preparation and whole cell recordings

Parasagittal neocortical brain slices (300 μm thick) from C57/B16J mice (P12–14, *n* = 8 and P20, *n* = 4) were prepared using a vibratome (Leica VT1000 S, Leica Microsystems GmbH, Germany) in ice-cold extracellular medium (ACSF; containing in mM: for P12–14 mice and incubation chamber, 125 NaCl, 25 NaHCO_{3}, 25 glucose, 2.5 KCl, 1.25 NaH_{2}PO_{4}, 2 CaCl_{2}, and 1 MgCl_{2}; for P20 mice slice cutting, 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}). Slices were transferred to a chamber containing ACSF bubbled with 95% O_{2}-5% CO_{2} and incubated at 35°C for 30 min and then at room temperature (20–22°C) until use. All experiments were performed at 35°C. Pyramidal neurons were identified by video-enhanced infrared microscopy using the 20×0.95NA water-immersion lens of an upright microscope (BX51WI; Olympus, Tokyo, Japan). Double somatic whole cell recordings were obtained from layer-5 pyramidal neurons. One pipette injected the current while the voltage was measured simultaneously from both. However, a two-electrode protocol for the measurement of intracellular voltage during current injection is not necessary if an active electrode compensation (AEC) technique (Brette et al. 2005, 2007a,b) is used that digitally removes the artifacts introduced by the electrode filter using a Wiener-Hopf optimal filtering. A review of this method, and the particular optimization principle used here, together with the comparison between single and dual electrode recordings can be found in appendix b. The injected current waveforms were constructed from two summed Ornstein-Uhlenbeck processes (see appendix a for details) with time constants τ_{fast} = 3 ms and τ_{slow} = 10 ms. There were two sets of variances (low, σ_{fast} = σ_{slow} = 0.18 and high, σ_{fast} = 0.36, σ_{slow} = 0.25) and four sets of DC biases 0.0, 0.02, 0.03, and 0.06 in relative units, yielding a total of 8 combinations of waveforms. Each waveform had a duration of 40 s and was preceded and followed by a 3-s null (0 current) stimulus used to measure background noise levels. The set of eight waveforms was injected twice at an interval of 15 min to test for the intrinsic reliability of the cells as well as to check the stability of the cellular properties over the duration of the recordings. To produce an average firing rate in the desired range (1–15 Hz), the waveforms were multiplied by a factor (in the range 250–750 pA) to give the current injected into the neurons. All current waveforms contained square pulses at the beginning and end of experiment, allowing for standard measures of input resistance and capacitance. The 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; 270 mosmol/l). Biocytin (3 mg/ml) was routinely added to allow post hoc staining of the recorded cells. Pipette capacitance was compensated, and whole cell recordings with access resistance <20MΩ were obtained using Multiclamp 700A amplifiers (Molecular Devices, Foster City, CA). The membrane potential was filtered at 10 kHz and digitized at 20 kHz using an ITC-18 (InstruTech, Port Washington, NY). All experiments were carried out following protocols approved by the Swiss Federal Veterinary Office.

### Mathematical model

The basic model used is the exponential integrate-and-fire (EIF) model (Fourcaud-Trocmé et al. 2003), given here by *Eqs. 7* and *9*, which was also extended to include refractory properties by making the three parameters *E*_{m}, 1/τ_{m}, and *V*_{T} dependent on the time since the last spike (called the rEIF model). For example, given that the last spike occurred at time *t*_{k} the value of *E*_{m} was calculated as *E*_{m}(*t*) = *E*_{m}^{0} + *A*1*e*−(t−tk)/τ1 + *A*2*e*−(t−tk)/τ2, where *E*_{m}^{0} denotes the value of *E*_{m} obtained away from a spike. The amplitudes *A*_{1}, *A*_{2} and time constants τ_{1}, τ_{2} were obtained by fitting the set of *I-V* curves measured in time slices after each output spike (see ⇓⇓Fig. 3). For 1/τ_{m} and *V*_{T}, a single exponential term was sufficient, and a dynamics for Δ_{T} was not required.

*Equation 7* with the definition of *F*(*V*) given in *Eq. 9* was integrated using an Euler scheme with a time step corresponding to the sampling rate of the recordings (20 kHz). The action potentials of the model appear as a rapid rise in voltage and the integration was stopped when *V* reached 0 mV. Because the downswing of the spike is not explicitly described in the model, the integration was interrupted for 2 ms after detection of an action potential (typical duration of the action potential downswing) and restarted at a reset *V*_{re}, which was in general close to or above the prespike threshold *V*_{T}^{0}. The shifting of the spike initiation threshold *V*_{T} toward a more depolarized value after a spike prevents the immediate generation of a second action potential. For the nonrefractory EIF model, it was necessary to use a reset value that is sufficiently subthreshold to prevent repetitive firing. For this model, *V*_{re} = −55 mV followed by a 10-ms refractory period in all the simulations was used.

The performance index (Jolivet et al. 2004; Kistler et al. 1997) used to compare spike trains generated by the model to the experimental recordings was taken to be the coincidence factor Γ (1)

*N*_{coinc} is the number of coincidences with precision Δ = 5 ms, 〈*N*_{coinc}〉 = 2*f*Δ*N*_{neuron} is the number of expected chance coincidences generated by a Poisson process with the same firing rate *f* as the neuron, *N*_{neuron} and *N*_{model} are the number of spikes in the spike trains of the neuron and the model, respectively (𝒩 is a normalization factor that is of no importance in the following as only ratios of Γ are considered). In Fig. 4*B*, the ratio Γ/Γ′ is plotted, where Γ evaluates the overlap between the prediction of the model and a target experimental spike train, and Γ′ is calculated between the target spike train and a second experimental recording obtained with the same driving current. Only pairs of trials with an experimental reliability Γ′ > 0.75 between recording sessions spaced by 15 min were used. Setting the threshold at this level still included neurons that fire at low rates (around 1 Hz) but excluded those in which a significant drift in cellular properties (significantly increasing or decreasing voltage baselines or firing rates) was seen over the time between the repeated measurements.

## RESULTS

To produce fluctuating voltage traces, a stochastic current was injected into somatosensory pyramidal cells via a whole cell somatic patch-clamp electrode. Basic electrophysiology provides an equation that relates the capacitive charging current and the summed effect *I*_{m}(*V*, *t*) of the transmembrane currents to this injected current *I*_{in}(*t*) (2)

In this formulation, *I*_{m} comprises both transmembrane currents and equilibrating currents flowing between the soma (point of current injection) and the dendrites and axon. The noise term *I*_{noise} comprises effects from weak background synaptic activity and other sources of high-frequency variability. *Equation 2* may be re-arranged to yield the dynamic transmembrane current (3)

The injected current waveform *I*_{in} is known a priori and the derivative d*V*/d*t* may be calculated from the experimentally measured voltage. Therefore once the capacitance *C* has been determined (conveniently achieved by examining the variance of *I*_{m}; see following text), the transmembrane current *I*_{m} can be calculated as a function of time *t* (Fig. 1 *A*).

### Dynamic I-V curve

Both *I*_{m}(*t*) and *V*(*t*) are parameterized by time, whereas for an *I-V* curve, a direct relation between the two is required. A scatter plot of *V*(*t*) on the *x* axis versus *I*_{m}(*t*) on the *y* axis is given in Fig. 1*B* where only data in the subthreshold voltage range and run up to the spike are included at this stage (by excluding all data that fall in a 200-ms window following each spike peak). For a given voltage, a Gaussian distribution is seen for the values of the current (Fig. 1*B*, *inset*), which is largely due to the intrinsic noise *I*_{noise} (83% of the SD). The average current, as a function of voltage, defines the dynamic *I-V* curve *I*_{d}(*V*) (4)

Practically, this quantity may be calculated by collecting all the points in the trajectory with a voltage in the range *V* ± Δ/2, where Δ is a bin of sufficiently small width. The dynamic *I-V* curve for a layer-5 pyramidal cell is given in Fig. 1*B* (red squares) and comprises linear and strong negative components.

### Determination of the capacitance

The calculation of the membrane current via *Eq. 3* requires the capacitance *C* to be known. Using an estimate *C*_{e} for the capacitance, this equation can be cast in the form (5)

The variance (see Fig. 1*B*, *inset*) at a voltage *V* takes the form (6)

The right-hand side of *Eq. 6* is minimized when *C*_{e} = *C* (Fig. 1, *C* and *D*), yielding the correct cellular capacitance. In calculating *Eq. 6*, the covariance of *I*_{in} and *I*_{m} at a voltage *V* is assumed to vanish, which is justifiable for the variance calculated over the subthreshold *I-V* curve. This approach yields a capacitance that is consistent with standard methods (confirmed by responses to square-pulse currents) and is justified by the low variability around the mean identified in the previous section.

### One-variable neuron models

An important family of neuronal models in widespread use (for recent reviews, see Burkitt 2006; Gerstner and Kistler 2002; Lindner et al. 2004) has a voltage dynamics of the form (7) where in general *F*(*V*) is a nonlinear function of voltage. These dynamics are supplemented by a threshold *V*_{th} above which an action potential is registered, followed by a reset to a subthreshold voltage *V*_{re}. In the absence of injected current, the rate of change of voltage is equal to the function *F*(*V*), so its zeros can be used to identify stable fixed points of the voltage such as the resting potential, and unstable points such as the spike-initiation threshold.

Different models for *F*(*V*) have been proposed with the most common choice being the linear, leaky integrate-and-fire model for which *F*(*V*) = (*E*_{m} − *V*)/τ_{m} where *E*_{m} is the resting potential and τ_{m} the membrane time constant. Nonlinear forms for *F*(*V*) have also been proposed that explicitly describe spike initiation, such as the quadratic integrate-and-fire model (Gutkin and Ermentrout 1998; Izhikevich 2004; Latham et al. 2000), which is also the canonical form of type I models (Ermentrout 1996; Ermentrout and Kopell 1986). An alternative proposal for the nonlinearity is the EIF model (Fourcaud-Trocmé et al. 2003) that includes the activation of the spike-generating sodium current. It has been shown (Fourcaud-Trocmé and Brunel 2005) that the functional form used to model the spike-initiation is a crucial determinant of the neuronal response to rapid synaptic signaling.

The dynamic *I-V* curve provides an experimental method for measuring *F*(*V*) directly. On comparing *Eqs. 2* and *7*, a natural interpretation is (8) where *I*_{d} was defined in *Eq. 4*. With this choice, the intrinsic current of the model will be the same as in the experimental data on average. Although the temporal effects of voltage-gated currents are averaged out in this one-variable (somatic voltage) model, it was shown in the previous section that after the intrinsic variability is accounted for, the unexplained variability around the mean is low, so it can be expected that *Eq. 8* would have the potential to provide an accurate description of neuronal dynamics.

In Fig. 2 *A*, the function *F*(*V*) (expansion of data in Fig. 1*B*) measured from a pyramidal cell is plotted (symbols). There are two fixed points where *F*(*V*) = 0 for the resting potential and spike-initiation threshold, a linear region where the cellular response is passive and a rapid nonlinear component at higher voltages. A semi-log plot (see *inset*) shows that the nonlinear component is exponential (*y* axis), and therefore these measurements provide convincing empirical evidence for the EIF model (Fourcaud-Trocmé et al. 2003) for pyramidal cells for which (9)

The EIF model has four parameters: the membrane time constant τ_{m}, the resting potential *E*_{m}, the spike-initiation threshold *V*_{T}, and the spike width Δ_{T}, which controls the sharpness of the initial phase of the spike. A least-squares fit of the full form of *Eq. 9* to the pyramidal-cell *I-V* curve data (Fig. 2*A*, solid line) is seen to be in excellent agreement with the pyramidal-cell dynamic *I-V* curve. It should be noted that the spike initiation for pyramidals is sharper than that of the hippocampal interneuron model (Wang and Buzsáki 1996) used in the original EIF paper (Fourcaud-Trocmé et al. 2003) in common with the “kink”-like spikes identified in recent analyzes of pyramidal-cell voltage recordings (McCormick et al. 2007; Naundorf et al. 2006).

### Variation across the population

By comparing the results from different recordings (*n* = 12), we found that the qualitative shape of the dynamic *I-V* curve was remarkably stable across layer-5 pyramidals with the EIF model always providing an excellent match. After fitting the model to data on a neuron-by-neuron basis, we examined the histograms of model parameters (Fig. 2*B*). The coefficients of variation (CV) for the membrane time constant, distance to threshold and spike width were 32, 25, and 33%, respectively, demonstrating a significant degree of heterogeneity in the population, with implications for network modeling.

### Postspike curve

The dynamic *I-V* method can also be applied to transient changes in response properties, such as the postspike recovery period excluded from the treatment illustrated in Figs. 1 and 2. To this end, we separated the postspike data into the time slices shown in Fig. 3*A*, and for each time slice, the dynamic *I-V* curve was again calculated. The resulting postspike *I-V* curves are quantitatively different and relax to the prespike form over a period of many tens of milliseconds. The EIF functional form (*Eq. 9*) also provided a good fit to the postspike data, which was unexpected, allowing for the refractory properties to be quantified in terms of the dynamic effects on the EIF parameters (Fig. 3*B*). The dynamics of the parameters could be fitted by a single exponential term except for the resting potential *E*_{m} (requiring a biexponential form), which showed a prominent sag indicative of the transient activation of a hyperpolarizing current or of the deactivation of a depolarizing current, such as *I*_{h}.

### Verification of the models

These pre- and postspike response properties were used to construct two models which were tested against further experiments. The first model comprises *Eqs. 7* and *9* with a “hard” subthreshold postspike reset and is identical to the EIF model (Fourcaud-Trocmé et al. 2003). The second is a novel refractory EIF model (rEIF) for which the “soft” refractory period relaxes with the dynamics identified in Fig. 3 (see methods).

For each cell investigated, a range of test currents, with different means and variances, were injected to produce voltage traces with firing rates in the range of 1–15 Hz. A unique set of model parameters was extracted individually for each cell via a fit to the dynamic *I-V* curve from a single 40-s voltage trace of intermediate firing rate (for the nonrefractory model, the voltage data in a 200-ms window following each spike peak were excluded from the analysis). The parameters of this model were then fixed to be tested on the other current waveforms and compared with the response from the same cell. Each test current was injected into the cell twice so that the intrinsic reliability of the cell could also be quantified. The measures of the quality of prediction of the rEIF and EIF models and the cellular reliability, comprised: the average firing rate, the number of predicted spikes (within ±5 ms of the spike peak) and the accuracy of the subthreshold voltage. For the latter, the root-mean-square (RMS) difference between the model voltage and the experimental voltage was examined using data at least ±50 ms away from spikes. For the voltage distributions, only data between 0.5 ms before and 3 ms after the spike were excluded. The results are plotted in Fig. 4.

Both models predict the positions of isolated spikes accurately and have a subthreshold voltage in close agreement with the experimental data. However, the hard reset of the EIF is a poor model of the postspike behavior of layer-5 pyramidal cells. In Fig. 4*A*, *bottom*, the EIF is seen to miss a spike in a closely spaced pair (*left*) and to add a spurious spike due to an underestimation of the threshold (*right*). At these firing rates, such effects were relatively rare and both models matched the firing rates well (Fig. 4, *B* and *C*) although with the EIF predicting a weak excess of spikes at higher rates. The positions of the spikes were also accurately predicted (see methods) where for the rEIF in Fig. 4*D*, it is seen that the model rarely predicts <75% of the spikes correctly with a mean ± SD of 83 ± 8% and where the EIF model has a mean ± SD of 72 ± 14% spikes predicted relative to the intrinsic reliability. Finally, the typical RMS difference of the subthreshold voltage for the rEIF model was 1.7 mV and for the EIF 2.0 mV, both of which again compare favorably with the RMS for repeated stimuli at 1.4 mV as can be seen in Fig. 4*E*, as well as in *F*, which compares the voltage distributions themselves.

### Efficiency of the method

The results presented so far were obtained from 40-s voltage traces. To compare the results against those from shorter recordings, we compared the *I-V* curves obtained by using only the first 10 s of data with those obtained using the whole 40-s trace. We found little difference between the model fits for the two recording durations (Fig. 5). It can be noted that the limiting factor for a clean fit is the collection of a sufficient number of spikes (in this case around 50 for the 10 s measurement) during the recording.

## DISCUSSION

We have demonstrated a novel method for measuring an *I-V* curve that characterizes the neuronal response under conditions of ongoing, naturalistic voltage fluctuations. For cortical pyramidal cells the prespike *I-V* curve was very well fitted by the EIF model (Fourcaud-Trocmé et al. 2003). The distribution of fitting parameters for a population of cells was also given with CV values demonstrating a significant degree of heterogeneity. The dynamic *I-V* curve was also measured in postspike time slices and the refractory properties quantified. An unexpected result was that the refractory *I-V* curves were still well fitted by the EIF functional form, allowing for the derivation of a refractory extension of the basic EIF model (rEIF).

Both the EIF and rEIF models were tested against further experiments and shown to predict spike times accurately with the rEIF showing particularly close agreement with the experimental data. In comparison with the neuronal response to repeated stimuli, it was shown that the predictive ability of models generated from the dynamic *I-V* method is largely limited by the intrinsic reliability of the neuron.

The method is highly efficient, requiring a few tens of seconds of experimental time. The low-noise data gathered here (see Fig. 2*A*) required 40-s voltage traces. However, if the neuron is firing at around a few hertz, ∼10–20 s is more than enough to categorize the pre- and postspike response properties (Fig. 5). This efficiency makes the method suitable for the high-throughput routine neuronal classification and model fitting to experimental data that is required for the generation of network models that include cellular heterogeneities.

### Applications and extensions

A number of elaborations of the method suggest themselves, specifically with a view to analyze the role of subthreshold voltage-activated currents. Apart from the postspike sag response of *E*_{m} (Fig. 3*B*), such effects were not explicit in the data: the subthreshold *I-V* curve was linear (Fig. 2*A*), and the additional dynamics from any voltage-activated channels were not required to accurately fit experimental data for traces with firing rates in the range 1–15 Hz. This is not entirely unexpected as the rapidly fluctuating voltage tonically activates voltage-gated channels with slower dynamics, and so their effect will largely be seen as a change of the total membrane conductance. The strength of this effect should vary with the different means and variances of the input current, though this was not significant over the range of currents used here. Nevertheless, more refined measurements that could combine fluctuating current with ramps (Swensen and Marder 2000) might be used to examine the potential effects of the h-current, persistent sodium current, or other currents that might be present. As a related point, it would be simple to extend the method to the case of fluctuating conductance injection (Destexhe et al. 2001). Although both current and conductance fluctuations lead to very similar voltage distributions (Burkitt 2001; Richardson 2004; Richardson and Gerstner 2005), the frequencies present in the spectrum of the fluctuations increases with the level of conductance (Destexhe and Rudolph 2004). The faster temporal structure of conductance-based input could potentially interact differently with the dynamics of intrinsic currents to alter the measured response properties. Such an effect might be measurable in neurons with subthreshold currents with fast dynamics (with time scales similar to the membrane time constant), but for currents such as *I*_{h} (with relatively slow dynamics), this effect could be expected to be weak. Finally, coupled to dendritic recordings or by stimulating presynaptic units simultaneously, the dynamic *I-V* curve may also be employed to study how the activation of dendritic conductances affects the processing of information in the soma.

### Implications for neural modeling

During the last couple of years, there has been much debate in the theoretical neuroscience community concerning the appropriate minimal description of subthreshold response and spike generation. In particular, the choice of a strict voltage threshold, as in the classical integrate-and-fire model (Lapicque 1907; Stein 1967), has been critically analyzed (Fourcaud-Trocmé et al. 2003), and various generalizations have been proposed, including canonical phase models (Ermentrout 1996; Ermentrout and Kopell 1986; Gutkin and Ermentrout 1998), quadratic, (Latham et al. 2000) and exponential (Fourcaud-Trocmé et al. 2003) integrate-and-fire models. Refractory properties have also been studied in several forms (Jolivet et al. 2004; Kistler et al. 1997; Lindner and Longtin 2005), and the inclusion of a second, slow variable has proven effective in reproducing more complex subthreshold dynamics such as resonance (Brunel et al. 2003; Izhikevich 2004; Richardson et al. 2003) and adaptation (Brette and Gerstner 2005; Clopath et al. 2007; Gigante et al. 2007; Jolivet et al. 2004). The dynamical *I-V* curve data for cortical pyramidals, and associated modeling, presented here provides clear empirical evidence for the EIF model (Fourcaud-Trocmé et al. 2003), which, even with a hard reset, was able to give an accurate prediction of the spike times. It is interesting to note that the measured spike width Δ_{T} for pyramidal cells is considerably sharper than that predicted from the original fits to an interneuron model with Hodgkin-Huxley spike-generating currents (Wang and Buzsáki 1996). Given that the leaky integrate-and-fire (IF) model is recovered from the EIF model in the limit Δ_{T} → 0, this empirical result goes some way in explaining the success of leaky IF models in predicting firing rates (Rauch et al. 2003) and spike times (Jolivet et al. 2006; Paninski et al. 2004).

Our findings also clearly highlight the importance of refractory properties for the correct modeling of cortical pyramidals. The measured refractory properties comprised conductance increase and, most importantly, a significantly increased threshold for spike initiation allowing for a postspike reset that is above the prespike threshold. Although the idea of exponentially relaxing refractory properties is far from new (Fuortes and Mantegazzini 1962; Geisler and Goldberg 1966), it has until recently (Chacron et al. 2007; Gerstner and van Hemmen 1993; Lindner and Longtin 2005) received little analytical attention. Although the refractory dynamics of τ_{m}, *E*_{m}, and *V*_{T} were used for the rEIF model (Fig. 4), a two-variable model comprising only voltage and an exponential relaxation for the threshold *V*_{T} captures the essential features of the pyramidal refractory period and still provides an accurate description of the postspike dynamics (model not shown); expending further analytical effort to understand this effect would be worthwhile.

Finally, as part of our analyses, the distribution of fitting parameters was measured across a sample population of cells. These revealed a degree of heterogeneity that could well have implications for the modeling of network stability and the sharpness of transitions between dynamic states, such as rhythmogenesis. The majority of analyzes of recurrent networks of spiking neurons has been performed on networks of homogeneous populations. The experimentally measured distributions provided here can be used in future studies that investigate the effects of heterogeneity on emergent states in neuronal tissue.

## APPENDIX A: ORNSTEIN-UHLENBECK PROCESS

The Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein 1930) is the Gaussian stochastic process defined by the equation (A1)

Here, ξ is a Gaussian white noise process with the properties 〈ξ(*t*)〉 = 0, 〈ξ(*t*)ξ(*t*′)〉 = δ(*t* − *t*′), where δ is the Dirac δ function. The process (*Eq. A1*) has stationary mean 〈*x*(*t*)〉 = μ, variance 〈*x*^{2}(*t*)〉 − 〈*x*(*t*)〉^{2} = σ^{2}, and exponential autocorrelation (A2)

To produce a numerical realization of the Ornstein-Uhlenbeck process *x*(*t*), we discretize the time interval into equally spaced points *t _{k}* =

*k*Δ

*t*, where 1/Δ

*t*is the desired sampling frequency. By integrating

*Eq. A1*, the change in the value of

*x*(

*t*) between

*t*and

_{k}*t*

_{k}_{+1}is found to be (A3) where Ψ

_{k}∼ 𝒩(0, 1) are independent Gaussian random numbers with zero mean and unit variance. For our recordings at 20 kHz, we generated the currents by applying

*Eq. A3*iteratively with a time step Δ

*t*= 0.05 ms. To produce a single current template, this was performed twice (with different random-number seeds) with the fast and slow time constants τ

_{fast}and τ

_{slow}and their amplitudes σ

_{fast}and σ

_{slow}, and the two processes summed together with an additional DC bias added to the current as required. This procedure was carried out over the set of biases and variances (each time with a different random-number seed) to yield the eight waveforms described in methods.

## APPENDIX B: SINGLE-ELECTRODE RECORDINGS

The faithful extraction of cell response properties requires accurate measurements of the membrane voltage. However, when injecting current and monitoring voltage simultaneously from a single electrode, the filtering properties of the electrode may affect the measurements to a considerable extent. Although standard techniques are available to compensate electrode capacitance and access resistance, in the case of strong, high-frequency stimuli, this correction is not sufficient (Fig. 6), and some additional compensation is required.

Recently an active electrode compensation (AEC) technique was introduced (Brette et al. 2005, 2007a,b) that solves this problem, allowing for the simultaneous injection of current and measurement of voltage using a single electrode. To verify this compensation technique for the dynamic *I-V* curve protocol, we performed recordings in double-patch configuration with one electrode (referred to as electrode 1) simultaneously injecting current and recording voltage, while the other electrode (control electrode 2) was used only to monitor the voltage. In general, the uncompensated voltage measured from electrode 1 showed a significant deviation from the control (Fig. 6). Following the AEC methodology, we will now explain how these artifacts can be removed from the voltage trace of electrode 1 through the estimation of the electrode filter without knowledge of the voltage from the control electrode, so as to allow for the measurement of dynamic *I-V* curves from single-electrode recordings.

### Electrode filter

We assume that the voltage across the electrode is a linearly filtered version of the input current *I*(*t*). The recorded voltage *V*^{rec} is the sum of the true membrane voltage *V* and the electrode voltage *V*^{el}, the latter of which we write as a convolution integral of the current with the (unknown) electrode filter *f*(*s*) (B1)

On the other hand, for a linear membrane the voltage can also be written in terms of a filter *g*(*s*) such that *V*(*t*) = *V*_{0} + ∫_{0}^{∞} *g*(*s*) *I*(*t* − *s*)d*s*. Therefore we cannot expect to determine the electrode filter directly but rather the combined filter *f̃*(*s*) = *f*(*s*) + *g*(*s*) of the system consisting of the electrode *and* the neuron. It is this filter that will be determined in the next section.

### Optimization procedure

We have seen in the main text that the membrane capacitance can be estimated from a variance minimization procedure. In our single-electrode approach, the variance must be calculated using the true, corrected membrane voltage, *V*(*t*) = *V*^{rec}(*t*) − ∫_{0}^{∞} *f*(*s*)*I*(*t* − *s*)ds. In analogy with *Eq. 6*, *left*, of the main text, we therefore minimize the function (B2) with respect to both the filter *f*(*s*) and the capacitance *C*_{e}, where the dot denotes the time derivative. Practically, the discrete form of this equation is used (B3)

Here, *V*_{i}^{rec} and *I _{i}* are the uncompensated voltage and current data at time step

*i*, and the derivatives are calculated as forward differences,

*V̇*

_{i}

^{rec}= (

*V*

_{i+1}

^{rec}−

*V*

_{i}

^{rec})/Δ

*t*, and

*İ*= (

_{i}*I*

_{i}_{+1}−

*I*)/Δ

_{i}*t*, where Δ

*t*denotes the time step of the recording (1/Δ

*t*= ν is the sampling frequency in kHz).

*f*is the

_{k}*k*th value of the filter and represents how much of the recorded voltage at a given time is due to the current that was injected

*k*time steps earlier.

*M*is the total number of points in the filter (a way to efficiently choose this parameter is discussed in the following text). The subscript

*i*∈ Ω indicates that the variance is calculated over the set Ω of those data points that are both (

*i*) close to the resting potential (|

*V*

^{rec}−

*V*

_{rest}| <0.5 mV), so that we can expect the membrane to behave linearly, and (

*ii*) lie > 200 ms after the peak of a spike, so that any refractory effects that would alter the membrane response are avoided. We emphasize that the cost function in

*Eq. B3*differs slightly from that in the standard AEC method (Brette et al. 2007a,b) because it acts on the derivatives of the voltage, making it independent of a potential voltage offset.

To minimize *Eq. 6*, *right*, we first introduce the notations (B4) (B5) (B6) for functions *x*(*t*) and *y*(*t*) that have been discretized in time and where the normalization factor *N* is the number of elements in the sums (i.e., the number of data points where *V*^{rec} is within 0.5 mV of the resting potential). It should be noted that in the following expressions, the indices *j* and *k* range from 0 to *M*, so that *i* − *j* or *i* − *k* can effectively be negative. This can be dealt with by padding the voltage data with a constant value for negative times. The minimization proceeds in two steps: we first cancel the derivative of *Eq. B3* with respect to *C*_{e}, which gives an expression for *C*_{e} as a function of the *f _{k}* (B7)

Then this expression is inserted back into *Eq. B3*, which is then minimized over the filter variables *f _{k}*. It is convenient to write the resulting filter (which describes the linear response of the combined electrode and neuron system) as a single vector, f⃗ = (

*f*

_{1},

*f*

_{2}, … ,

*f*), so that the result of the minimization can be written in compact form as (B8) where

_{M}*b⃗*is a vector of covariance between the uncorrected membrane current and the derivative of the injected current, and

*A*is a matrix of covariance involving only the injected current (

*A*

^{−}

^{1}denotes the matrix inverse to

*A*). Using the

*Eqs. B4*—

*B6*, the quantities

*A*and b⃗ are given by (B9) (B10) where the subscript

*F*in

*Eq. B10*denotes the uncorrected membrane current,

*F*=

*V̇*

^{rec}− (σ

_{V̇rec,I}/σ

_{I,I})

*I*.

In practice, the value of the cost function (*Eq. B3*) at the minimum diminishes with increasing filter length *M* but increases again for large *M* due to numerical errors; thus there is an optimal value for *M* that can be found by standard optimization methods. We observed that in general, a value of *M* such that the total length of the filter is approximately one membrane time constant (10–20 ms) provides a good choice.

### Correction for the membrane filtering

The optimal filter derived in the previous section contains a slow component that is due to the membrane filtering. This component must be identified and removed from the overall filter to obtain the correct electrode filter. As can be seen in Fig. 6, the filter comprises a fast, high-amplitude rise and fall followed by a slow relaxation. To decouple the component of the electrode from that of the membrane, the slow relaxation is fitted with the sum of two exponentials, one with a fast time constant (of the order of 0.2–0.5 ms) that is due to the electrode, the other with a slower time constant (of the order of tens of milliseconds) that is due to the membrane. Then the slow exponential is subtracted from the combined filter to yield the electrode filter (Fig. 6). Using this electrode filter, the component of the voltage due to the electrode can be calculated and subtracted from the voltage trace to obtain the correct membrane voltage (Fig. 6, *right*). The entire method described in the main text can then be applied to the corrected voltage trace. Finally, it can be noted that a more sophisticated electrode compensation method can also be used (Brette et al. 2007a,b) based on a deconvolution procedure.

## GRANTS

C.C.H. Petersen, W. Gerstner, and L. Badel acknowledge support from the Swiss National Science foundation. R. Brette acknowledges support from a French Agence Nationale de la Recherche grant (HR-CORTEX). M.J.E. Richardson was partially supported by a European grant Fast Analog Computing with Emergent Transient States during his stay in Lausanne and currently holds a Research Councils United Kingdom Academic Fellowship.

## DISCLOSURE

The authors state that there are no conflicts of interest.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society