Made-to-order spiking neuron model equipped with a multi-timescale adaptive threshold

Identical fluctuating currents elicit similar spike trains from the same neuron (Bryant and Segundo, 1976 ; Mainen and Sejnowski, 1995 ), but there are significant differences in the spike trains of different neurons to the same currents. The brain comprises a variety of neurons possessing such different response characteristics (Shinomoto et al., 2009 ). A model that is capable of precisely reproducing spike response of any neuron to any input current is an essential first step to both building a brain simulator (Brette et al., 2007 ; Diesmann et al., 1999 ; Gewaltig and Diesmann, 2007 ; Izhikevich and Edelman, 2008 ; Markram, 2006 ; McIntyre et al., 2004 ; Plesser and Diesmann, 2009 ), and understanding the computational mechanisms of neurons (Abeles, 1991 ; Beggs and Plenz, 2003 ; Diesmann et al., 1999 ; Ikegaya et al., 2004 ).

The Hodgkin–Huxley (HH) model is a standard model for describing neuronal firing, and has been continually revised by including ionic channels to reproduce standard electrophysiological measurements of cortical neurons (Hodgkin and Huxley, 1952 ; Koch, 1999 ). It has recently become possible to use elaborate simulation platforms, such as NEURON (Hines and Carnevale, 1997 ) and GENESIS (Bower and Beeman, 1995 ), for reproducing experimental data. Because of nonlinearity and complexity, however, parameter optimization of the HH type models is a notoriously difficult problem (Achard and De Schutter, 2006 ; Goldman et al., 2001 ; Huys et al., 2006 ), and these models require a high computational cost, which hinders performing the simulation of a massively interconnected network.

For these reasons, simplified phenomenological models such as the leaky integrate-and-fire (LIF) model (Stein, 1965 ) have been proposed, and these models have been useful for studying learning, memory, and the dynamics of neural networks (Gerstner and Kistler, 2002 ; Hansel and Sompolinsky, 1998 ). Due to this simplification, many details of the biophysical properties of neurons are abandoned. However, simplified models that extend beyond the qualitative description of neuronal firing have been developed (Fourcaud-Trocmé et al., 2003 ; Gerstner and Kistler, 2002 ; Izhikevich, 2004 ). Nonlinear integrate-and-fire models have been successfully applied to some typical neurons (Badel et al., 2008 ; Brette and Gerstner, 2005 ). The spike response model (SRM) exhibits fairly good predictive performances (Jolivet et al., 2004 , 2006 ; Kobayashi and Shinomoto, 2007 ). It is feasible to fit the LIF model to neurophysiological data, in particular, to reproduce the statistical features in interspike intervals (Inoue et al., 1995 ; Lansky and Ditlevsen, 2008 ; Lansky et al., 2006 ; Tuckwell and Richter, 1978 ), but parameter optimization in the SRM or nonlinear models for predicting spike times is a difficult problem because of their nonlinearity or a large number of parameters.

In this study, we propose a minimal model for accurately predicting a rich variety of spike responses in vitro and develop a procedure for optimizing model parameters. The model consists of the simple leaky integrator and a multiple-timescale adaptive threshold, which is referred to as the MAT model. The MAT model has low computational cost and high predictive performances for both stationary and nonstationary fluctuating currents. The model is adaptable to various types of neurons, including regular spiking (RS), intrinsic bursting (IB), and fast spiking (FS) neurons, by simply fitting three parameters, and is also capable of describing the spike patterns of chattering (CH) neurons.

The dynamics of a nonresetting leaky integrator proposed in this study is given as,

where τ_m, V(t), R, and I(t) are the membrane time constant, model potential, membrane resistance, and input current, respectively. When the model potential reaches the spike threshold θ(t), a spike is generated. In our model, V(t) is not reset even if the model assumes spikes. Instead, the spike threshold increases when the model assumes a spike and then decays toward the resting value. In the standard adapting threshold model, the spike threshold θ(t) decays exponentially (Brandman and Nelson, 2002 ; Chacron et al., 2003 ; Geisler and Goldberg, 1966 ). In the present study, we extended this adaptive threshold rule into a MAT rule.

