IEEE TRANSACTIONS ON NEURAL NETWORKS, VOL. 14, NO. 6, NOVEMBER 2003

1569

Simple Model of Spiking Neurons
Eugene M. Izhikevich
Abstract—A model is presented that reproduces spiking and bursting
behavior of known types of cortical neurons. The model combines the biologically plausibility of Hodgkin–Huxley-type dynamics and the computational efficiency of integrate-and-fire neurons. Using this model, one cansimulate tens of thousands of spiking cortical neurons in real time (1 ms
resolution) using a desktop PC.
Index Terms—Bursting, cortex, Hodgkin–Huxley, PCNN, quadratic integrate-and-fire, spiking, thalamus.

I. INTRODUCTION
To understand how the brain works, we need to combine experimental studies of animal and human nervous systems with numericalsimulation of large-scale brain models. As we develop such large-scale
brain models consisting of spiking neurons, we must find compromises
between two seemingly mutually exclusive requirements: The model
for a single neuron must be: 1) computationally simple, yet 2) capableof producing rich firing patterns exhibited by real biological neurons.
Using biophysically accurate Hodgkin–Huxley-type models is computationally prohibitive, since we can simulate only a handful of neurons
in real time. In contrast, using an integrate-and-fire model is computationally effective, but the model is unrealistically simple and incapableof producing rich spiking and bursting dynamics exhibited by corticalneurons.
In this paper, a simple spiking model (1), (2) is presented that is as
biologically plausible as the Hodgkin–Huxley model, yet as computationally efficient as the integrate-and-fire model. Depending on four parameters, the model reproduces spiking and bursting behavior of known
types of cortical neurons, as we illustrate in Fig. 1 and summarize in
Fig. 2.
Mathematical analysis of the model will be published in the monograph by Izhikevich [8]. The derivation of the first (1) is based on bifurcation theory and normal form reduction [2], [5], and the part v 0 =2
v + I is sometimes referred to as being a quadratic integrate-and-fireneuron. The full model was first published in [10, eqns. (4) and (5)with voltage reset discussed in Sect. 2.3.1 ] in a trigonometric form
more suitable for mathematical analysis. The form presented here ismore suitable for large-scale simulations.
II. THE MODEL
Bifurcation methodologies [8] enable us to reduce many biophysically accurate Hodgkin–Huxley-type neuronal models to a two-dimensional (2-D) system of ordinary differential equations of the form

20

= 0:04v + 5v + 1400= a(bv

v
u

0 )

0 +
u

I

(1)
(2)

u

with the auxiliary after-spike resetting
if v

 30 mV then
;

v

c

u

u + d:(3)

Here, v and u are dimensionless variables, and a, b, c, and d are dimensionless parameters, and 0 = d=dt, where t is the time. The variable
Manuscript received February 2003; revised April 2003. This work was supported in part by the Neurosciences Research Foundation and by a Grant fromthe Alafi Family Foundation.
The author is with The Neurosciences Institute, San Diego, CA 92121 USA(e-mail: Eugene.Izhikevich@nsi.edu).
Digital Object Identifier 10.1109/TNN.2003.820440

Fig. 1. The simple model (1), (2) can reproduce firing patters of neurons
recorded from the rat’s motor cortex. Data are kindly shared by N. Desai,model parameters as in Fig. 2.

represents the membrane potential of the neuron and u represents amembrane recovery variable, which accounts for the activation of K+
ionic currents and inactivation of Na+ ionic currents, and it provides
negative feedback to v . After the spike reaches its apex (+30 mV), themembrane voltage and the recovery variable are reset according to the(3). Synaptic currents or injected dc-currents are delivered via the variable I .
The part 0:04v 2 +5v +140 was obtained by fitting the spike initiation
dynamics of a cortical neuron (other choices also feasible) so that themembrane potential v has mV scale and the time t has ms scale. The
resting potential in the model is between 070 and 060 mV dependingon the value of b. As most real neurons, the model does not have a fixed
threshold; Depending on the history of the membrane potential prior to
the spike, the threshold potential can be as low as 055 mV or as highas 040 mV.
• The parameter a describes the time scale of the recovery variable
u. Smaller values result in slower recovery. A typical value isa = 0:02.
• The parameter b describes the sensitivity of the recovery variable
u to the subthreshold fluctuations of the membrane potential v .
Greater values couple v and u more strongly resulting in possible
subthreshold oscillations and low-threshold spiking dynamics. A
typical value is b = 0:2. The case b < a(b > a) corresponds
to saddle-node (Andronov–Hopf) bifurcation of the resting state[10].
• The parameter c describes the after-spike reset value of the membrane potential v caused by the fast high-threshold K+ conductances. A typical value is c = 065 mV.
v

