Julia で Hodgkin-Huxley モデルのシミュレーションをしてみる

NOTE

本記事で使用しているコードは, 「はじめての神経回路シミュレーション」¹における C コードを Julia で書き換えたものです.

Hodgkin-Huxley モデル

Hodgkin-Huxley モデルとは, A. L. Hodgkin, A. F. Huxley によって定式化された²世界初の単一ニューロンの数理モデルであり, あらゆるニューロンモデルの基礎となっている.

HH モデルは次式で表される：

\begin{equation*} C\dfrac{dV}{dt}=-\bar{g}_{\text{leak}}(V(t)-E_{\text{leak}})-g_{\text{Na}}(V,t)(V(t)-E_{\text{Na}})-g_{\text{K}}(V,t)(V(t)-E_{\text{K}})+I_{\text{ext}}(t) \end{equation*}

ここで, $C$ は膜のキャパシタンス $[\mathrm{\mu F/cm^{2}}]$ , $t$ は時間 $[\mathrm{ms}]$ , $V(t)$ は膜電位 $[\mathrm{mV}]$ , $\bar{g}_{\text{leak}}$ は定数でリークコンダクタンス $[\mathrm{mS/cm^{2}}]$ , $E_{\text{leak}}$ は主に $\mathrm{Cl^{-}}$ イオンの反転電位 $[\mathrm{mV}]$ , $g_{\text{Na}}(V,t)$ は電位依存の $\mathrm{Na^{+}}$ チャネルコンダクタンス $[\mathrm{mS/cm^{2}}]$ , $E_{\text{Na}}$ は $\mathrm{Na^{+}}$ イオンの反転電位 $[\mathrm{mV}]$ , $g_{\text{K}}(V,t)$ は電位依存の $\mathrm{K^{+}}$ チャネルコンダクタンス $[\mathrm{mS/cm^{2}}]$ , $E_{\text{K}}$ は $\mathrm{K^{+}}$ イオンの反転電位 $[\mathrm{mV}]$ , $I_{\text{ext}}(t)$ は細胞の外部から注入する電流 $[\mathrm{\mu A/cm^{2}}]$ である. キャパシタの総容量やコンダクタンス値を決定するチャネルの個数は膜の面積に比例するため, これらの値は単位面積あたりに正規化されている.

$g_{\text{Na}}(V,t)$ , $g_{\text{K}}(V,t)$ は, それぞれ次式で表される：

\begin{align*} &g_{\text{Na}}(V,t)=\bar{g}_{\text{Na}}m^{3}(V,t)h(V,t)\\ &g_{\text{K}}(V,t)=\bar{g}_{\text{K}}n^{4}(V,t) \end{align*}

ここで, $\bar{g}_{\text{Na}}$ , $\bar{g}_{\text{K}}$ はそれぞれ定数で最大コンダクタンス $[\mathrm{mS/cm^{2}}]$ , $m(V,t)$ , $h(V,t)$ , $n(V,t)$ は膜電位 $V$ に依存するゲート変数である. ゲート変数はイオンチャネルの開口率を表し, 次式で更新される：

\begin{align*} &\dfrac{dm}{dt}=\alpha_{m}(V)(1-m(V,t))-\beta_{m}(V)m(V,t)\\ &\dfrac{dh}{dt}=\alpha_{h}(V)(1-h(V,t))-\beta_{h}(V)h(V,t)\\ &\dfrac{dn}{dt}=\alpha_{n}(V)(1-n(V,t))-\beta_{n}(V)n(V,t) \end{align*}

$\alpha_{x}(V)$ , $\beta_{x}(V)$ は, それぞれ次式で定義される：

\begin{align*} &\alpha_{m}(V)=\dfrac{2.5-0.1V}{\exp(2.5-0.1V)-1}\\ &\beta_{m}(V)=4\exp\left(-\dfrac{V}{18}\right)\\ &\alpha_{h}(V)=0.07\exp\left(-\dfrac{V}{20}\right)\\ &\beta_{h}(V)=\dfrac{1}{\exp(3-0.1V)+1}\\ &\alpha_{n}(V)=\dfrac{0.1-0.01V}{\exp(1-0.1V)-1}\\ &\beta_{n}(V)=0.125\exp\left(-\dfrac{V}{80}\right) \end{align*}

