Python implementation of Hodgkin and Huxley's squid giant axon model using scipy ode solver.

hh_squid.py

# hh_squid.py --- 
# Author: Subhasis  Ray
# Created: Fri Oct 26 10:21:21 2018 (-0400)
# Last-Updated: Mon Oct 29 12:32:49 2018 (-0400)
#           By: Subhasis  Ray
# Version: $Id$

# Code:
"""A python implementation of Hodgkin-Huxley`s Squid Giant Axon model."""


import numpy as np
from scipy import interpolate
from scipy import integrate
from matplotlib import pyplot as plt


veq = -65.0
v = np.arange(-100, 25, 0.1)

a_n = 0.01 * (-(v - veq) + 10)/ (np.exp((-(v - veq) + 10)/10.0) - 1)
b_n = 0.125 * np.exp(-(v - veq) / 80.0)

a_m = 0.1 * (-(v - veq) + 25.0) / (np.exp((-(v -veq) + 25) / 10.0) - 1)
b_m = 4 * np.exp(-(v - veq) / 18.0)
a_h = 0.07 * np.exp(-(v - veq) / 20.0)
b_h = 1.0 / (np.exp((-(v - veq) + 30.0) / 10.0) + 1)

n_inf = a_n / (a_n + b_n)
tau_n = 1 / (a_n + b_n)
m_inf = a_m / (a_m + b_m)
tau_m = 1 / (a_m + b_m)
h_inf = a_h / (a_h + b_h)
tau_h = 1 / (a_h + b_h)

######### Start plotting state variables 

grid = plt.GridSpec(1, 3)
grid.update(left=0.08, right=0.98, bottom=0.55, top=0.9, hspace=0.3, wspace=0.5)
fig = plt.figure()
ax_m = fig.add_subplot(grid[0, 0])
ax_m.plot(v, a_m, label=r'$\alpha_{m}$')
ax_m.plot(v, b_m, label=r'$\beta_{m}$')
ax_h = fig.add_subplot(grid[0, 1])
ax_h.plot(v, a_h, label=r'$\alpha_{h}$')
ax_h.plot(v, b_h, label=r'$\beta_{h}$')
ax_n = fig.add_subplot(grid[0, 2])
ax_n.plot(v, a_n, label=r'$\alpha_{n}$')
ax_n.plot(v, b_n, label=r'$\beta_{n}$')
for ax in fig.axes:
    ax.set_ylabel(r'rate (ms$^{-1}$)')
grid_m_tau = plt.GridSpec(1, 2)
grid_m_tau.update(left=0.08, right=0.98, bottom=0.1, top=0.45, hspace=0.1, wspace=0.2)
ax_inf = fig.add_subplot(grid_m_tau[0, 0])
ax_inf.plot(v, m_inf, label=r'$m_\infty$')
ax_inf.plot(v, h_inf, label=r'$h_\infty$')
ax_inf.plot(v, n_inf, label=r'$n_\infty$')
ax_tau = fig.add_subplot(grid_m_tau[0, 1])
ax_tau.plot(v, tau_m, label=r'$\tau_m$')
ax_tau.plot(v, tau_h, label=r'$\tau_h$')
ax_tau.plot(v, tau_n, label=r'$\tau_n$')
# fig.tight_layout()
for ax in fig.axes:
    ax.set_xlabel('Vm (mV)')
    ax.legend()
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    # ax.set_frame_on(False)
ax_tau.set_ylabel(r'$\tau$ (ms)')
fig.suptitle('State variables in Hodgkin-Huxley model')
fig.set_size_inches(6, 4)
# plt.show()

######### End plotting state variables 

## These are not used as we use the ode solver to compute m, h and n
## as functions of time and voltage
n_inf_v = interpolate.interp1d(v, n_inf, bounds_error=False, fill_value=(n_inf[0], n_inf[-1]))
# tau_n_v = interpolate.interp1d(v, tau_n, bounds_error=False, fill_value=(tau_n[0], tau_n[-1]))
m_inf_v = interpolate.interp1d(v, m_inf, bounds_error=False, fill_value=(m_inf[0], m_inf[-1]))
# tau_m_v = interpolate.interp1d(v, tau_m, bounds_error=False, fill_value=(tau_m[0], tau_m[-1]))
h_inf_v = interpolate.interp1d(v, h_inf, bounds_error=False, fill_value=(h_inf[0], h_inf[-1]))
# tau_h_v = interpolate.interp1d(v, tau_h, bounds_error=False, fill_value=(tau_h[0], tau_h[-1]))

## Simulation time and sampling time steps
tend = 20.0
dt = 0.025
trange = np.arange(0, tend, dt)

