Representing the function and behavior of neurons in software is one of the core activities of computational neuroscience. As neurons communicate via electrochemical currents, this is typically accomplished through modeling the dynamical nature of the neuron's electrical properties. Several models treat the neuron as an equivalent electrical circuit, with its membrane potential described by one or more differential equations. In order to simulate the response of the neuron to various stimuli, these equations are numerically solved over some time interval for a given pattern of input current. There are several methods for numerically solving differential equations, and for the purpose of this series, I'm going to use the forward Euler method because it's one of the easiest to implement and understand while providing sufficient stability.

So, without further ado, let's start simulating some neurons in Python!

**Leaky Integrate-and-Fire**

The leaky integrate-and-fire (LIF) neuron model is one of the most basic formalisms of the spiking behavior of neurons. It is effectively a simple RC-circuit model (leaky integrator) combined with a delta-function (fire) tacked on when the membrane potential reaches a given threshold followed by a reset and optional refractory period. While the LIF model has limited basis physiologically, its strength lies in its computationally simple ability to generate spike events with an input-driven frequency. The LIF model is expressed as the single differential equation

,

where V is the membrane potential, tau_m is the membrane time constant equal to R_m*C_m, the equivalent resistance and capacitance of the neural membrane, respectively, and I is the input current. Adding a refractory period simply turns this into the piecewise equation

,

where t_ref is the time at which the neuron's refractory period ends. An example trace of an LIF neuron's membrane potential over time for a constant input it shown below.

And now what you've been waiting for: the Python code to generate this figure! As Python is white space sensitive, I highly recommend downloading the code as copy-pasting from this page may introduce errors.

Scripts: lif.py

`from numpy import *`

from pylab import *

## setup parameters and state variables

T = 50 # total time to simulate (msec)

dt = 0.125 # simulation time step (msec)

time = arange(0, T+dt, dt) # time array

t_rest = 0 # initial refractory time

## LIF properties

Vm = zeros(len(time)) # potential (V) trace over time

Rm = 1 # resistance (kOhm)

Cm = 10 # capacitance (uF)

tau_m = Rm*Cm # time constant (msec)

tau_ref = 4 # refractory period (msec)

Vth = 1 # spike threshold (V)

V_spike = 0.5 # spike delta (V)

## Stimulus

I = 1.5 # input current (A)

## iterate over each time step

for i, t in enumerate(time):

if t > t_rest:

Vm[i] = Vm[i-1] + (-Vm[i-1] + I*Rm) / tau_m * dt

if t >= Vth:

Vm[i] += V_spike

t_rest = t + tau_ref

## plot membrane potential trace

plot(time, Vm)

title('Leaky Integrate-and-Fire Example')

ylabel('Membrane Potential (V)')

xlabel('Time (msec)')

ylim([0,2])

show()

Hopefully most of this example is self-explanatory. However, there are a few key points I'd like to go over in more depth.

**NumPy functions: arange, zeros**

The arange(start, stop, step) function on line 7 generates an ordered array of numbers starting with start and incrementing by step until stop is reached. It is important to note that arange does not typically include stop in the final array. In this example, we want to generate a series of time points for the simulation, starting at 0 and going until 50 msec with a timestep of 0.125 msec. Since we want to include 50 in the final array, we simply make our stop value the desired time plus the timestep. *MATLAB equivalent: 0:dt:T;*

The zeros(shape) function on line 11 generates an n-dimensional array of zeros. There are two main differences between this and its MATLAB counterpart. First, if shape is a single integer, the array is 1-by-shape instead of a square matrix. Second, if a higher-order array is desired, zeros requires that shape be specified as a sequence e.g. zeros([2,3]). Those extra braces can be easy to forget when coming from a MATLAB background! In this example, we want to initialize an array for keeping track of our neuron's membrane potential at each time step. We use the len function to get the number of elements in the time array so the two will line up during the simulation phase. *MATLAB equivalent: zeros(size(time));*