where t_k is the kth spike time, L is the number of threshold time constants, τ_j (j = 1,…,L) are the jth time constants, α_j (j = 1,…,L) are the weights of the jth time constants, and ω is the resting value. To avoid singular bursting, we introduced an absolute refractory period of 2 ms during which the model is not allowed to fire even if the membrane potential exceeds the threshold. If the membrane potential still remains above threshold for the 2 ms, the model assumes another spike. We refer to the model as the MAT(L) model, including a single timescale model as the MAT(1). Here, we integrate the equation with the Euler method at a step size of 0.001 ms and verify results with backward Euler.

The MAT(L) model is specified by the 2L + 3 parameters {τ_m, R, τ_j, α_j (j = 1, 2,…,L), ω}. The model comprises two dynamic components: the membrane potential and spike threshold (Figure 1 A). We first determined the parameters of the membrane potential dynamics {τ_m, R} by fitting the voltage response to fluctuating currents with a nonresetting leaky integrator. These parameters were optimized for minimizing the squared error,

We determined the membrane time constant τ_m for 34 neurons (see Section “Electrophysiological Recordings”). The membrane time constant ranged from 2.0 to 12.0 ms, with a mean ± SD of 6.2 ± 2.2 ms. By confirming that predictive performance is not sensitive to this parameter, we adopted a common membrane time constant for all neurons, τ_m = 5 ms. The resistance R simply scales the potential and is unrelated to the spike time prediction of the model; therefore, we adopted a common resistance R = 50 MΩ for all neurons.

Second, we determined the parameters for the spike threshold dynamics {τ_j, α_j, ω}. In this study, we examined three cases, single timescale (L = 1), two timescales (L = 2), and three timescales (L = 3). We selected a set of time constants τ_j from trial parameters (10, 50, and 200 ms), so that the model performance defined below is maximized on average for all neurons. Given a set of time constants, the firing properties of individual neurons are solely represented by {α_j, ω}. These tailor-made parameters were optimized so that the coincidence factor Γ defined below was maximized for the set of sample data comprising three input currents (mean ± SD, 0.42 ± 0.14, 0.37 ± 0.25, and 0.29 ± 0.19 nA) and the corresponding spike trains. We used the Nelder–Mead method (Nelder and Mead, 1965 ) for the maximization procedure.

The performance of the spike time prediction is assessed with a coefficient measuring the degree of coincidence of a model spike train with a real spike train. Among various methods for measuring the similarity between two spike trains (Tsubo et al., 2004 ; Victor and Purpura, 1996 ), we adopted here the coincidence factor Γ (Jolivet et al., 2004 , 2006 , 2008a ; Kobayashi and Shinomoto, 2007 ), defined by

where N_data and N_model, respectively, represent the numbers of spikes in the data and the model, N_coinc is the number of coincident spikes with an allowable range of time Δ. Here, N_coinc is subtracted by the expected number of coincidences between the data and the random Poisson spikes of the given rate, 〈N_coinc〉 = 2vN_data, so that the coefficient takes a value close to 0 for mutually independent spike trains. The coefficient is divided by 1 − 2vΔ so that it takes a value 1 if all the spikes coincided within Δ. In this study, Δ is chosen as 2 ms.

Each neuron model was tailored to every neuron using sample current–voltage data, by maximizing the coincidence score Γ. The model was then tested with a new current to see how well it predicted spike times. The predictive performances was assessed based on the coincidence scores normalized with the trial-to-trial variation of the experiments, defined by

Various models were compared based on the predictive scores Γ_A evaluated for individual neurons under conditions in which the neurons were firing at rates close to 10 and 30 Hz (in case of FS, 70 and 140 Hz) respectively, with different degrees of fluctuation (SD/mean, 0.35 and 0.70). We took into account the predictive scores Γ_A, provided that spikes were elicited within a range from 5 to 200 Hz, and with a degree of trial-to-trial reproducibility, Γ_{experiment-experiment} > 0.1. As a result, 119 data among 4 × 34 = 136 data were available, and each neuron was evaluated by the average score for the three or four data.

The current injection experiments for 34 neurons in layers 2/3 and 5 of the rat motor cortex were performed exclusively for the present study, whereas the conductance injection experiments for nine neurons were performed previously for the study of Miura et al. (2007) . We abbreviate the animal preparation details that are common with the methods published in the paper, and describe the essentials needed in the present study. All experiments were performed in accordance with animal protocols approved by the Experimental Animal Committee of the RIKEN Institute.

