## Abstract

I review a class of hybrid models of neurons that combine continuous spike-generation mechanisms and a discontinuous ‘after-spike’ reset of state variables. Unlike Hodgkin–Huxley-type conductance-based models, the hybrid spiking models have a few parameters derived from the bifurcation theory; instead of matching neuronal electrophysiology, they match neuronal dynamics. I present a method of after-spike resetting suitable for hardware implementation of such models, and a hybrid numerical method for simulations of large-scale biological spiking networks.

## 1. Introduction

Some of the most promising models suitable for simulation of spiking dynamics of neurons are hybrid in nature: they combine a smooth excitable behaviour leading to the generation of the action potential, often called a spike, and a subsequent discontinuous reset of state variables due to the spike.

The simplest and the best-known example of such a model is the *leaky integrate-and-fire neuron*
where *v*∈*R* is the membrane potential of the neuron, *C* is the membrane capacitance, *g*_{leak} is the leak ohmic conductance and *E*_{leak} is the leak reverse potential. The input current *I*(*t*) brings the membrane potential of this linear model to the threshold *v*_{threshold}. Once *v* crosses the threshold, a spike is said to be fired, and the membrane potential is reset to a new (subthreshold) value *c*. Though technically not a spiking model (it lacks an intrinsic spike-generation mechanism, and hence is just a ‘threshold’ model), the leaky integrate-and-fire neuron was a popular model to simulate spiking networks.

Ermentrout (1996) introduced the *quadratic integrate-and-fire neuron*, which we write here in the form
where *v*_{rest} and *v*_{thresh} are the resting and the instantaneous threshold potentials, respectively (when *I* = 0), and *k* > 0 is a parameter. Because the right-hand side of this differential equation behaves as *kv*^{2} for large *v*, the variable escapes (blows up) to infinity in a finite time. To avoid simulating ‘infinity’, one clips the voltage trajectory at some sufficiently large value, *v*_{peak}, and resets it to a new value *c*. In the original paper, Ermentrout (1996) avoided infinity by transforming the model into a trigonometric form with the state variable *θ* defined on a circle (Ermentrout & Kopell 1986), which led many people to call it *θ*-*neuron*; his formulation is equivalent to the one above when and *c* = − ∞.

It follows from the bifurcation theory that any dynamical system near a saddle-node bifurcation on an invariant circle can be converted to this model by an appropriate change of variables (Ermentrout 1996). Since many cortical pyramidal neurons exhibit such a bifurcation, the model became an instant success among modellers. Escaping to infinity corresponds to firing a spike. In this sense, the model has an intrinsic nonlinear ‘autocatalytic’ (positive feedback) process responsible for the generation of action potentials, and the parameter *v*_{peak} is not a threshold value as in the leaky integrate-and-fire neuron, but the peak of the spike. Indeed, the spike has already been initiated and *v* is already approaching infinity when *v* reaches *v*_{peak} in figure 1*b* (contrast it with the integrate-and-fire model in figure 1*a*).

This model jump-started a whole new approach towards modelling the spiking behaviour of neurons—hybrid spiking models—as it combines a smooth spike-generation mechanism with autocatalytic upstroke of the spike () and a hard ‘after-spike’ reset. This is in stark contrast to traditional conductance-based models (Skinner 2006) such as the Hodgkin–Huxley model, which have separate conductance variables that are responsible for upstroke and downstroke of action potentials. All the variables and parameters in conductance-based models have definite electrophysiological meaning, and they all could be measured experimentally, at least in principle. In practice, though, it is quite challenging to estimate the parameters, as any electrophysiologist would attest, so many arbitrary choices have to be made based on the assumption of ‘reasonability’.

In contrast, the quadratic integrate-and-fire model above and a simple spiking model (2.1)–(2.2) below have few parameters and they can be tuned to match quantitatively neuronal responses to various inputs, thereby leading to simulated behaviour that is indistinguishable from *in vitro* and *in vivo* neuronal responses (Izhikevich 2007). Changing the parameters, one can reproduce the firing patterns, spiking and bursting behaviour of many neuronal types (figure 2).

