This is the second blog post in a series - you can find the previous blog post here

Integrate and Fire Models

Throughout this blog post we will focus on integrate and fire models. This class of model has been around for a long time, in fact longer than the Hodgkin-Huxley model. The first model was presented by Lapicque in 1907. Since then many alternative formulations have been presented. We can express the models in the form: $$ \tau \frac{dV}{dt} = -(V - E) + \psi(V) + R_m I(t) $$

Where:
$\tau$ represents membrane time constant ($C_m / g_L$ in notation used previously)
$E$ represents the rest potential
$\psi(V)$ is the spike generating current term
$R_m$ represents membrane resistance ($1/g_L$)
$I(t)$ is a function representing an applied current (in the Hodgkin-Huxley example we took a Gaussian white noise)

With integrate and fire models we have the issue that (typically) the action potential will shoot off to infinity. In order to stop this we implement a threshold ($V_{Th}$), when the process reaches this value it is reset (to $V_{Re}$) and the dynamics start again. This is not a big issue because for our purposes we are noly interested in the dynamics of an onset of an action potential, the mechanism of returning to normal levels is not of (much) interest. When modelling we may wish to hold the voltage at $V_{Re}$ following a spike for a short time ($\Delta_{T_{Rf}}$) to reflect the refractory period of a neuron. The refractory period can be modelled stochastically but usually a static value is succficient.

Integrate and fire models are thus defined by the form of the function $\psi(V)$ - this function may be linear or non-linear. Examples include the leaky-integrate and fire model (linear), Fitzhugh-Nagumo (polynomial) and the exponential integrate and fire (non-linear).

From Hodgkin-Huxley to Integrate and Fire

We now take a quick de-tour to justify the use of the integrate and fire model as an approximation to the Hodgkin-Huxley dynamics.

First we notice that gate m operates on a much faster time scale than gates n or h (and similarly much faster than the leak channel which controls the potential dynamics with all gates closed.) Given it is so much faster we can apply an instantaneous approximation, namely: $m(t) = \hat{m}(V_m(t))$ that is: the dynamics are defined by the membrane voltage. From plotting gate dynamics we can also observe that gates n and h are approximately translated reflections of each other. As an approximation we can create an adaptation variable $w$ with $n = aw$ and $h = b - w$ for constants $a, b$. We can then write down the equation: $$ C_m \frac{dV_m}{dt} = I_p - \overline{g_K}(aw)^4(V_m-V_K) - \overline{g_{Na}}(\hat{m}(V_m))^3(b - w)(V_m-V_{Na}) - \overline{g_L}(V_m-V_L) $$

There is a corresponding equation for the adaptation variable $w$ which we shall not concern ourselves with. We are only interested in the onset of spiking not the refractory period dynamics so we will take $w = w_{rest}$ to be a fixed value. We can therefore express the voltage dynamics as:
$$ C_m \frac{dV_m}{dt} = I_p - \overline{g_{eff}}(V_m-V_{eff}) - \lambda (\hat{m}(V_m))^3(V_m-V_{Na}) $$

By collecting terms and some re-arrangement. $g_{eff}, V_{eff}$ and $\lambda$ are all constant values. Therefore via the approximations outlined above we are left with an equation of the form:
$$ \tau \frac{dV}{dt} = -(V - E) + \psi(V) + R_m I(t) $$

Namely an integrate and fire model.

Exponential Integrate and Fire

Continuing with the line of reasoning above we shall consider the function: $\hat{m}(V_m)$. We assume the dynamics are so rapid that they are essentially in equilibrium:
$$ \frac{dm}{dt} = \alpha_m(V_m)(1-m) - \beta_m(V_m)m = 0 $$
So:
$$ \hat{m}(V_m) = \frac{\alpha_m(V_m)}{\alpha_m(V_m) + \beta_m(V_m)}$$

Since Hodgkin-Huxley suggested the following forms of these equations:
$\alpha_m(V_m) = \frac{0.1(25-V_m)}{e^{(2.5-0.1V_m)}-1} $
$\beta_m(V_m) = 4 e^{-V_m / 18} $

