Elemental spiking neuron model for reproducing diverse firing patterns and predicting precise firing times.

In simulating realistic neuronal circuitry composed of diverse types of neurons, we need an elemental spiking neuron model that is capable of not only quantitatively reproducing spike times of biological neurons given in vivo-like fluctuating inputs, but also qualitatively representing a variety of firing responses to transient current inputs. Simplistic models based on leaky integrate-and-fire mechanisms have demonstrated the ability to adapt to biological neurons. In particular, the multi-timescale adaptive threshold (MAT) model reproduces and predicts precise spike times of regular-spiking, intrinsic-bursting, and fast-spiking neurons, under any fluctuating current; however, this model is incapable of reproducing such specific firing responses as inhibitory rebound spiking and resonate spiking. In this paper, we augment the MAT model by adding a voltage dependency term to the adaptive threshold so that the model can exhibit the full variety of firing responses to various transient current pulses while maintaining the high adaptability inherent in the original MAT model. Furthermore, with this addition, our model is actually able to better predict spike times. Despite the augmentation, the model has only four free parameters and is implementable in an efficient algorithm for large-scale simulation due to its linearity, serving as an element neuron model in the simulation of realistic neuronal circuitry.

Introduction

A mathematical model of a single-neuron that is capable of accurately reproducing diverse spiking behaviors of biological neurons is required for simulating real neuronal circuitry (Diesmann et al., 1999; Izhikevich, 2004; McIntyre et al., 2004; Markram, 2006; Brette et al., 2007; Gewaltig and Diesmann, 2007; Izhikevich and Edelman, 2008; Plesser and Diesmann, 2009; Rossant et al., 2011). The leaky integrate-and-fire (LIF) neuron models, which had been regarded as a simple toy model of a single-neuron, have significantly developed in the direction of quantitatively analyzing electrical properties of real neurons by enhancing their adaptability to data (for a review of the LIF model, see Burkitt, 2006). Spike response models (SRMs; Gerstner and van Hemmen, 1992; Gerstner, 1995) based on linear integration of input currents have been successful in not only reproducing the data but also predicting the precise spike times for new inputs (Kistler et al., 1997; Gerstner and Kistler, 2002; Jolivet et al., 2004, 2006; Kobayashi and Shinomoto, 2007). Furthermore, non-linearity has been introduced to the model in a systematic manner and its suitability has been tested in typical neurons (Fourcaud-Trocmé et al., 2003; Izhikevich, 2003; Brette and Gerstner, 2005; Badel et al., 2008).

Given the idea that a good spiking neuron model can predict neuronal firings for new fluctuating inputs, the International Competition on Quantitative Single-Neuron Modeling was organized [International Neuroinformatics Coordinating Facility (INCF), 2009]. In the competition, spike neuron models were assessed based on their quantitative performance in accurately predicting spike times (Mainen and Sejnowski, 1995; Jolivet et al., 2008; Gerstner and Naud, 2009). The simplistic integrate-and-fire models described above displayed higher performance than the complicated biophysical neuron models of Hodgkin–Huxley type (Hodgkin and Huxley, 1952). Among the winning models, the multi-timescale adaptive threshold (MAT) model (Kobayashi et al., 2009; Shinomoto, 2010) performed best in prediction. This model was one of the simplest models and is based on the linear leaky integration of input currents with a moving threshold.

The criterion for assessing the spiking neuron model is not unique; in addition to quantitative reproducibility, neuron models are expected to possess the ability to manifest various firing patterns of biological neurons, as represented by the 20 different firing responses to transient currents demonstrated by Izhikevich (2004). In this respect, the complicated biophysical Hodgkin–Huxley models are capable of representing such a variety of firing responses, whereas the simple LIF models (including the MAT model) are incapable of reproducing diverse responses. Izhikevich (2003) proposed a simplistic non-linear model that can express this diversity. More recently, Mihalas and Niebur (2009) proposed a linear integrate-and-fire model that can produce a fairly rich variety of firing responses.

In this study, we apply improvements to the MAT model so that it is able to handle the 20 qualitative firing responses introduced by Izhikevich (collectively known as Izhikevich’s table), including inhibitory rebound spiking and resonate spiking. Our improved MAT model must not negatively impact the strong adaptability and predictability that the original MAT model possesses. Here the difficulty is that any model tends to be intractable when large degrees of freedom for diversity are added. In accordance with the wisdom of Occam’s razor, we refrain from defining a large number of parameters; instead, we augment only one parameter, and then observe whether the model becomes both capable of covering various firing patterns while maintaining its strong quantitative predictive performance.