1045-9227/03$17.00 © 2003 IEEE

1570

IEEE TRANSACTIONS ON NEURAL NETWORKS, VOL. 14, NO. 6, NOVEMBER 2003

Fig. 2. Known types of neurons correspond to different values of the parameters a, b, c, d in the model described by the (1), (2). RS, IB, and CH are corticalexcitatory neurons. FS and LTS are cortical inhibitory interneurons. Each inset shows a voltage response of the model neuron to a step of dc-current I = 10(bottom). Time resolution is 0.1 ms. This figure is reproduced with permission from www.izhikevich.com. (Electronic version of the figure and reproduction
permissions are freely available at www.izhikevich.com.)• The parameter d describes after-spike reset of the recovery variable u caused by slow high-threshold Na+ and K+ conductances.
A typical value is d = 2.
Various choices of the parameters result in various intrinsic firing patterns, including those exhibited by the known types of neocortical [1],[3], [4] and thalamic neurons as summarized in Fig. 2. A possible extension of the model (1), (2) is to treat u, a and b as vectors, and use u
instead of u in the voltage (1). This accounts for slow conductanceswith multiple time scales, but we find such an extension unnecessarily
for cortical neurons.

• IB (intrinsically bursting) neurons fire a stereotypical burst of
spikes followed by repetitive single spikes (Fig. 2IB). In themodel, this corresponds to c = 055 mV (high voltage reset)and d = 4 (large after-spike jump of u). During the initial burst,
variable u builds up and eventually switches the dynamics frombursting to spiking.
• CH (chattering) neurons can fire stereotypical bursts of closely
spaced spikes. The inter-burst frequency can be as high as 40 Hz.
In the model, this corresponds to c = 050 mV (very high voltage
reset) and d = 2 (moderate after-spike jump of u).
All inhibitory cortical cells are divided into the following two classes[4]:III. DIFFERENT TYPES OF DYNAMICS
Neocortical neurons in the mammalian brain can be classified into
several types according to the pattern of spiking and bursting seen in
intracellular recordings. All excitatory cortical cells are divided into thefollowing four classes [1], [3]:• RS (regular spiking) neurons are the most typical neurons in thecortex. When presented with a prolonged stimulus (injected stepof dc-current in Fig. 2RS, bottom) the neurons fire a few spikes
with short interspike period and then the period increases. This iscalled the spike frequency adaptation. Increasing the strength of
the injected dc-current increases the interspike frequency, though
it never becomes too fast because of large spike-afterhyperpolarizations. In the model, this corresponds to c = 065 mV (deepvoltage reset) and d = 8 (large after-spike jump of u).

• FS (fast spiking) neurons can fire periodic trains of action potentials with extremely high frequency practically without any adaptation (slowing down), as one can see in Fig. 2FS. In the model,
this corresponds to a = 0:1 (fast recovery).
• LTS (low-threshold spiking) neurons can also fire high-frequency
trains of action potentials (Fig. 2LTS), but with a noticeable spikefrequency adaptation. These neurons have low firing thresholds,
which is accounted for by b = 0:25 in the model. To achieve a
better quantitative fit with real LTS neurons, other parameters ofthe model need to be changed as well.
In addition, our model can easily reproduce behavior of thalamo-cortical neurons, which provide the major input to the cortex• TC (thalamo-cortical) neurons have two firing regimes: Whenat rest (v is around 060 mV) and then depolarized, they exhibit