## 2. Derivation of simple model of spiking neurons

Many neuronal models, including conductance-based models, have phase portraits with nullclines as in figure 3*a*. The decision to fire or not to fire a spike is made in the shaded region, which is magnified in figure 3*b*. Notice that the fast *V* -nullcline looks like a square parabola and the slow *u*-nullcline can be approximated by a straight line in the shaded region. Moreover, there is a continuous change of variables that can transform the model in the region into the form with a quadratic fast *V* -nullcline and linear slow *u*-nullcline (Izhikevich 2003, 2007)
2.1
and
2.2
with the after-spike resetting
2.3
where *u* is a recovery variable, and *a*,*b*,*c* and *d* are independent parameters (the other parameters can be removed by an appropriate rescaling of variables and time). As in the case of the quadratic integrate-and-fire neuron, the voltage variable *v* escapes to infinity in a finite time; this corresponds to *V* in figure 3*a* leaving the shaded square and approaching the right-hand side branch of the N-shaped *V* -nullcline, i.e. generating the upstroke of an action potential. Instead of following the trajectory all the way until it returns to the shaded square, i.e. instead of simulating the peak and the downstroke of the action potential, the state variables of the simple spiking model are reset back to the shaded square instantaneously via equation (2.3); see figure 1. Such an after-spike reset simplifies the model significantly, yet retains the nonlinear spike-generation mechanism, excitability, spiking and bursting behaviour. Izhikevich (2007, ch. 8) provides examples of parameters of the model that simulate many cortical, thalamic and hippocampal neurons.

## 3. Other hybrid models

The simple spiking model can be transformed into the form
and
by an appropriate rescaling of variables. It belongs to a family of hybrid dynamical systems of the form
3.1
and
3.2
with the after-spike reset (2.3), though the parameters *c* and *d* have to be rescaled too (similarly, the parameters *a*,*b*,*c* and *d* in Izhikevich (2003) have the same meaning as those in Izhikevich (2007) but differ in values because of the rescaling).

Here, the function *f*(*v*) describes the current–voltage characteristic of the membrane potential near the threshold (Izhikevich 2007), and it typically looks like a parabola. Other choices besides *f*(*v*) = *v*^{2} are possible. For example, one can take *f*(*v*) = |*v*|^{3} or *f*(*v*) = 1/(1 − *v*)^{n} − *v* or *f*(*v*) = |*v*|^{n}_{+} − *v* with the linear rectifier function |*v*|_{+} = 0 when *v* ≤ 0 and |*v*|_{+} = *v* when *v* > 0. Brette & Gerstner (2005) suggested the *adaptive exponential integrate-and-fire model* with *f*(*v*) = *e*^{v} − *v*, whereas Touboul (2009) suggested the *quartic model* with *f*(*v*) = *v*^{4} + 2*αv*.

A whole new class of hybrid spiking models is given by the system of the form
3.3
and
3.4
where *u* plays the role of a conductance and *E* is its reverse potential, which could be assumed to take values ±1 or 0 after appropriate rescaling. Sometimes, it makes sense to use a nonlinear equation for (Izhikevich 2007). As long as *f*(*v*) scales as *v*^{1+ε} for some *ε* > 0 when , the voltage variable escapes to infinity in a finite time and all the hybrid models above can serve to simulate a spike-generation mechanism with the reset equation (2.3).

## 4. Hardware implementation: slow variable reset (parameter *d*)

It is relatively straightforward to implement equations (2.1)–(2.2) in a digital or analogue circuit. The after-spike reset of *u* by a constant *d* in equation (2.3) takes into account the action of high-threshold outward conductances that are activated during the spike and affect the timing of the following spikes. The reset, however, poses an engineering problem: implementing such a hard reset requires a disproportionately large number of transistors relative to the rest of the circuit, which is a limiting factor for large-scale neuromorphic hardware. Fortunately, the model allows a shortcut that removes the need for such a reset completely, thereby resulting in slimmer and faster hardware implementations.

