Hodgkin Huxley 神经元模型及实现

最近在接触脉冲神经网络相关知识，如有错误，请多多指教！

一 原理

1. 电路图

Hodgkin-Huxley(HH)模型可以看作是rc电路模型的扩展，该模型除了泄漏通道外还包含钠(Na)和钾(K)通道。离子通道被建模为与电池串联的电阻器。电池表示特定离子的平衡电势，电阻器反映通道对特定离子的渗透性。与泄漏通道的固定电阻相比，Na和K通道的电阻有自己的动态特性，这取决于跨膜的电压。
更精确地说，对于给定的膜电位，每个通道都有相应的开启和关闭速率。这些通道的打开还是关闭是由霍奇金根据经验建立的函数和函数来描述的。这些和函数可以用在微分方程中来描述开放通道的比例如何随时间变化：
$$\frac{dn}{dt}={\alpha }_n(V)(1-n)-{\beta }_n(V)n（公式一）$$
alpha_n和 beta_n描述打开通道的比例如何随时间变化
$${\tau }_n=\frac{1}{{\alpha }_n+{\beta }_n}(公式二)$$
表示多久通道才会打开
$$n_{\infty }=\frac{{\alpha }_n}{{\alpha }_n+{\beta }_n}(公式三)$$
表示通道打开的概率

2. 常微分方程

由于HH模型需要使用4个常微分方程进行表示模型，故计算较为复杂，很难进行大规模仿真。但其精确地描绘出膜电压的生物特性，能够很好地与生物神经元的电生理实验结果相吻合

二 代码实现-python版本

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

1. 模拟离子通道

# 离子通道被建模为与电池串联的电阻器。电池表示特定离子的平衡电势，电阻器反映通道对特定离子的渗透性
V = np.arange(-80.0,80,0.01) # [mV]#n用来限制k离子的通透性，m和h用来限制na离子，因为na离子含有失活门，即使通道本身是打开的，如果失活门阻塞通道，离子电流也不能流动。
alpha_n = 0.01 *(10-V)/(np.exp((10-V)/10)-1)
alpha_m = 0.1*(25-V)/((np.exp((25-V)/10))-1)
alpha_h = 0.07 * np.exp(-(V/20))beta_m = 4* (np.exp(-V/18))
beta_h = 1/(np.exp((30-V)/10)+1)
beta_n = 0.125 * np.exp(-(V/80))#公式二、公式三实现
tau_n = 1/(alpha_n + beta_n)
inf_n = alpha_n*tau_ntau_m = 1/(alpha_m + beta_m)
inf_m = alpha_m*tau_m tau_h = 1/(alpha_h + beta_h)
inf_h = alpha_h*tau_h#绘图
plt.clf()
plt.subplot(1,2,1)
plt.plot(V,inf_n)
plt.plot(V,inf_m)
plt.plot(V,inf_h)
plt.title('Steady state values') 
plt.xlabel('Voltage (mV)')
plt.subplot(1,2,2)
plt.plot(V,tau_n)
plt.plot(V,tau_m)
plt.plot(V,tau_h)
plt.title('Time constants') 
plt.xlabel('Voltage (mV)')
plt.legend(['n','m','h'])

运行结果如下图：由图可知，K离子选通变量$n$表明在平衡电位，大多数K趋于关闭。随着膜电位的增加，通道开始打开，但有一定的延迟。对于钠通道，当膜电位增加时，m值也增加，时间常数表明变化几乎是瞬时的（值很小）。然而，$h$趋于下降，这说明，在相对高的电压下，许多钠通道会打开，但大多数通道是失活的，因此几乎不会有离子电流。

2. 模拟膜电位

$$-C_m\frac{dV}{dt}=g_L(V-E_L)+g_K{n}^4(V-E_K)+g_{Na}{m}^{3}h(V-E_{N}a)-I_{extern}（公式四）$$