**Simulation method**

The most interesting part of this example is the actual method of simulation so I'm going to look at it closely.

`for i, t in enumerate(time):`

This is the beginning of the simulation loop. In English, this is telling the Python interpreter "for every element in the array time, run through the proceeding code block setting i to the current index and t to the current value." Python has good built-in support for iterators so you can avoid dealing with array indices if their only purpose is to immediately perform a lookup for the desired element. For instance, for c in "characters" would set c to each letter in the string "characters" in turn. Sometimes, however, you do have a use for the index, in which case you can use the enumerate function depicted here.

I'm going to take time now to remind any reader who needs a refresher that a quirk of Python is that instead of using keywords e.g. end or curly braces to indicate code blocks, it uses a colon and white space. The colon at the end of the for line signifies the beginning of the loop code block, while every line within the block must be indented the same number of spaces. The Python interpreter knows the block ends when the next line of code starts with the indentation as the line with the initiating colon. Nested blocks use further levels of indentation.

` if t > t_rest:`

A simple conditional checking to ensure the neuron currently isn't in a refactory period. If it is -- the current simulation time is less than the end of the most recent refractory period -- then we save some time and skip the numerical integration step, otherwise we continue.

` Vm[i] = Vm[i-1] + (-Vm[i-1] + I*Rm) / tau_m * dt`

This is the forward Euler method for numerical integration applied to the LIF differential equation.

` if Vm[i] >= Vth:`

Vm[i] += V_spike

t_rest = t + tau_ref

Finally, if the current value of the membrane potential exceeds the threshold value of 1V, we add on the spike delta (0.5V) and set the refractory period to end 4 msec from now.

**Plotting**

The plotting functions provided by Matplotlib generally follow those of MATLAB. The key thing to note here is the last line of code. While MATLAB will display a figure by default the moment it is created, Matplotlib draws plots in the background and needs to be explicitly instructed to make the figure visible on screen. This is accomplished by the show function.

**Conclusion**

This was the first example of using Python for neural modeling. The LIF neuron model was showcased for its simplicity and ubiquity, though there are more interesting and useful ones, such as the Hodgkin-Huxley model which I will demonstrate next. Hopefully you are intrigued by what Python has to offer in this area, as there is plenty more to come!

What are the major advantages/reasons for working with Python, instead of (say) Matlab or Octave?