定数は, $C=1\, \mathrm{\mu F/cm^{2}}$ と正規化したとき, $\bar{g}_{\text{leak}}=0.3\, \mathrm{mS/cm^{2}}$ , $E_{\text{leak}}=10.6\, \mathrm{mV}$ , $g_{\text{Na}}=120\, \mathrm{mS/cm^{2}}$ , $E_{\text{Na}}=115\, \mathrm{mV}$ , $g_{\text{K}}=36\, \mathrm{mS/cm^{2}}$ , $E_{\text{K}}=-12\, \mathrm{mV}$ である.

ここで, ゲート変数について補足する. 先ほど示したゲート変数の式は,

\begin{equation*} \tau_{x}(V)\dfrac{dx}{dt}=-(x-x_{0}(V)) \end{equation*}

と変形でき,

\begin{align*} &\tau_{x}(V)=\dfrac{1.0}{(\alpha_{x}(V)+\beta_{x}(V))}\\ &x_{0}(V)=\dfrac{\alpha_{x}(V)}{(\alpha_{x}(V)+\beta_{x}(V))} \end{align*}

が得られる.

つまり, 定常状態で $x$ の値は $x_{0}(V)$ にあり, 時定数 $\tau_{x}(V)$ で変化する.

以上より, HH モデルは4変数 $(V,n,m,h)$ からなる微分方程式であり, これを数値的に解くことで, ニューロンの挙動, 特にスパイク発射の挙動を再現することが出来る.

HH モデルのシミュレーション

HH モデルの $1$ ニューロンのシミュレーションを Julia で実装する.

using Plots
using LaTeXStrings

plot_font = "Computer Modern"
default(fontfamily=plot_font, linewidth=1, label=nothing, grid=true)

const E_REST = -65.0           # mV
const C = 1.0                  # micro F / cm^2
const G_LEAK = 0.3             # mS / cm^2
const E_LEAK = 10.6 + E_REST   # mV
const G_NA = 120.0             # mS / cm^2
const E_NA = 115.0 + E_REST    # mV
const G_K = 36.0               # mS / cm^2
const E_K = -12.0 + E_REST     # mV

const DT = 0.01   # 10 micro s
# const T = 1000    # 1000 ms
# const NT = T / DT

@inline alpha_m(v) = (2.5 - 0.1 * (v - E_REST)) / (exp(2.5 - 0.1 * (v - E_REST)) - 1.0)
@inline beta_m(v) = 4.0 * exp(-(v - E_REST) / 18.0)
@inline alpha_h(v) = 0.07 * exp(-(v - E_REST) / 20.0)
@inline beta_h(v) = 1.0 / (exp(3.0 - 0.1 * (v - E_REST)) + 1.0)
@inline alpha_n(v) = (0.1 - 0.01 * (v - E_REST)) / (exp(1 - 0.1 * (v - E_REST)) - 1.0)
@inline beta_n(v) = 0.125 * exp(-(v - E_REST) / 80.0)

@inline m0(v) = alpha_m(v) / (alpha_m(v) + beta_m(v))
@inline h0(v) = alpha_h(v) / (alpha_h(v) + beta_h(v))
@inline n0(v) = alpha_n(v) / (alpha_n(v) + beta_n(v))

@inline tau_m(v) = 1.0 / (alpha_m(v) + beta_m(v))
@inline tau_h(v) = 1.0 / (alpha_h(v) + beta_h(v))
@inline tau_n(v) = 1.0 / (alpha_n(v) + beta_n(v))

@inline dmdt(v, m) = (1.0 / tau_m(v)) * (-m + m0(v))
@inline dhdt(v, h) = (1.0 / tau_h(v)) * (-h + h0(v))
@inline dndt(v, n) = (1.0 / tau_n(v)) * (-n + n0(v))
@inline dvdt(v, m, h, n, i_ext) = (-G_LEAK * (v - E_LEAK) - G_NA * m^3 * h * (v - E_NA) - G_K * n^4 * (v - E_K) + i_ext) / C