Wistar rats (postnatal days, 17–27) were deeply anesthetized with diethyl ether gas and then decapitated. Data used for examining the spike time prediction was taken from rat cortical slices (400-μm thick) sectioned using a microslicer (PRO-7; Dosaka, Kyoto, Japan). After incubation for 30 min at 31 °C and at least 1 h of recovery at room temperature, each slice was transferred to a submerged-type recording chamber continuously circulating with normal artificial cerebrospinal fluid (ACSF; 32°C), which consisted of (in mM) 124 NaCl, 2.5 KCl, 1.2 KH₂PO₄, 26 NaHCO₃, 1.2 MgSO₄, 2.5 CaCl₂, and 25 D-glucose and was saturated with 95% O₂ and 5% CO₂ gas. We recorded in the current-clamp or conductance-clamp whole-cell configuration from the soma using patch pipettes filled with (in mM) 140 K-gluconate, 2 NaCl, 1 MgCl₂, 10 HEPES, 0.2 EGTA, 2 5-ATPNa₂, 0.5 GTPNa₂, and 10 biocytin, pH 7.4 (Tsubo et al., 2007 ). The neuron membrane potentials were recorded using a current clamp amplifier (Axoclamp 2B; Molecular Devices, Union City, CA, USA) in the conventional bridge mode. Recorded signals were digitized using an analog–digital interface (Digidata 1322A; Molecular Devices) at 40 kHz in current injection experiments and 50 kHz in conductance injection experiments. To block the spontaneous input via AMPA, NMDA, and GABA_A receptors, a cocktail of 6-cyano-7-nitroquinoxaline-2,3-dione (10 μM; Sigma, St Louis, MO, USA), DL-2-amino-5-phosphonopentanoic acid (25 μM; Sigma), and bicuculline methiodide (30 μM; Sigma) was added to ACSF and bath-applied to the cortical slices. Based on the criteria suggested by Nowak et al. (2003) , the 34 recorded cells were classified into 22 RS, 8 IB, and 4 FS neurons.

For the current injection experiments, the input current

was prepared by mimicking natural inputs comprising the synaptic input from thousands of unitary excitatory and inhibitory post-synaptic currents (EPSCs and IPSCs) with amplitudes I_exc = 0.1 nA and I_inh = 0.033 nA, and decay-time constants τ_exc = 1 ms and τ_inh = 3 ms, respectively (Tsubo et al., 2004 ). A = 0.1–1.2 is a scale factor, and the voltage dependence of EPSCs and IPSCs were neglected (Stevens and Zador, 1998 ). Arrival time of excitatory and inhibitory spikes {t_k} and {t_j} were drawn from the Poisson processes with statistical rates r_exc and r_inh, respectively. During the stationary fluctuating current trials, the statistical rates were fixed so that the mean ± SD of the input currents were 0.40 ± 0.14 nA (r_exc = 6.88 kHz, r_inh = 2.88 kHz) or 0.40 ± 0.28 nA (r_exc = 24.52 kHz, r_inh = 20.52 kHz). Then, these injected currents were scaled by a factor A; for instance, the current ID “0.20 ± 0.14 nA” indicates the current made from “0.40 ± 0.28 nA” scaled by 0.5. For nonstationary current trials, the statistical rates were sinusoidally time-varied so that the mean ± SD of the input currents were [0.30 + 0.10sin (2πt)] ± [0.21 + 0.07sin (2πt/2)] nA. One trial lasted 5 s and the identical currents were repeated twice with a 15-s intertrial interval.

For conductance injection experiments, the excitatory and inhibitory input conductances g_exc(t) and g_inh(t) were described by the Ornstein–Uhlenbeck processes (OUP; Uhlenbeck and Ornstein, 1930 ) with time constants τ_exc = 2.7 ms and τ_inh = 10.5 ms, mean conductances g_exc0 and g_inh0 = Bg_exc0, and SDs σ_exc = g_exc0/2 and σ_inh = Bg_exc0/4, respectively. The excitatory mean conductance g_exc0 was set to 0.018 and 0.024 μS for parameter determination, and 0.012 and 0.030 μS for model assessment. The ratio of excitatory and inhibitory mean conductance B = g_exc0/g_inh0 was set to −0.348 and −0.370, respectively. Using a conductance clamp method (Robinson and Kawai, 1993 ; Sharp et al., 1993 ), we injected these conductances as a fluctuating current I(t):

where V(t) is the recorded membrane potential, and E_exc = 0 mV and E_inh = −75 mV are the reversal potential of excitatory and inhibitory synaptic currents, respectively. One trial lasted 10 s. In each trial, a set of time constants, mean values, and SDs of excitatory and inhibitory OUP conductances were fixed. The identical conductances were repeated twice with a 30-s intertrial interval.

