神经元模型 Hodgkin Huxley model

最近在接触脉冲神经网络的基础知识，于是准备先了解几个常见的神经元。首先便是最经典的Hodgkin Huxley模型了。

一、介绍

二、电路图

三、离子电流表征

四、代码实现

一、介绍

Hodgkin Huxley模型是Alan Hodgkin和Andrew Huxley在1952年描述的。两人通过测量鱿鱼巨型轴突中动作电位的对不同电流的反应，在没有计算机的情况下拟合出了一个模型。这个模型由一组四个常微分方程组成，虽然描述较为准确，但是计算复杂，不适合大范围仿真。

二、电路图

Hodgkin Huxley模型将细胞建立为如上图所示电路。其中脂质双层被表示为电容 $C_{m}$ ；钠、钾离子通道分别用一个电源和电导来表示（图中统一建立为电源 $E_{n}$ 和电导 $g_{n}$ ）；泄露通道也用一个电源 $g_{L}$ 和电导 $E_{n}$ 来表示；最后产生的离子泵由电流 $I_{p}$ 表示；膜电位由电压 $V_{m}$ 表示。值得一提的是，其中钠、钾离子通道的电导与时间、电压相关，而泄漏通道的电导与他们不同，是恒定值。

由电路公式我们可以简单得知，流过脂质双层的电流大小为：

$I_{c} = C_{m}\frac{dV_{m}}{dt}$

又由欧姆定律，对于钠、钾离子通道，其电流都可以表达为：

$I_{Na} = g_{Na} \cdot (V_{m}-V_{Na})$

$I_{K} = g_{K} \cdot (V_{m}-V_{K})$

当然，因为膜内外离子分布不同，细胞外钠离子多，细胞内钾离子多，所以这里的两种电流的方向一定是相反的，二者的数值在符号上一定相反。泄漏电流的公式也与两者一样，所以不再表出。

由此，可以得到流过脂质双层的电流大小表达为：

$I = C_{m}\frac{dV_{m}}{dt}+ g_{Na}\cdot (V_{m} - V_{Na}) +g_{K}\cdot (V_{m} - V_{K})+g_{L}\cdot (V_{m} - V_{L})$

三、离子电流表征

通过上一小节，我们可以了解到了流过脂质双层的电流大小的表达式。而Hodgkin Huxley通过一系列的实验，开发了一个模型，其中可兴奋细胞的特性由一组四个常微分方程表达。

$I = C_{m}\frac{dV_{m}}{dt} + \overline{g_{K}} \cdot n^{4}(V_{m}-V_{K})+\overline{g_{Na}} \cdot m^{3}\cdot h(V_{m}-V_{Na})+\overline{g_{L} }\cdot (V_{m}-V_{L})$

$\frac{dn}{dt} = \alpha_{n} (V_{m}) (1-n) -\beta _{n}(V_{m})n$

$\frac{dm}{dt} = \alpha_{m} (V_{m}) (1-m) -\beta _{m}(V_{m})m$

$\frac{dh}{dt} = \alpha_{h} (V_{m}) (1-h) -\beta _{h}(V_{m})h$

其中 $I$ 为流经单位面积的电流。在该神经元模型中，每一个离子通道都由多个亚基控制，公式中的n、m、h分别对应钾离子通道亚基激活概率、钠离子通道亚基激活概率、钠离子通道亚基失活概率，因为在鱿鱼轴突中，例如钾离子通道由4个亚基控制，所以公式中要乘上3次方，其他几个乘方也是同理。而 $\overline{g_{i}}$ 代表对应离子通道的最大电导率，乘上对应的概率之后，就得到了对应该时刻的电导率。

在下三个公式中， $\alpha _{i}$ 代表亚基从失活到激活的速率， $\beta _{i}$ 代表亚基从激活到失活的速率，它们与这一时刻的膜电位相关。乘上对应激活亚基和失活亚基的比例，就可以得到这一时刻的变化量。

通过这四个公式，我们就可以较准确地表达出膜电位的生物特性。而在作者原论文中，给出了 $\alpha$ 和 $\beta$ 的表达式。

$\alpha_{n} = \frac{0.01(10-V)}{e^{\frac{10-V}{10}}-1}$

$\alpha_{m} = \frac{0.1(25-V)}{e^{\frac{25-V}{10}}-1}$

$\alpha _{h} = 0.07e^{-\frac{V}{20}}$

$\beta _{n} = 0.125e^{-\frac{V}{80}}$

$\beta _{m} = 4e^{-\frac{V}{18}}$

$\beta _{h} = \frac{1}{e^{\frac{30-V}{10}}+1}$

其中， $V = V_{rest}-V_{m}$

而关于形如 $\frac{dx}{dt} = \alpha_{x} (1-x) -\beta _{x}x$ 的方程，我们可以得到其解为：

