Integrate and fire models with XPP

Integrate-and-fire models are straightforward to simulate using XPP (or MATLAB, for that matter.) The difficulty with such models is that they have a discontinuous reset. Here is a typical integrate and fire model:

tau dV/dt = Vrest - V + Rm Ie(t)

along with the reset condition:

V(t+)=Vreset when V(t-)=Vthresh

Here is the commented XPP file for this model with a general applied current which is a mixture of sinusoids and a constant bias current iaf.ode 

# integrate and fire
# here is the ODE
v'=(-v+vrest + rm*ie(t))/tau
init v=-65
#
# here is the reset condition - when v-vthreh=0 in the increasing 
# direction, reset the voltage
global 1 v-vthresh {v=vreset}
# here is the stimulus
ie(t)=a0+a1*sin(w1*t)+a2*sin(w2*t)
# here are parameters
# voltage - mV  time -  msec current - nanoamp R - MOhm 
par vrest=-65,vthresh=-50,vreset=-70,tau=15
par rm=10 
# stimulus parameters
par a0=1.5,a1=.75,w1=.05,a2=.75,w2=.12345
# plot the stimulus
aux stim=rm*ie(t)-vrest
# set up simulation parameters
@ total=500,dt=.1
# set up plot parameters
@ nplot=2,xlo=0,xhi=500,ylo=-75,yhi=-35, yp2=stim
# extend the maximal storage - XPP stores only 5000 pts by default
@ maxstor=20000
done

Run this file and look at the output. You will see that the I&F model tends to fire on the upstroke of the stimulus rather than the peak. Try a pure sine stimulus: change a1=1.5,a2=0 and vary the frequency w1 Count the spikes in a 1000 msec stimulation for a variety of frequencies, say 2 > w1 > .01

The following is a similar file but this time the stimulus is low-pass filtered white noise iafnoise.ode 

 
# integrate and fire with noise
# here is the ODE
v'=(-v+vrest + rm*ie)/tau
init v=-65
#
# here is the reset condition - when v-vthreh=0 in the increasing 
# direction, reset the voltage
global 1 v-vthresh {v=vreset}
# here is the stimulus
ie'=(-ie+a1*normal(0,1))*w1
# here are parameters
# voltage - mV  time -  msec current - nanoamp R - MOhm 
par vrest=-65,vthresh=-50,vreset=-70,tau=15
par rm=10 
# stimulus parameters
par a1=30,w1=.05
# plot the stimulus
aux stim=rm*ie+vrest
# set up simulation parameters
@ total=500,dt=.1,meth=euler,bound=1000
@ nplot=2,xlo=0,xhi=500,ylo=-75,yhi=-35
@ yp2=stim,maxstor=20000
done

The stimulus satisfies:

dIe/dt = (N(t) - Ie)w1

where w1 is the filter frequency and N(t) is delta-correlated noise with amplitude a1 . Try different values of w1. Whenever you have noise, you should use Euler's method. Here are some notes numerical integration of differential equations.

The integrate and fire with adaptation

Cortical neurons have spike frequency adaptation (SFA). This is an outward (hyperpolarizing) current which makes the neuron less likely to fire once it has previously fired. It takes some time to build up. If a current step is applied, then the first couple spikes occur rapidly followed by spikes with larger ISIs. We can modify the I&F model to incorporate SFA. Here is a general I&F model with SFA:

cm dV/dt = -gL (V-Vrest) -g(V-EK) + I(t)/A

dg/dt = (-g + a delta(t))/taua

If we multiply this by the total membrane resistance and rescale g, we get:

tau dV/dt = Vrest - V + Rm I(t) -g(V-EK)

dg/dt = (-g + a/gL delta(t))/taua

Each time the cell fires, g is incremented by dg=a/gL and decays with a time constant taua . Here is the XPP file iafadap.ode : 