In this paper, we first provide a full analysis of the dynamics of the original MAT model, and describe our enhanced model by introducing the voltage dependent term. Next, we examine the model to determine whether it is capable of reproducing the entire set of qualitative features. Finally, we examine the augmented model to determine whether it possesses the ability to precisely predict spike times of biological neurons.

Dynamic Characteristics of the Original MAT Model

Before extending the MAT model, we first analyze the dynamic characteristics of the original MAT model, which consists of two parts: the non-resetting leaky integrator and the MAT (Kobayashi et al., 2009). Figure 1 illustrates the dynamics of the MAT model.

Figure 1

Dynamics of the multi-timescale adaptive threshold (MAT) model. The model potential V(t) (blue) is obtained from input current I(t) (green) with the non-resetting leaky integration given in Eq. 1. If the model potential hits threshold θ(t) (red), the model assumes a spike and increases the threshold. The adaptive threshold θ(t) decays toward resting value ω with multiple timescales according to Eqs 2 and 3.

The dynamics of the non-resetting leaky integrator is given by

where V(t) and I(t) are the model potential and input current, and τ and R are parameters representing the leak time constant and input resistance, respectively. The model assumes a spike when V(t) reaches a threshold value indicated by θ(t). In this model, V(t) is not reset as in the standard LIF model (Stein, 1965, 1967), but instead, threshold θ(t) is increased, and then decays toward its resting value. In standard adaptive threshold models, the spike threshold θ(t) decays with a single exponential function (Geisler and Goldberg, 1966; Brandman and Nelson, 2002; Chacron et al., 2003), but the present model contains multiple exponential decays, given by

where t is the kth spike time, L is the number of threshold time constants, τ is the jth time constant, α is the weight of the jth time constant, and ω is the resting value. In addition, to avoid singular bursting, the neuron is not allowed to fire within an absolute refractory period τ even when the membrane potential exceeds the threshold. If the membrane potential lies above the threshold after the period τ, the model assumes another spike.

The model is therefore fully specified by the parameters {τ, R, τ, τ, α, (j = 1, 2,…, L), ω}. We inherit many of the parameter values adopted in the original MAT model of L = 2, as R = 50 MΩ, τ = 2 ms, τ1 = 10 ms, and τ2 = 200 ms, which were determined using experimental data obtained by current injection experiments on cortical neurons, including regular-spiking (RS), intrinsic-bursting (IB), and fast-spiking (FS) neurons (Kobayashi et al., 2009). Here, the multiple timescales are regarded as respectively expressing different biological ionic currents, such that 10 ms represents fast transient Na+ current and delayed rectifier K+ current, and 200 ms represents non-inactivating K+ current, hyperpolarization-activated cation current, and Ca2+-dependent K+ current (Koch, 1999; Hille, 2001). The only one parameter that we modified from the original MAT model is the membrane time constant; we changed it from τ = 5 to 10 ms. This is because the time constant fitted to the present data turned out to be 11.7 ms, the original time constant 5 ms is too small compared to the widely accepted range of membrane time constant 10–20 ms (McCormick et al., 1985; Yang et al., 1996), and the spike time prediction performed with τ = 10 ms gave rise to the highest predictive performance among the choices of 5, 10, 11.7, 15, and 20 ms.

The MAT model is capable of reproducing a variety of spike responses manifested by a broad class of cortical neurons, including RS, IB, and FS neurons. Furthermore, the MAT model is capable of demonstrating repeated burst firing called “chattering” (CH; Gray and McCormick, 1996). To examine the basic capability of the model in expressing a wide variety of qualitative firing features, we first use bifurcation theory (Hoppensteadt and Izhikevich, 1997) to investigate the firing response of the original MAT model given a constant current injection.

Using bifurcation analysis, we find that the model parameter space can be divided into domains in which the model exhibits qualitatively different dynamic behavior, as represented by type I or II excitability, which is defined according to whether the frequency–current (f–I) relationship is continuous or discontinuous at the threshold current, respectively (Hodgkin, 1948; Rinzel and Ermentrout, 1989), and bursting (Izhikevich, 2007). Figure 2 depicts the allocation of the three domains on a plane of parameters α1 and α2, which is obtained by sectioning a three-dimensional parameter space at a given residual parameter ω.

Figure 2

Phase diagram in the α. The MAT model is capable of exhibiting not only type I (yellow region) and type II (blue region) excitabilities, but also repeated burst firing (red region). Mathematically unfeasible region is indicated in gray. The respective boundaries between different phases are given by α2 = 0 and Eqs 13–15 (see the text for details).

Because the MAT model is equipped with the resetting operation in addition to an ordinary differential equation, a standard bifurcation analysis cannot be directly applied to the model; however, an apt analogy can be found for its bifurcation phenomena. Figure 3A depicts the features of phase dynamics in a plane spanned by ω and α2 at a given positive value of α1. The geometric representation of the model dynamics, called the phase plane analysis, is given in the figure.