# Current injection
inject = 20.0   #  uA/cm2
delay = 5.0    # ms
duration = 2.0  # ms

irange = np.zeros_like(trange)
irange[(trange > delay) & (trange < (delay + duration))] = inject
irange_t = interpolate.interp1d(trange, irange, bounds_error=False, fill_value=(0, inject))

## Channel conductances
g_na = 120    # mS/cm2
g_k = 36     # mS/cm2
gleak = 0.3  # mS/cm2
cm = 1.0        # 1uF/cm2
## Reversal potentials
eleak = -54.3  # mV
ena = 50.0     # mV
ek = -77.5     # mV

## Initial conditions
vm0 = veq
m0 = m_inf_v(vm0)
h0 = h_inf_v(vm0)
n0 = n_inf_v(vm0)

## Look-up tables for the transition rates
a_n_v = interpolate.interp1d(v, a_n, kind='cubic', bounds_error=False, fill_value=(a_n[0], a_n[-1]))
b_n_v = interpolate.interp1d(v, b_n, kind='cubic', bounds_error=False, fill_value=(b_n[0], b_n[-1]))
a_h_v = interpolate.interp1d(v, a_h, kind='cubic', bounds_error=False, fill_value=(a_h[0], a_h[-1]))
b_h_v = interpolate.interp1d(v, b_h, kind='cubic', bounds_error=False, fill_value=(b_h[0], b_h[-1]))
a_m_v = interpolate.interp1d(v, a_m, kind='cubic', bounds_error=False, fill_value=(a_m[0], a_m[-1]))
b_m_v = interpolate.interp1d(v, b_m, kind='cubic', bounds_error=False, fill_value=(b_m[0], b_m[-1]))

def dv_dt(t, y):
    """The function for computing the derivatives. This is passed to the
    ODE integrator for initial value problem (ivp). Instead of
    computing an exponential at each step for m, h and n, we pass the
    computation to the ode integrator as well.

    t: float - time

    y : array-like of length 4 containing [0]m, [1]h, [2]n, [3]vm

    NOTE: m, h, n and Vm are the only variables that require solving
    with ODE integrator. To plot the other variables, gna,
    gk, ina, ik, store them in global variables.

    """
    m, h, n, vm = y
    dm_dt = a_m_v(vm) * (1 - m) - b_m_v(vm) * m
    dh_dt = a_h_v(vm) * (1 - h) - b_h_v(vm) * h
    dn_dt = a_n_v(vm) * (1 - n) - b_n_v(vm) * n
    gna = g_na * m**3 * h
    ina = gna * (vm - ena)
    gk = g_k * n**4
    ik = gk * (vm - ek)
    ileak = (vm- eleak) * gleak
    i = irange_t(t)
    dv_dt = (i - (ina + ik + ileak)) / cm
    # print(f't:{t} m:{m}, h:{h}, n:{n}, Vm:{vm}, i:{i}, gk:{gk}, gna:{gna}, dv_dt:{dv_dt}')
    return np.array([dm_dt, dh_dt, dn_dt, dv_dt])


ret = integrate.solve_ivp(dv_dt, (trange[0], trange[-1]), [m0, h0, n0, vm0], t_eval=trange, method='RK45', atol=1e-9, rtol=1e-6)

gna_arr = g_na * ret.y[0, :]**3 * ret.y[1, :]
ina_arr = gna_arr * (ret.y[3, :] - ena)
gk_arr = g_k * ret.y[2, :]**4
ik_arr = gk_arr * (ret.y[3, :] - ek)

fig, axes = plt.subplots(nrows=5, sharex='all')
axes[0].plot(ret.t, ret.y[0, :], label='m')
axes[0].plot(ret.t, ret.y[1, :], label='h')
axes[0].plot(ret.t, ret.y[2, :], label='n')
axes[1].plot(ret.t, gna_arr, label=r'$g_{Na}$')
axes[1].plot(ret.t, gk_arr, label=r'$g_{K}$')
axes[2].plot(ret.t, ina_arr, label=r'$i_{Na}$')
axes[2].plot(ret.t, ik_arr, label=r'$i_{K}$')
axes[3].plot(ret.t, ret.y[3, :], label='Vm (mV)')
axes[4].plot(trange, irange, label=r'injected current ($\mu$A/cm$^{2}$)')
axes[4].set_xlabel('Time (ms)')

for ax in axes.flat:
    ax.legend()
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
fig.suptitle('Time course of membrane voltage\n and state variables in Hodgkin-Huxley model')
fig.set_size_inches(6, 6)
plt.show()


    
# 
# hh_squid.py ends here