Code
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
class HHM:
def __init__(self, Iapp):
self.C = 1.
self.Vn, self.Vk, self.Vl = 115., -12., 10.613
self.gn, self.gk, self.gl = 120., 36., .3
self.am = lambda V: .1 * (V - 25) / (1 - np.exp(-(V - 25) / 10))
self.bm = lambda V: 4. * np.exp(-V / 18)
self.ah = lambda V: .07 * np.exp(-V / 20)
self.bh = lambda V: 1. / (1 + np.exp(-(V - 30) / 10))
self.an = lambda V: .01 * (V - 10) / (1 - np.exp(-(V - 10) / 10))
self.bn = lambda V: .125 * np.exp(-V / 80)
self.Iapp = Iapp
def f(self, y, t):
V, n, m, h = y
return np.array([
-(self.gk * (n ** 4) * (V - self.Vk)
+ self.gn * (m ** 3) * h * (V - self.Vn)
+ self.gl * (V - self.Vl) + self.Iapp(t)) / self.C,
self.an(V) * (1 - n) - self.bn(V) * n,
self.am(V) * (1 - m) - self.bm(V) * m,
self.ah(V) * (1 - h) - self.bh(V) * h])
def visualize(T, Iapp = lambda t: -12 if t % 100 < 2.5 else 0, g = False):
m = HHM(Iapp); I = [m.Iapp(t) for t in T]
V, N, M, H = odeint(m.f, np.array([-0., .32, .05, .6]), T).T
plt.plot(T, V, label = 'V'); plt.plot(T, I, label = 'I'); plt.legend(); plt.show()
if g:
plt.plot(T, N, label = 'n'); plt.plot(T, M, label = 'm'); plt.plot(T, H, label = 'h'); plt.legend(); plt.show()
Result
visualize(np.linspace(0, 500, 1001), g = True)
visualize(np.linspace(0, 5000, 10001))
visualize(np.linspace(0, 500, 1001), Iapp = lambda t: -12 if t % 10 < 2.5 else 0)