Hodgkin Huxley Model

参考

[1] Neuron Dynamics

[2] Spiking Neurons

[3] The Hodgkin Huxley Model

[4] 动作电位的形成机制

[5] 霍奇金赫胥黎枪乌贼巨型轴突

1952年发表，1963年获诺贝尔生理学奖

通过对乌贼的巨型轴突进行实验，运用数学方法分析以及三种离子描述了细胞膜电势的动态变化，是其他简化神经元模型的基础。

为什么选择乌贼的巨型轴突呢？

鱿鱼巨型轴突以其异常大的直（〜0.5毫米）而著称。鱿鱼神经系统和其他神经系统中的大多数轴突通常至少薄100倍。鱿鱼巨型轴突的大尺寸是一种专门用于迅速传导动作电位的动作电位，当其从掠食者逃脱时会触发鱿鱼的地幔收缩。除了对鱿鱼有益之外，大直径轴突对实验也有利，因为它允许进行迄今为止在生物物理研究中使用的较小轴突在技术上不可行的操纵。

细胞膜两侧离子浓度差是产生了电势差，在神经元动力学中有重要作用。

电势差导致离子浓度差，离子浓度差产生电势差（奈奎斯特电势）Δu=Ktqln⁡n1n2\Delta u = \frac{Kt}{q} \ln \frac{n_1}{n_2}Δu=qKtlnn2n1，其中n1n_1n1为膜外离子浓度，n2n_2n2为膜内离子浓度。
反电势：当电势差小于ENaE_{Na}ENa时，钠离子进入，减小这个差异；反之，当电势差大于ENaE_{Na}ENa时，钠离子流出，减小这个差异。即当Δu\Delta uΔu从ENaE_{Na}ENa开始变化时，电流的方向就与其变化方向相反，使得其重新回到平衡电势，又称反电势。

这里解释一下为什么膜内钾离子浓度高，而其平衡电势确实膜外为正：因为每个钾离子携带一个正电荷，当钾离子通过通道离开细胞时，同时带出了正电荷，把与其配对的负电荷留在膜内.。正负电荷在空间上分离，不就会产生电场吗？电荷是同性相斥，异性相吸。这个横跨细胞膜的电场会形成一种阻止钾离子外流的力场。随着钾离子不断流出，电荷在细胞膜两侧不断积累。阻止钾离子外流的力场也不断增加,很快达到与驱动离子外流的扩散力势均力敌，钾离子外流就会停止。

静息电位urest≈−65mVu_{rest} \approx -65mVurest≈−65mV

HH Model

对巨型章鱼轴突进行实验，主要考察钠、钾和漏（以cl−cl^-cl−为代表）三种通道。

通过实验，Hodgkin和Huxley发现这三个通道的并行关系（控制变量进行实验），宏观电流可以被认为是大量嵌入在细胞膜上的微观离子通道综合作用引起的，所以有了上面的等效电路图，外界电流（电荷）的通路有四条：一条用于给电容充电，另外三条用于流过电阻器。

细胞膜是一个半透膜，将膜内外液体分离，它能够存储电荷，起到了电容的作用，记作CCC
钠钾离子通道可以通过电流，其允许离子通过的能力可以用电阻RNaR_{Na}RNa、RKR_{K}RK，或电导ggg表示（为了表达的简便），g=1/Rg=1/Rg=1/R，或者说离子流过的每一个状态都对应一个电阻，因为细胞膜是有对应的离子泵，其通过相应离子的能力可以改变，因此是一个可变电阻器
因为细胞膜内外离子浓度不同，与之构建起的平衡的每种离子的电场也不同，根据奈奎斯特公式，每个离子对应一个平衡电势，用ENaE_{Na}ENa、EKE_{K}EK表示
对于Leaky通道，以cl−cl^-cl−为代表，其等效为RRR，为常数，平衡电势用ELE_LEL表示，该通道的电流相较而言较小

根据电学公式C=q/UC=q/UC=q/U、IC=CdudtI_C=C\frac{du}{dt}IC=Cdtdu以及并联电流公式，可得
I(t)=IC(t)+∑kIk(t)+IL(t)Cdudt=−∑kIk(t)+I(t)I(t) = I_C(t) + \sum \limits_k I_k(t) + I_L(t)\\ C\frac{du}{dt} = -\sum\limits_k I_k(t) + I(t) I(t)=IC(t)+k∑Ik(t)+IL(t)Cdtdu=−k∑Ik(t)+I(t)
式中，∑kIk\sum\limits_kI_kk∑Ik表示通过细胞膜的离子电流的和，例如对漏通道，考虑到平衡电势，IL=gL(u−EL)I_L = g_L(u-E_L)IL=gL(u−EL)。