We first study the timescales essential for predicting spike times of cortical neuron in vitro. We design the minimal MAT model for predicting spike times. Second, we show that the MAT model accurately predicts the spike responses of a variety of cortical neurons and those for various input conditions by adjusting only three parameters. We compare the predictive performances of this model with three standard models: a HH model, LIF model, and SRM. The parameters of each model are fitted to data, using widely accepted methods (Appendix).

To study the timescales essential for predicting spike times, we compare the predictive performances of MAT(L) models. We examine three cases, single timescale (L = 1), two timescales (L = 2), and three timescales (L = 3). Timescales of the adaptive threshold dynamics τ_j are selected from 10, 50, and 200 ms. The predictive performances achieved by all possible combinations of the three timescales are summarized in Table 1 . The optimal timescale for the MAT(1) model is τ₁ = 50 ms. The optimal timescales for the MAT(2) model are τ₁ = 10 ms and τ₂ = 200 ms, and this model resulted in the best performance for all combinations of the three timescales. It is noteworthy that the predictive performance of this MAT(2) model was even better than that of the MAT(3) model. Though adding parameters may improve fitting, a wider class model does not necessarily acquire the better predictive ability by the learning from a finite set of samples, because a model with overabundant parameters may exhibit overfitting for finite samples. We refer to the MAT(2) model with τ₁ = 10 ms and τ₂ = 200 ms as the MAT* model. In the following subsections, we adopt the MAT* model for predicting spike times.

Identical fluctuating currents induce different responses in different neurons (Figure 1 B). Neurons are categorized into three types based on conventional categorical criteria given by McCormick et al. (1985) and Nowak et al. (2003) : RS, IB, and FS. A variety of spiking patterns among these neurons elicited by the same input current are predicted simply by selecting three parameters of the MAT* model. It should be noted that the model potentials are computed with the same leaky integrator and therefore are identical for all neurons. In the MAT* model, different firing patterns are manifested solely due to the difference in the threshold dynamics (Figure 1 C).

For a set of 34 neurons, three parameters α₁, α₂, and ω of the MAT* model are determined so that the fraction of common spikes between the model and experiment is maximized for the sample dataset. RS, IB, and FS neurons that are categorized based on the conventional categorical criteria are separately localized in a component plane of the fast and slow threshold dynamics α₁ and α₂ (Figure 2 A).

H(t) (Eq. 3) are depicted with respect to three prototypical neurons (Figure 2 B): RS neurons contain both large amounts of fast α₁ and slow α₂ components of the threshold dynamics, whereas IB neurons contain a small amount of the fast component α₁ and a large amount of the slow component α₂, and FS neurons contain a small amount of the fast component α₁ and a negligible amount of the slow component α₂.

We compare the MAT* model with a HH model, LIF model, and SRM for their ability to predict spike times elicited from fluctuating currents by neurons in a rat cortical slice preparation (Figure 3 B). In each model, the model parameters are optimized from sample data consisting of injected currents and the induced membrane potential of an individual neuron, and then test with a new current for the ability to predict spike times of the same neuron. The predictive ability is evaluated with a coincidence factor Γ between predicted spikes and target spikes within a range of 2 ms (see Section “Materials and Methods”). Identical fluctuating currents elicit similar spike trains from the same neuron (Figure 3 A). The score is normalized by variation among experimental trials so that the best prediction could be evaluated as Γ_A close to 1. The mean ± SD of the predictive scores Γ_A are evaluated for 34 neurons (Figure 3 C).

The HH model is fitted to experimental conductance data, but is weak in its ability to predict spike times, resulting in a predictive scores Γ_A = 0.51 ± 0.26. The LIF model and SRM with an adaptive threshold yield Γ_A = 0.66 ± 0.26, and Γ_A = 0.70 ± 0.26, respectively. The MAT* model achieves the best predictive score of Γ_A = 0.89 ± 0.21.

In addition to the stationary situations in which synaptic input currents fluctuate finely, we also mimic the nonstationary situation in which behavioral stimuli or motor commands greatly modulate the input currents. In this study, the temporal mean of the fluctuating current is modulated sinusoidally with a period of 1 s, while the SD of fluctuation is slowly modulated with a period of 2 s (Figure 4 A). The predictive scores Γ_A evaluated for an RS neuron for the entire 5 s interval are 0.49, 0.78, 0.64, and 0.84 for the HH, LIF, SRM, and MAT* models, respectively.

