Hodgkin–Huxleyモデルの勘所 + pythonによるシミュレーション
Hodgkin-Huxleyモデルとは何か
Hodgkin-Huxleyモデルは、ニューロンの膜特性を記述した現象論的なモデルである*1。1952年にヤリイカの巨大な軸索の研究を基にして作られた。電気生理学において最も重要な方程式の1つであり、ニューロンの電気的な振る舞いを定量的に理解するために欠かすことの出来ないモデルである。
Hodgkin-Huxleyモデルはイカのように表される
またはそれぞれイカの微分方程式で表される。
Hodgkin-Huxleyモデルの勘所
このモデルを理解する上での勘所は、イオンチャネルをどのようにモデル化しているかを理解することである。Hodgkin-Huxleyモデルは、電位依存性イオンチャネルを、電位に対してイオンチャネルがどの程度活性化するかと、活性化するまでにどの程度時間が掛かるか、という2つのパラメータによってモデル化している。
数式を用いて説明する。はイオンチャネルがどの程度活性化しているかを表す変数であり、活性化パラメータと呼ぶ。この変数がある電位になって十分時間が経過した時に、取る値が極限値であり、どの程度の速さで極限値に向かうかが時定数で表される。
式(2)は、が一定の時、解析的に解けて、はカキのようになる。
すわなちは, で初期値を取り、時定数でに収束する式である。極限値と時定数はの関数とすることで、イオンチャネルの電位依存性を表現している。
極限値と時定数の電位依存性
Hodgkin-Huxleyモデルの原著論文では、極限値と時定数を遷移確率によって表しているがここでは説明を省き、時定数と, 極限値の電位依存性のみをカキに示す。ちなみに原著論文では、静止膜電位は0 mVである。
まずは極限値から見ていく。図1青はカリウムチャネルの極限値であり、静止膜電位から電位が上昇していく毎に増加していく。図1オレンジはナトリウムチャネルの極限値であり、静止膜電位付近では、傾きは緩やかであるが、10 mVを超えた付近から傾きが急激に上昇する。これがスパイクの閾値と関係してくる。また図1緑は、ナトリウムチャネルの極限値であり、40 mV程度を超えるとほとんど0となる。
続いて時定数を見ていく。図2青はカリウムチャネルの活性化パラメータの時定数である。電位が上昇する毎に速くなる。図2オレンジはナトリウムチャネルの活性化パラメータの時定数である。カリウムチャネルの時定数とくらべて非常に速く、全ての電位で1 ms以下の値である。図2緑はナトリウムチャネルの不活性化パラメータの時定数である。静止膜電位付近では、8 ms程度とかなり遅いが、電位が上がる毎に急激に早くなり、40 mVを超えると1 ms近くなる。
長くなったので今回はここまで。
最後に今回シミュレーションに使ったコードを貼り付けておく。使用の際は自己責任でお願いします。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import time class HH(): def __init__(self): # 膜電位 Vm [mV] self.V = np.arange(-12,120,0.1) # 膜容量 Cm [uF/cm2] self.Cm = 1. # 平衡電位 E [mV] self.E_Na = 120. self.E_K = -12. self.E_L = 10.6 # 最大コンダクタンス gbar [mS/cm2] self.gbar_Na = 120. self.gbar_K = 36. self.g_L = 0.3 def Sim_Init(self,dt=0.01,TotalTime=7): # サンプリング時間 [ms] self.dt = dt # シミュレーション時間 [ms] self.TotalTime = TotalTime # 時間 [ms] self.t = np.arange(0,TotalTime,dt) # サンプリング数 n = self.t.size self.n = n # ------------------------------------------ # Initialize # ------------------------------------------ self.I = np.zeros(n) self.Vm = np.zeros(n) self.Na_I = np.zeros(n) self.Na_g = np.zeros(n) self.Na_m = np.zeros(n) self.Na_h = np.zeros(n) self.Na_m_inf = np.zeros(n) self.Na_h_inf = np.zeros(n) self.Na_m_tau = np.zeros(n) self.Na_h_tau = np.zeros(n) self.Na_m_alpha= np.zeros(n) self.Na_h_alpha= np.zeros(n) self.Na_m_beta = np.zeros(n) self.Na_h_beta = np.zeros(n) self.Na_DF = np.zeros(n) self.K_I = np.zeros(n) self.K_g = np.zeros(n) self.K_n = np.zeros(n) self.K_n_inf = np.zeros(n) self.K_n_tau = np.zeros(n) self.K_n_alpha = np.zeros(n) self.K_n_beta = np.zeros(n) self.K_DF = np.zeros(n) self.L_I = np.zeros(n) self.L_DF = np.zeros(n) # 初期値 self.Na_m[0] = self.Calc_m_inf(0) self.Na_h[0] = self.Calc_h_inf(0) self.K_n[0] = self.Calc_n_inf(0) def Sim_Stim_Step(self,Stim_del=1.,Stim_dur=5.,Stim_amp=25): self.Stim_del = Stim_del # 刺激開始 [ms] self.Stim_dur = Stim_dur # 刺激幅 [ms] self.Stim_amp = Stim_amp # 刺激振幅 [uA/cm2] self.I[ int(round(self.Stim_del/self.dt)) : int(round((self.Stim_del+self.Stim_dur)/self.dt)) ] = self.Stim_amp def Sim_Run(self): st = time.time() for i in range(self.n-1): # n self.K_n_alpha[i] = self.Calc_n_alpha(self.Vm[i]) self.K_n_beta[i] = self.Calc_n_beta(self.Vm[i]) # m self.Na_m_alpha[i] = self.Calc_m_alpha(self.Vm[i]) self.Na_m_beta[i] = self.Calc_m_beta(self.Vm[i]) # h self.Na_h_alpha[i] = self.Calc_h_alpha(self.Vm[i]) self.Na_h_beta[i] = self.Calc_h_beta(self.Vm[i]) # 定常状態 n∞ self.Na_m_inf[i] = self.Calc_inf(self.Na_m_alpha[i], self.Na_m_beta[i]) self.Na_h_inf[i] = self.Calc_inf(self.Na_h_alpha[i], self.Na_h_beta[i]) self.K_n_inf[i] = self.Calc_inf(self.K_n_alpha[i], self.K_n_beta[i]) # 時定数 τ self.Na_m_tau[i] = self.Calc_tau(self.Na_m_alpha[i],self.Na_m_beta[i]) self.Na_h_tau[i] = self.Calc_tau(self.Na_h_alpha[i],self.Na_h_beta[i]) self.K_n_tau[i] = self.Calc_tau(self.K_n_alpha[i], self.K_n_beta[i]) self.Na_m[i+1] = self.Na_m[i] + self.dt*( self.Na_m_inf[i]-self.Na_m[i] ) / self.Na_m_tau[i] self.Na_h[i+1] = self.Na_h[i] + self.dt*( self.Na_h_inf[i]-self.Na_h[i] ) / self.Na_h_tau[i] self.K_n[i+1] = self.K_n[i] + self.dt*( self.K_n_inf[i]- self.K_n[i] ) / self.K_n_tau[i] # コンダクタンス [mS/cm2] self.Na_g[i+1] = self.Calc_g_Na(self.Na_m[i+1],self.Na_h[i+1]) self.K_g[i+1] = self.Calc_g_K(self.K_n[i+1]) # Driving Force [mV] self.Na_DF[i+1] = self.Calc_DF(self.Vm[i],self.E_Na) self.K_DF[i+1] = self.Calc_DF(self.Vm[i],self.E_K) self.L_DF[i+1] = self.Calc_DF(self.Vm[i],self.E_L) # イオン電流 [uA/cm2] self.Na_I[i+1] = self.Calc_I(self.Vm[i],self.E_Na,self.Na_g[i+1]) self.K_I[i+1] = self.Calc_I(self.Vm[i],self.E_K, self.K_g[i+1]) self.L_I[i+1] = self.Calc_I(self.Vm[i],self.E_L, self.g_L) #膜電位 Vm [mV] self.Vm[i+1] = self.Vm[i] +self.dt/self.Cm*(-self.Na_I[i+1] - self.K_I[i+1] - self.L_I[i+1] + self.I[i+1]) elapsed_time = time.time() - st print ("Elapsed Time : %0.3f" % elapsed_time) def Plot_Sim(self): N = 7 NN= 1 fig=plt.figure() plt.subplots_adjust(wspace=0, hspace=0.1) ax=plt.subplot(N,1,NN) plt.plot(self.t,self.I) plt.ylabel('I [uA/cm2]') plt.grid(ls=':') NN +=1 plt.subplot(N,1,NN, sharex=ax) plt.plot(self.t,self.Vm) plt.ylabel('Vm [mV]') plt.grid(ls=':') NN +=1 plt.subplot(N,1,NN, sharex=ax) ax3=plt.plot(self.t,self.K_n_inf ,':') ax1=plt.plot(self.t,self.Na_m_inf,':') ax2=plt.plot(self.t,self.Na_h_inf,':') plt.plot(self.t,self.K_n ,ax3[0].get_color(),label='n') plt.plot(self.t,self.Na_m,ax1[0].get_color(),label='m') plt.plot(self.t,self.Na_h,ax2[0].get_color(),label='h') plt.legend(loc='lower right') plt.grid(ls=':') plt.ylim([0,1]) NN +=1 plt.subplot(N,1,NN, sharex=ax) plt.plot(self.t,self.K_n_tau,label='n') plt.plot(self.t,self.Na_m_tau,label='m') plt.plot(self.t,self.Na_h_tau,label='h') plt.ylabel('τ [ms]') plt.legend(loc='lower right') plt.grid(ls=':') NN +=1 plt.subplot(N,1,NN, sharex=ax) plt.plot(self.t,self.K_g,label='K') plt.plot(self.t,self.Na_g,label='Na') plt.ylabel('g [mS/cm2]') plt.legend(loc='lower right') plt.grid(ls=':') NN +=1 plt.subplot(N,1,NN, sharex=ax) plt.plot(self.t,self.K_DF,label='K') plt.plot(self.t,self.Na_DF,label='Na') plt.ylabel('V [mV]') plt.legend(loc='lower right') plt.grid(ls=':') NN +=1 plt.subplot(N,1,NN, sharex=ax) plt.plot(self.t,self.K_I,label='K') plt.plot(self.t,self.Na_I,label='Na') plt.plot(self.t,self.L_I,label='L') plt.plot(self.t,self.Na_I+self.K_I+self.L_I,label='Ion') plt.ylabel('I [uA/cm2]') plt.legend(loc='lower right') plt.grid(ls=':') plt.xlim([0,self.TotalTime]) plt.xlabel('Time [ms]') axs=fig.axes # for ax in axs[:-2]: # ax.get_xaxis().set_ticklabels([]) for ax in axs: ax.set_ylim(ax.get_ylim()) return fig def Calc_n_alpha(self,Vm): n_alpha = 0.01 * (10.-Vm) / ( np.exp( (10.-Vm)/10. ) - 1. ) return n_alpha def Calc_n_beta(self,Vm): n_beta = 0.125 * np.exp( - Vm/80. ) return n_beta def Calc_m_alpha(self,Vm): m_alpha = 0.1 * (25.-Vm) / ( np.exp( (25.-Vm)/10. ) - 1. ) return m_alpha def Calc_m_beta(self,Vm): m_beta = 4. * np.exp( - Vm/18. ) return m_beta def Calc_h_alpha(self,Vm): h_alpha = 0.07 * np.exp( - Vm/20. ) return h_alpha def Calc_h_beta(self,Vm): h_beta = 1. /( np.exp( (30.-Vm)/10. ) + 1. ) return h_beta def Calc_inf(self,alpha,beta): inf = alpha / (alpha + beta) return inf def Calc_n_inf(self,Vm): n_inf = self.Calc_inf(self.Calc_n_alpha(Vm),self.Calc_n_beta(Vm)) return n_inf def Calc_m_inf(self,Vm): m_inf = self.Calc_inf(self.Calc_m_alpha(Vm),self.Calc_m_beta(Vm)) return m_inf def Calc_h_inf(self,Vm): h_inf = self.Calc_inf(self.Calc_h_alpha(Vm),self.Calc_h_beta(Vm)) return h_inf def Calc_tau(self,alpha,beta): tau = 1 / (alpha + beta) return tau def Calc_g_Na(self,m,h): g_Na = self.gbar_Na * m ** 3 * h return g_Na def Calc_g_K(self,n): g_K = self.gbar_K * n ** 4 return g_K def Calc_DF(self,Vm,E): DF = Vm - E return DF def Calc_I(self,Vm,E,g): I = g * (Vm - E) return I def Calc_Vm_Charac(self): # n self.n_alpha = self.Calc_n_alpha(self.V) self.n_beta = self.Calc_n_beta(self.V) # m self.m_alpha = self.Calc_m_alpha(self.V) self.m_beta = self.Calc_m_beta(self.V) # h self.h_alpha = self.Calc_h_alpha(self.V) self.h_beta = self.Calc_h_beta(self.V) # 定常状態 n∞ self.m_inf = self.Calc_inf(self.m_alpha, self.m_beta) self.h_inf = self.Calc_inf(self.h_alpha, self.h_beta) self.n_inf = self.Calc_inf(self.n_alpha,self.n_beta) # 時定数 τ self.m_tau = self.Calc_tau(self.m_alpha, self.m_beta) self.h_tau = self.Calc_tau(self.h_alpha, self.h_beta) self.n_tau = self.Calc_tau(self.n_alpha, self.n_beta) def Plot_inf_Charac(self): fig = plt.figure() plt.plot(self.V,self.n_inf,label='n') plt.plot(self.V,self.m_inf,label='m') plt.plot(self.V,self.h_inf,label='h') plt.xlabel('Vm [mV]') plt.ylabel('$n_∞$') plt.xlim([self.V[0],self.V[-1]]) plt.legend(loc='lower right') plt.ylim([0,1]) return fig def Plot_tau_Charac(self): fig=plt.figure() plt.plot(self.V,self.n_tau,label='n') plt.plot(self.V,self.m_tau,label='m') plt.plot(self.V,self.h_tau,label='h') plt.xlabel('Vm [mV]') plt.ylabel('τ [ms]') plt.xlim([self.V[0],self.V[-1]]) plt.legend(loc='lower right') ax = fig.axes ax[0].set_ylim(ax[0].get_ylim()) return fig if __name__ == '__main__': hh= HH() hh.Calc_Vm_Charac() hh.Plot_inf_Charac() hh.Plot_tau_Charac() hh.Sim_Init() hh.Sim_Stim_Step() hh.Sim_Run() hh.Plot_Sim() plt.show()
*1:
Hodgkin, A. L., & Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117(4), 500–544.