关于III的方向问题：惯例是电流由膜内流向膜外为正。

平衡电势EEE可测，膜电势uuu以及电流III都可以测量，但各离子通道的电导ggg该如何测量呢？Hodgkin和Huxley成功测量了各通道的电阻（电导）是如何随时间和膜电势变化的，他们引入门控变量mmm，nnn，hhh来对给定时间通道开启的概率进行建模，得出如下关系式。这里mmm和nnn共同负责调控NaNaNa离子通道，hhh负责调控KKK离子通道。mmm，nnn，hhh是根据相应的时间常数不断变化的。

动作电位的形成

到这里，谈一下动作电位的形成过程：静息状态时，NaNaNa离子激活态门处于关闭状态而失活态门处于开放状态，细胞膜受刺激时对NaNaNa离子的通透性增加，即去极化阶段，表现在NaNaNa离子的两个门都处于开放阶段，当膜电势达到阈值时，膜对NaNaNa离子的通透性显著增加，膜电位迅速由负变正接近Na离子的平衡电位，但未真正达到，此时NaNaNa离子失活态门开始关闭，Na离子的通透性开始下降到静息水平，K离子的通道门缓慢开放，导致动作电位从峰值返回到静息状态水平，称为再极化，变化的膜电势时NaNaNa离子通道完全关闭，由于KKK离子通道关闭的速度较慢，使得稍微过量的KKK离子流出使得细胞膜电势比静息电位还低，形成超极化电位。

可见，NaNaNa离子的流通与两个门有关，也就是mmm和hhh，KKK离子的流通与一个门有关，也就是nnn。

门控变量

每个离子通道都有两种状态，允许和不允许状态（permissive or non-permissive）。用pi(t)p_i(t)pi(t)表示某离子通道某时刻的概率，值的大小从0到1变化，由不允许到允许转换的速率为αi(t)\alpha_i(t)αi(t)，由允许到不允许状态转换的速率为βi(t)\beta_i(t)βi(t)。于是对打开概率进行建模为dpidt=αi(u)(1−pi)+βi(u)pi\frac{dp_i}{dt} = \alpha_i(u)(1-p_i) + \beta_i(u)p_idtdpi=αi(u)(1−pi)+βi(u)pi。当某个离子通道打开后，相应的宏观电流也发生变化，因此有Gk=gk∏ipiG_k = g_k \prod \limits _i p_iGk=gki∏pi，其中gkg_kgk表示当所有离子通道打开时的归一化系数。
dmdt=αm(u)(1−m)+βm(u)mdhdt=αh(u)(1−h)+βh(u)hdndt=αn(u)(1−n)+βn(u)n\frac{dm}{dt}=\alpha_m(u)(1-m) + \beta_m(u)m\\ \frac{dh}{dt}=\alpha_h(u)(1-h) + \beta_h(u)h\\ \frac{dn}{dt}=\alpha_n(u)(1-n) + \beta_n(u)n dtdm=αm(u)(1−m)+βm(u)mdtdh=αh(u)(1−h)+βh(u)hdtdn=αn(u)(1−n)+βn(u)n

允许与不允许状态转换速率是与膜电势有关的，实验测得六个参数分别是

xxx	αx(u/mV)\alpha_x(u/mV)αx(u/mV)	βx(u/mV)\beta_x(u/mV)βx(u/mV)
nnn	(0.1−0.01u)/[exp⁡(1−0.1u)−1](0.1-0.01u)/[\exp(1-0.1u)-1](0.1−0.01u)/[exp(1−0.1u)−1]	0.125exp⁡(−u/80)0.125\exp(-u/80)0.125exp(−u/80)
mmm	(2.5−0.1u)/[exp⁡(2.5−0.1u)−1](2.5-0.1u)/[\exp(2.5-0.1u)-1](2.5−0.1u)/[exp(2.5−0.1u)−1]	4exp⁡(−u/18)4\exp(-u/18)4exp(−u/18)
hhh	0.07exp⁡(−u/20)0.07\exp(-u/20)0.07exp(−u/20)	1/[exp⁡(3−0.1u)+1]1/[\exp(3-0.1u)+1]1/[exp(3−0.1u)+1]