global t_val = Float64[]
global v_val = Float64[]
global m_val = Float64[]
global h_val = Float64[]
global n_val = Float64[]

function hh!(x, y, T)
    global t_val = Float64[]
    global v_val = Float64[]
    global m_val = Float64[]
    global h_val = Float64[]
    global n_val = Float64[]
    v = E_REST
    m = m0(v)
    h = h0(v)
    n = n0(v)

    i_ext = 9.0 # micro A / cm^2

    for nt in 1:T/DT
        t = DT * nt
        push!(t_val, t)
        push!(v_val, v)
        push!(m_val, m)
        push!(h_val, h)
        push!(n_val, n)

        dmdt1 = dmdt(v, m)
        dhdt1 = dhdt(v, h)
        dndt1 = dndt(v, n)
        dvdt1 = dvdt(v, m, h, n, i_ext)

        dmdt2 = dmdt(v + 0.5 * DT * dvdt1, m + 0.5 * DT * dmdt1)
        dhdt2 = dhdt(v + 0.5 * DT * dvdt1, h + 0.5 * DT * dhdt1)
        dndt2 = dndt(v + 0.5 * DT * dvdt1, n + 0.5 * DT * dndt1)
        dvdt2 = dvdt(v + 0.5 * DT * dvdt1, m + 0.5 * DT * dmdt1, h + 0.5 * DT * dhdt1, n + 0.5 * DT * dndt1, i_ext)

        dmdt3 = dmdt(v + 0.5 * DT * dvdt2, m + 0.5 * DT * dmdt2)
        dhdt3 = dhdt(v + 0.5 * DT * dvdt2, h + 0.5 * DT * dhdt2)
        dndt3 = dndt(v + 0.5 * DT * dvdt2, n + 0.5 * DT * dndt2)
        dvdt3 = dvdt(v + 0.5 * DT * dvdt2, m + 0.5 * DT * dmdt2, h + 0.5 * DT * dhdt2, n + 0.5 * DT * dndt2, i_ext)

        dmdt4 = dmdt(v + DT * dvdt3, m + DT * dmdt3)
        dhdt4 = dhdt(v + DT * dvdt3, h + DT * dhdt3)
        dndt4 = dndt(v + DT * dvdt3, n + DT * dndt3)
        dvdt4 = dvdt(v + DT * dvdt3, m + DT * dmdt3, h + DT * dhdt3, n + DT * dndt3, i_ext)

        m += DT * (dmdt1 + 2 * dmdt2 + 2 * dmdt3 + dmdt4) / 6.0
        h += DT * (dhdt1 + 2 * dhdt2 + 2 * dhdt3 + dhdt4) / 6.0
        n += DT * (dndt1 + 2 * dndt2 + 2 * dndt3 + dndt4) / 6.0
        v += DT * (dvdt1 + 2 * dvdt2 + 2 * dvdt3 + dvdt4) / 6.0
    end
    plot(getfield(Main, Symbol(string(x) * "_val")), getfield(Main, Symbol(string(y) * "_val")), dpi=300)
end

実際に, シミュレーションを可視化してみる.

hh!(:t, :v, 1000)
xlabel!("Time [ms]")
ylabel!("Membrane Potential [mV]")

上記の実行により, 以下が得られる.

図1: $1000\, \mathrm{ms}$ 間の計算結果

図1では, スパイク発射の様子が分かりにくいため, $T=100$ としてシミュレーションする.

hh!(:t, :v, 100)
xlabel!("Time [ms]")
ylabel!("Membrane Potential [mV]")

上記の実行により, 以下が得られる.

図2: $100\, \mathrm{ms}$ 間の計算結果

スパイク発射のメカニズムを理解するために, 各変数の振る舞いを並列表示する.