Figure 3

Bifurcation diagram of type I/II excitabilities. (A) Phase diagram in the ω–α2 plane for arbitrary positive value of α1 and geometric representation of the model dynamics. The MAT model exhibits saddle-node bifurcation and saddle-node off invariant circle bifurcation in the upper and lower parts of the phase diagram, α2 ≥ 0 and α2 < 0, respectively. (B) The corresponding dynamics of adaptive threshold θ(t). (C) The f–I curve. The firing frequency induced by a constant input current increases from zero in the type I phase α2 ≥ 0, while it exhibits a hysteresis in the type II phase, α2 < 0. The model parameters for type I: α1 = 15, α2 = 3, and ω = 5, and for type II: α1 = 15, α2 = −0.05, and ω = 5.

Type I excitability

In the upper part of the phase diagram of Figure 3A, α2 ≥ 0, H(t) is a monotonic decreasing function (Figure 3B), and the periodic firing is initiated with the saddle-node bifurcation at IR = ω. Accordingly, the frequency of the oscillation increases continuously from zero, as demonstrated by the continuous f–I curve shown in Figure 3C.

We can analyze the model dynamics by decomposing θ(t)–ω into two variables, and , and drawing the trajectory in the plane of (θ1, θ2). Since these variables satisfy differential equations dθ1/dt = −θ1/τ1 and dθ2/dt = −θ2/τ2, we can obtain the differential equation for the trajectory of (θ1, θ2) by eliminating time t as

which leads to

where τ1/τ2 = 10/200 = 0.05 for the original choice of two timescales, and c is an integration constant.

Given a current I, the MAT model generates a spike if the state (θ1, θ2) is in the “firing region,”

and if an absolute refractory period τ has expired from the last spike. When the model generates a spike, the state is shifted according to

Figure 4A shows the trajectory for the case of type I excitability. The origin (θ1, θ2) = (0, 0) remains as a stable fixed point, or a stable node, if it lies outside the firing region. If the firing region contains (θ1, θ2) = (0, 0), the stable node disappears and the repetitive firing, or the limit cycle of oscillation, starts.

Figure 4

Dynamics of type I and type II firing. Geometric representation of the firing dynamics in the plane of . (A) Type I. (B) Type II. Once the state point meets the boundary of the firing region θ1 + θ2 ≤ IR − ω (red solid lines), the state (θ1, θ2) is reset to (θ1 + α1, θ2 + α2) (green dashed lines). The state returns to the initial position along the trajectory (green solid line) given by Eq. 5(gray dashed line). The origin (θ1, θ2) = (0, 0) remains as a stable fixed point [a purple closed circle in (B)] if it is outside the firing region. The model parameters for (A): α1 = 10, α2 = 1 and IR − ω = 0.5, and for (B): α1 = 10, α2 = −0.2 and IR − ω = −0.5.

The period of limit cycle oscillation T is determined by the condition that the shifted state (θ1 + α1, θ2 + α2) will come back to the initial state (θ1, θ2) after the relaxation for a period T. This condition is represented by a set of equations

By eliminating θ1 and θ2 from these equations, we obtain the equation that determines the period of oscillation T,

The f–I curve can be obtained by solving this equation for given current I. In a limit of large T, the second term in the left-hand-side dominates, and the equation is approximated as and thus giving

Type II excitability

In the lower part of ω − α2 phase diagram of Figure 3A, α2 < 0, the system exhibits a bifurcation that corresponds to the saddle-node off invariant circle bifurcation for an ordinary differential equation system. In this parameter range, H(t) expresses a dent in response to a single spike (Figure 3B), which makes the system bistable, exhibiting both the resting state and autonomous repetitive firing. The bistability induces a hysteresis in the f–I curve (Figure 3C), with a finite frequency gap induced by the dent in H(t).

The equation for determining the frequency–current relation, Eq. 9, also holds for this case of α2 < 0. The points of difference from the case of α2 > 0 are that the trajectory giving a repetitive firing cycle can coexist with the stable fixed point (θ1, θ2) = (0, 0) (Figure 4B), and that the limit cycle trajectory disappears before the frequency f = 1/T vanishes.

Burst firing

For the parameter range of α1 < 0 and α2 > 0, Kobayashi et al. (2009) have shown that the MAT model is capable of exhibiting the burst firing. Here, we analyze the bursting in terms of the trajectory of the state (θ1, θ2). Figure 5 depicts a closed loop trajectory in which the burst firing is repeated with intermissions of the relaxation periods. In the case of α1 + α2 < 0, the occurrence of a spike rather lowers the threshold. If the state (θ1, θ2) still remains in the firing region after the absolute refractory period τ, the neuron generates another spike. The neuron repeats this resetting until the state escapes from firing region.