Touboul (2009) pointed out that the behaviour of *u* is sensitive to the value of the parameter *v*_{peak}. Such a sensitivity arises owing to the fact that the slow variable *u* blows up to infinity too during the spike, i.e. while the fast variable *v* approaches infinity. This is not the case for the exponential integrate-and-fire model, where *u* stays bounded when . In fact (Touboul 2009), the slow variable is bounded for any other spiking model where the right-hand side of the voltage equation scales as *v*^{2+ε} for some *ε* > 0 when , i.e. it is faster than the square. The simple model (2.1)–(2.2) is thus unique.

Let us take advantage of this unique feature of the model: instead of stopping *v* at *v*_{peak} and then resetting *u* by a constant *d*, we can allow *v* to increase beyond *v*_{peak} until *u* is increased by *d* owing to the dynamics of equation (2.2). That is, we can move the parameter *v*_{peak} up so as to allow the slow variable to grow to the desired ‘after-spike’ value on its own.

The amount of adjustment of *v*_{peak} depends on the other parameters of the model. If one wants to remove the resetting of *u* by *d*, the new peak value is approximately *v*_{peak}*e*^{dk/(abC)} (derived elsewhere), but its precise value could be obtained numerically for any given parameter set. In fact, the simple spiking model (2.1)–(2.2) can be rewritten without the recovery variable reset entirely (i.e. *d* = 0), but treating *v*_{peak} as the fourth independent parameter (after *a*, *b* and *c*). Such a model would result in a slim hardware implementation.

## 5. Simulation methods

In this section, we present a few useful methods for effective simulation of the spiking model (2.1)–(2.2); the methods would also work for other types of hybrid spiking models. Typically, large-scale simulations are implemented using the explicit forward Euler method of the form
5.1
and
5.2
with the time step, *τ*, often taken to be *τ* = 1 ms. One can also use *v*(*t* + *τ*) instead of *v*(*t*) in the right-hand side of equation (5.2), as the new value of *v* is already known from the first equation.

There are two major issues with this method: detecting the crossing of *v*_{peak} corresponding to the peak of the spike and dealing with it, and numerical instability due to strong synaptic conductances.

### (a) Detecting *v*_{peak} crossing

As seen in §4, it is important to detect the crossing of the value *v*_{peak} so that the slow variable *u* is adequately modified and the timing of the spike is registered appropriately. A straightforward though highly inefficient way is to use a really small adaptive time step. Instead, we use a large time-step simulation, but do extra work when the neuron fires a spike.

Once the voltage value *v*(*t* + *τ*) is above *v*_{peak}, detected by a simple if–then statement, we use linear interpolation to determine the intermediate time (between *t* and *t* + *τ*) where the crossover actually occurred (figure 4)
(Nonlinear interpolation methods would also work, though they would require more computations per spike, and possibly extra memory requirement to store previous values of *v*.) The slow variable is then updated according to equation (5.2), but instead of *τ*, we use the smaller value *t*_{peak} − *t*. Thus, the slow variable *u* endures only partial update. One can use the value of *v* from the previous time step, *v*(*t*), the peak value *v*_{peak} or a combination of these values in equation (5.2).

### (b) Avoiding numerical instability