# 参数定义 可参考
# 平衡电位
E_Na = 115.0   # [mV]
E_K  = -12.0   # [mV]
E_L  = 10.6    # [mV]# 最大电导
g_Na = 120.0   # [mS]
g_K  = 36.0    # [mS]
g_L  = 0.3     # [mS]dt = 0.01      # [ms]
T = 40         # [ms]
t = np.arange(0,T,dt)V = np.zeros(len(t))
n = np.zeros(len(t))
m = np.zeros(len(t))
h = np.zeros(len(t))I_E = 0.0
V[0] = 0.0
h[0] = 0.59
m[0] = 0.05
n[0] = 0.31# 在10ms是注入10mA的电流，然后15ms再次置电流为0
for i in range(1,len(t)):if i == 1000:I_E = 10.0if i == 1500:I_E = 0.0# Calculate the alpha and beta functions （代码同上）alpha_n = (0.1 - 0.01*V[i-1]) / ( np.exp(1   - 0.1*V[i-1]) - 1)alpha_m = (2.5 - 0.1 *V[i-1]) / ( np.exp(2.5 - 0.1*V[i-1]) - 1)alpha_h = 0.07*np.exp(-V[i-1]/20.0)beta_n = 0.125*np.exp(-V[i-1]/80.0)beta_m = 4.0 * np.exp(-V[i-1]/18.0)beta_h = 1 / ( np.exp(3 - 0.1*V[i-1]) + 1)# Calculate the time constants and steady state valuestau_n = 1.0/(alpha_n + beta_n)inf_n = alpha_n*tau_ntau_m = 1.0/(alpha_m + beta_m)inf_m = alpha_m*tau_mtau_h = 1.0/(alpha_h + beta_h)inf_h = alpha_h*tau_h# 更新n,m,h  n[i]  = (1-dt/tau_n)*n[i-1] + (dt/tau_n)*inf_nm[i]  = (1-dt/tau_m)*m[i-1] + (dt/tau_m)*inf_mh[i]  = (1-dt/tau_h)*h[i-1] + (dt/tau_h)*inf_h# 公式四# 更新膜电位方程I_Na = g_Na*(m[i]**3)*h[i]  * (V[i-1]-E_Na)I_K  = g_K *(n[i]**4)       * (V[i-1]-E_K)I_L  = g_L                  * (V[i-1]-E_L)#dv = -(I_Na + I_K + I_L - I_E)dv = I_E - (I_Na + I_K + I_L)V[i]  = V[i-1] + dv*dtplt.clf()
plt.subplot(1,3,1)
plt.plot(t,V)
# 膜电位变化
plt.title('Membrane potential')
plt.subplot(1,3,2)
plt.plot(t,n)
plt.plot(t,m)
plt.plot(t,h)
plt.title('Gating variables')
plt.legend(['n','m','h'])
plt.subplot(1,3,3)
plt.plot(t,g_Na*h * m**3)
plt.plot(t,g_K*n**4)
plt.title('Ionic currents')
plt.legend(['Na','K'])

运行结果如下图：

三 代码实现-matlab版本

%% Basic Hodgkin-Huxley Implementation
% Equations from Hodgkin, Huxley, J. Physiol (1952) 117, 500-544clear;
figure(1);clf;dt = 0.001;
t = 0:dt:40;
% Nernst Potentials
Ena = 115; Ek = -12; El = 10.6;% Maximum Conductances
gna = 120; gk = 36; gl = 0.3;% Membrane Capacitance
C = 1;% 离子通道打开
an = @(u) (0.1-0.01*u)/(exp(1-0.1*u)-1);
am = @(u) (2.5-0.1*u)/(exp(2.5-0.1*u)-1);
ah = @(u) 0.07*exp(-u/20);% 离子通道关闭
bn = @(u) 0.125*exp(-u/80);
bm = @(u) 4*exp(-u/18);
bh = @(u) 1/(exp(3-0.1*u)+1);% DEFINE FORMULAE FOR STEADY STATE GATE ACTIVATIONSm_inf = @(u) am(u) / ( am(u) + bm(u) );
n_inf = @(u) an(u) / ( an(u) + bn(u) );
h_inf = @(u) ah(u) / ( ah(u) + bh(u) );% INITIALIZE STATE VARIABLESm = 0*t + m_inf(0);
n = 0*t + n_inf(0);
h = 0*t + h_inf(0);
v = 0*t;% DEFINE STIMULUS STRENGTH, DURATION, & DELAYtSTIM_START = 10;
tSTIM_DUR = 1;
STIM_STRENGTH = 10;% DEFINE MISCinRange = @(x,a,b) (x>=a) & (x<b);% MAIN LOOPfor i = 1:length(t)-1% extract membrane voltageu = v(i);% solve for membrane currentsIk  = gk  * n(i)^4      * (u - Ek);Ina = gna * m(i)^3*h(i) * (u - Ena);Il  = gl  *               (u - El);   Imem = Ik + Ina + Il;% determine stimulus current, if anyIstim = 0;if inRange( t(i) , tSTIM_START , tSTIM_START+tSTIM_DUR)Istim = STIM_STRENGTH;end% define the state variable derivativesdv = (Istim - Imem)/C;dm = am(u) * (1-m(i)) - bm(u) * m(i);dh = ah(u) * (1-h(i)) - bh(u) * h(i);dn = an(u) * (1-n(i)) - bn(u) * n(i);% use forward euler to increment the state variablesv(i+1) = v(i) + dv*dt;m(i+1) = m(i) + dm*dt;h(i+1) = h(i) + dh*dt;n(i+1) = n(i) + dn*dt;end% plot the results
plot(t,v);
axis([t(1) t(end) -20 120]);
xlabel('time (ms)');
ylabel('membrane voltage (mV)');

运行结果

四 参考

1. Hodgkin Huxley
2. matlab参考代码