Figure 5

Dynamics of bursting. (A) The time dependence of dynamic threshold θ(t) and the resulting burst firing. (B) Geometric representation of the bursting dynamics in the plane of . Once the state point meets the boundary of the firing region (a red solid line), the system starts to repeat spiking (black crosses), resetting (red dashed lines) and short relaxation for τ = 2 ms (green curves). After escaping from the firing region (an open square), the system enters relaxation period (a gray dashed line). The model parameters: α1 = −0.8, α2 = 0.5 and IR − ω = 0.5.

The burst firing of N spikes brings the state point (θ1, θ2) to

For instance, the burst firing of two spikes starts with the condition of

This condition determines the boundary between the type I and burst firing regions. In the case of IR − ω = 0 as in Figure 2, this is

Unfeasible region

In the burst firing region, the number of spikes contained in each burst increases by decreasing α1. By decreasing α1 further, the system enters an unfeasible region in which the system cannot escape from the firing region. The boundary of this unfeasible region is given by taking the limit N → ∞ in Eq. 11and comparing it with firing condition of Eq. 6. This is given by

In the type II firing region, the firing frequency increases by decreasing α2. Because the firing frequency should be bounded by 1/τ due to the refractory period, the condition of the firing frequency reaching 1/τ determines the boundary between the type II and unfeasible region. In the case of IR − ω = 0, the condition becomes

In the limit of τ2 > τ1 ≫ τ, the boundary between the unfeasible region and the burst firing region (Eq. 14) and the boundary between the unfeasible region and the type II region (Eq. 15) asymptotically approach to a single straight line, α2 = −(τ1/τ2)α1.

Augmentation of the MAT Model with Voltage Dependency

We have shown above that the MAT model is capable of exhibiting basic dynamic characteristics, such as type I/II excitability and burst firing, but we show in the next section that the original model is not sufficiently flexible to exhibit a richer variety of responses of biological neurons to transient input currents, summarized as Izhikevich’s table (Izhikevich, 2004).

The reason that the original MAT model cannot represent a variety of firing responses to dynamic input current is that the adaptive threshold does not depend on the present state of voltage. Therefore, we suggest extending the adaptive threshold dynamics of the MAT model by adding a term representing the voltage dependency of the form

where β is a parameter for controlling the voltage dependency, and K(s) is an α-function kernel of timescale τ,

This dependency of the spiking threshold on the derivative of membrane voltage was actually observed in biological neurons (Azouz and Gray, 2000) and is considered to have resulted from activation/inactivation dynamics of voltage-gated sodium channels (Naundorf et al., 2006; Platkiewicz and Brette, 2011), or the backpropagation of axonal spike to the soma (Yu et al., 2008). The timescale τ of the influence of dV/dt dependency is estimated as a few ms (Azouz and Gray, 2000); we fixed this time constant at 5 ms. Thus, the tunable parameter for this voltage dependency is only the coefficient β.

Due to the linear kernel integration with the alpha function kernel, the numerical integration of this model can be expedited by rewriting the evolution equation into a set of differential equations. By rewriting the evolution equation in this way, it is possible to implement the algorithm of exact subthreshold integration (Morrison et al., 2007; Plesser and Diesmann, 2009), which we discuss in the Section “Appendix.”

Variety of Firing Responses Realized by the Original and Augmented MAT Models

In this section, we consider numerous transient input current types, including step currents, ramp currents, and a set of short current pulses, and we examine whether the original and augmented MAT models are capable of handling the 20 qualitative firing responses defined by Izhikevich (2004).

Original MAT model

We found that the original MAT model is capable of reproducing 9 of the 20 basic firing responses, as shown in Figures 6A–I. More specifically, the model produces basic tonic spiking in response to a step current (Figure 6A), frequency adaptation (Figure 6B), integrator (Figure 6C), class 1 (type I) excitability (Figure 6D), class 2 (type II) excitability (Figure 6E), and bistability (Figure 6F). The terms class 1/2 are used in the original paper by Hodgkin (1948), whereas they are also referred to as type I/II in terms of the connection to bifurcations (Rinzel and Ermentrout, 1989). For the type II parameter region shown in Figure 2, the threshold may exhibit a dent after a spike, which can be interpreted as depolarizing after-potential that helps the occurrence of a successive spike (Figure 6G). In the region of bursting shown in Figure 2, the MAT model exhibits tonic bursting to a step current injection (Figure 6H). Given parameters near the boundary between the bursting and type I regions of Figure 2, the model shows burst firing only at the onset of a step current injection, and then starts tonic spiking (Figure 6I).

Figure 6