When simulating realistic neuronal networks, the term *I* in the simple model denotes the sum of all input currents. In particular, it includes the synaptic current
where *g*_{i}(*t*) is the time-varying conductance and *E*_{i} is the reverse potential for a particular synaptic current *i*, e.g. *i* = NMDA, AMPA, GABA_{A} and GABA_{B}. It is convenient to rewrite the synaptic term in the form
where is the total conductance and *E*(*t*) = ∑ (*g*_{i}(*t*)*E*_{i})/*g*(*t*) is the total reverse potential.

For the sake of clarity and to emphasize the synaptic term, we will rewrite the voltage equation in the model (2.1)–(2.2) in the more general form
5.3
To simulate this equation efficiently, one often employs the simplest numerical method—the explicit forward Euler method—with the time step *τ* = 1 ms:
5.4
This, however, may result in numerical instabilities often seen as zig-zag solutions illustrated in figure 5*a*. Indeed, positive values of *g*(*t*) push the membrane voltage towards the reverse potential *E*(*t*). However, when *g*(*t*) is large, the term *τg*(*t*)[*E*(*t*) − *v*(*t*)] in equation (5.4) becomes large, resulting in overshoot and divergence from *E*(*t*). Notice that, no matter how small the simulation step *τ* is, the conductance variable could always become so large as to create the instability. The standard way to avoid the numerical instability is to use either very small adaptive time step *τ* or more accurate numerical methods, such as implicit Euler methods or high-order Runge–Kutta methods. These, however, slow down the simulation considerably.

Since the numerical instability is caused by the linear term, one can use the hybrid numerical method
that combines the simplicity and efficiency of explicit methods and the numerical stability of implicit methods. Notice that the right-hand side contains the present (already known) value of the voltage variable, *v*(*t*), as well as its future (as yet unknown) value *v*(*t* + *τ*). However, the dependence on *v*(*t* + *τ*) is linear, and the equation can be solved for *v*(*t* + *τ*):
5.5
resulting in an efficient yet stable method illustrated in figure 5*b*.

## 6. Discussion

Hybrid models of spiking neurons turned out to be an invaluable tool for large-scale realistic simulations of brain models. For example, in 2005 Izhikevich simulated a thalamo-cortical system with six cortical layers, three thalamic nuclei and 22 different neuronal types, having the size of the human brain, i.e. 10^{11} neurons and almost one quadrillion synapses (http://www.izhikevich.org). A smaller network of only one million neurons with multi-compartment dendritic trees and human white-matter anatomy is described in Izhikevich & Edelman (2008). Each dendritic compartment was simulated using the same simple spiking model with an appropriate choice of parameters to reproduce slow dendritic excitability, forward- and back-propagating spikes. Such simulations would not be possible with conductance-based Hodgkin–Huxley models unless one has access to supercomputing resources.

In such large-scale brain models, the emphasis is on faithful reproduction of neuronal dynamics with the view to understand neural computations. Nevertheless, one can still use hybrid spiking models to understand the contribution of various ionic channels and conductances to neural computations and brain dynamics. For example, to model the effect of pharmacological manipulation of, say, delayed rectifier K channels, one needs to fit the parameters of the simple spiking model using recordings of control (normal) neurons and neurons manipulated by the pharmacological agent. These would result in two parameter sets, which then should be used in large-scale simulations to study the effect of pharmacological agents on brain dynamics.

While simulating large-scale brain models, there is a natural simulation step of 1 ms for all communications between neurons. This step is sufficiently small in comparison with the membrane time constant (a few milliseconds for most neurons) and with the axonal conduction delays (up to tens of milliseconds for cortico-thalamic and cortico-cortical (contralateral) connections), yet it is quite large from the numerical method point of view. Technically, the system (5.1)–(5.2) should not even be treated as being a discrete-time approximation of model (2.1)–(2.2), but rather as being an independent map-based neuronal spiking model on its own. No wonder such mappings are prone to become unstable when synaptic conductances are too strong. The hybrid simulation method proposed in §5*b* is probably the simplest and the most effective method resulting in stable dynamics. Combined with the procedure of dealing with crossing the peak of the spike (§5*a*), the simulation would not need to use small adaptive time steps or higher-order numerical methods.

Analogue hardware implementations, on the other hand, have their own constraints. An unexpectedly severe one is the requirement to reset the slow variable *u* by a constant *d* after each spike. The method described in §4 removes this requirement, thereby cutting down on the number of transistors needed to implement the simple spiking model.

## Footnotes

One contribution of 12 to a Theme Issue ‘Theory of hybrid dynamical systems and its applications to biological and medical systems’.

- © 2010 The Royal Society