We can approximate $\hat{m}(V_m)$ with a logistic function:
$$ \hat{m}(V_m) \approx (1 + e^{-\beta(V_m - \theta)})^{-1}$$
We can then express the Taylor expansion of this as:
$$ \hat{m}(V_m) = \sum (-1)^k e^{\beta(V_m - \theta)(1+k)}$$
Which we can see from the expansion of $1/(1+y)$ with $y = e^{-\beta(V_m - \theta)}$. So a first order approximation is:
$$ \hat{m}(V_m) \approx e^{\beta(V_m - \theta)}$$

Then the current corresponding to the sodium channel can be expressed approximately:
$$I_{Na} = g_{Na}(b - w_{rest})(V_m - V_{Na}) e^{3\beta(V_m - \theta)} $$

If we take $(V_m - V_{Na}) \approx (V_{rest} - V_{Na}) < 0$ as an approximation we get approximate voltage dynamics as:
$$ C_m \frac{dV_m}{dt} = I_p - \overline{g_{eff}}(V_m-V_{eff}) + \hat{\lambda} e^{\beta (V_m - \theta)} $$

This is known as the exponential integrate and fire (EIF) model.This model has been shown to fit experimental data (and the Hodgkin-Huxley model) very well in practice. Typically we use a parameterization of the spiking term:
$$ \psi(V) = \Delta_T e^{\left( \frac{V - V_T}{\Delta_T} \right)} $$ We can see this model is highly non-linear and without applying the threshold mechanics the membrane potential would shoot to infinity. The 2 new parameters in this model are:
$V_T$ which represents the voltage scale at which the exponential term becomes significant in the dynamics
$\Delta_T$ representing the sharpness of the spike

With a Gaussian white noise term $\xi_t$ we can fully specify the dynamics using: $$\tau \frac{dV_m}{dt} = (V_L - V_m) + \Delta_T e^{\left( \frac{V_m - V_T}{\Delta_T} \right)} + \sigma \sqrt{2 \tau}\xi_t $$

As before we can solve this using a foward Euler scheme (although the forcing term is non-linear this still provides a reasonable solution for small enough time steps.) An implementation of this can be seen below:

# Implementation of a noisy EIF neuron using a forward Euler scheme

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# Set seed for repeatability
np.random.seed(123)

# Set time step dt (ms) and number of steps N
dt = 0.001
N = 50000

# Set model parameters
# Membrane time constant tau (ms) and leak reversal potential VL (mV)
tau = 30
VL = -70

# Spike sharpness DelT (mV) and exponential potential threshold VT (mV)
DelT = 3
VT = -60

# Variation in gaussian noise sig
sig = 25

# Set voltage spike threshold Vth (mV), reset voltage Vr (mV) and refractory period Tref (ms)
Vth = 30
Vr = -70
Tref = 5

# Set up arrays for time T (ms) and voltage V (mV)
T = np.arange(N) * dt
V = np.zeros(N)
V[0] = Vr

# Set up refractory period counter Tc (ms)
Tc = 0

for i in range(1, N):
    if Tc > 0:
        V[i] = Vr
        Tc -= 1
    else:    
        Vtemp = V[i-1] + dt/tau*(VL - V[i-1]) + DelT*dt/tau*np.exp((V[i-1] - VT)/DelT) + sig*np.sqrt(2*dt/tau)*np.random.normal(0,1,1)
        if Vtemp > Vth:
            V[i] = Vr
            Tc = np.ceil(Tref/dt)
        else:
            V[i] = Vtemp

# Plot voltage trajectory
plt.plot(T, V)
plt.xlabel("Time (ms)")
plt.ylabel("Membrane Potential (mV)")
plt.title("Membrane Voltage Trajectory")
plt.show()

Conclusion

The voltage trajectory displays realistic neuron dynamics for the onset of spiking. However as expected through the use of the refractory period implementation the depolarizing phase is not captured well. This is ok since we consider the information to be carried by the spike itself not the behaviour shortly afterwards. We can also see that this implementation is considerably simpler than that of Hodgkin-Huxely since there is only one ODE. This is of benefit when modelling large networks of neurons where time/computational constraints become a consideration.

It has been shown that the EIF model can predict spiking behaviour very well in practice, despite the somewhat cavalier assumptions made during its derivation from the Hodgkin-Huxley.

In a future blog post we will make use of the mathematical tractability of the EIF model to analyse neurons in more depth.

References

This blog post is the second part of a series - you can find the next blog post here