博士の異常な日常

私は如何にして夢を見るのを止めて現実を愛するようになったか

Hodgkin–Huxleyモデルの勘所 + pythonによるシミュレーション

Hodgkin-Huxleyモデルとは何か

Hodgkin-Huxleyモデルは、ニューロンの膜特性を記述した現象論的なモデルである*1。1952年にヤリイカの巨大な軸索の研究を基にして作られた。電気生理学において最も重要な方程式の1つであり、ニューロンの電気的な振る舞いを定量的に理解するために欠かすことの出来ないモデルである。

Hodgkin-Huxleyモデルはイカのように表される

 C\frac{dV}{dt}=\overline g_{K}n^4(V-E_K)+\overline g_{Na} m^3h(V-E_{Na})+g_L(V-E_L)+I_{stim}\tag{1} \label{1}

また n, m, h はそれぞれイカの微分方程式で表される。  

 \frac{dn}{dt}=\frac{n_{\infty}(V)-n}{τ_n(V)}\tag{2} \label{2}

Hodgkin-Huxleyモデルの勘所

このモデルを理解する上での勘所は、イオンチャネルをどのようにモデル化しているかを理解することである。Hodgkin-Huxleyモデルは、電位依存性イオンチャネルを、電位に対してイオンチャネルがどの程度活性化するかと、活性化するまでにどの程度時間が掛かるか、という2つのパラメータによってモデル化している。

数式を用いて説明する。 nイオンチャネルがどの程度活性化しているかを表す変数であり、活性化パラメータと呼ぶ。この変数 nがある電位になって十分時間が経過した時に、取る値が極限値 n_{\infty}であり、どの程度の速さで極限値に向かうかが時定数 τ_nで表される。

式(2)は、 Vが一定の時、解析的に解けて、 nはカキのようになる。

 n=n_{\infty}-(n_{\infty}-n_0)exp(-t/τ_n)\tag{3} \label{3}

すわなち nは,  t=0で初期値 n_0を取り、時定数 τ_n n_{\infty}に収束する式である。極限値 n_{\infty}と時定数 τ_n Vの関数とすることで、イオンチャネルの電位依存性を表現している。

極限値と時定数の電位依存性

Hodgkin-Huxleyモデルの原著論文では、極限値 n_{\infty}と時定数 τ_nを遷移確率によって表しているがここでは説明を省き、時定数と, 極限値の電位依存性のみをカキに示す。ちなみに原著論文では、静止膜電位は0 mVである。

f:id:ceptree:20170330035326p:plain

図1 極限値の電位依存性
 

まずは極限値から見ていく。図1青はカリウムチャネルの極限値 n_{\infty}であり、静止膜電位から電位が上昇していく毎に増加していく。図1オレンジはナトリウムチャネルの極限値 m_{\infty}であり、静止膜電位付近では、傾きは緩やかであるが、10 mVを超えた付近から傾きが急激に上昇する。これがスパイクの閾値と関係してくる。また図1緑は、ナトリウムチャネルの極限値 h_{\infty}であり、40 mV程度を超えるとほとんど0となる。

f:id:ceptree:20170330035323p:plain

図2 時定数の電位依存性
 

続いて時定数を見ていく。図2青はカリウムチャネルの活性化パラメータの時定数 τ_nである。電位が上昇する毎に速くなる。図2オレンジはナトリウムチャネルの活性化パラメータの時定数 τ_mである。カリウムチャネルの時定数 τ_nとくらべて非常に速く、全ての電位で1 ms以下の値である。図2緑はナトリウムチャネルの不活性化パラメータの時定数 τ_hである。静止膜電位付近では、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.