# iafadap.ode
v'=(vrest-v + rm*ie(t)-g*(v-ek))/tau
g'=-g/taua
init v=-65,g=0
# I increment adaptation whenever the spike fires -- g=g+dg  
global 1 v-vthresh {v=vreset;g=g+dg}
# heav(t) is the step function
ie(t)=amp*heav(t-ton)*heav(toff-t)
par vrest=-65,vthresh=-50,vreset=-65,tau=15,rm=10
par ek=-85,dg=.1,taua=100
par amp=4,ton=50,toff=200
aux stim=rm*ie(t)+vrest
@ total=500,dt=.1,maxstor=10000
@ xhi=500,yhi=-45,ylo=-85,nplot=2,yp2=stim
done

Try this. Then set the adaptation to zero dg=0 and see the difference. Try different parameters such as the magnitude of the stimulus, amp and the time constant of adaptation taua of the degree of adaptation dg .

Math geek stuff

We could actually solve this analytically to get the steady-state firing rate as a function of the applied current. However, there will be a number of integrals for which there is no closed form. So instead, we can find the steady-state firing rate numerically. One way is to keep the current on and look at the ISI after a long periodi of time. But there is a much cooler way to do this. Suppose that we are at steady-state and spiking with unknown period P. Then at t=0 , V(0)=Vreset and g(0)=g(P)+dg since right after a spike, g is incremented. At t=P , V(P)=Vthresh . So, we have a boundary value problem with two ODEs and three boundary conditions. However, we also have the free parameter, P so we rescale time to the interval [0,1] and must solve:

tau dV/dt = (Vrest - V + Rm I -g(V-EK))P

dg/dt = -Pg/taua

V(0)=Vreset

V(1)=Vthresh

g(0)=g(1)+dg

Here is the XPP file iafrbvp.ode if you want to try it. Use the Bndryval Show option to see it solve the BVP. Change parameters a litte at a time and resolve it if you want.

Other simple models

The integrate-and-fire model is the simplest possible model for a spiking neuron. However, the FI curve does not match inhibitory neurons very well and there are other examples of cortical neurons whose first few spikes are poorly matched. The I&F model produces an almost linear FI curve while many cortical neurons have a more nonlinear FI curve. Of course, we can get these kinds of FI curves from biophysical models. But, there is a simple model which arises near certain critical parameter values (called bifurcation points) which is called the "theta-model" or the quadratic integrate-and-fire model. The equations for this model are:

tau dV/dt = a (V-Vrest)(V-Vthr) + R I(t)

Note that a is a positive parameter with dimensions of 1/volt that governs the "excitability" of the model. If V exceeds the threshold, Vthr, then the solutions to this equation blow up in a finite amount of time. This is the spike. The voltage is then reset to -infinity. If I(t) is constant and large enough, the neuron will blow up periodically. I will leav it as an exercise for you to find the period:

P = tau int[ dV/{a (V-Vrest)(V-Vthr) + R I},V=-infinity,V=+infinity]

Since this is really not possible to handle numerically (since infinity is not a computer-friendly concept), we can make a change of variables to obtain the so-called "theta" model. Let

V(t) = (Vthr+Vrest)/2+btan(x(t)/2)

where

b = (Vthr-Vrest)/2

Then x(t) satisfies:

tau b dx/dt = ab2[1-cos(x)] + [1+cos(x)][RI -a b2]

Whenever x(t) crosses an odd multiple of pi, a spike is emitted so it is like a less tractable version of the integrate-and-fire model. However, with a constant current, it has a very nice square-root FI curve which matches many cortical neurons much better than the I&F. The frequency or firing rate is :

F = Sqrt[ a R I - (ab)2]/(Pi tau)

Here is an example of the first ISI of a cortical pyramidal cell with a square-root fit: F=70 Sqrt[I-.82].

Find parameters in the theta model which match this assuming that R=50 MOhm and tau=30 msec. The only free parameters are a,b. Note the firing rate in the figure is in Hz but the FI curve given in the formula above is in spikes per millisecond!