A variety of firing responses manifested by the original MAT model. (A) basic tonic spiking; (B) frequency adaptation; (C) integrator; (D) class 1 excitability; (E) class 2 excitability; (F) bistability; (G) depolarizing after-potential; (H) tonic bursting; and (I) mixed mode. The blue bars indicate the generated spikes. The red, blue and green traces represent the adaptive threshold, the model potential, and the input current, respectively.

Augmented MAT model

The original MAT model cannot complete the entire set of firing responses, because the dynamics of the adaptive threshold is only dependent on past spikes and does not depend on the present state of the membrane voltage. To enable the MAT model to respond variably to time-dependent input current, we have augmented the adaptive threshold dynamics of the model by adding a term representing the voltage dependency (using Eqs 16 and 17 above). Figures 7A–K shows 11 firing patterns that are successfully realized by our augmented MAT model; each firing pattern is listed below. Due to the sensitivity to dV/dt, the augmented model can exhibit phasic spiking or bursting; more specifically, the neuron emits a spike or a burst of spikes only at the onset of a step current injection (Figures 7A,B), exhibits delayed response to short pluses (Figure 7C), shows rebound spiking or bursting (Figures 7D,E), and exhibits threshold variability (Figure 7F).

Figure 7

A variety of firing responses manifested by the augmented MAT model. (A) phasic spiking; (B) phasic bursting; (C) latency; (D) rebound spiking; (E) rebound bursting; (F) threshold variability; (G) subthreshold oscillations; (H) resonator; (I) accommodation; (J) inhibition-induced spiking; and (K) inhibition-induced bursting. Note that this entire set of firing responses follows definitions in Izhikevich (2004). The blue bars indicate the generated spikes. The red, blue, and green traces represent the adaptive threshold, the model potential, and the input current, respectively. The gray and orange traces are the value of dV/dt and the voltage dependent term, respectively.

Furthermore, our augmented MAT model can mimic subthreshold oscillation in which the dynamic threshold follows the increases and decreases in voltage (Figure 7G). According to the transient oscillation mimicked by the dynamic threshold, the model may behave as a resonator, firing in response to a pair of input pulses with a particular interval (Figure 7H); for detailed conditions of such a resonator, please see the Section “Appendix.” The dV/dt dependency induces an effect of accommodation (Figure 7I). With our augmented model, the neuron is also able to elicit spikes or bursts in response to inhibitory input, an approach called inhibition-induced spiking or bursting (Figures 7J,K).

The model parameters and input conditions that were used for reproducing these responses are summarized in Table 1. The locations of these model parameters in the space of α1, α2, and β are depicted in Figure 8. Some of the firing responses, such as class 1 vs class 2, or integrator vs resonator are mutually exclusive by definition. However, there are cases in which multiple firing response types can be realized with the same set of model parameters solely by changing the input conditions. Such degeneracy occurs, for instance, between phasic spiking, latency, threshold variability, accommodation, and rebound spiking, or between phasic bursting, rebound bursting, and tonic bursting.

Table 1

Model parameters and input currents used for reproducing firing patterns in Figures .

Firing patterns	Model parameters				Input currents
	α1	α2	β	ω
Tonic spiking	10	0		5	Ic = 0.15 nA
Adaptation	10	1		5	Ic = 0.15 nA
Integrator	10	0		5	Ip = 0.28 nA, Δp = 2 ms, IPIs = 8, 118, 13 ms
Class 1 excitable	15	3		5	dI/dt = 2.5 pA/s
Class 2 excitable	15	−0.05		5	dI/dt = 2.5 pA/s
Bistability	20	−0.4		5	Ic = 0.95 nA, Ip = + 0.20, −0.04 nA, Δp = 3, 15 ms, IPI = 22 ms
Depolarizing after-potential	25	−1		5	Ip = 0.2 nA, Δp = 5 ms
Tonic bursting	−0.5	0.35		5	Ic = 0.15 nA
Mixed mode	−0.8	0.7		5	Ic = 0.15 nA
Phasic spiking	10	0	−0.3	5	Ic = 0.08 nA
Phasic bursting	−0.5	0.35	−0.3	5	Ic = 0.08 nA
Latency	10	0	−1	5	Ip = 0.58 nA, Δp = 0.5 ms
Rebound spiking	10	0	−2.5	5	Ip = −0.6 nA, Δp = 1 ms
Rebound bursting	−0.5	0.35	−2.5	5	Ip = −0.6 nA, Δp = 1 ms
Threshold variability	10	0	−0.1	5	Ip = 0.2, −0.2, 0.2 nA, Δp = 2 ms, IPIs = 118, 13 ms
Subthreshold oscillations	10	0	0.1	5	Ip = 0.2 nA, Δp = 2 ms
Resonator	10	0	0.1	5	Ip = 0.36 nA, Δp = 2 ms, IPIs = 8, 118, 13 ms
Accommodation	10	0	−0.5	5	dI/dt = 1, 4.5 nA/s, Δramp = 90, 20 ms
Inhibition-induced spiking	20	0	2	5	Ip = −0.30 nA, Δp = 40 ms
Inhibition-induced bursting	−0.5	0.35	2	5	Ip = −0.16 nA, Δp = 60 ms

