Frontiers in Computational Neuroscience Computational Neuroscience Made-to-order Spiking Neuron Model Equipped with a Multi-timescale Adaptive Threshold

Information is transmitted in the brain through various kinds of neurons that respond differently to the same signal. Full characteristics including cognitive functions of the brain should ultimately be comprehended by building simulators capable of precisely mirroring spike responses of a variety of neurons. Neuronal modeling that had remained on a qualitative level has recently advanced to a quantitative level, but is still incapable of accurately predicting biological data and requires high computational cost. In this study, we devised a simple, fast computational model that can be tailored to any cortical neuron not only for reproducing but also for predicting a variety of spike responses to greatly fl uctuating currents. The key features of this model are a multi-timescale adaptive threshold predictor and a nonresetting leaky integrator. This model is capable of reproducing a rich variety of neuronal spike responses, including regular spiking, intrinsic bursting, fast spiking, and chattering, by adjusting only three adaptive threshold parameters. This model can express a continuous variety of the fi ring characteristics in a three-dimensional parameter space rather than just those identifi ed in the conventional discrete categorization. Both high fl exibility and low computational cost would help to model the real brain function faithfully and examine how network properties may be infl uenced by the distributed characteristics of component neurons. For these reasons, simplifi ed phenomenological models such as the leaky integrate-and-fi re (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 sim-plifi cation, many details of the biophysical properties of neurons are abandoned. However, simplifi ed models that extend beyond the qualitative description of neuronal fi ring have been developed (Fourcaud-Trocmé et al. Nonlinear integrate-and-fi re 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 fi t the LIF model to neurophysiological data, in particular, to reproduce the statistical features in interspike intervals (Inoue et al. but parameter optimization in the SRM or nonlinear models for predicting spike times is a diffi cult problem because of their non-linearity or a large number of parameters. In this study, we propose a minimal model for accurately predicting …


INTRODUCTION
Identical fl uctuating currents elicit similar spike trains from the same neuron (Bryant and Segundo, 1976;Mainen and Sejnowski, 1995), but there are signifi cant 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 fi rst 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 fi ring, 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 diffi cult 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. various types of neurons, including regular spiking (RS), intrinsic bursting (IB), and fast spiking (FS) neurons, by simply fi tting three parameters, and is also capable of describing the spike patterns of chattering (CH) neurons.

MULTI-TIMESCALE ADAPTIVE THRESHOLD MODEL
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 fi re 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.

FITTING PROCEDURE
The MAT(L) model is specifi ed 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 1A). We fi rst determined the parameters of the membrane potential dynamics {τ m , R} by fi tting the voltage response to fl uctuating currents with a nonresetting leaky integrator. These parameters were optimized for minimizing the squared error,ˆ( (4) 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 confi rming 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 defi ned below is maximized on average for all neurons. Given a set of time constants, the fi ring properties of individual neurons are solely represented by {α j , ω}. These tailor-made parameters were optimized so that the coincidence factor Γ defi ned 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.

ASSESSMENT OF PREDICTIVE PERFORMANCE
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(Jolivet et al., , 2006(Jolivet et al., , 2008aKobayashi and Shinomoto, 2007), defi ned 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 coeffi cient takes a value close to 0 for mutually independent spike trains. The coeffi cient 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, defi ned by Various models were compared based on the predictive scores Γ A evaluated for individual neurons under conditions in which the neurons were fi ring at rates close to 10 and 30 Hz (in case of FS, 70 and 140 Hz) respectively, with different degrees of fl uctuation (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.

ELECTROPHYSIOLOGICAL RECORDINGS
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.
For the current injection experiments, the input current where g t was prepared by mimicking natural inputs comprising the synaptic input from thousands of unitary excitatory and inhibitory postsynaptic 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 fl uctuating current trials, the statistical rates were fi xed 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 fl uctuating 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 fi xed. The identical conductances were repeated twice with a 30-s intertrial interval.

RESULTS
We fi rst 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 fi tted to data, using widely accepted methods (Appendix).

TIMESCALES ESSENTIAL FOR PREDICTING SPIKES
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 τ 1 = 50 ms. The optimal timescales for the MAT(2) model are τ 1 = 10 ms and τ 2 = 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 fi tting, a wider class model does not necessarily acquire the better predictive ability by the learning from a fi nite set of samples, because a model with overabundant parameters may exhibit overfi tting for fi nite samples. We refer to the MAT(2) model with τ 1 = 10 ms and τ 2 = 200 ms as the MAT* model. In the following subsections, we adopt the MAT* model for predicting spike times.

FREE PARAMETERS OF THE MAT* MODEL REPRESENTING DIFFERENT KINDS OF NEURONS
Identical fl uctuating currents induce different responses in different neurons ( Figure 1B). 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 fi ring patterns are manifested solely due to the difference in the threshold dynamics ( Figure 1C). For a set of 34 neurons, three parameters α 1 , α 2 , 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 α 1 and α 2 (Figure 2A).
H(t) (Eq. 3) are depicted with respect to three prototypical neurons ( Figure 2B): RS neurons contain both large amounts of fast

MODEL PERFORMANCE
We compare the MAT* model with a HH model, LIF model, and SRM for their ability to predict spike times elicited from fl uctuating currents by neurons in a rat cortical slice preparation ( Figure 3B). 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 fl uctuating currents elicit similar spike trains from the same neuron ( Figure 3A). 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 3C). The HH model is fi tted 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.

PREDICTION OF SPIKE RESPONSES TO NONSTATIONARY CURRENT INPUTS
In addition to the stationary situations in which synaptic input currents fl uctuate fi nely, 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 fl uctuating current is modulated sinusoidally with a period of 1 s, while the SD of fl uctuation is slowly modulated with a period of 2 s ( Figure 4A). 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 fi ring 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 4B). 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 fl uctuating currents, slowly modulated nonstationary currents, and the onset of current injection.

PREDICTION OF SPIKE RESPONSES TO CONDUCTANCE INPUTS
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.

A VARIETY OF RESPONSES OF MAT* MODEL TO A RECTANGULAR CURRENT
We also explore a space of the parameters α 1 , α 2 , and ω to examine if MAT* model can produce the typical fi ring responses to a rectangular current. As we confi rmed 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 5A,B). It is noteworthy that the specifi c fi ring pattern of the "chattering" (CH) neuron (Gray and McCormick, 1996) can be reproduced in the parameter range α 1 < 0 and α 2 > 0: A negative factor for the fast timescale, α 1 < 0, lowers the fi ring threshold and initiates bursting with the given absolute refractory period. However, as bursting continues, the positive factor for the longer timescale, α 2 > 0, accumulates and eventually terminates bursting. This cycle is repeated to produce "chattering".

DISCUSSION
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 fi ring 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 fl uctuating inputs (Figure 4). It also has a potential for generating the specifi c fi ring pattern of the CH neurons (Figure 5).

BIOLOGICAL ENTITIES RELATED TO ADAPTIVE THRESHOLD DYNAMICS
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 rectifi er K + current, [50 ms]: hyperpolarization-activated cation current, and K + A-current, and [200 ms]: noninactivating K + current, hyperpolarization-activated cation current, and Ca +2 -dependent K + current.
Eventually, 10 and 200 ms were selected for the MAT* model based on its predictive performance. Different types of neuronal fi ring 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 fi ring characteristics are modifi ed by some external operations on ion channel densities, it would be worthwhile to fi t the MAT model and examine how the weight parameters vary in such a process.
The effectiveness of our MAT model implies that neuronal fi ring 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-fi t 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 identifi cation is fi nitely available, the model should be as concise as possible.

TIME CONSTANT OF THE MEMBRANE POTENTIAL DYNAMICS
The membrane time constant τ m fi tted 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 signifi cant 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 fi ring 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.

LIMITATION OF MAT MODELS AND POSSIBLE MODEL EXTENSION
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 fi ring 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-fi re models are effi cient (Badel et al., 2008;Brette and Gerstner, 2005;Fourcaud-Trocmé et al., 2003). However, the parameter optimization of these models is a very diffi cult 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 identifi cation method cannot always fi nd 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 fi rst 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 effi ciency 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.

ADVANTAGE OF USING THE MAT MODEL FOR SIMULATING THE CORTICAL CIRCUIT
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 benefi cial 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.

APPENDIX
In this study, we describe three standard neuron models namely a Hodgkin-Huxley (HH) model, leaky integrated-and-fi re (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 suffi ciently small time step of 0.001 ms. Accuracy of the numerical integration is ascertained using the backward Euler method.

Model equations
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 rectifi er" K + current is described by, where E K = −90 mV is the reversal potential of K + . A noninactivating K + current is described by, Due to the current I M the model exhibits spike frequency adaptation.

Fitting procedure
The Destexhe model is specifi ed 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 fi xed 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 defi ned the spike time of the HH model neuron as the time when the derivative of the membrane potential exceeded 10 mV/ms.

Model equations
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), If V(t) > θ, then V(t + 2 ms) = θ − 6 mV (28)

Fitting procedure
The LIF model is specifi ed by three parameters {τ m , R, θ}. As with the MAT model, we fi rst determined {τ m , R} by fi tting the voltage response to fl uctuating currents with the leaky integrator. By confi rming 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.

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

V t s I t s ds t t f
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 fi red and decays exponentially (Brandman and Nelson, 2002;Chacron et al., 2003;Geisler and Goldberg, 1966) as,