The four models are also compared for their ability to predict transient firing rate adaptation of an RS neuron to the initiation of current injections. All preexisting models fail to reproduce the transient adaptation, but the MAT* model is able to predict the adaptive phenomena (Figure 4 B). The predictive scores Γ_A evaluated for the initial 1 s interval are 0.49, 0.52, 0.38, and 0.68 for the four models.

In this way, the MAT* model is superior to the three standard models in predicting spikes for stationary fluctuating currents, slowly modulated nonstationary currents, and the onset of current injection.

In addition to the current injection experiments, we also examine the predictive ability of the four models for the conductance injection experimental data. We analyze the data assuming that τ_m = 5 ms and R = 50 MΩ. Predictive scores Γ_A for the HH, LIF, SRM, and MAT* models are 0.52 ± 0.23, 0.64 ± 0.17, 0.45 ± 0.21, and 0.84 ± 0.18; thus, the MAT* model consistently exhibits a better performance for predicting spike times.

We also explore a space of the parameters α₁, α₂, and ω to examine if MAT* model can produce the typical firing responses to a rectangular current. As we confirmed in the application to biological data, the MAT* model can represent response characteristics of the RS, IB, and FS neurons. Numerical simulations were performed with various parameter values (Figures 5 A,B).

It is noteworthy that the specific firing pattern of the “chattering” (CH) neuron (Gray and McCormick, 1996 ) can be reproduced in the parameter range α₁ < 0 and α₂ > 0: A negative factor for the fast timescale, α₁ < 0, lowers the firing threshold and initiates bursting with the given absolute refractory period. However, as bursting continues, the positive factor for the longer timescale, α₂ > 0, accumulates and eventually terminates bursting. This cycle is repeated to produce “chattering”.

We have proposed a new model named “MAT” that is equipped with a multi-timescale adaptive threshold predictor and a nonresetting leaky integrator. The minimal MAT* model is capable of reproducing a rich variety of neuronal spike responses from RS, IB, and FS neurons by simply adjusting three adaptive threshold parameters (Figure 1 ). The model expresses a variety of the firing characteristics rather than a conventional discrete categorization (Figure 2 ). The MAT* model provides the highest predictive performance among four models including the conventional HH, LIF, and SRM (Figure 3 ), and is superior for predicting spikes for greatly fluctuating inputs (Figure 4 ). It also has a potential for generating the specific firing pattern of the CH neurons (Figure 5 ).

An essential component of the MAT* model is the introduction of multiple time constants into the threshold dynamics. Adaptive threshold time constants were selected from among trial parameters of 10, 50, and 200 ms. These time constants would be related to biological ionic currents (Hille, 2001 ; Koch, 1999 ): [10 ms]: Fast transient Na⁺ current and delayed rectifier K⁺ current, [50 ms]: hyperpolarization-activated cation current, and K⁺ A-current, and [200 ms]: noninactivating K⁺ current, hyperpolarization-activated cation current, and Ca⁺²-dependent K⁺ current.

Eventually, 10 and 200 ms were selected for the MAT* model based on its predictive performance. Different types of neuronal firing can be reproduced and predicted by simply adapting components of shorter and longer timescale dynamics to experimental data. RS, IB, and FS neurons are separately localized in a component plane of the two timescale dynamics. RS neurons have a large amount of both shorter and longer timescale threshold dynamics; IB neurons have a smaller amount of shorter one; and FS neurons have a small amount of shorter one and a negligible amount of longer one. If the firing characteristics are modified by some external operations on ion channel densities, it would be worthwhile to fit the MAT model and examine how the weight parameters vary in such a process.

The effectiveness of our MAT model implies that neuronal firing dynamics are mainly governed by multiple refractory effects with different timescales, and the characteristics are mainly determined by their combination. The time constants (10, 50, and 200 ms) were chosen as trial parameters, and turned out to be good. When targeting individual neurons, there is still room for improvement in performance, by adding more timescales or selecting the more suitable set of time constants for individual neurons.

We examined a three timescale MAT(3) model that employs all three timescales (10, 50, and 200 ms) for its ability to predict spike times. The results showed that the MAT(3) model performance was slightly inferior to that of the two timescale MAT(2) model (Table 1 ). Mathematically, the MAT(2) model is a special case of the MAT(3) model. These would have occurred because a model with too many parameters tends to over-fit to sample data and its predictive performance is poorer than the simpler models. In devising a model of reality, we have to remember the principle of Occam’s razor. As long as the experimental time for parameter identification is finitely available, the model should be as concise as possible.

