はじめに
近年，機械学習分野ではニューラルネットワーク(ANN)がよく利用されます．
通常，ニューラルネットワークでは，脳神経細胞の簡単な数理モデル「形式ニューロンモデル」というものが用いられていますが，この形式ニューロンモデルは簡約化しすぎており，実際の神経細胞とは大きくかけ離れています．
何が違うのかというと，形式ニューロンモデルでは，実数値を情報として扱うのに対し，本物の神経細胞は{0,1}のデジタル信号(スパイク)を情報として扱います．
この本来の特徴を精緻に再現した数理モデルを，スパイキングニューロンモデル，スパイキングニューラルネットワーク(SNN)と言います．
- Spiking Neural Networkについて
ここではSNNについて深く説明はしません！
以下の記事に詳細はまとめているので，そちらを一読の上で，この記事を読むとより深く理解できるかと思います．
【参考：Spiking Neural Networkとは何なのか】
Hodgkin-Huxleyニューロンモデル
さて，本題です．
スパイキングニューロンモデルにはいくつか種類があります．
違いは，どれだけ本物に近づけるか，です．
スパイキングニューラルネットワーク分野では，計算量が少ないLIF(Lealy integrate-and-fire)モデルがよく使われます．(LIFモデルについてもいつか紹介記事書きますね)
神経解析学の分野では，表題にもあるHodgkin-Huxleyニューロンモデル (HHモデル)がよく使われるようです．(読み方はホジキン・ハクスレーモデル)
今回はこのHHモデルの紹介と，Pythonでの実装を行なっていきます．
- 注意点
そういえば機械学習タグをつけましたが，実際機械学習には適用するのは超非現実的です．
(そのタグで来た方ごめんなさい)
なので，「へえ，そんな分野もあるんだな」程度に読んでもらえれば幸いです．
また私も勉強中の身ですので，解釈の間違い等あるかもしれません，ご了承ください．
- そもそもHHモデルとは何か
(コードを早く見せい！って方はこちら)
HHモデルは，イギリスの電気生理学者のホジキンさんとハクスレーさんが，イカ(🦑)の神経細胞を観察・測定して神経細胞(正確には興奮性のもの)の電気的特性を方程式に落とし込んだものを指します．
神経細胞の内外では，多く存在するイオンの種類の違いにより電位差が生じており，それを膜電位 (Membrane Potential)と呼ぶのですが，
この膜電位の計算方法をいかに精密に計算するか，または簡単かつなるべく正確に計算するか，がスパイキングニューロンモデルでは鍵になります．
このHHモデルは，現在定義されているスパイキングニューロンモデルの中では，かなり正確かつ精密に再現されているモデルの一つです．
- HH方程式
HHの膜電位方程式は以下の複数の微分方程式からなります．
たくさんあって難しそうに見えますが，実際読んでみるとそうでもありません．
$$ I=C_m\frac{dV}{dt}+g_Kn^4(V-E_K)+g_{Na}m^3h(V-E_{Na})+g_l(V-E_l) \ \ \ (1)$$
これがメインの微分方程式です．
各種変数は，$ I $が入力電流，$ V $が膜電位，$ C_m $が膜容量，$ g_* $がイオンチャネルコンダクタンス，$ n,m,h $は各チャネルの(不)活性パラメータでチャネルが開いている度合いを示します．また，$ E_* $はそのチャネルの平衡電位です．
ちなみに下添字については，$ K $がカリウムチャネル，$ Na $がナトリウムチャネル，$ l $がその他のチャネル(Leakの頭文字)です．
(1)式がなぜ，膜電流イコールになっているのかはわかりませんが，膜電位を観察したいので以下のように変形して考えます．
$$ \frac{dV}{dt}=\frac{1}{C_m}[I-g_Kn^4(V-E_K)-g_{Na}m^3h(V-E_{Na})-g_l(V-E_l)] \ \ \ (1')$$
メインの微分方程式はこれで考えるとして，他の方程式も見てみましょう．
$$\frac{dn}{dt}=\alpha_n(V)(1-n)-\beta_n(V)n \ \ \ (2)$$
$$\frac{dm}{dt}=\alpha_m(V)(1-m)-\beta_m(V)m \ \ \ (3)$$
$$\frac{dh}{dt}=\alpha_h(V)(1-h)-\beta_h(V)h \ \ \ (4)$$
$$\alpha_n(V)=\frac{0.01(10-(V-E_{rest}))}{\exp(\frac{10-(V-E_{rest})}{10})-1} \ \ \ (5)$$
$$\alpha_m(V)=\frac{0.1(25-(V-E_{rest}))}{\exp(\frac{25-(V-E_{rest})}{10})-1} \ \ \ (6)$$
$$\alpha_h(V)=0.07\exp(\frac{-(V-E_{rest})}{20}) \ \ \ (7)$$
$$\beta_n(V)=0.125\exp(\frac{-(V-E_{rest})}{80}) \ \ \ (8)$$
$$\beta_m(V)=4\exp(\frac{-(V-E_{rest})}{18}) \ \ \ (9)$$
$$\beta_h(V)=\frac{1}{\exp(\frac{30-(V-E_{rest})}{10})+1} \ \ \ (10)$$
です．
めちゃくちゃ式が多いですね．
何度も言いますが，複雑そうに見えてそんなに複雑ではないです．
どれも似たような式ばかりですしね．
さて式(2)~(10)では何度も$ (V-E_{rest}) $が登場しますが，これは膜電位初期値$ E_{rest} $を実験によって変える場合の式です．
実際には$ E_{rest}=65 $あたりで定義して先に計算した式がよく見られます．
Pythonで実装する
HH方程式をなんとなく理解したら，実装していきます．
こういった時間微分方程式を解くときは，タイムステップを決めてコードに落とすことがよくあります．
式(1')は，
$$ dV=\frac{1}{C_m}[I-g_Kn^4(V-E_K)-g_{Na}m^3h(V-E_{Na})-g_l(V-E_l)]\cdot dt \ \ \ (1'')$$
としてしまう，ということです．
この手法で方程式を全て実装し，実験してみましょう！
"""
Hodgkin-Huxley Neuron model
Copyright(c) HiroshiARAKI
"""
import numpy as np
import matplotlib.pyplot as plt
class HodgkinHuxley:
def __init__(self, time, dt, rest=-65., Cm=1.0, gNa=120., gK=36., gl=0.3, ENa=50., EK=-77., El=-54.387):
"""
Initialize Neuron parameters
:param time: experimental time
:param dt: time step
:param rest: resting potential
:param Cm: membrane capacity
:param gNa: Na+ channel conductance
:param gK: K+ channel conductance
:param gl: other (Cl) channel conductance
:param ENa: Na+ equilibrium potential
:param EK: K+ equilibrium potential
:param El: other (Cl) equilibrium potentials
"""
self.time = time
self.dt = dt
self.rest = rest
self.Cm = Cm
self.gNa = gNa
self.gK = gK
self.gl = gl
self.ENa = ENa
self.EK = EK
self.El = El
def calc_v(self, i):
""" compute membrane potential """
# initialize parameters
v = self.rest
n = 0.32
m = 0.05
h = 0.6
v_monitor = []
n_monitor = []
m_monitor = []
h_monitor = []
time = int(self.time / self.dt)
# update time
for t in range(time):
# calc channel gating kinetics
n += self.dn(v, n)
m += self.dm(v, m)
h += self.dh(v, h)
# calc tiny membrane potential
dv = (i[t] -
self.gK * n**4 * (v - self.EK) - # K+ current
self.gNa * m**3 * h * (v - self.ENa) - # Na+ current
self.gl * (v - self.El)) / self.Cm # other current
# calc new membrane potential
v += dv * self.dt
# record
v_monitor.append(v)
n_monitor.append(n)
m_monitor.append(m)
h_monitor.append(h)
return v_monitor, n_monitor, m_monitor, h_monitor
def dn(self, v, n):
return (self.alpha_n(v) * (1 - n) - self.beta_n(v) * n) * self.dt
def dm(self, v, m):
return (self.alpha_m(v) * (1 - m) - self.beta_m(v) * m) * self.dt
def dh(self, v, h):
return (self.alpha_h(v) * (1 - h) - self.beta_h(v) * h) * self.dt
def alpha_n(self, v):
return 0.01 * (10 - (v - self.rest)) / (np.exp((10 - (v - self.rest))/10) - 1)
def alpha_m(self, v):
return 0.1 * (25 - (v - self.rest)) / (np.exp((25 - (v - self.rest))/10) - 1)
def alpha_h(self, v):
return 0.07 * np.exp(-(v - self.rest) / 20)
def beta_n(self, v):
return 0.125 * np.exp(-(v - self.rest) / 80)
def beta_m(self, v):
return 4 * np.exp(-(v - self.rest) / 18)
def beta_h(self, v):
return 1 / (np.exp((30 - (v - self.rest))/10) + 1)
if __name__ == '__main__':
# init experimental time and time-step
time = 100
dt = 0.05
# create Hodgkin-Huxley Neuron
neu = HodgkinHuxley(time, dt)
# Input data (sin curve)
input_data = np.sin(0.5 * np.arange(0, time, dt))
input_data = np.where(input_data > 0, 20, 0) + 10 * np.random.rand(int(time/dt))
input_data_2 = np.cos(0.3 * np.arange(0, time, dt) + 0.5)
input_data_2 = np.where(input_data_2 > 0, 10, 0)
input_data += input_data_2
# calc neuron status
v, m, n, h = neu.calc_v(input_data)
# plot
x = np.arange(0, time, dt)
plt.subplot(3, 1, 1)
plt.title('Hodgkin-Huxley Neuron model Simulation')
plt.plot(x, input_data)
plt.ylabel('I [μA/cm2]')
plt.subplot(3, 1, 2)
plt.plot(x, v)
plt.ylabel('V [mV]')
plt.subplot(3, 1, 3)
plt.plot(x, n, label='n')
plt.plot(x, m, label='m')
plt.plot(x, h, label='h')
plt.xlabel('time [ms]')
plt.ylabel('Conductance param')
plt.legend()
plt.show()
今回は，ノイズを載せた矩形サインカーブ二つを入力電流としてみました．
実際に上記コードを実行してみると，
このようになりました．
LIFモデルのように不応期を別で定義したり，発火閾値を定義しなくとも膜電位の脱分極や過分極が再現できていることがわかります．
おわりに
今回は自分の勉強がてらにHodgkin-Huxleyニューロンモデルを実装してみました．
勉強してわかりましたが，やはりこのモデルは機械学習には向いてはいないかもしれませんね．
LIFと比べてどれほど計算時間に差が生まれるかはわかりませんが，とにかくハイパーパラメータが多いのは機械学習にとって厄介です．
SNNでは結局デジタル信号として情報を扱うので，膜電位がいくら正確でも無駄になってしまう部分も多いかもしれませんね．