PythonでHodgkin–Huxleyモデルを実装してみた
はじめに
ふだんはパッチクランプと行動実験をひたすらラボにこもってやっている筆者ですが、ふと思い立って違うことをやってみた記録です。
ホジキン・ハクスレー（Hodgkin–Huxley）モデルは、神経細胞の「活動電位（action potential）」を定量的に説明した超有名なモデルです。
1952年、英国ケンブリッジ大学のA. L. HodgkinとA. F. Huxleyはヤリイカの巨大軸索（直径1mmもある巨大な神経線維です。めちゃくちゃ太い。）を用いて生体電気信号を電極で計測する実験を行いました。内容としては膜電位を一定の条件に固定（ボルテージクランプ法）することで細胞膜のイオンチャネルを通過するイオン電流を記録し、それを定量的に数式でモデル化するというものです。すごい。
この記事では、そのHodgkin–HuxleyモデルをPythonで実装し、
神経細胞の発火を数値的にシミュレーションしてみたいと思います。
今回も例によってchatGPT大先生の多大なお力添えをいただいております。
Hodgkin–Huxleyモデル？
🔌 生理学的背景
神経細胞（ニューロン）は、電気的な信号によって情報を伝えます。
その信号の単位が「活動電位（スパイク）」です。
通常、ニューロンの膜電位は約−65 mVの「静止状態」にありますが、
一定以上の刺激（閾値）を受けると、膜が急激に脱分極して＋電位に転じます。
この瞬間的な電位変化が「発火」です。
背後では次のようなイオンの流れが起こっています：
- Na⁺チャネルが開く → Na⁺が流入し脱分極
- わずかに遅れてK⁺チャネルが開く → K⁺が流出し再分極
- チャネルが閉じて膜電位が元に戻る
これらの時定数の異なるイオン電流が、スパイク波形を形成します。
数理モデルとしての構造
Hodgkin–Huxleyモデルでは、膜電位の時間変化を次の微分方程式で記述します：
ここで：
| 記号 | 意味 |
|---|---|
| 膜電位 (mV) | |
| 膜容量 | |
| 外部からの注入電流 | |
| 各イオンチャネルの最大コンダクタンス | |
| 各イオンの平衡電位 | |
| チャネルの開閉確率（ゲーティング変数） |
Na⁺チャネルは**m（開く確率）とh（閉じる確率）の2つのゲート、
K⁺チャネルはn（開く確率）**の1つのゲートで制御されます。
これらのゲートは膜電位依存的に変化(Voltege gated)し、複雑なイオン流を生み出します。
直感的な理解：電気回路モデル
Hodgkin–Huxleyモデルは、細胞膜を電気回路として捉えることができます。
- 細胞膜自身はコンデンサー（電荷を蓄える）
- イオンチャネルは可変抵抗（電位依存的に変化）
- Na⁺やK⁺の電流は電流源
この電気回路を微分方程式として表すことで、神経活動を再現できます。
Pythonによる実装
以下のコードでは、scipy.integrate.solve_ivp を用いてHodgkin–Huxley方程式を数値的に解きます。
chatGPT先生ありがとう。僕はclaude codeに浮気しません。
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
def hh_model(t, y):
V, m, h, n = y # V in mV
Cm = 1.0 # μF/cm^2
gNa, gK, gL = 120.0, 36.0, 0.3 # mS/cm^2
ENa, EK, EL = 50.0, -77.0, -54.4 # mV
# step current: 10 μA/cm^2 from 10 to 40 ms
Iext = 10.0 if 10.0 <= t <= 40.0 else 0.0
# Shift membrane potential so that -65 mV -> 0 mV
v = V + 65.0
# Rate functions with numerical guards (use expm1 for stability)
def a_m(v):
denom = np.expm1((25.0 - v)/10.0) # exp(x)-1
return 0.1*(25.0 - v)/denom if denom != 0 else 1.0 # limit as v->25 is 1.0
def b_m(v): return 4.0*np.exp(-v/18.0)
def a_h(v): return 0.07*np.exp(-v/20.0)
def b_h(v): return 1.0/(np.exp((30.0 - v)/10.0) + 1.0)
def a_n(v):
denom = np.expm1((10.0 - v)/10.0)
return 0.01*(10.0 - v)/denom if denom != 0 else 0.1 # limit as v->10 is 0.1
def b_n(v): return 0.125*np.exp(-v/80.0)
am, bm = a_m(v), b_m(v)
ah, bh = a_h(v), b_h(v)
an, bn = a_n(v), b_n(v)
dmdt = am*(1.0 - m) - bm*m
dhdt = ah*(1.0 - h) - bh*h
dndt = an*(1.0 - n) - bn*n
INa = gNa * (m**3) * h * (V - ENa)
IK = gK * (n**4) * (V - EK)
IL = gL * (V - EL)
dVdt = (Iext - INa - IK - IL) / Cm
return [dVdt, dmdt, dhdt, dndt]
# Initial conditions (classic HH steady ~ -65 mV)
y0 = [-65.0, 0.0529, 0.596, 0.317] # m, h, n near steady at -65 mV
t_span = (0.0, 50.0)
t_eval = np.linspace(*t_span, 2000)
sol = solve_ivp(hh_model, t_span, y0, t_eval=t_eval, max_step=0.05, rtol=1e-6, atol=1e-9)
plt.plot(sol.t, sol.y[0])
plt.xlabel("Time (ms)")
plt.ylabel("Membrane potential (mV)")
plt.title("Hodgkin–Huxley with step current (10–40 ms)")
plt.tight_layout()
plt.show()
まとめ
スパイクがなんとかかんとか表示されましたね。
なかなか発火せずコードを調整して色々と格闘してとりあえずは出力できました。
Pythonはやはりライブラリが強力で色々できるので便利です。例えばpyABFでABFファイルの電流データの解析もできるし、GUIでデータをエクスプローラから直で引っ張ってこれるような解析パイプライン構築も可能ですね。
以下参考文献
-
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.
https://doi.org/10.1113/jphysiol.1952.sp004764 -
杉 晴夫 (2015).
『神経とシナプスの科学 現代脳研究の源流』
講談社ブルーバックス. -
山本 拓都.
Juliaで学ぶ計算論的神経科学
https://compneuro-julia.github.io/preface.html
Discussion