The membrane time constant τ_m fitted to the 34 neurons were distributed as 6.2 ± 2.2 ms, which was much shorter than the most common membrane time constant estimates of 10–20 ms (McCormick et al., 1985 ; Yang et al., 1996 ). The significant difference between these estimates was due to the difference in conditions; traditional experiments measure the neuronal reaction from the resting state, whereas we were measuring the ongoing neuronal reaction as neurons were firing in a frequency close to the in vivo condition. Therefore, our membrane time constant represented an effective leak that sums up the effects of all ion channels.

Although the MAT model devised in this study outperformed other models for predicting spikes of cortical neurons, it still failed to predict approximately 10% of the spikes compared to the trial-to-trial variation in experiments. The apparent imperfection of our model is associated with a lack of nonlinearity in the model membrane dynamics. There are a number of nonlinear models that are successful in representing a variety of firing patterns in response to rectangular and ramp currents (Izhikevich, 2004 ). Nonlinear models are also capable of predicting spike times, and some models such as the adaptive nonlinear integrate-and-fire models are efficient (Badel et al., 2008 ; Brette and Gerstner, 2005 ; Fourcaud-Trocmé et al., 2003 ). However, the parameter optimization of these models is a very difficult problem; the different optimization methods and the different initial conditions can lead to vastly different results (Achard and De Schutter, 2006 ; Goldman et al., 2001 ). It should be noted that even if we have a complete set of models, the model identification method cannot always find the true parameters, but rather it often determines parameters that render the model worse than simpler models. The single timescale adaptive threshold model MAT(1) is less powerful than nonlinear models in describing the mathematical details of spiking mechanism, but it nevertheless won the first place at the “International Competition on Quantitative Neuron Modeling” in both 2007 and 2008 (Jolivet et al., 2008a ,b ).

Another possible extension of the MAT model is the introduction of soma shape and dendritic spatial structure (Pinsky and Rinzel, 1994 ). In another competition for predicting somatic spike times in response to current injection in both the soma and apical dendrite (Jolivet et al., 2008b ; Larkum et al., 2004 ), we extended the leaky integrator in the MAT(1) model into a two compartment model and succeeded in predicting spike times with the highest score. This fact indicates the efficiency of the compartment model.

There is also room for revising adaptive threshold dynamics. In our MAT model, the dynamic threshold depends simply on previous spike times. The level of membrane potential for the initiation of an action potential depends not only on the preceding spike times but also on the instantaneous rising trend of the potential level (Azouz and Gray, 2000 ; Kobayashi and Shinomoto, 2007 ).

Recently, Mihalas and Niebur (2009) proposed a dynamic threshold model that is equipped with not only multiple timescales but also the above-mentioned voltage dependence. The Mihalas and Niebur model and our MAT model have contrasting advantages that are also shortcomings, mostly due to the difference in the model complexity; the Mihalas and Niebur model can reproduce rich phenomena such as the inhibitory rebound due to elaborate setups including the voltage-dependent threshold, while the MAT model can readily achieve a good predictive performance due to its simplicity and the small number of free parameters.

Modeling neuronal microcircuitry has been initiated based on either the HH models (Markram, 2006 ; Traub et al., 2005 ) or LIF model (Diesmann et al., 1999 ; Gewaltig and Diesmann, 2007 ). The HH models are biologically plausible in its biophysical features, but have disadvantages in their hard parameter optimization and high computational cost. In contrast, the LIF model need much lower computational cost, but it cannot account for the spiking characteristics in vivo (Shinomoto et al., 1999 ).

A nonlinear model proposed by Izhikevich (2004) is worth notice, because this is not only simple but also versatile enough to represent a variety of neurons. In contrast to the Izhikevich model, our MAT model is essentially linear. Nonlinear changes of the state of the neuron are only required at the time of input or output spikes. Accordingly, it may be possible to speed up the numerical simulation of the neuronal circuitry by exactly integrating the subthreshold dynamics (Plesser and Diesmann, 2009 ). Thus the linearity may be beneficial in considering large-scale simulation.