初始化时，x0=αi/(αi+βi)x_0 = \alpha_i / (\alpha_i+\beta_i)x0=αi/(αi+βi)

综上，离子通道的电流之和可以用下式表述：
∑kIk=gNam3h(u−ENa)+gKn4(u−EK)+gL(u−EL)\sum\limits_k I_k=g_{Na}m^3h(u-E_{Na}) + g_Kn^4(u-E_K) + g_L(u-E_L) k∑Ik=gNam3h(u−ENa)+gKn4(u−EK)+gL(u−EL)
通过上述公式，我们可以对膜电势进行建模，使用恒定电流刺激，仿真结果如下：

上图为单脉冲，当恒定输入很强时，神经元会产生脉冲序列，如下图所示：

脉冲发放率

可以发现，在某强度下恒定电流刺激时，神经元发放一定数量的脉冲，则单位时间下神经元的脉冲发放数量即为脉冲发放率。对各电流强度下进行实验，即可得到F-I图：

激活阈值：当神经元膜电势超过阈值，神经元发放脉冲

当刺激电流比较小的时候，神经元不会发放脉冲信号，其膜电势会不断leaky，直到回落到静息电位。那么这能说明存在一个电流阈值使得其刺激神经元时介于发放与不发放之间吗？我们做一个实验，首先实验一个电流强度和持续刺激时间合适，使得神经元不发放脉冲（比放好不发放脉冲的程度小一点），然后将电流强度提升一倍，刺激时间减小一倍，再去刺激神经元发现结果和之前一样，这说明不是电流强度决定是否发放脉冲，而是累计的电荷量，也就是膜电势阈值决定脉冲发放。但这在HH模型中并未显式提出。

问题来了，是不是当电流强到一定程度，脉冲发放率会无穷的大？这就要提到动作电位的不应期（Refractory）。

不应期：大多数 Na+ 通道失活时被称为绝对不应期

神经元动作电位先后经过去极化、再极化、超极化阶段，有峰值电位回落时回到比静息电位还要低，此时离子泵要经过调整才能回到静息状态，准备再次发放脉冲，也就是说，在不应期中，无论加多大刺激，神经元都不会再发放脉冲了。

MATLAB实现

具体步骤如下：

1、设置参数，包括平衡电势、电容、电阻、静息电位、最大电导、dt、以及仿真时长

2、循环开始

更新α\alphaαβ\betaβ

若i为1，利用上述参数得到初始化mnh

根据mnh更新电导

计算各通道电流，得到−∑kIk(t)+I(t)-\sum\limits_k I_k(t) + I(t)−k∑Ik(t)+I(t)

更新膜电势

根据膜电势更新mnh

%this a Matlab implement of HH Model%=======================================
% Function name: run_model
% Inputs (2): -I_inj, the injection current array in uA/cm^2
%             -dt, the time step in ms
%
% Function description:
% This function runs the Hodgkin Huxley model for a given injection current
% and time step over 100 ms. The membrane voltage, K+ conductance, and Na+
% conductance are plotted afterwards.
%=======================================
function run_model(I_inj, dt)
global num_spike
num_spike = 0;
% Initialize time vector and time step
t = 0:dt:100; %Time (ms)% Initialize m, n, and h, which represent K+ channel activation, Na+
% channel activation, and Na+ channel inactivation, respectively.
m = zeros(1, length(t));
n = zeros(1, length(t));
h = zeros(1, length(t));% Initialize alphas and betas, the rate constants
alpha_m = zeros(1,length(t));
beta_m = zeros(1,length(t));
alpha_n = zeros(1,length(t));
beta_n = zeros(1,length(t));
alpha_h = zeros(1,length(t));
beta_h = zeros(1,length(t));% Initialize membrane voltage and conductance vectors
V_m = zeros(1,length(t));
g_K = zeros(1,length(t));
g_Na = zeros(1, length(t));% Constants
g_K_max = 36;       %K+ channel max conductance (mS/cm^2)
g_Na_max = 120;     %Na+ channel max conductance (mS/cm^2)
g_L = 0.3;      %Leakage current conductance (mS/cm^2)
E_K = -12;      %K+ channel Nernst potential (mV)
E_Na = 115;     %Na+ channel Nernst potential (mV)
E_L = 10.6;     %Leakage channel Nernst potential (mV)
V_rest = -70;   %Resting voltage (mV)
C_m = 1.0;      %Membrance capacitance (uF/cm^2)for i = 1:(length(t)-1)
%     V_m(1) = 60;% Calculate all alpha and beta valuesalpha_m(i) = 0.1*((25 - V_m(i))/(exp((25-V_m(i))/10)-1));beta_m(i) = 4*exp(-1*V_m(i)/18);alpha_n(i) = 0.01*((10-V_m(i))/(exp((10-V_m(i))/10)-1));beta_n(i) = 0.125*exp(-1*V_m(i)/80);alpha_h(i) = 0.07*exp(-1*V_m(i)/20);beta_h(i) = 1/(exp((30-V_m(i))/10)+1);% Initial conditionsif(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));end% Calculate conductancesg_Na(i) = m(i)^3*g_Na_max*h(i);g_K(i) = n(i)^4*g_K_max;% Calculate currentsI_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_inj(i) - I_K - I_Na - I_L;% Calculate membrane voltage, m, n, and h using Euler's methodV_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;%check spikeif V_m(i) > 0 - V_restif (V_m(i)>V_m(i-1))  && (V_m(i)>V_m(i+1))num_spike  = num_spike + 1;endend
endV_m = V_m + V_rest;

