Hodgkin-Huxley 模型简介
目录
0. 写在前面
1. 前言
2. 预备知识
2.1. 神经元结构
2.2. 静息膜电位
2.3. 等效电路
2.4. 其他
3. 霍奇金-赫胥黎(HH)模型的基本构成
3.1. 数学表达
3.2. 基于BrainPy的HH Model的实现
4. 结论
0. 写在前面
最近在北京大学吴思老师课题组进行了为期一个月的实习,收获颇丰,现将自己的学习成果进行记录。
记录的第一个模型就是HH Model。
内容参考《神经计算建模实战》吴思等 著 以及网上有有关资料
1. 前言
Hodgkin-Huxley (HH) 模型是由英国生理学家 Alan Lloyd Hodgkin 和 Andrew Huxley 在1952年提出的一种用于解释和模拟神经元如何生成动作电位的数学模型。这项工作为理解神经系统的基本功能提供了理论基础,并因此获得了诺贝尔生理学或医学奖。
有意思的是,艾伦·霍奇金(Allen Hodgkin)和安德鲁·赫胥黎(Andrew Huxley)两人在1939年首次合作,但仅几个月后,第二次世界大战爆发,两人在战争期间均转向研究武器。 战争结束后,两人继续合作,并于1950年提出HH模型。 1963年两人因此共同获得了诺贝尔生理学或医学奖。 可见战争对科技的阻碍以及和平的不易!
2. 预备知识
2.1. 神经元结构
为了了解HH模型的构成,我们需要首先了解神经元的电导模型。
神经元是大脑信息处理的基本单位,其计算功能由形态、电学特性、连接方式和位置决定。神经元主要由细胞体(包含细胞核和细胞器)和两种突起组成:长而细的轴突和粗而密集的树突。轴突末端与下一神经元的树突相连,形成传递信息的突触结构。
2.2. 静息膜电位
神经元的细胞膜是由磷脂双分子层构成的结构,能够有效分隔细胞内外环境。由于磷脂双分子层不允许水溶性或带电分子自由通过,因此需要特定的运输蛋白来帮助这些物质穿过细胞膜。运输蛋白主要包括载体蛋白、通道蛋白和离子泵三类。离子通道作为通道蛋白的一种,能够让如钠离子(Na⁺)、钾离子(K⁺)和氯离子(Cl⁻)等离子通过。根据其是否受到闸门控制,离子通道可以是无闸门的或是由电、化学或物理刺激控制的。在静息状态下,神经元细胞膜对钾离子(K⁺)的通透性高于钠离子(Na⁺),这是由于无闸门控制的钾离子通道数量多于钠离子通道。
为了维持细胞内外离子浓度的差异,神经元依赖钠钾泵进行离子交换,这种泵使用ATP能量将钠离子运出细胞并将钾离子运入细胞内,形成并维持了细胞内外的离子浓度梯度。具体来说,钠钾泵每消耗一个ATP分子就能将2个钾离子移入细胞内,并将3个钠离子移出细胞外。这一过程导致细胞外钠离子(Na⁺)和氯离子(Cl⁻)浓度显著高于细胞内部,而细胞内的钾离子(K⁺)浓度则远高于细胞外部。正是由于这种离子浓度梯度的存在,在静息状态下,细胞膜内外建立了不同的离子分布,为动作电位的发生奠定了基础。
细胞膜内外的离子浓度差产生了一种化学梯度(chemical gradient),这种梯度使得离子倾向于从高浓度区域流向低浓度区域。这部分不是文章的重点,因此不赘述。想进一步了解的可以查询Nernst 方程式和戈德曼-霍奇金-卡茨方程(Goldman-Hodgkin-Katz equation,GHK 方程)以便进一步的了解。
2.3. 等效电路
细胞的电学性质由穿过细胞膜的离子决定,
G
H
K
GHK
GHK方程定量描述了细胞膜电位,但未解释通透性参数
P
x
P_x
Px如何随时间变化。由于离子通道的通透性不断变化,神经细胞展现出复杂的电学特性。为了更好地理解这一动态过程,可以通过绘制神经元的等效电路来描述和分析这些变化。
模型具体包含以下三个部分:
- 代表离子通道的电阻器(resistor)
- 代表离子浓度梯度的电源(battery)
- 代表细胞膜储存电荷能力的电容器(capacitor)
考虑较为实际的神经元的情况,细胞膜对
K
+
K^+
K+渗透,还对
N
a
+
Na^+
Na+和
C
l
−
Cl^−
Cl−渗透,并且此时电路还接收一个外界输入电流
I
(
t
)
I(t)
I(t),如下图所示:
此时细胞膜单位面积的离子电流等于:
i
ion
=
−
g
Cl
(
V
M
−
E
Cl
)
−
g
K
(
V
M
−
E
K
)
−
g
Na
(
V
M
−
E
Na
)
i_{\text{ion}} = -g_{\text{Cl}}(V_{\text{M}} - E_{\text{Cl}}) - g_{\text{K}}(V_{\text{M}} - E_{\text{K}}) - g_{\text{Na}}(V_{\text{M}} - E_{\text{Na}})
iion=−gCl(VM−ECl)−gK(VM−EK)−gNa(VM−ENa)
为表示单位膜面积接收到的电流,此我们将 𝐼(𝑡) 除以神经元的总表面积 𝐴,得到:
根据给定的关系,我们可以写出:
I
(
t
)
A
=
c
M
d
V
M
d
t
+
i
ion
\frac{I(t)}{A} = c_{\text{M}} \frac{dV_{\text{M}}}{dt} + i_{\text{ion}}
AI(t)=cMdtdVM+iion
将离子电流
i
ion
i_{\text{ion}}
iion 代入上述方程,得到:
c
M
d
V
M
d
t
=
−
g
Cl
(
V
M
−
E
Cl
)
−
g
K
(
V
M
−
E
K
)
−
g
Na
(
V
M
−
E
Na
)
+
I
(
t
)
A
c_{\text{M}} \frac{dV_{\text{M}}}{dt} = -g_{\text{Cl}}(V_{\text{M}} - E_{\text{Cl}}) - g_{\text{K}}(V_{\text{M}} - E_{\text{K}}) - g_{\text{Na}}(V_{\text{M}} - E_{\text{Na}}) + \frac{I(t)}{A}
cMdtdVM=−gCl(VM−ECl)−gK(VM−EK)−gNa(VM−ENa)+AI(t)
进一步展开并整理:
c
M
d
V
M
d
t
=
−
(
g
Cl
+
g
K
+
g
Na
)
V
M
+
g
Cl
E
Cl
+
g
K
E
K
+
g
Na
E
Na
+
I
(
t
)
A
c_{\text{M}} \frac{dV_{\text{M}}}{dt} = -(g_{\text{Cl}} + g_{\text{K}} + g_{\text{Na}})V_{\text{M}} + g_{\text{Cl}}E_{\text{Cl}} + g_{\text{K}}E_{\text{K}} + g_{\text{Na}}E_{\text{Na}} + \frac{I(t)}{A}
cMdtdVM=−(gCl+gK+gNa)VM+gClECl+gKEK+gNaENa+AI(t)
为便于处理,上述重写为:
c
M
d
V
M
d
t
=
−
V
M
−
E
R
r
M
+
I
(
t
)
A
c_{\text{M}} \frac{dV_{\text{M}}}{dt} = -\frac{V_{\text{M}} - E_{\text{R}}}{r_{\text{M}}} + \frac{I(t)}{A}
cMdtdVM=−rMVM−ER+AI(t)其中
r
M
=
1
g
Cl
+
g
K
+
g
Na
,
E
R
=
(
g
Cl
E
Cl
+
g
K
E
K
+
g
Na
E
Na
)
r
M
r_{\text{M}} = \frac{1}{g_{\text{Cl}} + g_{\text{K}} + g_{\text{Na}}}, \quad E_{\text{R}} = (g_{\text{Cl}}E_{\text{Cl}} + g_{\text{K}}E_{\text{K}} + g_{\text{Na}}E_{\text{Na}})r_{\text{M}}
rM=gCl+gK+gNa1,ER=(gClECl+gKEK+gNaENa)rM
式中,
r
M
r_M
rM表示单位膜电阻,
E
R
E_R
ER表示细胞的静息膜电位。
对于一个电学特性不变化的被动细胞膜来说,当电导和外界电流恒定时,细胞膜电位最终将达到一个稳态:
V
ss
=
g
Cl
E
Cl
+
g
K
E
K
+
g
Na
E
Na
+
I
A
g
Cl
+
g
K
+
g
Na
V_{\text{ss}} = \frac{g_{\text{Cl}}E_{\text{Cl}} + g_{\text{K}}E_{\text{K}} + g_{\text{Na}}E_{\text{Na}} + \frac{I}{A}}{g_{\text{Cl}} + g_{\text{K}} + g_{\text{Na}}}
Vss=gCl+gK+gNagClECl+gKEK+gNaENa+AI
在没有外界电流的情况下,稳态膜电位是三个离子电流的平衡电位的加权和:
V
ss
,
I
=
0
=
E
R
=
g
Cl
E
Cl
+
g
K
E
K
+
g
Na
E
Na
g
Cl
+
g
K
+
g
Na
V_{\text{ss}, I=0} = E_{\text{R}} = \frac{g_{\text{Cl}}E_{\text{Cl}} + g_{\text{K}}E_{\text{K}} + g_{\text{Na}}E_{\text{Na}}}{g_{\text{Cl}} + g_{\text{K}} + g_{\text{Na}}}
Vss,I=0=ER=gCl+gK+gNagClECl+gKEK+gNaENa
方程中所需要的电导和平衡电位数值可以很容易通过实验获得。
2.4. 其他
笔者是计算机专业的,对于生物学的知识不是很清楚,因此只是为了了解HH Model的动力学以及如何编程即可。关于离子电流的测定:电压钳 、泄漏电流的测定、 I N a I_{Na} INa 和 I K I_K IK 的测定。感兴趣的读者可以自行百度或 G o o g l e Google Google。
3. 霍奇金-赫胥黎(HH)模型的基本构成
3.1 数学表达
经过各种生理实验和数学推导Hodgkin 和 Huxley 最终成功搭建了 HH 模型。总结起来,HH 模型是一个四变量模型,每个变量对应一个常微分方程,其中一个是关于膜电位 𝑉 的方程,另外三个是关于离子通道门控变量 𝑚、ℎ 和 𝑛 的方程。HH 模型的表达式为
{
c
d
V
d
t
=
−
g
ˉ
Na
m
3
h
(
V
−
E
Na
)
−
g
ˉ
K
n
4
(
V
−
E
K
)
−
g
ˉ
L
(
V
−
E
L
)
+
I
ext
d
n
d
t
=
ϕ
[
α
n
(
V
)
(
1
−
n
)
−
β
n
(
V
)
n
]
d
m
d
t
=
ϕ
[
α
m
(
V
)
(
1
−
m
)
−
β
m
(
V
)
m
]
d
h
d
t
=
ϕ
[
α
h
(
V
)
(
1
−
h
)
−
β
h
(
V
)
h
]
\left\{ \begin{aligned} c \frac{dV}{dt} &= -\bar{g}_{\text{Na}} m^3 h (V - E_{\text{Na}}) - \bar{g}_{\text{K}} n^4 (V - E_{\text{K}}) - \bar{g}_{\text{L}} (V - E_{\text{L}}) + I_{\text{ext}} \\ \frac{dn}{dt} &= \phi [\alpha_n(V) (1 - n) - \beta_n(V)n] \\ \frac{dm}{dt} &= \phi [\alpha_m(V) (1 - m) - \beta_m(V)m] \\ \frac{dh}{dt} &= \phi [\alpha_h(V) (1 - h) - \beta_h(V)h] \end{aligned} \right.
⎩
⎨
⎧cdtdVdtdndtdmdtdh=−gˉNam3h(V−ENa)−gˉKn4(V−EK)−gˉL(V−EL)+Iext=ϕ[αn(V)(1−n)−βn(V)n]=ϕ[αm(V)(1−m)−βm(V)m]=ϕ[αh(V)(1−h)−βh(V)h]
上式中,各个门控变量的转换速率表达式如下:
α
n
(
V
)
=
0.01
(
V
+
55
)
1
−
exp
(
−
V
+
55
10
)
,
β
n
(
V
)
=
0.125
exp
(
−
V
+
65
80
)
\alpha_n(V) = \frac{0.01(V + 55)}{1 - \exp\left(-\frac{V + 55}{10}\right)}, \quad \beta_n(V) = 0.125 \exp\left(-\frac{V + 65}{80}\right)
αn(V)=1−exp(−10V+55)0.01(V+55),βn(V)=0.125exp(−80V+65)
α
h
(
V
)
=
0.07
exp
(
−
V
+
65
20
)
,
β
h
(
V
)
=
1
exp
(
−
V
+
35
10
)
+
1
\alpha_h(V) = 0.07 \exp\left(-\frac{V + 65}{20}\right), \quad \beta_h(V) = \frac{1}{\exp\left(-\frac{V + 35}{10}\right) + 1}
αh(V)=0.07exp(−20V+65),βh(V)=exp(−10V+35)+11
α
m
(
V
)
=
0.1
(
V
+
40
)
1
−
exp
(
−
V
+
40
10
)
,
β
m
(
V
)
=
4
exp
(
−
V
+
65
18
)
\alpha_m(V) = \frac{0.1(V + 40)}{1 - \exp\left(-\frac{V + 40}{10}\right)}, \quad \beta_m(V) = 4 \exp\left(-\frac{V + 65}{18}\right)
αm(V)=1−exp(−10V+40)0.1(V+40),βm(V)=4exp(−18V+65)
式中,
I
ext
I_{\text{ext}}
Iext 是外界输入的电流,
ϕ
\phi
ϕ 是温度因子(temperature factor)。进行实验的温度非常重要——由于离子通道的开放关闭本质上是随机的,它们对温度很敏感,因此开关状态的转换速率与温度成指数关系。温度越高,开放关闭的转换速度越快。特别地,温度系数与温度的关系表达为:
ϕ
=
Q
10
(
T
−
T
base
)
/
10
\phi = Q_{10}^{(T - T_{\text{base}})/10}
ϕ=Q10(T−Tbase)/10式中,
Q
10
Q_{10}
Q10 是温度升高 10°C 的速率比。对于乌贼的巨轴突,
T
base
=
6.
3
∘
C
T_{\text{base}} = 6.3^\circ C
Tbase=6.3∘C 并且
Q
10
=
3
Q_{10} = 3
Q10=3
3.2 基于BrainPy的HH Model的实现
理论知识介绍完毕后,我们接下来探讨如何通过编程实现Hodgkin-Huxley (HH) 模型。由于吴思老师课题组开发了基于JAX的BrainPy的脑动力学编程框架,因此笔者主要在北大使用基于BrainPy的编程环境进行学习和研究,以下代码示例将基于BrainPy框架进行展示。尽管其他神经计算框架(如NEURON、Brian2等)也能实现类似的模型,但限于笔者的能力和经验,本文仅展示基于BrainPy的具体实现。
关于BrainPy的相关操作详见:BrainPy文档
我们首先定义相关的HH 类:
import brainpy as bp
import brainpy.math as bm
import matplotlib.pyplot as plt
# from scipy.special import kwargs
# 生成中文字体
plt.rcParams.update({"font.size": 15})
plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False
class HH(bp.dyn.NeuDyn):
def __init__(self, size, ENa=50., gNa=120., EK=-77., gK=36.,
EL=-54.387, gL=0.03, V_th=20., C=1.0, T=6.3):
# 初始化
super(HH, self).__init__(size=size)
# 定义神经元参数
self.ENa = ENa
self.EK = EK
self.EL = EL
self.gNa = gNa
self.gK = gK # 钾离子通道的最大电导值
self.gL = gL
self.C = C
self.V_th = V_th # 阈值电压
self.Q10 = 3. # 温度身高10度的速率比
self.T_base = 6.3
self.phi = self.Q10 ** ((T - self.T_base) / 10)
# 定义神经元变量,一定要定义变量,只有将其声明为变量,在做JIT的时候,量才会改变,只限于BrainPy
self.V = bm.Variable(-70.68 * bm.ones(self.num)) # 膜电位
# 这里的self.num实际上是根据上面的size来进行给定的,表示的是总体的size是多大 ,eg (10×10,100),其num就是10000
# 膜电位的初始值
self.m = bm.Variable(0.0266 * bm.ones(self.num)) # 离子通道m
self.h = bm.Variable(0.772 * bm.ones(self.num)) # 离子通道h
self.n = bm.Variable(0.235 * bm.ones(self.num)) # 离子通道n
# 神经元接收到的输入电流,外部输入电流
self.input = bm.Variable(bm.zeros(self.num))
# 神经元发放状态
self.spike = bm.Variable(bm.zeros(self.num, dtype=bool))
# 神经元上次发放的时刻
self.t_last_spike = bm.Variable(bm.ones(self.num) * -1e7)
# 定义积分函数,指数欧拉
self.integral = bp.odeint(f=self.derivative, method='exp_auto')
@property # 使用这个将类的方法转换为只读属性
def derivative(self): # 将上述公式联立在一起,不联立的话精度就会差一点
return bp.JointEq([self.dV, self.dm, self.dh, self.dn])
# 算的是HH模型的第三个式子,就是对着对应的数学式子编码即可
def dm(self, m, t, V):
alpha = 0.1 * (V + 40) / (1 - bm.exp(-(V + 40) / 10))
beta = 4.0 * bm.exp(-(V + 65) / 18)
dmdt = alpha * (1 - m) - beta * m
return self.phi * dmdt
# 算的是HH模型的第四个式子
def dh(self, h, t, V):
alpha = 0.07 * bm.exp(-(V + 65) / 20.)
beta = 1 / (1 + bm.exp(-(V + 35) / 10))
dhdt = alpha * (1 - h) - beta * h
return self.phi * dhdt
def dn(self, n, t, V):
alpha = 0.01 * (V + 55) / (1 - bm.exp(-(V + 55) / 10))
beta = 0.125 * bm.exp(-(V + 65) / 80)
dndt = alpha * (1 - n) - beta * n
return self.phi * dndt
def f_gNa(self, m, h):
return (self.gNa * m ** 3.0 * h)
def f_gK(self, n):
return (self.gK * n ** 4.0)
# 定义的第一个式子,t前面的就是变量,后面的就是参数
def dV(self, V, t, m, h, n):
I_Na = (self.gNa * m ** 3.0 * h) * (V - self.ENa)
I_K = (self.gK * n ** 4.0) * (V - self.EK)
I_leak = self.gL * (V - self.EL) # 泄露电流
dVdt = (- I_Na - I_K - I_leak + self.input) / \
self.C # self.input是外部输入电流
return dVdt
# 更新函数:每个时间步都会运行此函数完成变量更新
def update(self):
t = bp.share.load('t')
dt = bp.share.load('dt')
# 更新下一时刻变量的值
V, m, h, n = self.integral(self.V, self.m, self.h, self.n, t, dt=dt)
# 判断神经元是否产生膜电位,上一个时刻小于这个阈值,下一个时刻大于这个阈值
self.spike.value = bm.logical_and(self.V < self.V_th, V >= self.V_th)
# 更新神经元发放的时间
self.t_last_spike.value = bm.where(self.spike, t, self.t_last_spike)
self.V.value = V
self.m.value = m
self.h.value = h
self.n.value = n
self.input[:] = 0. # 重置神经元接收到的输入
当我们定义好一个HH Model类后,我们通过模拟并可视化不同大小的外部输入电流对HH模型神经元膜电位的影响,展示神经元在不同刺激条件下的动态行为:
# 设置函数时长为2ms的不同大小的外部输入电流
def simple_run():
current_sizes = [1., 2., 4., 8., 10., 15.] # 定义了不同的电流大小
currents, length = bp.inputs.section_input( # 生成一个电流输入序列
values=[0., bm.asarray(current_sizes), 0.], # 三个阶段的值
durations=[10, 2, 25], # 三个阶段的持续时间
return_length=True # 返回仿真的时长
)
hh = HH(currents.shape[1]) # shape[1]表示的是神经元的个数
runner = bp.DSRunner(hh, # 传入HH模型的实例
monitors=['V', 'm', 'h', 'n'], # 监控四个变量
inputs=['input', currents, 'iter'])
runner.run(length)
# 可视化
# line_styles = [(0, (1, 1)), (0, (5, 5)), (0, (5, 1)), 'dotted', 'dashed', 'solid', ]
line_styles = ['solid', 'dashed', 'dotted', ] # 实线,虚线,点线
colors = ['tab:gray', 'black'] # 灰色,黑色
for i in range(hh.num):
plt.plot(runner.mon.ts, runner.mon.V[:, i], # 提取所有行和第一列的数据
label=f'I={current_sizes[i]} mA',
linestyle=line_styles[i % 3],
color=colors[i // 3]
)
# plt.xlabel('t (ms)')
plt.ylabel('V (mV)')
# 将电流变化画在膜电位变化的下方
plt.plot(runner.mon.ts,
bm.where(currents[:, -1] > 0, 10., 0.).to_numpy() - 90,
label='Current',
linestyle=(0, (1, 1)))
plt.tight_layout() # 调整子图参数,使图形布局紧凑
plt.legend()
plt.show()
运行结果如下:
注意到 HH 模型存在这“全或无”的特征:当施加的电流低于某一阈值时,膜电位迅速恢复到静息状态;当电流超过某个阈值时,就会产生动作电位。
下面探究 HH 模型对不同大小的外界输入的响应,使用 b r a i n p y . i n p u t s . s e c t i o n i n p u t ( ) brainpy.inputs.section_input() brainpy.inputs.sectioninput()函数来设置时长 2 毫秒的不同大小的外界输入,
# 和simple_run逻辑一致
def visualize_hh_responses_in_book():
# currents
current_sizes = [1., 2., 4., 8., 10., 15.]
currents, length = bp.inputs.section_input(
values=[0., bm.asarray(current_sizes), 0.],
# 电流在前10ms的时候,全是0,之后的2ms对应的是一个array,然后之后的25ms还是0
durations=[10, 2, 25],
return_length=True
)
# simulation
hh = HH(currents.shape[1])
runner = bp.DSRunner(hh, monitors=['V', 'm', 'h', 'n'],
inputs=['input', currents, 'iter'])
# hh传入的HH模型的实例,需要监控四个变量,指定仿真输入方式为迭代方式
runner.run(length)
# visualization
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
# fig用于容纳绘图内容,gs用于管理子图布局
ax = fig.add_subplot(gs[0, 0]) # 向fig中添加一个子图,gs[0,0]表示网络布局中的第一行、第一列的单元格
plt.plot(runner.mon.ts, runner.mon.V[:, 0]) # 记录膜电位V的时间序列,选择的是第一个神经元的膜电位序列
plt.annotate(f'I={current_sizes[0]} mA', xy=(12.2, -69), xytext=(14.2, -83),
arrowprops=dict(arrowstyle="->"))
# 在图中添加注释,选择的是带箭头的直线,从注释文本指向目标点
plt.plot(runner.mon.ts, runner.mon.V[:, 1])
plt.annotate(f'I={current_sizes[1]} mA', xy=(11.7, -68.5), xytext=(1.5, -63),
arrowprops=dict(arrowstyle='fancy', color='green', lw=1))
# lw设置线宽的大小
plt.plot(runner.mon.ts, runner.mon.V[:, 2])
plt.annotate(f'I={current_sizes[2]} mA', xy=(16.3, -61), xytext=(20.5, -55),
arrowprops=dict(arrowstyle='fancy', color='green', lw=1))
plt.plot(runner.mon.ts, runner.mon.V[:, 3])
plt.annotate(f'I={current_sizes[3]} mA', xy=(12.9, -4), xytext=(2.9, -1),
arrowprops=dict(arrowstyle='fancy', color='green', lw=1))
plt.plot(runner.mon.ts, runner.mon.V[:, 4])
plt.annotate(f'I={current_sizes[4]} mA', xy=(12.5, 8), xytext=(2.5, 11),
arrowprops=dict(arrowstyle='fancy', color='green', lw=1))
plt.plot(runner.mon.ts, runner.mon.V[:, 5])
plt.annotate(f'I={current_sizes[5]} mA', xy=(12, 19), xytext=(2, 22),
arrowprops=dict(arrowstyle='fancy', color='green', lw=1))
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.ylabel(r'$V$ (mV)')
plt.xlabel(r'$t$ (ms)')
# 将电流变化画在膜电位变化的下方
plt.plot(runner.mon.ts, bm.where(
currents[:, -1] > 0, 10., 0.).to_numpy() - 90)
plt.annotate('Current', xy=(30, -87), xytext=(30, -87), )
# plt.savefig('HH_responses_with_different_currents.pdf',transparent=True, dpi=500)
plt.show()
这张图片展示的与上一张没啥区别,只是更加的直观
def visualize_constant_current_response_in_book():
# currents
currents, length = bp.inputs.section_input(values=[0., 10., 0.],
durations=[10, 50, 10],
return_length=True)
# simulation
hh = HH(1)
runner = bp.DSRunner(hh, monitors=['V', 'm', 'h', 'n'],
inputs=['input', currents, 'iter'])
runner.run(length)
# visualization
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
ax = fig.add_subplot(gs[0, 0])
plt.plot(runner.mon.ts, runner.mon.V, label='V')
plt.plot(runner.mon.ts, currents - 90, label='Current (I=5mA)')
plt.xlabel(r'$t$ (ms)')
plt.ylabel(r'$V$ (mV)')
plt.text(0, -67, r'$V$')
plt.text(0, -88, r'$I$=5mA')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# plt.savefig('HH_responses_with_constant_current.pdf', transparent=True, dpi=500)
# transparent 用来 控制是否使图像背景透明
plt.show()
下图可以发现:当输入电流足够大并保持足够长的时间时,模型将产生周期性响应
下面探究 HH 神经元对一个持续时长的输入的响应
def conductance_during_action_potential():
# currents
currents, length = bp.inputs.section_input(values=[0., 10., 0.], durations=[
10, 5, 20], return_length=True)
# simulation
hh = HH(1)
runner = bp.DSRunner(hh, monitors=['V', 'm', 'h', 'n'],
inputs=['input', currents, 'iter'])
runner.run(length)
# visualization
fig, gs = bp.visualize.get_figure(3, 1, 2, 8)
ax = fig.add_subplot(gs[0, 0])
plt.plot(runner.mon.ts, runner.mon.V, lw=2) # lw设置线条宽度
plt.ylabel(r'$V$ (mV)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax = fig.add_subplot(gs[1, 0])
plt.plot(runner.mon.ts, hh.f_gNa(
runner.mon.m, runner.mon.h), label='gNa', lw=2)
plt.plot(runner.mon.ts, hh.f_gK(runner.mon.n), label='gK', lw=2)
plt.text(13.8, 33, 'gNa')
plt.text(16.5, 10, 'gK')
plt.ylabel('Conductance')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax = fig.add_subplot(gs[2, 0])
plt.plot(runner.mon.ts, runner.mon.m, label='m', lw=2)
plt.plot(runner.mon.ts, runner.mon.h, label='h', lw=2)
plt.plot(runner.mon.ts, runner.mon.n, label='n', lw=2)
plt.text(14.5, 0.94, 'm')
plt.text(26.1, 0.60, 'h')
plt.text(26.1, 0.31, 'n')
plt.ylabel('Channel')
plt.xlabel(r'$t$ (ms)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# 保存为pdf文件
# plt.savefig('HH_conductance_during_AP.pdf', transparent=True, dpi=500)
plt.show()
上图展示了在 HH 神经元一次动作电位爆发期间钠离子和钾离子通道的电导变化情况以及各个门控变量的变化情况。可以发现在动作电位产生的初始阶段,𝑚 和 𝑛 逐渐增加,ℎ 逐渐减小。并且,𝑚 增长最快,导致钠离子通道的电导首先增大,大量
N
a
+
Na^+
Na+ 涌入神经元细胞。由
I
Na
=
g
Na
(
V
−
E
Na
)
I_{\text{Na}} = g_{\text{Na}} (V - E_{\text{Na}})
INa=gNa(V−ENa)可知,在这一阶段,增大的
g
Na
g_{\text{Na}}
gNa 会使神经元膜电位被驱动力
(
V
−
E
Na
)
(V - E_{\text{Na}})
(V−ENa)拉向钠离子通道的反转电位
E
Na
≈
55
mV
E_{\text{Na}} \approx 55 \, \text{mV}
ENa≈55mV。这一阶段对应着神经元动作电位的上升过程。随后,钾离子通道 𝑛 逐渐打开,使得钾离子通道的电导
g
K
g_{\text{K}}
gK 逐渐增大,越来越多的
K
+
K^+
K+ 流出细胞外;与此同时,钠离子通道的失活闸门 ℎ 逐渐发挥作用,钠离子通道被阻塞,
g
Na
g_{\text{Na}}
gNa减小。随着
g
K
g_{\text{K}}
gK 占主导,神经元膜电位被驱动力
(
V
−
E
K
)
(V - E_{\text{K}})
(V−EK)拉向
K
+
K+
K+ 的反转电位
E
K
≈
−
77
mV
E_{\text{K}} \approx -77 \, \text{mV}
EK≈−77mV。这对应着动作电位的下降过程。
if __name__ == '__main__':
pass
simple_run()
visualize_hh_responses_in_book()
visualize_constant_current_response_in_book()
conductance_during_action_potential()
4. 结论
到此,我们就已经全部介绍完了HH Model的相关知识,当然,限于笔者的能力,不能全部涉及到,同时文章中有些许错误之处,还望之处,并多多包涵。
HH Model确实很好的模拟了神经元的动作电位,又能解释其背后的生物学原理,但是由于较为复杂的原因,后续引入了大量的简化神经元模型,详细的模型后续进行介绍。