Neural Modeling with Python (Part 3)

By Byron Galbraith | February 2, 2011

So far we've looked at how to simulate a simple LIF model neuron and a complex Hodgkin-Huxley model neuron. The LIF neuron is computationally simple but physiologically implausible, while Hodgkin-Huxley gives us a very good representation of actual neural dynamics but is parameter-heavy and computationally expensive. An intriguing compromise between the two exists -- one that can generate a wide variety of observed neural spiking behavior while doing so with limited computational demand. It is called the quadratic integrate-and-fire model neuron, or simply Izhikevich neuron.

Izhikevich Model Neuron

In 2003, Eugene Izhikevich published a paper entitled "Simple Model of Spiking Neurons" (Izhikevich, 2003) in which he describes how the quadratic integrate-and-fire neuron can be used to efficiently replicate the observed spiking dynamics of several different classes of cortical neurons. Surprisingly, this can be accomplished using just two differential equations and four parameters, a slight increase from the LIF model but a drastic reduction from Hodgkin-Huxley. The equations are
with a reset similar to the LIF neuron defined by
The variables and parameters in the model are all dimensionless values that can be thought of as representing various states or properties of the neuron. In these equations, v is membrane potential, u is membrane recovery, a is the time scale of u, b relates subthreshold sensitivity of u to v, c is the reset value of v after spiking, and d is the effect spiking has on u at reset.

This may all seem rather abstract to you and I'd agree. My main criticism of the Izhikevich model is that it seems to be plagued with so-called magic numbers -- constants and parameter values that appear without justification or intuitive basis. It is not immediately obvious how to choose the values of a, b, c, and d to achieve the desired behavior, nor is why the multipliers in the equation for dv/dt are what they are. However, the point of the Izhikevich neuron is not to gain an understanding of how a neuron works, but to be able to replicate a wide range of experimentally observed neural spiking behavior cheaply for use in large-scale simulations.

Object-Oriented Python

The last two examples were straightforward procedural scripts simulating a single set of parameters, so the MATLAB equivalents wouldn't be that much different. For replicating Izhikevich's 20 different neural spiking behaviors as demonstrated in "Which Model to Use for Cortical Spiking Neurons?" (Izhikevich, 2004) I decided to go with a simple object-oriented approach. I'm not going to present the entire script on this page as in the previous two posts for brevity's sake, though as before the full script will be available for download.

The IzhNeuron class

The first class I define is IzhNeuron, a representation of an Izhikevich model neuron. This class will be used to initialize and keep track of the internal state of a neuron during simulation.

class IzhNeuron:
  def __init__(self, label, a, b, c, d, v0, u0=None):
    self.label = label

    self.a = a
    self.b = b
    self.c = c
    self.d = d

    self.v = v0
    self.u = u0 if u0 is not None else b*v0