p1 = hh!(:t, :v, 15)
ylabel!(L"V\quad [\mathrm{mV}]")
p2 = hh!(:t, :m, 15)
ylabel!(L"m")
p3 = hh!(:t, :h, 15)
ylabel!(L"h")
p4 = hh!(:t, :n, 15)
ylabel!(L"n")
xlabel!(L"t\quad [\mathrm{ms}]")

plot(p1, p2, p3, p4, layout=(4, 1))

上記の実行により, 以下が得られる.

図3: スパイク発射のメカニズム

定常状態において, ごく短い時間, $1\, \mathrm{ms}$ だけ電流を加えると, $1$ 発だけスパイクが発射される (図3). 電流を加えることによって膜電位が脱分極すると, それに伴って $\mathrm{Na^{+}}$ チャネルのゲート変数 $m$ の値が大きくなる. すると $\mathrm{Na^{+}}$ イオンの電流が発生し, 膜電位はさらに脱分極することになる. このポジティブフィードバックにより, 膜電位が急に $100\, \mathrm{mV}$ 程度上昇する.

膜電位の値が十分大きくなると, 今度はゲート変数 $h$ の値が $0$ になるため, 発生した $\mathrm{Na^{+}}$ イオンの電流は停止し, 膜電位は下降する. これにより膜電位が急速上昇・下降し, スパイクが発射される. それと並行して, $h$ の変化と同程度の時間スケールで $\mathrm{K^{+}}$ イオンの電流が発生し, 膜電位を静止電位よりも過分極させる. この期間は不応期とよばれ, $\mathrm{Na^{+}}$ チャネルが不活性状態にあるため, 次のスパイクを発射できない. 最終的に $m$ , $h$ , $n$ の値は定常状態に戻り, 次にスパイクを発射する準備が整う.

このように, スパイクは $\mathrm{Na^{+}}$ チャネルと $\mathrm{K^{+}}$ チャネルのコンダクタンスの非常に素早いダイナミクスによって生成される. また, それ以外の膜電位のダイナミクスは, 時定数 $\tau=RC=C/\bar{g}_{\text{leak}}$ に従う相対的に遅いダイナミクスである. ここで, $R$ は膜抵抗であり, $\bar{g}_{\text{leak}}$ の逆数である.

発火頻度

先ほどのシミュレーションでは外部電流として $I_{\text{ext}}=9.0\, \mathrm{\mu A/cm^{2}}$ と設定していたが, $I_{\text{ext}}$ の値を変化させるとスパイク発射はどのように変化するのだろうか.

外部電流によるスパイク発射の変化をシミュレーションするために, hh! の第4引数で i_ext を設定できるようにし, 以下を実行する.

@gif for i in 5.0:0.01:10.0
    hh!(:t, :v, 1000, i)
    title!(L"I_{\mathrm{ext}} = " * string(i))
    xlabel!("Time [ms]")
    ylabel!("Membrane Potential [mV]")
end

上記の実行により, 以下が得られる.

NOTE

GIF ファイルのサイズが大きいため, 読み込みに時間がかかります. 参考までに, 私の環境では30秒弱かかりました. 気長にお待ちください.

図4: 外部電流の変化に伴うスパイク発射シミュレーション

図4から分かるように, $I_{\text{ext}}=6.2\, \mathrm{\mu A/cm^{2}}$ のときはスパイクを発射しないが, $I_{\text{ext}}=6.3\, \mathrm{\mu A/cm^{2}}$ のときにはスパイクを発射する. このように, 外部電流がある閾値を超えると突然高頻度でスパイクを発射するニューロンを Type II ニューロンとよぶ.

一方, 大脳皮質のニューロンのように, 入力される外部電流の強度に比例して連続的に発火頻度が上昇するニューロンを Type I ニューロンとよぶ.

Footnotes

山﨑匡, 五十嵐潤, "はじめての神経回路シミュレーション: 1ニューロンからヒト全脳モデルまで", 講談社, 2021, ISBN: 978-4-627-85631-8. ↩
A. L. Hodgkin, A. F. Huxley, "A quantitative description of membrane current and its application to conduction and excitation in nerve," The Journal of Physiology, First published: 28 August 1952, https://doi.org/10.1113/jphysiol.1952.sp004764. ↩