Hogkin-Huxley minimal sim

hh.py

# Adapted from https://hodgkin-huxley-tutorial.readthedocs.io/en/latest/_static/Hodgkin%20Huxley.html

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

C_m  =   1.0
g_Na = 120.0
g_K  =  36.0
g_L  =   0.3
E_Na =  50.0
E_K  = -77.0
E_L  = -54.387

def alpha_m(V):       return 0.1*(V+40.0)/(1.0 - np.exp(-(V+40.0) / 10.0))
def beta_m(V):        return 4.0*np.exp(-(V+65.0) / 18.0)
def alpha_h(V):       return 0.07*np.exp(-(V+65.0) / 20.0)
def beta_h(V):        return 1.0/(1.0 + np.exp(-(V+35.0) / 10.0))
def alpha_n(V):       return 0.01*(V+55.0)/(1.0 - np.exp(-(V+55.0) / 10.0))
def beta_n(V):        return 0.125*np.exp(-(V+65) / 80.0)

def I_Na(V, m, h):    return g_Na * m**3 * h * (V - E_Na)
def I_K(V, n):        return g_K  * n**4 * (V - E_K)
def I_L(V):           return g_L * (V - E_L)

def dydt(t, y, iapp):
    V, m, h, n = y
    dVdt = (iapp - I_Na(V, m, h) - I_K(V, n) - I_L(V)) / C_m
    dmdt = alpha_m(V)*(1.0-m) - beta_m(V)*m
    dhdt = alpha_h(V)*(1.0-h) - beta_h(V)*h
    dndt = alpha_n(V)*(1.0-n) - beta_n(V)*n
    return dVdt, dmdt, dhdt, dndt

iapp = 0.5
res = solve_ivp(dydt, (0, 450), [-65, 0.05, 0.6, 0.32], args=(iapp,), vectorized=True)
t, (V, m, h, n) = res.t, res.y

plt.plot(t, V)