The top 9 and the bottom 11 firing patterns are the ones represented by the original MAT model (β = 0) and augmented MAT model (β ≠ 0), respectively. .

Figure 8

Distribution of model parameters that realized a variety of firing patterns. (A) Parameter values (α1, α2) of the original MAT model, which are used to reproduce nine firing patterns in Figure 6. The colors of the firing type regions are the same as Figure 2. (B) Parameter values (α1, α2, β) of the augmented MAT model, which are used to reproduce 11 firing patterns in Figure 7. The numerical values of the respective model parameters are summarized in Table 1.

Reproducing the quantitative features of biological neuronal firing

It was reported that the original MAT model can reproduce the qualitative features of four representative firing types, FS, RS, IB, and CH neurons (Kobayashi et al., 2009). However, the model does not necessarily account for them quantitatively. To be precise, the original MAT model is capable of reproducing FS, RS, and CH firing, but the model should be augmented with the voltage dependent term in order to reproduce IB firing in a range of biologically acceptable timescales. Here we demonstrate the typical parameter setting of the model in reproducing the biological neuronal firing patterns.

FS neurons

They fire high-frequency tonic spikes with little or no frequency adaptation. The original MAT model (β = 0) is capable of reproducing the tonic firing of 200 Hz in response to a rectangular current of 0.6 nA (McCormick et al., 1985) with the model parameters of α1 = 10, α2 = 0, and ω = 15. Note that even the multiple timescales are not necessarily needed for this case, and the simple LIF model is enough to reproduce the phenomena.

RS neurons

They fire tonic spikes with adapting frequency in response to rectangular current, and have class 1 excitability. This can also be represented by the original MAT model. For instance, the model is capable of reproducing the situation in which the neuron is injected a rectangular current of 0.6 nA and exhibits initial firing frequency of 120 Hz as defined from the first interspike interval and adapts the firing frequency to 30 Hz (McCormick et al., 1985) with the parameters of α1 = 20, α2 = 2, and ω = 20. Having two timescales is necessary for reproducing the frequency adaptation.

CH neurons

They fire high-frequency bursts of a few (<5) spikes with relatively constant inter-burst interval of 15–50 ms (Gray and McCormick, 1996). When reproducing CH firing by the original MAT model, the inter-burst interval in the tonic bursting is determined by the relaxation of the slow component. The MAT model can reproduce tonic bursting of inter-burst frequency of 20 Hz and intra-burst frequency of 500 Hz in response to a rectangular current of 0.6 nA, respectively (Nowak et al., 2003), with the parameters of α1 = −2.5, α2 = 2, and ω = 28.

IB neurons

They generate a burst of three to five spikes at the beginning of a strong depolarizing pulse of current, and then switch to tonic spiking. Intra-burst frequency is low compared to CH neurons (Gray and McCormick, 1996; Nowak et al., 2003). Though we demonstrated that the original MAT model is capable of representing the mixed mode (Figure 6I), the intra-burst frequency in this case becomes 1/τ = 500 Hz, which is outside the biologically plausible range (≤350 Hz; Nowak et al., 2003). Thus the MAT model has to be augmented with the voltage dependent term to generate the depolarizing wave evoked by the current injection (McCormick et al., 1985). The augmented MAT model can reproduce the IB firing with frequencies of the burst and tonic modes of 300 Hz and 30 Hz, respectively (Nowak et al., 2003) with the model parameters of α1 = 9, α2 = 0.3, β = −0.3, and ω = 28.

Other types of neurons

In addition to those four representative firing types, there are more abundant firing responses (Kawaguchi, 1995; Kawaguchi and Kondo, 2002; Markram et al., 2004; Izhikevich, 2007; Ascoli et al., 2008). In the following, we show the parameter values with which the augmented MAT model is capable of reproducing them.

The low-threshold spiking (LTS) neurons respond in a manner identical to RS neurons to a prolonged suprarheobasic current injection at depolarized potentials, and respond to current injections just greater than rheobase with a transient bursting of two or more spikes with short (<20 ms) ISIs at hyperpolarized potentials. The augmented MAT model is capable of reproducing the response with the parameters of α1 = 10, α2 = 1, β = −0.2, and ω = 25.