Since this is the first class definition I've presented, I'll go over the basics in detail. First, a class is defined by using the class keyword followed by the name of the class and a colon to indicate the start of the class definition code block. This class defines the single special method __init__ which is called by the Python interpreter whenever a new instance of the class is created. This is equivalent to constructors in Java, for instance. The __init__ method can take any number of arguments, though, as is true for all class methods, the first argument passed in will be a reference to the object itself. Another point to mention here is that we can assign default values to arguments, as seen with u0=None at the end of the argument list. This means that if the u0 argument is omitted from the initialization call it will be assigned the value of None (Python's version of null or undefined).

The __init__ method's purpose here is to initialize the neuron's state variables and parameters, as well as assign a label to the neuron used to identify it when plotting the results. We also need to make sure we are assigning the values to the particular instance of the class and not to the general program itself, thus the need for the self reference. Finally, the u0 parameter is being assigned using Python's ternary operator. This is convenient shorthand for an if-else statement, and can be read in this example as "use the value of u0 if u0 is not None else return the value of b*v0."

The IzhSim class

In addition to having 20 differently defined neurons, Izhikevich also specifies different simulation setups for each, which I represent with the IzhSim class. This class keeps track of the neuron it's simulating, the time and time step used for integration, the input stimulus, and a class method for performing the integration.

class IzhSim:
  def __init__(self, n, T, dt=0.25):
    self.neuron = n
    self.dt     = dt
    self.t      = t = arange(0, T+dt, dt)
    self.stim   = zeros(len(t))
    self.x      = 5
    self.y      = 140
    self.du     = lambda a, b, v, u: a*(b*v - u)

  def integrate(self, n=None):
    if n is None: n = self.neuron
    trace = zeros((2,len(self.t)))
    for i, j in enumerate(self.stim):
      n.v += self.dt * (0.04*n.v**2 + self.x*n.v + self.y - n.u + self.stim[i])
      n.u += self.dt * self.du(n.a, n.b, n.v, n.u)
      if n.v > 30:
        trace[0,i] = 30
        n.v        = n.c
        n.u       += n.d
        trace[0,i] = n.v
        trace[1,i] = n.u
    return trace

This is similar to the IzhNeuron class in that it has an __init__ method for initializing the various instance properties needed for the simulation. The x and y properties are there because it turns out that in a few of Izhikevich's examples he changes up the values of some of constants in the dv/dt equation. There is also a particular example where he uses a different equation for du/dt, thus the du lambda expression. In those cases where these parameters need to be modified, I can override them after creating the IzhSim instance.

Setting up neurons

Now that we have these two classes defined, we can easily setup the 20 example neurons before simulating and plotting them. Following are just a few to show how it can be done.

sims = []

## (A) tonic spiking
n = IzhNeuron("(A) tonic spiking", a=0.02, b=0.2, c=-65, d=6, v0=-70)
s = IzhSim(n, T=100)
for i, t in enumerate(s.t):
  s.stim[i] = 14 if t > 10 else 0

## (G) Class 1 exc.
n = IzhNeuron("(G) Class 1 exc.", a=0.02, b=-0.1, c=-55, d=6, v0=-60)
s = IzhSim(n, T=300)
s.x = 4.1
s.y = 108
for i, t in enumerate(s.t):
  s.stim[i] = 0.075*(t-30) if t > 30 else 0

## (R) accomodation
n = IzhNeuron("(R) accomodation", a=0.02, b=1, c=-55, d=4, v0=-65, u0=-16)
s = IzhSim(n, T=400, dt=0.5)
s.du = lambda a, b, v, u: a*(b*(v + 65))
for i, t in enumerate(s.t):
  if t < 200:     s.stim[i] = t/25
  elif t < 300:   s.stim[i] = 0
  elif t < 312.5: s.stim[i] = (t-300)/12.5*4
  else:           s.stim[i] = 0

The three chosen models demonstrate a typical setup, one that requires different constants in the dv equation, and one that requires a different integration time step and du equation, respectively. I also use explicitly named arguments when initializing the IzhNeuron objects, though this is for clarity only as I could just as effectively used only the values as long as they were in the correct order.

Simulation and plotting

Now that we've initialized all 20 neural simulations, we can use a simple loop to evaluate and plot the results.

for i,s in enumerate(sims):
  res = s.integrate()
  ax  = subplot(5,4,i+1)

  ax.plot(s.t, res[0], s.t, -95 + ((s.stim - min(s.stim))/(max(s.stim) - min(s.stim)))*10)

  ax.set_ylim([-100, 35])
  ax.set_title(s.neuron.label, size="small")

Here we iterate over each simulation in the sims list, first calling its integrate method, which returns the traces of the v and u variables. The rest is plot formatting to mimic the figure presented by Izhikevich, though it's worth noting a couple points about this. First, while Python starts array indexing at 0, subplot indexing starts at 1 like MATLAB, thus the i+1. Second, we remove the x- and y-axis labels by setting them to empty lists. That's it!


I demonstrated how to use simple object-oriented programming to greatly simplify the simulation and rendering of 20 distinct Izhikevich neurons. Not only does this setup now allow for additional neurons and simulation setups to be added or removed with ease, it does so in one third the lines of code compared to the example MATLAB code. The Izhikevich neuron model itself is, while perhaps on the abstract side, a powerful tool for replicating many complex spiking patterns. Furthermore it can do this with the most complicated term being a quadratic, which is a significant advantage when it comes to computational efficiency.

One Response to Neural Modeling with Python (Part 3)

  1. Richard Tomsett says:

    Enjoying these python articles, thanks! Just a pedantic point – the Quadratic IAF model is similar to Izhikevich’s neuron, but doesn’t have the slow recovery variable term (u), so cannot display as wide a range of behaviour (it is closely related to the theta neuron: ).

Leave a Reply

Your email is never published nor shared. Required fields are marked *


You may use these HTML tags and attributes:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>