I address that question in my previous post (http://www.neurdon.com/2010/12/20/python-in-computational-neuroscience/). My main claims are that Python is free, does not suffer from license restrictions (a concern when thinking about scalability), and is a general programming language. While the first two points apply to Octave as well, the latter is the biggest standout. In Python, I get many of the benefits of MATLAB combined with the power of C++, Java, or Perl in one environment.

Excellent tutorial!!!!

Here some comments:

I personally like using the namespaces of the python modules before using the methods of the module. For example:

>>> import numpy as np

>>> V = np.zeros(len(time))

You avoid loading the whole module in your current namespace.

To add more physiological realism, you could start the simulation with the Vm vector in the resting membrane potential. This can make the trick:

Vrest = -70 # in mV

Vm = ones(len(time))*Vrest

Now you could set the action potential threshold -50 mV

Vtrh = =50 # in mV

I like expressing the membrane input resistance (Rm) in terms of ‘leak’ conductance (1/gl), to let the current term (Iapp) alone in the equation. In this way, we can set Iapp easily as an argument to the function. Moreover, expressing the voltage in terms of dirving force( V-El), where El is the reversal potential of the resistor (-65 mV), would allow you to see the effect of the voltage to current injections in two directions (positive, depolarizing, and negative, hyperpolarizing).

The equation results:

Cm*dV/dt = -gl*(V-El) + Iapp

I implemented a object-oriented version of the integrate-and-fire model. You can see it here:

“””

The integrate-and-fire neuron model is a simple model of spiking

behaviour that sacrificies biophysical realism for mathematical

simplicity. Bellow threshold, the membrane potential V obeys the

differential equation:

Cm*dV/dt = Iapp – gl*(V-El)

If V reaches a threshold (Vthr), then neuron is said to spike, and

V is instantaneously reset to a value (usually to the resting

membrane potential). The term Iapp represent the current (eg. current

injection through an electrode, a synapse, etc…). The equation is

described in Eq 5.7 in the Dayan and Abbot book (page 163)

“””

import numpy as np

from matplotlib.pyplot import figure, show

class SimpleIF(object):

“”” Simple Integrated-and-fire model of the form:

Cm*dV/dt + gl*(V-El) = Iapp

Cm is the membrane capacitante

gl is the leak conductance

El is the resting membrane potential

Ipp is the applied current

“””

def __init__(self, Cm, gl, El, Vinit=-70.):

self._Cm = float(Cm) # in ms

self._gl = float(gl) # in uS

self._El = float(El) # in mV

self._Vinit = float(Vinit) # in mV

# fixed threshold

self._Vthr = -50 # in mV

def dVdt(self, Vmb, Iapp):

“”” returns the instantaneous voltage change

as a function of Cm, gl, El and Iapp

“””

Cm = self._Cm

gl = self._gl

El = self._El

return (Iapp – gl*(Vmb – El) )/Cm

def timecourse(self, current):

“”” returns the time course of the voltage as

a function of the current injected when solved

by the implicit Euler method:

f(x) = f(x-1) + dt*f'(x)

we assume dt is the sampling interval of the current

“””

voltage = np.empty(len(current))

# initial condition

voltage[0] = self._Vinit

dt = 0.02 # step size integration

# Eulers solver

for i in range(1, len(current)):

dVdt = self.dVdt(Vmb = voltage[i-1], Iapp = current[i-1])

voltage[i] = voltage[i-1] + dt*dVdt

# action potential threshold

if voltage[i] > self._Vthr:

voltage[i-1] = 0.0 # overshooting

voltage[i] = self._El # reset voltage

return voltage

if __name__ == ‘__main__’:

# current injection

dt = 0.02 # in ms

current = int(250/dt)*[0] + int(1000/dt)*[2.5] + int(250/dt)*[0]

mycell = SimpleIF(Cm =4.9, gl = .16, El=-65., Vinit = -75) # spike threshold = -50 mV!!!

voltage = mycell.timecourse(current)

# figure

fig1 = figure()

ax1 = fig1.add_subplot(111)

time = np.linspace(0, len(current)*(dt/1000), len(current)) # transform in sec

ax1.plot(time, voltage)

ax1.set_ylabel(‘Voltage (mV)’)

ax1.set_xlabel(‘Time (ms)’)

show()

Thanks Nin! This is a much more Pythonic take on this example.

Thanks Byron for your excellent tutorials, i’m learning Theoretical/Computational Neuroscience and your tutorials are helping me to understand and write neuron models in python.

Cheers!!

I’m sorry if there is a reason behind this, or it has been mentioned before and not corrected on purpose, but there is a difference between the code you have posted in the tutorial and the code in the lif.py file:

In the tutorial, the condition for the second if statement states

if t >= Vth, whilst in the file it statesif Vm[i] >= Vth, which is the correct version.TUTORIAL CODE:

for i, t in enumerate(time):

if t > t_rest:

Vm[i] = Vm[i-1] + (-Vm[i-1] + I*Rm) / tau_m * dt

if t >= Vth:Vm[i] += V_spike

t_rest = t + tau_ref

FILE CODE:

for i, t in enumerate(time):

if t > t_rest:

Vm[i] = Vm[i-1] + (-Vm[i-1] + I*Rm) / tau_m * dt

if Vm[i] >= Vth:Vm[i] += V_spike

t_rest = t + tau_ref

Otherwise, thank you for the great tutorial, even someone like me that knows nothing about programming and the mathematical foundation of this model found it easy to understand and very accessible!