IEEE TRANSACTIONS ON NEURAL NETWORKS, VOL. 14, NO. 6, NOVEMBER 2003

1571

we choose the ratio of excitatory to inhibitory neurons to be 4 to 1, andwe make inhibitory synaptic connections stronger. Besides the synaptic
input, each neuron receives a noisy thalamic input.
In principle, one can use RS cells to model all excitatory neurons
and FS cells to model all inhibitory neurons. The best way to achieve
heterogeneity (so that different neurons have different dynamics), is
to assign each excitatory cell (ai ; bi ) = (0:02; 0:2) and (ci ; di ) =2(065; 8) + (15; 06)ri , where ri is a random variable uniformly distributed on the interval [0,1], and i is the neuron index. Thus, ri = 0
corresponds to regular spiking (RS) cell, and ri = 1 corresponds tothe chattering (CH) cell. We use ri2 to bias the distribution toward RS
cells. Similarly, each inhibitory cell has (ai ; bi ) = (0:02; 0:25) +(0:08; 00:05)ri and (ci ; di ) = (065; 2).
The model belongs to the class of pulse-coupled neural networks(PCNN): The synaptic connection weights between the neurons are
given by the matrix S = (sij ), so that firing of the j th neuron instantaneously changes variable vi by sij .

Fig. 3. Simulation of a network of 1000 randomly coupled spiking neurons.
Top: spike raster shows episodes of alpha and gamma band rhythms (verticallines). Bottom: typical spiking activity of an excitatory neuron. All spikes wereequalized at 30 mV by resetting v first to 30 mV and then to c.

++

tonic firing as in Fig. 2TC, left voltage trace. However, if a negative current step is delivered so that the membrane potential is
hyperpolarized (v is around 090 mV), the neurons fire a reboundburst of action potentials, as in Fig. 2TC, right voltage trace.
The model can exhibit other interesting types of dynamics.
• RZ (resonator) neurons have damped or sustained subthresholdoscillations, as in Fig. 2RZ. They resonate to rhythmic inputs
having appropriate frequency (as the resonate-and-fire model[9]). This behavior corresponds to a = 0:1 and b = 0:26.
Notice that there is a bistability of resting and repetitive spiking
states: The neuron can be switched between the states by anappropriately timed brief stimuli.
Dynamics of other neuronal types, including those in brainstem, hippocampus, basal ganglia, and olfactory bulb, can also be described byour model.
Our “one-fits-all” choice of the function 0:04v2 + 5v + 140 in (1) isjustified when large-scale networks of spiking neurons are simulated,
as we discuss below. However, if one is interested in the behavior ofa single neuron, then other choices of the function are available, andsometimes more preferable. For example, the function 0:04v 2 +4:1v +108 with b = 00:1 is a better choice for the RS neuron, since it leads tothe saddle-node on invariant circle bifurcation and Class 1 excitability[10].
IV. PULSE-COUPLED IMPLEMENTATION
We have used this model to simulate a sparse network of 10 000
spiking cortical neurons with 1 000 000 synaptic connections in realtime (resolution 1 ms) using a 1 GHz desktop PC and C++ programming language. The following MATLAB program (also available on
author’s webpage) simulates a network of randomly connected 1000
neurons in real time. Motivated by the anatomy of a mammalian cortex,% Created by Eugene M.
% Excitatory neurons
Ne=800;re=rand(Ne,1);
a=[0.02*ones(Ne,1);b=[0.2*ones(Ne,1);c=[-65+15*re.^2;
d=[8-6*re.^2;
S=[0.5*rand(Ne+Ni,Ne),

Izhikevich, February 25, 2003
Inhibitory neurons
Ni=200;ri=rand(Ni,1);
0.02+0.08*ri];
0.25-0.05*ri];-65*ones(Ni,1)];
2*ones(Ni,1)];-rand(Ne+Ni,Ni)];