In conclusion, the multi-timescale adaptive threshold model can be tailored to any cortical neuron to reproduce and predict precise spike times for any input current. This MAT model is a step toward building a brain simulator that comprises a variety of neurons with broadly distributed characteristics. The next step will be to include interactions between model spike neurons.

In this study, we describe three standard neuron models namely a Hodgkin–Huxley (HH) model, leaky integrated-and-fire (LIF) model, and spike response model (SRM), which are compared with the MAT model in their ability to predict spike times. As with the MAT model, all models are numerically integrated with the Euler method with a sufficiently small time step of 0.001 ms. Accuracy of the numerical integration is ascertained using the backward Euler method.

As a HH model, we adopted the Destexhe model (Destexhe and Paré, 1999 ),

where C is the membrane capacitance, V is the membrane potential, E_L = −80 mV is the reversal potential of the leak current, and I(t) is the input current.

The voltage-dependent Na⁺ current is described by,

where E_Na = 50 mV is the reversal potential of Na⁺ and V_S = 10 mV is the inactivation shift. V_T was optimized to each neuron from the sample data.

The “delayed rectifier” K⁺ current is described by,

Due to the current I_M the model exhibits spike frequency adaptation.

The Destexhe model is specified by the six parameters {C, g_Na, g_Kd, g_M, g_L, V_T}. First, we determined the conductance parameters a = {C, g_Na, g_Kd, g_M, g_L} with a fixed value of the threshold V_T by the method of Huys et al. (2006) , which has widely been used for reproducing the membrane dynamics. These parameters were optimized for minimizing the squared error with non-negative constraints,

Second, the threshold V_T was optimized so that the coincidence factor Γ (see Section “Materials and Methods) was maximized for the sample data set. We defined the spike time of the HH model neuron as the time when the derivative of the membrane potential exceeded 10 mV/ms.

In the LIF model, the dynamics of membrane potential are given by,

where τ_m, V(t), R, and I(t) are the leak time constant, model potential, membrane resistance, and input current, respectively. The model assumes a spike if the membrane potential V exceeds a threshold θ. In this study, we adopt a partial resetting rule proposed by Bugmann et al. (1997) and Troyer and Miller (1997) ,

The LIF model is specified by three parameters {τ_m, R, θ}. As with the MAT model, we first determined {τ_m, R} by fitting the voltage response to fluctuating currents with the leaky integrator. By confirming that predictive performance is not sensitive to these parameters, we adopted a common time constant and resistance for all neurons, τ_m = 5 ms and R = 50 MΩ. Second, we optimized θ for maximizing the coincidence factor Γ (see Section “Materials and Methods”) for the sample data set.

In SRM (Gerstner and Kistler, 2002 ), the membrane potential is given by,

where t_f is the last spike time, and κ(t) and η(t) represent the integration kernel and a spike shape, respectively. Among the various types of SRM, we adopted the adaptive threshold type (Jolivet et al., 2006 ) in which the spike threshold θ(t) is lifted when fired and decays exponentially (Brandman and Nelson, 2002 ; Chacron et al., 2003 ; Geisler and Goldberg, 1966 ) as,

In this equation, t_k is the kth spike time, α is the weight, τ_α is the time constant of the threshold, and ω is the resting value of the threshold. As with the MAT model, we introduced an absolute refractory period of 2 ms.

The SRM with adaptive threshold is specified by two functions and three parameters {κ(t), η(t), τ_α, α, ω}. We first fitted {κ(t), η(t)} by fitting the voltage response to fluctuating currents. These functions were optimized for minimizing the squared error (Kobayashi and Shinomoto, 2007 ),

Second, we determined the spike threshold parameters {τ_α, α, ω}. The time constant τ_α was selected from 10, 50, and 200 ms, and we adopted a common time constant for all neurons. Then, we optimized {α, ω} to maximize the coincidence factor Γ for the sample data set. The best performance was achieved when τ_α = 50 ms.

Mathematically, the MAT model is very similar to a special case of SRM. We consider the MAT model (Eqs. 1–3) and take the new variable U: to translate into an SR-like model. U is not the membrane potential, but U is given by the SRM-like equation as,

In this representation, the threshold is constant. When the potential U(t) exceeds the threshold ω, a spike is generated.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

We are grateful to T. Fukai for helpful comments. This study is supported in part by Grants-in-Aid for Scientific Research to SS (20300083, 20020012), and to YT (20700280) from the MEXT, Japan. RK was supported by the Research Fellowship of the Japan Society for the Promotion of Science for Young Scientists.