Phasic spiking or phasic bursting was observed from inhibitory interneurons in the rat frontal cortex (Kawaguchi and Kondo, 2002). Our augmented MAT model can reproduce phasic spiking and bursting in response to a rectangular current of 0.6 nA with the parameters of α1 = 20, α2 = 1, β = −0.3, and ω = 35 and α1 = −2.2, α2 = 2, β = −0.3, and ω = 35, respectively. By increasing current intensity, the model with these sets of parameters can generate tonic spiking or bursting.

The thalamocortical (TC) neurons are known to possess firing regimes of tonic spiking and rebound bursting as in LTS neurons, but their adaptation is smaller than LTS neurons (Llinás and Jahnsen, 1982; Destexhe and Sejnowski, 2001). For instance, our augmented MAT model can reproduce such a response with the parameters of α1 = 10, α2 = 0, β = −0.2, and ω = 25. Although inhibition-induced spiking or bursting (Figures 7J,K) is a specific feature of TC neurons (Izhikevich, 2003), we have not found a set of model parameters with which the model can realize both the inhibition-induced firing and rebound bursting (see Table 1).

In contrast to TC neurons, a certain type of thalamic interneuron does not have a burst regime, though they can fire rebound spikes (Pape and McCormick, 1995). Except for rebound spiking, this type of thalamic interneuron can fire high-frequency tonic spikes without pronounced frequency adaptation, like FS neurons. Our augmented MAT model can reproduce such a response with the parameters of α1 = 5, α2 = 0, β = −0.1, and ω = 25. Note that there is a different type of interneuron in thalamus that generates robust action potential bursts after an inhibitory input (Bal and McCormick, 1993), which our model cannot reproduce.

Reproducing and Predicting Biological Neuronal Responses to Fluctuating Currents

In the preceding section, we showed that our augmented MAT model is capable of reproducing rebound spiking and a variety of responses induced by threshold variability. As the number of parameters increases, in general, it is likely to become less capable of accommodating the numerical data. In this section, we examine whether the augmented MAT model still has the ability to predict biological spike times through the augmentation of voltage dependency.

In our examination, we use publicly available data from the International Competition on Quantitative Single-Neuron Modeling 2009 [Gerstner and Naud, 2009; International Neuroinformatics Coordinating Facility (INCF), 2009] in which the original MAT model demonstrated the ability to predict spikes, and therefore won first place in the competition. The publicly available data consisted of in vivo-like fluctuating input current and the voltage responses of 13 trials recorded from a L5 pyramidal neuron.

We adjusted the model parameters to match the data of 10 randomly selected trials, and then validated the ability of our augmented MAT model to predict the remaining 3 trials. We evaluated average predictive ability by repeating this random sampling process 100 times. The adaptive capability of the model in predicting spike times was assessed using a coefficient measuring the fraction of spikes that coincided in 4 ms between a model spike train and a real biological spike train. Among various measures (including Victor and Purpura, 1996; Tsubo et al., 2004), we adopted coincidence factor Γ as introduced in the competition (Jolivet et al., 2004, 2006, 2008; Kobayashi and Shinomoto, 2007).

Based on the above procedure, the predictive performance of the original MAT model was 0.77, whereas the augmented MAT model was 0.84, as summarized in Table 2. We also compared the performance of these models using the different coincidence windows ranging from 2 to 10 ms, and confirmed the consistent superiority of the augmented MAT model to the original model. Clearly, the augmentation does not cause predictive ability to deteriorate, but rather improves it. Figure 9 demonstrates the spike prediction of the original and augmented MAT models. From the sampling shown in the figure, both models attain similar predictive performance in the stationary period of high input current (17.5–18 and 20–21 s); however, the augmented MAT model is more proficient than the original model in predicting spikes in the period of weak input current (18.8–19.5 and 21–21.5 s) by responding to small peaks of model potential.

Table 2

Optimized parameters of the original and augmented MAT models fitted to the biological data and their fitting and predictive performances.

	Original MAT	Augmented MAT
α1	180	180
α2	1	3
β	None	0.2
ω	21	15
Fitting performance Γ	0.80	0.84
Predictive performance Γ	0.77	0.84

Figure 9

