We introduce a two-dimensional integrate-and-fire model that combines an exponential spike mechanism with an adaptation equation, based on recent theoretical findings. We describe a systematic method to estimate its parameters with simple electrophysiological protocols (current-clamp injection of pulses and ramps) and apply it to a detailed conductance-based model of a regular spiking neuron. Our simple model predicts correctly the timing of 96% of the spikes (±2 ms) of the detailed model in response to injection of noisy synaptic conductances. The model is especially reliable in high-conductance states, typical of cortical activity in vivo, in which intrinsic conductances were found to have a reduced role in shaping spike trains. These results are promising because this simple model has enough expressive power to reproduce qualitatively several electrophysiological classes described in vitro.
While detailed conductance-based neuron models with tens of ion channels of the Hodgkin-Huxley type (Hille 2001; Hodgkin and Huxley 1952) can, in principle, capture the electrophysiological behavior of neurons in great detail, the tuning of up to a hundred or more parameters to specific neuron types has so far been limited to a few case studies (Bower and Beeman 1997; Roth and Häusser 2001; Segev et al. 1989), and systematic procedures for automatic extraction of parameters from experimental data are unknown. Besides, it has been shown that spike patterns remain unchanged along certain directions of the parameter space of detailed conductance-based models (Goldman et al. 2001; Prinz et al. 2003), which suggests that these models could be simplified without drastic consequences for the spiking behavior. Simplified phenomenological neuron models of the integrate-and-fire type are certainly less exact from a biophysical point of view, but because they have fewer parameters, automatic procedures for parameter extraction from a limited amount of experimental data become feasible (Jolivet et al. 2004; Keat et al. 2001; Paninski et al. 2004). A traditional leaky integrate-and-fire model (Lapicque 1907; Stein 1967; Tuckwell 1988), which combines linear filtering of input currents with a strict voltage threshold, is widely used as a simple model of a spiking neuron, but it has some obvious limitations. Therefore generalizations have been proposed that go in at least three directions. First, an extension to quadratic (Ermentrout 1996; Latham et al. 2000) or exponential (Fourcaud-Trocme et al. 2003) integrate-and-fire neurons allows replacement of the strict voltage threshold by a more realistic smooth spike initiation zone. Second, the addition of a second variable allows inclusion of subthreshold resonances or adaptation (Izhikevich 2003; Richardson et al. 2003). Third, a change in the stimulation paradigm from current injection to conductance injection allows moving the integrate-and-fire models closer to a situation that cortical neurons would experience in vivo (Destexhe et al. 2003).
In this report, we will introduce an adaptive exponential integrate-and-fire (aEIF) model that combines the three extensions mentioned above and shows that all model parameters can be systematically extracted from a series of standard stimulation paradigms. To show the feasibility of the approach, the method in this paper was applied to artificial data generated from a detailed conductance-based model of a regular spiking neuron.
Detailed neuron model
The reference model was a single-compartment model of a regular spiking pyramidal cell with voltage-dependent currents IM, INa, and IK (McCormick et al. 1993), with the parameter values used in Destexhe et al. (1998) (code for the Neuron simulator available online at http://senselab.med.yale.edu/senselab/modeldb/ShowModel.asp?model=3817). It includes spike-related conductances and adaptation; the model consists of five differential equations and comprises 31 independent parameters. In vivo–like synaptic activity was modeled by two fluctuating synaptic conductances, one excitatory (AMPA, reversal potential Ee = 0 mV) and one inhibitory (GABAA, reversal potential Ei = −75 mV), described by Ornstein-Uhlenbeck processes with respective time constants τe = 2.728 ms and τi = 10.49 ms (Destexhe et al. 2001). The simulations consisted of 15 stimulations lasting 20 s with various parameter values (Table 1), classified in three groups according to the level of total conductance (defined as the sum of synaptic and leak conductances): high conductance (HC; ratio total conductance to leak conductance 5:1), medium conductance (MC; 3:1), and low conductance (LC; 2:1). The conductance ratio of the LC state corresponds roughly to that in a quiet state in vivo, whereas the HC level corresponds to high synaptic activity, following the observations of Pare et al. (1998). The detailed model was simulated using NEURON simulation software (Hines and Carnevale 1997). Each of the 15 simulations was run twice with different seeds for the random number generator, one run being used for estimating the parameters and another one for testing the performance of the model. In additional simulations, we also used current injection and applied either slow ramps or short pulses.
Adapting the aEIF model
We consider an integrate-and-fire model with adaptation defined by (1) where C is the membrane capacitance, f(V) is a function that characterizes the passive properties and the spiking mechanism, w is an adaptation variable, and I is the synaptic current. Following Fourcaud-Trocme et al. (2003), we take for f(V) a combination of linear and exponential functions (2) where gL is the leak conductance, EL is the resting potential, ΔT is the slope factor, and VT is the threshold potential (Fig 1A). Formally, a model such as the quadratic or exponential integrate-and-fire model is said to generate a spike if the potential V grows rapidly toward infinity. In practice, a spike in our aEIF model is triggered when the voltage reaches a threshold Vpeak = 20 mV, but the exact value is not critical because it only shifts spike times by a fraction of millisecond. After a spike has been triggered, integration of the equation is restarted from a reset value Vr, with Vr = EL. The slope factor determines the sharpness of the threshold: in the limit ΔT → 0, the model becomes a standard integrate-and-fire model (Lapicque 1907; Stein 1967) with threshold VT.
The adaptation current w is defined by (3) where τw is the time constant and a represents the level of subthreshold adaptation. At each firing time, the variable w is increased by an amount b, which accounts for spike-triggered adaptation. Izhikevich (2003) showed that when the quadratic integrate-and-fire model is augmented by this adaptation current, it can be tuned to reproduce qualitatively all major classes of neurons, as defined electrophysiologically in vitro (e.g., regular spiking, intrinsic bursting, fast spiking). For example, high reset values (Vr > VT) induce bursting; large values of b give strong spike-frequency adaptation; high values of the parameter a yield subthreshold oscillations, and medium values cause overshoots in response to current pulses.
Simulation of the aEIF model and analysis of the results were done with MATLAB (The Mathworks, Natick, MA) on a portable PC with a 2.4-GHz Intel processor.
Parameter values for the aEIF model were extracted from data generated by the detailed model using a series of standard electrophysiological paradigms (injection of current pulses and current ramps as discussed below and random conductance injection as described above). The resulting values are shown in Table 1.
The passive membrane properties (C, gL, EL) were obtained from an exponential fit to the response of the detailed model to a current pulse (0.1 nA for 100 ms). To determine the value of parameter a in the second equation (subthreshold adaptation), we observe that when the potential V is fixed, the adaptation current w approaches a(V−EL). Hence far away from the threshold (i.e., V < VT − ΔT, so that the exponential term in Eq. 2 can be neglected), voltage and current satisfy the following linear relationship We injected a very slow current ramp (0.01 nA/s—one could alternatively use a voltage clamp) to estimate the I–V relationship between −70 and −53 mV (Fig 1B, solid line). This range was chosen to cover typical polarization levels encountered during stimulations of biological neurons. We found that the I–V curve was approximately linear across this range—an inflexion was seen if the range was extended above. We obtained the value of parameter a from the slope of the best linear fit to the I–V curve after subtraction of the slope gL expected from the passive properties.
To determine the value of b (spike-triggered adaptation) and the adaptation time constant τw, we depolarized the membrane potential to −60 mV (using constant current injection) to get close to the average potential during synaptic stimulation (Destexhe et al. 2001), and then we injected a periodic series of short current pulses (2 nA for 5 ms at frequency 5, 10, and 20 Hz), large enough to trigger spikes (Fig 1C). Considering that we are far away from threshold, we can deduce the level of adaptation just before pulse onset from the speed of membrane depolarization dV/dt according to the following formula The difference between this estimation of w and the contribution from subthreshold adaptation estimated earlier gives the contribution from spike-triggered adaptation. Fitting this series of values to an exponential gave an estimation of b and τw. The value of the time constant depended a little on the stimulating frequency (163, 144, and 120 ms at 5, 10, and 20 Hz, respectively), but the estimation of b was very robust (81, 80, and 78 pA at 5, 10, and 20 Hz, respectively). As the final value of these two parameters, we took b = 80.5 pA and τw = 144 ms and kept these values fixed thereafter.
Last, we determined the threshold parameters VT and ΔT in Eq. 2 by the following procedure. For a given value of the slope factor ΔT, we calculated for each of the input scenarios the value of the threshold VT so that the aEIF model and the detailed one have the same average firing rate (Fig 1D). We called this value the effective threshold. For zero slope factor, i.e., for the standard integrate-and-fire (IF) model, the effective threshold varied a lot across the 15 input scenarios (0.35 mV2); thus the standard IF model with a fixed threshold was unable to predict correctly the firing rate of the detailed model. Surprisingly, the effective threshold was very stable across different stimulation regimens (Fig 1D) if a slope factor of ΔT = 2 mV was chosen (variance of 0.02 mV2 around a mean of VT = −50.4 mV). We therefore chose VT = −50.4 mV and ΔT = 2 mV as the optimal values and kept these parameters fixed throughout the results section. We also compared these values with the ones obtained by an alternative method (perhaps less applicable with real noisy recordings): VT was estimated as the potential of the inflexion point in response to a constant current that triggered spikes, whereas the slope factor was estimated by comparing VT and VS, where VS is the larger solution to f(VS) = 0 (Fig 1A), obtained as the largest membrane potential that the detailed model could reach without spiking in response to short current pulses. The resulting values were consistent with the previous ones (VT = −50.7 mV, ΔT = 2.2 mV).
We compared the spike trains of the aEIF model to the ones of the reference neuron model in terms of percentage of missing spikes M (relative to the number of spikes in the reference model) and percentage of extra spikes E (relative to the number of spikes in the aEIF model). Both measures take values in the interval 0–100%. We considered that two spikes match if they lie within 2 ms of each other. For greater convenience and comparison with previously published results (Jolivet et al. 2004), we also used a single performance measure, the coincidence factor Γ, defined in Kistler et al. (1997), ranging from 0 to 1, which is, for small E and M, related to the previous two measures by Γ = 1 − (E + M)/2. Thus the coincidence factor takes extra and missing spikes equally into account. Note that we only want to match the spike trains, not the subthreshold voltage traces.
Performance of the aEIF model
After the aEIF model was calibrated as described in methods, we compared the voltage traces of the aEIF model with those of the reference neuron model using an identical stimulation for both models. During random conductance injection, the voltage traces of the two different models were nearly indistinguishable (Fig 2A), and most spikes occurred with the same timing (±2 ms). Averaged across the 15 different random conductance paradigms (see methods and Table 1 for details of the paradigms), the aEIF model emitted 3% extra spikes and missed 4% of the spikes of the reference model (Fig. 2B), which yields a coincidence factor of Γ = 0.96. A quantitative comparison showed that the degree of similarity between the two spike trains depended only weakly on the characteristics (LC, MC, or HC) of the conductance injection paradigm (Γ = 0.95 for LC; 0.95 for MC; and 0.96 for HC; Fig 2B). Under step current injection (Fig 2C), the voltage traces are indistinguishable in the subthreshold regimen, and both traces show overshoots caused by adaptation. For superthreshold step currents, the timing of the first spike is identical, but because of a slight mismatch in the mean firing rate after adaptation, the two spike trains drift apart later on.
In a second series of experiments, we checked the influence of several components of our model. When the exponential spike mechanism was replaced by the sharp threshold of the standard integrate-and-fire model, i.e., with f(V) = −gL(V − EL) and VT = −47.1 mV (average effective threshold at slope factor ΔT = 0 mV, Fig. 1D), the performance was significantly impaired: the model fired on average 12% extra spikes and missed 12% of the spikes of the detailed model (Γ = 0.87). In detail, the performance depended on the conductance level: the model fired more extra spikes in HC states and missed more spikes in LC states (Fig 2D). When we used the effective threshold (optimized for each input scenario) instead of the fixed average threshold, the performance was balanced but not significantly better (11% extra and missing spikes; Γ = 0.88). Thus the exponential spike mechanism brought a great improvement in prediction performance (Γ = 0.88 to Γ = 0.96) beyond the notion of effective threshold, which confirms earlier theoretical studies (Fourcaud-Trocme et al. 2003).
Relevance of model components
We analyzed the importance of the different mechanisms included in the aEIF model, depending on the total conductance level.
When spike-triggered adaptation was removed, i.e., with b = 0, the performance was seriously degraded (Γ = 0.67), and as expected, the model tended to fire too many spikes (38% extra spikes, 17% missing spikes). However the effect was much more pronounced in LC states (Γ = 0.60, 45% extra spikes) than in HC states (Γ = 0.76, 30% extra spikes). We compensated this removal by inserting a constant current equal to the average adaptation current at 15 Hz (which amounts to lowering the resting potential; Fig. 3A ), and the degradation was still significant but less important (Γ = 0.81) and more balanced (16% extra spikes, 15% missing spikes). Again, the degradation was much stronger in LC states (Γ = 0.73) than in HC states (Γ = 0.87). The absence of spike-triggered adaptation could not be compensated further by adjustment of the threshold parameters, because the variance of the effective threshold (see methods) increased significantly (best value >0.8 mV2 compared with 0.02 mV2 in the original aEIF model). The input dependence of the effective threshold reflects a problem of generalization, which is also apparent in Fig. 3A, where the performance degrades as soon as the firing rate departs from 15 Hz.
Similarly, we removed subthreshold adaptation by setting a = 0 (or, equivalently, using the slope of the I–V curve at rest, where adaptation was minimal) and calculated the optimal threshold parameters again, which were now VT = −49.3 mV and ΔT = 1.4 mV. As expected, the effective threshold was more variable (0.15 mV2 instead of 0.02 mV2). The performance in prediction was impaired (Γ = 0.94), although much less than when spike-triggered adaptation was removed. We found that we could compensate for the absence of subthreshold adaptation by inserting a constant current equal to the average adaptation current at −58 mV (the average potential of the detailed model for the MC input scenario at 15 Hz), with the corresponding optimal threshold parameters (VT = −50.6 mV and ΔT = 2 mV). In this case the coincidence factor was on average Γ = 0.952 instead of 0.957 with the original parameter values. We note that the replacement of adaptation by a constant current does not reduce the number of parameters because the amplitude of the current needs to be optimized. Moreover, without subthreshold adaptation, we lose the overshoot in Fig. 2C that is well reproduced by the full aEIF model.
We also tested the sensitivity of the model to the value of the adaptation time constant (Fig. 3B). When adaptation was made twice as fast (τw = 72 ms and b = 161 pA—b had to be doubled so as to conserve the same mean level of adaptation), the performance of the aEIF model was reduced (Γ = 0.91) and again was better in HC states (Γ = 0.93) than in LC states (Γ = 0.88). A similar tendency was seen when adaptation was made twice slower (τw = 288 ms and b = 40 pA), although the degradation was less obvious (from Γ = 0.94 in LC states to Γ = 0.95 in HC states). Changing the value of τw had no significant impact on the optimal values of the threshold parameters VT and ΔT.
Reliability of the model in a noisy setting
With real experimental recordings, the performance of the model could be impaired in three different ways: 1) various sources of experimental noise could lead to a misestimation of parameters; 2) the parameters of the real neuron could drift during the course of an experiment; and 3) when testing the model, noise in current injections could degrade the prediction performance.
With respect to measurement noise, we identified the parameters b and τw as the most critical ones because we extract them from the speed of membrane depolarization (see methods; estimation of passive and stationary properties is less sensitive to noise). To check the stability of our approach, we artificially added noise to the voltage trace and repeated parameter extraction. Even with large noise (Fig. 1C, gray trace; random Gaussian numbers with SD 0.5 mV added to each sample), parameters b and τw could be extracted with an error of <10%. This error could be further reduced by repetitions of the measurements. Other errors in parameter extraction could result from misestimation of the bridge resistance (in a single-electrode set-up). To evaluate the consequences of this and other parameter mismatches or simply of a drift of neuronal properties during the time-course of an experiment, we arbitrarily modified all parameter values by 5% (gL = 275 nS, C = 295 pF, τw = 151 mS, a = 4.2 nS, b = 0.076 nA) and reran the procedure for optimization of the threshold parameters, which yielded VT = −50.05 mV and ΔT = 1.6 mV. The aEIF model with these altered parameters was tested on the same recordings of the detailed model used before. We found that the performance in predicting the spike trains of the detailed model was reduced, but to a very reasonable extent: the average coincidence factor was Γ = 0.95 (instead of 0.96 with the correct parameter values), and there were on average 5% extra spikes and 4% missing spikes (instead of 3% and 4%). Interestingly, the nonoptimality of the parameter values could be noticed in the minimal variance of the effective threshold (see methods), which rose from 0.02 mV2 with the correct values to 0.17 mV2 with the incorrect ones.
Prediction performance could also be impaired by noise in the injection of synaptic conductances (experimentally, by the dynamic clamp protocol; Sharp et al. 1993). To test the potential impact of such errors on our procedure, we added noise to the synaptic conductances in the form of random Gaussian numbers added to each sample, with SD equal to 10% of the average conductance. First, it had no effect at all on the outcome of the optimization procedure for the threshold parameters. This is not surprising because that procedure relies only on the firing rates. Second, it caused only a minor reduction in prediction performance: the coincidence factor was Γ = 0.952 instead of 0.957. When the noise level was raised to 20%, Γ was reduced to 0.94. This high reliability can be explained by the fact that cortical neurons in vitro (Mainen and Sejnowski 1995) and noisy spiking neuron models (Brette and Guigon 2003) display reproducible spike trains in response to fluctuating inputs.
Our results show that an exponential integrate-and-fire model with adaptation can accurately predict the spike trains of a detailed regular spiking neuron model driven by realistic, conductance-based, synaptic inputs. The aEIF model combines an exponential spike mechanism (Fourcaud-Trocme et al. 2003) with an adaptation equation with reset (Izhikevich 2003). It predicted, on average, 96% of the spikes with 2-ms precision, improving significantly on the leaky integrate-and-fire model with adaptation (88%).
We also found that adaptation currents had an important role in shaping spike trains in LC states, but this role was reduced in HC states, typical of in vivo activity: the aEIF model was still able to reproduce the target spike trains with high accuracy when the adaptation mechanism was seriously modified. This finding is consistent with the fact that synaptic activity in vivo can considerably alter the electrophysiological neuron classes defined in vitro (Steriade 2004).
Interestingly, the Izhikevich model was designed to reproduce qualitatively the major electrophysiological neuron classes (e.g., regular spiking, bursting, chattering) by just changing a few parameters, based on a mathematical theory. The same properties hold for our model, which differs from the model of Izhikevich only by the spike mechanism: indeed, if the reset value Vr is changed from −70 to −47 mV, the aEIF model turns into a bursting neuron (Fig. 3C); for other parameter values, the model exhibits postinhibitory rebound spikes after a prolonged hyperpolarization (Fig. 3D). Our results show that the simple aEIF model can quantitatively reproduce the spike trains of a detailed model of a regular spiking neuron with realistic in vivo–like inputs, but the fact that it is versatile enough to reproduce qualitatively several electrophysiological classes opens promising perspectives—although fitting the model to bursting neurons might be a harder task. On the practical side, these results make the aEIF model a good candidate for large-scale simulations of realistic cortical networks; on the theoretical side, the analysis of networks of aEIF neurons could also help us understand how intrinsic properties defined electrophysiologically in vitro affect the neuron behavior in vivo.
We thank R. Jolivet for stimulating discussions.
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 © 2005 by the American Physiological Society