v=-65*ones(Ne+Ni,1);
u=b.*v;firings=[];% Initial values of v
% Initial values of u
% spike timings

for t=1:1000% simulation of 1000 ms
I=[5*randn(Ne,1);2*randn(Ni,1)]; % thalamic input
fired=find(v>=30);% indices of spikes
firings=[firings; t+0*fired,fired];v(fired)=c(fired);u(fired)=u(fired)+d(fired);
I=I+sum(S(:,fired),2);v=v+0.5*(0.04*v.^2+5*v+140-u+I); % step 0.5 ms
v=v+0.5*(0.04*v.^2+5*v+140-u+I); % for numerical
u=u+a.*(b.*v-u);% stabilityend;plot(firings(:,1),firings(:,2),’.’);

One can see in Fig. 3 that the network exhibits cortical-like asynchronous dynamics; that is, neurons fire Poisson spike trains with mean
firing rates around 8 Hz. Dark vertical lines indicate that there are occasional episodes of synchronized firings in the alpha and gamma frequency range (10 and 40 Hz, respectively). Although the network isconnected randomly and there is no synaptic plasticity, the neurons
self-organize into assemblies and they exhibit collective rhythmic behavior in the frequency range corresponding to that of the mammaliancortex in the awake state. Changing the relative strength of synaptic
connections and the strength of the thalamic drive can produce othertypes of collective behavior, including spindle waves and sleep oscillations. We can easily observe and study these cortical states because oursimple spiking model describes accurately dynamics of known typesof cortical neurons. Thus, there is no longer a contradiction between
biological plausibility and computational efficiency of model neuralnetworks.

1572

IEEE TRANSACTIONS ON NEURAL NETWORKS, VOL. 14, NO. 6, NOVEMBER 2003

V. CONCLUSION

Global and Partial Synchronism inPhase-Locked Loop Networks

In this paper, a simple model that reproduces the rich behavior of
biological neurons, including spiking, bursting, and mixed mode firing
patters, post-inhibitory (rebound) spikes and bursts, continuous spiking
with frequency adaptation, spike threshold variability, bistability ofresting and spiking states, and subthreshold oscillations and resonanceis presented (the latter are discussed in [7] and [9]).
Our model is the simplest possible model that can reproduce thesetypes of neuronal behavior: It consists of only two equations and has
only one nonlinear term, i.e., v 2 . Yet, the model is canonical in thesense that the difference between it and a whole class of biophysicallydetailed and accurate Hodgkin–Huxley-type models, including thoseconsisting of enormous number of equations and taking into account all
possible information about ionic currents, is just a matter of coordinatechange [6].
We show how to use the model to build networks of spiking neurons
capable of exhibiting collective dynamics and rhythms similar to thoseof the mammalian cortex. Due to the extreme computational simplicityof the model, we can simulate thalamo-cortical networks consisting oftens of thousands of spiking neurons in real time with 1 ms resolution
using an old 1-GHz desktop PC.

L. H. A. Monteiro, N. C. F. Canto, J. G. Chaui-Berlinck, F. M. Orsatti,
and J. R. C. Piqueira

ACKNOWLEDGMENT
The author thanks F. C. Hoppensteadt, D. Golomb, R. Pichevar, E.
W. Gall, A.l K. Seth, J. R. Iversen, and Y. Chen, who read the first draftof the manuscript and made a number of useful suggestions.