运行结果图如下：

SNN系列｜神经元模型篇(1) Hodgkin Huxley相关推荐

深入理解javascript函数系列第二篇——函数参数
前面的话 javascript函数的参数与大多数其他语言的函数的参数有所不同.函数不介意传递进来多少个参数,也不在乎传进来的参数是什么数据类型,甚至可以不传参数.本文是深入理解javascript函数 ...
Webpack系列-第一篇基础杂记
系列文章 Webpack系列-第一篇基础杂记 Webpack系列-第二篇插件机制杂记 Webpack系列-第三篇流程杂记前言公司的前端项目基本都是用Webpack来做工程化的,而Webpack虽然 ...
SQL Server调优系列进阶篇（如何维护数据库索引）
前言上一篇我们研究了如何利用索引在数据库里面调优,简要的介绍了索引的原理,更重要的分析了如何选择索引以及索引的利弊项,有兴趣的可以点击查看. 本篇延续上一篇的内容,继续分析索引这块,侧重索引项的日常 ...
pytorch 指定卡1_[原创][深度][PyTorch] DDP系列第一篇：入门教程
引言 DistributedDataParallel(DDP)是一个支持多机多卡.分布式训练的深度学习工程方法.PyTorch现已原生支持DDP,可以直接通过torch.distributed使用,超 ...
用番茄工作法提升工作效率（三）工作任务的管理（系列完结篇）
一.本文主题程咬金有三板斧,本文章是本系列文章的最后一篇. 将介绍工作任务的管理,主要涉及周报,月报,年报的等工作计划的管理.有人可能会有疑问:什么?这有什么可介绍的,按照公司要求就可以了. 不同的 ...
微信小程序购物商城系统开发系列-工具篇
微信小程序购物商城系统开发系列-工具篇微信小程序开放公测以来,一夜之间在各种技术社区中就火起来啦.对于它估计大家都不陌生了,对于它未来的价值就不再赘述,简单一句话:可以把小程序简单理解为一个新的操 ...
SQL Server调优系列基础篇（子查询运算总结）
前言前面我们的几篇文章介绍了一系列关于运算符的介绍,以及各个运算符的优化方式和技巧.其中涵盖:查看执行计划的方式.几种数据集常用的连接方式.联合运算符方式.并行运算符等一系列的我们常见的运算符.有兴 ...
（Hibernate进阶）Hibernate系列——总结篇(九)
这篇博文是hibernate系列的最后一篇,既然是最后一篇,我们就应该进行一下从头到尾,整体上的总结,将这个系列的内容融会贯通. 概念 Hibernate是一个对象关系映射框架,当然从分层的角度看,我 ...
SQL Server调优系列基础篇（联合运算符总结）
前言上两篇文章我们介绍了查看查询计划的方式,以及一些常用的连接运算符的优化技巧,本篇我们总结联合运算符的使用方式和优化技巧. 废话少说,直接进入本篇的主题. 技术准备基于SQL Server200 ...
return error怎么定义_SpringBoot 系列 web 篇之自定义返回 Http Code 的 n 种姿势
200105-SpringBoot 系列 web 篇之自定义返回 Http Code 的 n 种姿势虽然 http 的提供了一整套完整.定义明确的状态码,但实际的业务支持中,后端并不总会遵守这套规则 ...

SNN系列｜神经元模型篇(1) Hodgkin Huxley