HH神经元模型基础知识与实现

HH神经元模型的建立

1963年，Hodgkin和Huxley通过他们对神经元动作电位的解释与Eccles分享了诺贝尔奖，这类神经元模型也被称为HH模型。尽管HH模型年代久远，而且仍然有许多不足，但是这一模型在现在仍然在使用，并且有望在计算神经科学领域发挥作用

Hodgkin以及Huxley的工作

1.他们利用电压钳位放大器来研究电流与电压的关系。通过改变细胞外离子 的浓度（主要是钠离子和钾离子），推断离子所携带的电流。
2. 将实验结果与数学模型相拟合
3. 求解数学模型。通过比较模型的动作电位与实验记录的动作电位验证模型的正确性

基础知识

在开始介绍具体内容前，我们先简单复习一下神经元的基础知识。神经元产生动作电位主要有五个阶段，分别是极化，去极化，反极化，复极化，超极化。
1.极化：神经纤维膜外正电荷聚集，膜内负电荷聚集（外正内负）。即神经元的静息状态。
2.去极化： 神经受到刺激时，门控Na⁺通道打开，膜外的Na⁺短期内大量涌入膜内，造成了内膜电位升高。
3.反极化：门控Na⁺通道关闭，门控K⁺通道打卡，K⁺大量外流，膜电位恢复静息电位。
4.超极化：钾离子通道完全打开后，这时钾离子的大量外流使得质膜再度极化，以至于超过原来的静息电位，称为超极化；超极化后膜电位会恢复至静息电位。

模型构建

在HH模型中主要研究三种离子电流，分别是上文提到的钠离子电流 I_Na，钾离子电流 I_K，以及他们称之为leak current 的泄露电流 I_L。
神经元的总体电流方程如下：
$$I=I_C+I_i=C_m\frac{dV}{dt}+I_i.\text{\hspace{1em}(1.1)}$$
其中：
$$I_i=I_{N}a+I_K+I_L\text{\hspace{1em}(1.2)}$$
且C_m为常数。
每种粒子的电流大小都是由对应的电压乘以该离子的膜电导所计算出来:
$$\begin{gathered} I_{N}a=g_{N}a(V-E_{N}a) \\ {I}_K=g_K(V-E_k) \\ {I}_L={\bar{g}}_L(V-E_L)\text{\hspace{1em}(1.3)} \end{gathered}$$
由于漏电流中的电导为常数，因此我们主要介绍如何计算钾离子以及钠离子电流。

钾离子电流

Hodgkin和Huxley通过实验发现，对细胞膜上的钾离子通道进行了假设。他们假设细胞膜上含有多个离子通道，且每个离子通道都有多个独立的门控粒子控制，假设n是单个门控粒子打开状态的概率。在他们的实验中发现，每个钾离子通道都由四个独立的门控离子所控制，因此钾离子电导g_k和门控粒子打开状态概率n的关系为：
$$g_k={\bar{g}}_k{n}^4\text{\hspace{1em}(1.4)}$$
其中g_K bar 是一个常数，表示的是钾离子的最大电导值。
另外，如果处于开放状态的粒子占比为n，那么关闭状态的粒子就是1-n，并且他们的状态可以相互转化。我们假设从关闭状态变为开放状态粒子的概率为α_n， 从开放状态变为关闭状态的粒子的概率为β_n, 那么我们可以得到n随着时间变化的公式：
$$\frac{dn}{dt}={\alpha }_n(1-n)-{\beta }_{n}n\text{\hspace{1em}(1.5)}$$
通过实验数据可以得到α_n以及β_n的公式：
$$\begin{gathered} {\alpha }_n=0.01\frac{V+55}{1-exp(-(V+55)/10)} \\ {\beta }_n=0.125exp(-(V+65/80))\text{\hspace{1em}(1.6)} \end{gathered}$$
另外，我们设时间为无穷大时的n为：
$$n_{\infty }=\frac{{\alpha }_n}{{\alpha }_n+{\beta }_n}\text{\hspace{1em}(1.7)}$$
时间常数为：
$${\tau }_n=\frac{1}{{\alpha }_n+{\beta }_n}\text{\hspace{1em}(1.8)}$$

钠离子电流