REFERENCES[1] B. W. Connors and M. J. Gutnick, “Intrinsic firing patterns of diverse
neocortical neurons,” Trends in Neurosci., vol. 13, pp. 99–104, 1990.
[2] G. B. Ermentrout and N. Kopell, “Parabolic bursting in an excitablesystem coupled with a slow oscillation,” SIAM J. Appl. Math., vol. 46,
pp. 233–253, 1986.
[3] C. M. Gray and D. A. McCormick, “Chattering cells: superficial pyramidal neurons contributing to the generation of synchronous oscillationsin the visual cortex,” Science, vol. 274, no. 5284, pp. 109–113, 1996.
[4] J. R. Gibson, M. Belerlein, and B. W. Connors, “Two networks of electrically coupled inhibitory neurons in neocortex,” Nature, vol. 402, pp.
75–79, 1999.
[5] F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks. New York: Springer-Verlag, 1997.
[6], “Canonical neural models,” in Brain Theory and Neural Networks,
2nd ed, M. A. Arbib, Ed. Cambridge, MA: MIT, 2002.
[7] E. M. Izhikevich, N. S. Desai, E. C. Walcott, and F. C. Hoppensteadt,“Bursts as a unit of neural information: selective communication via
resonance,” Trends in Neurosci., vol. 26, pp. 161–167, 2003.
[8] E. M. Izhikevich, Dynamical Systems in Neuroscience: The Geometryof Excitability and Bursting, to be published.
[9], “Resonate-and-fire neurons,” Neural Networks, vol. 14, pp.
883–894, 2001.
[10], “Neural excitability, spiking, and bursting,” Int. J. Bifurc. Chaos,vol. 10, pp. 1171–1266, 2000.

Abstract—We analytically investigate the existence of global and partial
synchronism in neural networks where each node is represented by a phaseoscillator. Partial synchronism, which is important to pattern recognition,
can be caused by increasing the natural frequency of an oscillator and restricting the frequencies of others in certain ranges.
Index Terms—Phase-locked loop (PLL), phase oscillator, synchronism,
visual cortex.

I. INTRODUCTION
Hoppensteadt and Izhikevich [1] and Wang [2] have proposed network models where the neural activity is described by differential equations. Both architectures can be used for pattern recognition via associative memory, which occurs when a group of neurons fires synchronously. These models were inspired on findings and ideas about thefunctioning of the mammalian visual cortex.
It is known that different features of an object appearing in a visual scene, as its color, shape, velocity, and the direction of its motion, are processed in different cortical areas (e.g., [3]). Several authors have proposed that these features are linked through temporalcorrelations of neural activities (e.g., [4]). Thus, each feature is represented by a neural group oscillating in a synchronized way, and distinct
features of the same object are simultaneously represented by distinct
synchronous groups. Neural groups corresponding to different objects
must be desynchronized from each other, in order to avoid ambiguous
conjunctions. Experimental observations in mammals seem to support
this theory (e.g., [5]).
In the Wang model [2], the activity of a cortical column is describedby two nonlinear differential equations developed by Wilson and
Cowan [6] representing the interactions between two populations ofneurons that are distinguished by the fact that their synapses are either
excitatory or inhibitory.
In the Hoppensteadt-Izhikevich model [1], the activity of a corticalcolumn is represented by the nonlinear differential equation of a phaselocked loop (PLL) (e.g., [7], [8]) describing the temporal evolution of
the phase associated to the oscillation.
In both network models, the connectivity among nodes follows
simple laws. In this sense, these networks are classified as regular.
Despite this regularity, both networks are considered intractable
analytically, because it is a very hard task to find expressions relating
the parameter values of the oscillators in order to assure global orpartial synchronism.
We present analytical results about a simple version of the Hoppensteadt-Izhikevich network. We try to give some hints for answering theManuscript received January 31, 2003; revised April 4, 2003. The work of L.
H. A. Monteiro and J. R. C. Piqueira was supported by CNPq.
L. H. A. Monteiro and N.F.C. Canto are the with Pós-graduação em Engenharia Elétrica, Universidade Presbiteriana Mackenzie. They are also with the
Departamento de Engenharia de Telecomunicações e Controle, Escola Politécnica, Universidade de São Paulo, Brazil.
J. G. Chaui-Berlinck is with Departamento de Fisiologia, Instituto de Biocincias, Universidade de São Paulo, Brazil.
F. M. Orsatti and J. R. C. Piqueira are with Departamento de Engenharia
de Telecomunições e Controle, Escola Politécnica, Universidade de São Paulo,
Brazil (e-mail: luizm@mackenzie.com.br).
Digital Object Identifier 10.1109/TNN.2003.820441

1045-9227/03$17.00 © 2003 IEEE