Model prediction for spiking response of a L5 pyramidal neuron to fluctuating current. The top column is a voltage trace of the experimental data obtained from an L5 neuron (data #9, Challenge A in the International Competition on Quantitative Single-Neuron Modeling 2009). The middle and bottom columns are spike times predicted by the original and augmented MAT models. The model potential and adaptive threshold are depicted in blue and red, respectively. Spikes are asterisked in blue or red, according to whether they coincided with experimental data within 4 ms. Dotted lines represent experimental spike times.

Discussion

Numerical simulations of realistic neuronal microcircuitry have been primarily based on the Hodgkin–Huxley model (Traub et al., 2005; Markram, 2006) or the LIF model (Diesmann et al., 1999; Gewaltig and Diesmann, 2007). As noted by Izhikevich (2004), these models typically have a conflict between biological reality and computational efficiency. The non-linear model proposed by Izhikevich (2004) and the linear model proposed by Mihalas and Niebur (2009) would be exceptional models that evade such a conflict.

The original MAT model possesses the strongest ability to reproduce and predict precise spike times, but is unable to reproduce the full variety of firing responses demonstrated by Izhikevich (2004). In this paper, we have proposed a revision of the MAT model in which the full variety of firing responses given by Izhikevich (2004) is successfully reproduced. We thoroughly analyzed the dynamic characteristics of the original model and augmented model, the latter including a voltage dependency term. We also showed that our augmented model does not decrease its predictive ability in comparison to the original MAT model, but rather improves quantitative predictive ability.

Voltage dependency of spiking threshold

The variability of spike threshold has been observed in numerous physiological experiments (Azouz and Gray, 2000; Ferragamo and Oertel, 2002; de Polavieja et al., 2005; Wilent and Contreras, 2005; Naundorf et al., 2006), and Hodgkin–Huxley-style models that can reproduce the threshold variability were proposed (Azouz and Gray, 2003; Naundorf et al., 2006; McCormick et al., 2007; Colwell and Brenner, 2009). More recently, Platkiewicz and Brette (2010) deduced a threshold equation that explained the manner in which the spike threshold varies with the membrane potential, depending on ionic channel properties.

There are also simplified models that introduce voltage dependency into the threshold in the LIF model (originally in Hill, 1936). Dodla et al. (2006) proposed a model that explains post-inhibitory facilitation by introducing the threshold depending on the instantaneous value of voltage with mono-exponential decay. Mihalas and Niebur (2009) adopted this form of the voltage dependency and demonstrated that the model was able to reproduce a wide variety of firing responses to various transient current inputs. Experimental research by Azouz and Gray (2000) indicated that the spike threshold related inversely with the time derivative of membrane potential. We incorporated this relationship straightforwardly into the threshold dynamics of the original MAT model while preserving its linearity, succeeding in reproducing the variety of firing patterns in Izhikevich’s table.

It would be worth noticing that the origin of the threshold variability is still in controversy. A recent study suggested that the threshold variability seen in experiments could be an artifact (Yu et al., 2008); because spikes are initiated not in the soma but in the axon initial segment, the threshold variability of the firing seen in intracellular recordings can be largely attributed by its backpropagation to the soma. Contrariwise, Platkiewicz and Brette (2011) asserted that the backpropagation cannot be the dominant cause of the threshold variability; they suggested that sodium channel inactivation is a strong candidate for the mechanism responsible for the threshold variability.

Comparison with other simplistic spike neuron models

To end our discussion, we compare our augmented MAT model with other eminent simplistic spike neuron models. The Izhikevich model is a concise non-linear model that has succeeded in demonstrating a diversity of firing responses by controlling only a few parameters. The adaptive exponential integrate-and-fire (AdEx) model (Brette and Gerstner, 2005), which combines features of the exponential integrate-and-fire model (Fourcaud-Trocmé et al., 2003) with the two-variable Izhikevich model, is also able to account for a large variety of firing patterns (Naud et al., 2008; Touboul and Brette, 2009). In comparison with these models, a key strength of our augmented MAT model is its linearity, which makes it straightforward to efficiently simulate a large network of spiking neurons using the method of exact subthreshold integration (Morrison et al., 2007; Plesser and Diesmann, 2009).

The Mihalas–Niebur model is a tractable linear model that describes burst firing in terms of after-potentials in a more realistic manner than that of the MAT model, which mechanistically repeats the regular firing given by a prefixed absolute refractory period τ = 2 ms. The Mihalas–Niebur model is capable of demonstrating most firing responses, except for resonate firing in Izhikevich’s table. The reason that the augmented MAT model is able to reproduce resonate firing is that the alpha function kernel naturally induces a delay in the voltage dependency.

The original MAT model can be regarded as a subtype of the SRM, which may describe arbitrary dependency of neuronal excitability on previous spiking history within a linear framework (Gerstner and Kistler, 2002). With an augmentation of voltage dependency in the threshold equation, the augmented MAT model can be considered as one expression of the Mihalas–Niebur model. Indeed, the SRM or the Mihalas–Niebur model has a great ability to adapt to experimental data due to its large degree of freedom. However, such a high adaptability of the model sometimes causes over-fitting to the finite data, and as a result, reduces its predictive performance. Thus our model, which has a restricted number of free parameters, can be more competent in simulating neuronal circuitry faithfully, given practically available data.

Conflict of Interest Statement

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.