我们可以使用推到钾离子电流的类似方法推到钠离子电流，但是钠离子电流与钾离子电流不同的是，在反极化的过程中钠离子通道会关闭，也就是失活，而钾离子通道会保持一直开放。因此Hodgking和Huxley认为不能想钾离子一样只包含一个门控变量n，为此他们引入了另一个门控变量h来表示失活水平，类似的，h随着时间变化的公式为：
$$\frac{dh}{dt}={\alpha }_h(1-h)-{\beta }_{h}h\text{\hspace{1em}(1.9)}$$
并且钠离子通道中的门控粒子打开的概率设为m，m与n和h一样，与时间变化的公式为：
$$\frac{dm}{dt}={\alpha }_m(1-m)-{\beta }_{m}m\text{\hspace{1em}(1.10)}$$
最后，通过实验验证我们可以得到钠离子电导的公式：
$$gNa={\bar{g}}_{N}a{m}^{3}h\text{\hspace{1em}(1.11)}$$
m和h在时间无穷大的计算方式与n相同，只需将式子中的α和β换成对应下标的参数即可，在这里就不一一赘述。
α_m和β_m , α_h和β_h的计算公式：
$$\begin{gathered} {\alpha }_m=0.1\frac{V+40}{1-exp(-(V+40)/10)} \\ {\beta }_m=4exp(-(V+65)/18)\text{\hspace{1em}(1.12)} \end{gathered}$$
$$\begin{gathered} {\alpha }_h=0.07exp(-(V+65)/20) \\ {\beta }_h=\frac{1}{exp(-(V+35)/10)+1}\text{\hspace{1em}(1.13)} \end{gathered}$$

模型计算

简要介绍一下如何计算并实现HH神经元模型，具体参数可参考实现代码
1.通过初始膜电压（静息电位）带入α和β的相关公式将α和β的初始值计算出来。
2.将计算出来的α和β值代入1.7式计算出初始的n，m，h的值。这一步可以理解为n，m，h的初值就是静息态时间无穷大的时候的值
3. 将计算出来的m， n， h待会电流表达式，计算出神经元的总电流
4. 将总电流带入1.1式计算神经元膜电压

实现代码 -python

import numpy as np
import matplotlib.pyplot as plt

num_spike = 0
I_p = 10.0
dt = 0.01
T = 100
t = np.arange(0, T, dt)

m = np.zeros(len(t))
n = np.zeros(len(t))
h = np.zeros(len(t))

alpha_m = np.zeros(len(t))
beta_m = np.zeros(len(t))
alpha_n = np.zeros(len(t))
beta_n = np.zeros(len(t))
alpha_h = np.zeros(len(t))
beta_h = np.zeros(len(t))

V_m = np.zeros(len(t))
g_K = np.zeros(len(t))
g_Na = np.zeros(len(t))

g_K_max = 36
g_Na_max = 120
g_L = 0.3
E_K = -12
E_Na = 115
E_L = 10.6
V_rest = -70
C_m = 1.0

for i in range(1, len(t) - 1):
    # V_m[1] = -65
    alpha_m[i] = 0.1 * ((25 - V_m[i]) / (np.exp((25 - V_m[i]) / 10) - 1))
    beta_m[i] = 4 * np.exp(-1 * V_m[i] / 18)
    alpha_n[i] = 0.01 * ((10 - V_m[i]) / (np.exp((10 - V_m[i]) / 10) - 1))
    beta_n[i] = 0.125 * np.exp(-1 * V_m[i] / 80)
    alpha_h[i] = 0.07 * np.exp(-1 * V_m[i] / 20)
    beta_h[i] = 1 / (np.exp((30 - V_m[i]) / 10) + 1)

    if i == 1:
        m[i] = alpha_m[i] / (alpha_m[i] + beta_m[i])
        n[i] = alpha_n[i] / (alpha_n[i] + beta_n[i])
        h[i] = alpha_h[i] / (alpha_h[i] + beta_h[i])

    g_Na[i] = m[i] ** 3 * g_Na_max * h[i]
    g_K[i] = n[i] ** 4 * g_K_max

    I_Na = g_Na[i] * (V_m[i] - E_Na)
    I_K = g_K[i] * (V_m[i] - E_K)
    I_L = g_L * (V_m[i] - E_L)
    I_ion = I_p - I_K - I_Na - I_L

    V_m[i + 1] = V_m[i] + I_ion / C_m * dt
    m[i + 1] = m[i] + (alpha_m[i] * (1 - m[i]) - beta_m[i] * m[i]) * dt
    n[i + 1] = n[i] + (alpha_n[i] * (1 - n[i]) - beta_n[i] * n[i]) * dt
    h[i + 1] = h[i] + (alpha_h[i] * (1 - h[i]) - beta_h[i] * h[i]) * dt

    if V_m[i] > 0 - V_rest:
        if (V_m[i] > V_m[i - 1]) and (V_m[i] > V_m[i + 1]):
            num_spike = num_spike + 1

V_m = V_m + V_rest

plt.plot(t, V_m)
plt.show()

在这里使用的是欧拉法求解上文中的微分方程。

仿真结果

参考文献

[1]: Principles of computational modelling in neuroscience-Cambridge University Press