$x(t) = x_{0}-(x_{0}-x{\infty })(1-e^{\frac{t}{\tau _{x}}})$

其中 $x_{\infty } = \frac{\alpha_{x}}{\alpha_{x} + \beta_{x}}$ ， $\tau _{x} = \frac{1}{\alpha _{x} + \beta _{x}}$

易知，t=0时， $x_{0 } = \frac{\alpha_{x}}{\alpha_{x} + \beta_{x}}$

当然，为了能让这个模型顺利运行起来，我们还需要其他参数。原论文中也给出了作者对这些参数的测量值，如下表所示。分别对应脂质双层的电容（ $\mu F/cm^{2}$ ），三种离子通道的静息电位( $mV$ )，三种离子通道的最大电导率( $mS/cm^{2}$ )，以及静态膜电位（ $mV$ ）。

$C_{m}$	$E_{Na}$	$E_{K}$	$E_{L}$	$\overline{g_{Na}}$	$\overline{g_{K}}$	$\overline{g_{L}}$	$V_{m}$
1	50	-77	-54.4	120	36	0.3	-65

接着我们就可以用代码来仿真这个神经元了。

四、代码实现

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline#参数定义
#平衡电位 单位mV
E_Na = 50.0
E_K = -77.0
E_L = -54.4
V_m = -65.0#最大电导 单位 mS/cm^2
g_Na = 120.0
g_K = 36.0
g_L = 0.3#脂质双层电容
C_m = 1
#单位时间
dt = 0.01#总时间
T = 40#每一时刻
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 = 0.0V[0] = V_m#计算n,m,h初始值
V_tmp = V[0]-V_m
alpha_n = (0.1 - 0.01 * V_tmp) / (np.exp(1 - 0.1 * V_tmp) - 1)
alpha_m = (2.5 - 0.1 * V_tmp) / (np.exp(2.5 - 0.1 * V_tmp) - 1)
alpha_h = 0.07 * np.exp(-V_tmp/20.0)beta_n = 0.125 * np.exp(-V_tmp / 80.0)
beta_m = 4 * np.exp(-V_tmp / 18.0)
beta_h = 1/(np.exp(3 - 0.1 * V_tmp)+1)n[0] = alpha_n / (alpha_n + beta_n)
m[0] = alpha_m / (alpha_m + beta_m)
h[0] = alpha_h / (alpha_h + beta_h)for i in range(1, len(t)):#计算alpha和betaV_tmp = V[i-1] - V_malpha_n = (0.1 - 0.01 * V_tmp) / (np.exp(1 - 0.1 * V_tmp) - 1)alpha_m = (2.5 - 0.1 * V_tmp) / (np.exp(2.5 - 0.1 * V_tmp) - 1)alpha_h = 0.07 * np.exp(-V_tmp/20.0)beta_n = 0.125 * np.exp(-V_tmp / 80.0)beta_m = 4 * np.exp(-V_tmp / 18.0)beta_h = 1/(np.exp(3 - 0.1 * V_tmp)+1)#计算delta_n, delta_m, delta_htau_n = 1.0/(alpha_n + beta_n)tau_m = 1.0/(alpha_m + beta_m)tau_h = 1.0/(alpha_h + beta_h)inf_n = alpha_n *tau_ninf_m = alpha_m *tau_minf_h = alpha_h *tau_h#更新 n,m,h,使用欧拉一阶近似法n[i] = n[i-1] + dt * (alpha_n * (1-n[i-1]) - beta_n * n[i-1])m[i] = m[i-1] + dt * (alpha_m * (1-m[i-1]) - beta_m * m[i-1])h[i] = h[i-1] + dt * (alpha_h * (1-h[i-1]) - beta_h * h[i-1])#更新电流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 - (I_Na + I_K + I_L))*dt/C_mV[i] = V[i-1] + dvplt.clf()
plt.subplot(1,1,1)
plt.plot(t,V)
# 膜电位变化
plt.title('Membrane potential')
plt.legend(['V'])
plt.ylim([-100, 100])

接下来我们就可以设置不同大小、时间的电流，来给神经元刺激，观察神经元的膜电位变化情况。

将刺激电流设置为0，得到的图像如下：

可见，膜电位保持在-65mV静息电位。

当给予一个短而大的刺激电流时，膜电位会产生一个很明显的突变，图像如下：

如果持续一段时间地给予较大的电流刺激，膜电位就会在这段时间里持续地产生波形，图像如下：

但当持续一段时间地给予过大的电流刺激时，膜电位在产生一个突变之后便会保持一个比静息电位高的数值不变，直到刺激结束恢复静息电位，图像如下：

当给予神经元一次短而小的刺激时，并不能让膜电位突变到正值；但若短时间间隔内再次用小电流刺激后，膜电位反而会产生一次突变，图像如下：

通过不断改变电流的大小和时间，便可以探究神经元的不同特性了。