神经元内信号传递的计算模型-HH模型

神经元动力学可以被设想为一个总和过程(有时也称为“集成”过程)，并结合一种触发动作电位高于临界电压的机制。

在这里主要介绍Hodgkin Huxley模型模拟离子通道，膜电势的改变。

原理

基础知识

细胞膜两侧离子浓度差产生了电势差，该电势差被称之为电位，而膜电位的改变与 $Na^+$ 、 $K^+$ 、leaky(以 $cl^-$ 为代表）的三种通道的流入流出息息相关， $Na^+$ 通道有激活和失活两个门控开关， $K^+$ 通道只有激活门；每个离子通道分别存在允许和不允许两种状态，而状态的转变与膜电势有关，可以表示每一时刻每一门控开关打开的概率，换句话说也就是可以衡量各离子通道允许离子通过的能力，进而记录膜电势的变化，达到对于神经元内信号传递过程进行仿真的目的

膜作为电路的生物物理学

将上文的基础知识进行拆分可以将神经元内信号传递过程归纳为如下几步：
1.电刺激达到阈值，膜电势出现改变，各离子通道允许与不允许状态的概率发生改变：

每个离子通道都有两种状态，允许和不允许状态（permissive or non-permissive）。
用 $p_i(t)$ 表示某离子通道某时刻的概率，值的大小从0到1变化，由不允许到允许转换的速率为 $\alpha_i(t)$ ,由允许到不允许状态转换的速度为 $\beta_i(t)$

2.离子通道的状态改变就对应着相应的门控开关开放概率的改变：

引入门控变量m ，n，h（ $Na^+$ 通道有激活和失活两个门控开关， $K^+$ 通道只有激活门）来对给定时间通道开启的概率进行建模，这里m和n共同负责调控Na离子通道，h负责调控K离子通道。
允许与不允许状态转换速率是与膜电势有关的，不同门控开关对应不同参数：
根据上图六个转换速率的参数可以进一步推导出开放通道的比例如何随时间变化：

以及各通道打开的概率（以n为例）：

3.各离子通道允许离子通过的能力可以用电导 $g$ 来进行衡量：

离子流过的每一个状态都对应一个电阻，因为细胞膜是有对应的离子泵，其通过相应离子的能力可以改变，因此是一个可变电导
$g_{Na}=m^3\times g_{Na\_max}\times h$
$g_{K}=n^4\times g_{K\_max}$

4.由 $g = 1 / R$ , $I = U / R$ 我们可以计算出流经各离子通道的电流：

$I_{Na}=g_{Na}\times (u-E{_{Na}})$
$I_{K}=g_{K}\times (u-E{_{K}})$
$I_{Leaky}=g_{L}\times (u-E{_{L}})$

5.根据电学公式 $C=q/U,I_c=C\frac{du}{dt}$ 以及并联电流公式我们可以计算出此时的膜电势：
在这里插入图片描述
式中 $u$ 代表膜电势， $\sum\limits_{k}I_k(t)$ 代表通过细胞膜的离子电流的和， $I (t)$ 则表示输入的刺激电流

整个过程总结来说也就是 $\alpha$ 、 $\beta$ — $m$ 、 $h$ 、 $n$ — $g$ — $I_C$ — $E_C$ 这一系列变量的循环更新。

python实现

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

模拟离子通道

具体步骤如下：

1.设置参数，此处设置膜电势为-80至80mV这个区间
2.计算在不同膜电势情况下离子通道门控开关m，h，n对应的状态转换速率参数： $\alpha_m,\alpha_h,\alpha_n,\beta_m,\beta_h,\beta_n$
3.计算各通道打开概率 $p=\alpha /(\alpha+\beta)$

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
##1.设置参数
# 离子通道被建模为与电池串联的电阻器。电池表示特定离子的平衡电势，电阻器反映通道对特定离子的渗透性
V = np.arange(-80.0,80,0.01) # [mV]
##2.计算状态速率转换参数
#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))##3.计算概率
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)')

在这里插入图片描述

模拟膜电位

具体步骤如下：

1.设置参数，包括平衡电势、电容、电阻、静息电位、最大电导、dt、以及仿真时长
2.循环开始
更新 $\alpha,\beta$
若i为1，利用上述参数得到初始化 $m, n, h$
根据 $m, n, h$ 更新电导 $g$
计算各通道电流 $I$ 得到 $-\sum\limits_{k}I_k(t)+I(t)$
更新膜电势 $u$
根据膜电势 $u$ 更新 $m, n, h$

# 参数定义 可参考
# 平衡电位
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，m，n的值
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'])

在这里插入图片描述
下面的代码为持续恒定电流刺激的膜电位变化：

'''
=======================================
% 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.
#=======================================
'''
def run_model(I_inj, dt):num_spike = 0# Initialize time vector and time stept = np.arange(0,dt,100) #Time (ms)# Initialize m, n, and h, which represent K+ channel activation, Na+# channel activation, and Na+ channel inactivation, respectively.m = np.zeros(len(t))n = np.zeros(len(t))h = np.zeros(len(t))# Initialize alphas and betas, the rate constantsalpha_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))# Initialize membrane voltage and conductance vectorsV_m = np.zeros(len(t))V_m[0] = 60g_K = np.zeros(len(t))g_Na = np.zeros(len(t))# Constantsg_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 in range(0,len(t)):#Calculate all alpha and beta valuesalpha_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)# Initial conditionsif i == 0: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])# 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])*dtn[i+1] = n[i] + (alpha_n[i]*(1-n[i])-beta_n[i]*n[i])*dth[i+1] = h[i] + (alpha_h[i]*(1-h[i])-beta_h[i]*h[i])*dt#check spikeif 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 + 1V_m = V_m + V_restreturn V_m,t