2 一個のニューロン ── 膜電位のダイナミクスと発火率
前章で、人工ニューロンが生物のニューロンから何を捨てているかを見た。時間を捨て、スパイクを捨て、樹状突起の構造を捨てていた。この章では、捨てる前のものを見る。
なぜそんなことをするのか。理由は二つある。第一に、活性化関数の由来を知るためである。シグモイド関数はどこから来たのか。ReLU はなぜあの形なのか。深層学習の教科書では「経験的にうまくいくから」で済まされることが多いが、もともとは実在するニューロンの入出力特性の近似だった。その出所を見ておくと、後で活性化関数を選ぶときの感覚が違ってくる。
第二に、力学系の道具を手に入れるためである。微分方程式・相平面・分岐—これらは次章のウィルソン・コーワン方程式でそのまま使う。そして本章第4節では、本書を貫くことになるヤコビアンが最初に顔を出す。
ホジキンとハクスレーは1952年、ヤリイカの巨大軸索を使って活動電位の発生機構を解明した。彼らの方程式は、いまも計算論的神経科学の出発点である。まずはそこから始めよう。
各節は「なぜこれをやるか → 必要な数学 → 導出 → 数値で確かめる」の順に進む。
第1節の常微分方程式に見覚えがあれば飛ばしてよい。第2節で膜を電気回路として書き、第3節でホジキン・ハクスレー方程式に進む。第4〜5節の相平面と分岐は、第3章のネットワークでそのまま使う道具なので、ここで手に馴染ませておきたい。
第6節のスパイク統計は、第6章のデコーディングで使う。ポアソン過程の性質(とくに平均と分散が等しいこと)だけは押さえてほしい。
第7節は、この章を深層学習につなぐ節である。急ぐならここだけでもよい。
1. 必要な数学:常微分方程式と数値解法
一階線形の方程式
本章と次章で繰り返し現れるのは、いまの値が一定の目標値に近づいていく形の方程式である。時間を t、時間とともに変わる一つの量を x = x(t)、その目標値を定数 I と書く。ここでの I は x と同じ単位を持つ量で、まだ電流そのものを意味しない。正の定数 \tau は近づく速さを決める時間尺度で、t と同じ時間の単位を持つ。式にすると、次の形になる。
\tau \frac{dx}{dt} = -x + I
読み方はこうだ—x は I に向かって減衰していく。x > I なら右辺が負なので x は減り、x < I なら増える。x = I で止まる。I が定数なら、解は手で書ける。y = x - I と置くと \tau\, dy/dt = -y となり、
y(t) = y(0)\, e^{-t/\tau} \quad\Longrightarrow\quad x(t) = I + \big(x(0) - I\big)\, e^{-t/\tau}
\tau は時定数である。t = \tau で、初期値と最終値の差が 1/e \approx 0.37 倍に縮む。\tau が小さいほど速く応答する。この形の方程式は本書に何度も出る。膜電位の緩和(本章)、集団活動の緩和(第3章)、そして勾配流(第8章)。「今の値と目標値の差に比例して動く」という構造は、それだけ普遍的である。
数値解法
一般には解析解が得られないので、数値的に解く。もっとも単純なのが オイラー法である。g(x, t) を「値が x、時刻が t のときの変化率 dx/dt を返す関数」とし、\Delta t を一回に進める短い時間とする。いまの変化率に \Delta t を掛ければ、そのあいだの変化量を近似できる。ここでの g は関数の名前であって、次節に出てくるコンダクタンスとは別の記号である。
\frac{dx}{dt} = g(x, t) \quad\Longrightarrow\quad x(t + \Delta t) \approx x(t) + \Delta t\, g\big(x(t), t\big)
「いまの傾きで \Delta t だけ直進する」を繰り返すだけだ。上の方程式で x(0)=0、I=10、\tau=10、\Delta t=1 とすれば、最初の傾きは (10-0)/10=1 なので一歩後は x \approx 1。解析解は 10(1-e^{-0.1}) \approx 0.95 で、少し進みすぎている—一歩のあいだに傾きが小さくなるのを無視しているからだ。実装は数行で済む。
def euler(g, x0, T, dt):
n = int(T / dt)
xs = np.empty((n, *np.shape(x0)))
x = np.array(x0, dtype=float)
for i in range(n):
xs[i] = x
x = x + dt * g(x, i * dt)
return xs注意すべきは刻み幅である。\Delta t が時定数 \tau に比べて大きいと、数値解が振動したり発散したりする。目安は \Delta t \ll \tau。ホジキン・ハクスレー方程式(Hodgkin–Huxley、第3節)では \tau が 0.1 ms 程度まで小さくなるので、\Delta t = 0.01 ms 程度が要る。
より精度の高い方法(Runge–Kutta 法など)もあるが、本書ではオイラー法で足りる場面がほとんどである。
2. 膜を RC 回路として見る
RC 回路とは、抵抗(R)と電荷を蓄えるコンデンサ(C)を組み合わせた回路である。
膜電位とは何か
ニューロンの細胞膜は、内外でイオンの濃度が違う。細胞内はカリウムイオン \mathrm{K}^+ が多く、細胞外はナトリウムイオン \mathrm{Na}^+ が多い(上付きの + は正の電荷を持つ印である)。この濃度差があり、しかも膜が特定のイオンだけを通す(静止時は主に \mathrm{K}^+)ために、膜の内外に電位差が生じる。これが膜電位 V である。静止時には約 -70 mV(細胞内が負)。
膜は脂質でできていて電気を通さないが、イオンチャネルという穴が開いていて、特定のイオンだけを通す。この構造は、電気回路として書ける。
以下では、容量・コンダクタンス・電流をいずれも膜の単位面積あたりで表す。時間 t は ms、膜電位 V は細胞外を基準にした細胞内の電位で mV である。
- 膜そのもの — 電気を通さない絶縁体が二つの導体(細胞内液と細胞外液)を隔てている。これはコンデンサである。電荷を蓄える能力を表す膜容量を C とし、単位は μF/cm² とする
- イオンチャネル — イオンを通す経路。これは抵抗(コンダクタンス g)である。コンダクタンスは電気の通りやすさ、つまり抵抗の逆数で、単位は mS/cm² である
- 濃度差 — イオンを押し流す力。これは電池(平衡電位 E)である。濃度差による力と電気的な力が釣り合って、そのイオンの正味の電流がゼロになる膜電位が平衡電位で、単位は mV である
膜方程式
回路の法則(電流の保存)から式が出る。記号を決めておこう。g_L は漏れの経路のコンダクタンス(mS/cm²)、E_L は漏れ電流がゼロになる電位(mV)で、添字の L は漏れ(leak)を表す。I_{\text{ext}} は外から注入する電流密度(μA/cm²)で、正なら膜電位を上げる向きとする。以下では電流密度を単に電流と呼ぶ。外部から注入した電流は、コンデンサに溜まる分と、膜を漏れ出ていく分に分かれる—I_{\text{ext}} = C\,dV/dt + g_L(V - E_L) である。右辺第2項は外向きを正とした漏れ電流で、各項の単位はいずれも μA/cm² になる。C\,dV/dt について解けば、次の形になる。
C \frac{dV}{dt} = -g_L (V - E_L) + I_{\text{ext}}
右辺第一項が漏れ電流(leak)である。V が E_L より高ければ電流が外へ流れ、V を下げる。負号がついているのはそのためだ。この式を、第1節の形に整理してみよう。両辺を g_L で割る。
\underbrace{\frac{C}{g_L}}_{=\,\tau_m} \frac{dV}{dt} = -(V - E_L) + \frac{I_{\text{ext}}}{g_L}
第1節の \tau\, dx/dt = -x + I そのものである。膜は「入力に向かって時定数 \tau_m = C/g_L で緩和する系」だと分かる。典型的には \tau_m \approx 10〜20 ms。
リーク積分発火モデル
上の式には、スパイクが出てこない。ずっと緩和し続けるだけである。そこで、手で閾値を入れる。「V が閾値 V_{\text{th}} に達したらスパイクを出し、V を V_{\text{reset}} に戻す」という規則を付け加えるのだ。これがリーク積分発火モデル(leaky integrate-and-fire, LIF)である。R = 1/g_L はリークコンダクタンスの逆数、つまり膜抵抗である。
\tau_m \frac{dV}{dt} = -(V - E_L) + R I_{\text{ext}}, \qquad V \ge V_{\text{th}} \;\Rightarrow\; \text{スパイク発火},\; V \leftarrow V_{\text{reset}}
乱暴だが、驚くほど有用なモデルである。スパイクの波形は諦めるが、「いつ発火するか」だけなら実際のニューロンをかなりよく予測する。
f–I 曲線
LIF モデルで、入力電流を上げると発火率がどう変わるかを計算してみよう。これが後で活性化関数につながる。
先に数で見よう。\tau_m = 10 ms、E_L = V_{\text{reset}} = -70 mV、V_{\text{th}} = -50 mV、入力による上昇分を RI = 30 mV とすると、電位は -40 mV へ向かうので閾値に届く。-50 = -40 - 30e^{-T/10} を解くと T = 10\ln 3 \approx 11 ms、発火率は 1000/11 \approx 91 Hz である。
一般に書こう。定常入力 I を与える。V が V_{\text{reset}} から V_{\text{th}} まで上がるのにかかる時間を求めればよい。第1節の解を使うと、V(t) = E_L + RI + (V_{\text{reset}} - E_L - RI)e^{-t/\tau_m} である。V(T) = V_{\text{th}} と置いて T について解く。
T = \tau_m \ln \frac{E_L + RI - V_{\text{reset}}}{E_L + RI - V_{\text{th}}}
同じ時間 T ごとにスパイクが出るのだから、単位時間あたりのスパイク数である発火率は f = 1/T である。T を秒で表せば f の単位は Hz(1秒あたりの回数)になる。T を ms で計算したときは、1/T を1000倍して Hz に直す。入力電流 I と発火率 f の関係を f–I 曲線(frequency–current curve)と呼ぶ。
形を読もう。定常入力 I を与え続けたとき V が向かう先は V_\infty = E_L + RI である。発火するのは、この行き先が閾値を超えているとき—V_\infty > V_{\text{th}} のとき—に限る。超えていなければ、V は閾値の下で止まってしまい、いつまでも発火しない。f = 0 である。式のうえでは、V_\infty > V_{\text{th}} のときに対数の中身の分母 V_\infty - V_{\text{th}} が正になり、V_{\text{reset}} < V_{\text{th}} と合わせて分子 V_\infty - V_{\text{reset}} も正になる—そこで初めて T が意味を持つ。向きに注意してほしい。閾値をわずかに超えただけのときは T が大きく、f はごく小さい(対数が発散するからである)。入力をさらに上げると T が縮み、f が立ち上がる。入力をさらに上げると f は増え続けるが、増え方は鈍る(対数の中身が 1 に近づくため)。
「閾値までゼロ、超えたら急に立ち上がり、やがて増え方が緩やかになる」—この形が、活性化関数の原型である。ただし LIF そのものは飽和しない(電流を上げれば発火率はいくらでも上がる)。飽和が現れるのは、発火直後に次の発火が起きにくい期間(不応期)を入れたときで、第7節でそこに戻る。
3. ホジキン・ハクスレー方程式
LIF モデルは閾値を手で入れた。ホジキンとハクスレーがやったのは、閾値そのものを導くことだった。
イオンチャネルは電位で開閉する
鍵となる発見はこうである。\mathrm{Na}^+ チャネルと \mathrm{K}^+ チャネルのコンダクタンスは、膜電位に依存して変化する。
しかも変化の仕方が違う。膜電位が上がると、
- \mathrm{Na}^+ チャネルは速く開き、その後ゆっくり閉じる(不活性化する)
- \mathrm{K}^+ チャネルはゆっくり開く
この時間差が、活動電位という現象を生む。
方程式
ホジキンとハクスレーは、コンダクタンスをゲート変数で表した。記号を先に読んでおこう。\bar{g}_{\mathrm{Na}} と \bar{g}_{\mathrm{K}} は、それぞれナトリウムとカリウムのチャネルがすべて開いたときの最大コンダクタンス(上の棒が最大値の印。単位は mS/cm²)、E_{\mathrm{Na}} と E_{\mathrm{K}} は各イオンの平衡電位(mV)である。m はナトリウムの活性化ゲート、h はナトリウムの不活性化ゲート、n はカリウムの活性化ゲートが開いている確率で、いずれも単位を持たない 0 から 1 の数である(h が大きいほど、不活性化で塞がれていない)。ゲートが独立に開くと考えると、ナトリウムは三つの m ゲートと一つの h ゲートが同時に開く確率 m^3 h、カリウムは四つの n ゲートが同時に開く確率 n^4 になる。これを最大コンダクタンスに掛け、さらに電位差を掛けたものが各イオンの電流である。
C \frac{dV}{dt} = -\underbrace{\bar{g}_{\mathrm{Na}}\, m^3 h\, (V - E_{\mathrm{Na}})}_{\mathrm{Na}^+ \textsf{ 電流}} -\underbrace{\bar{g}_{\mathrm{K}}\, n^4\, (V - E_{\mathrm{K}})}_{\mathrm{K}^+ \textsf{ 電流}} -\underbrace{g_L (V - E_L)}_{\textsf{漏れ電流}} + I_{\text{ext}}
ゲート変数 m, h, n はそれぞれ [0,1] の値を取り、次の方程式に従う(x は m, h, n のいずれか)。ここで \alpha_x(V) は閉じたゲートが開く速さ、\beta_x(V) は開いたゲートが閉じる速さで、どちらも膜電位 V の関数であり、単位は ms^{-1} である。閉じている割合 1-x に開く速さを掛けた分だけ x が増え、開いている割合 x に閉じる速さを掛けた分だけ減る、と読める。
\frac{dx}{dt} = \alpha_x(V)\,(1 - x) - \beta_x(V)\, x
これも第1節の形である。実際、\tau_x(V) = 1/(\alpha_x + \beta_x)、x_\infty(V) = \alpha_x/(\alpha_x+\beta_x) と置くと、
\tau_x(V) \frac{dx}{dt} = -x + x_\infty(V)
「x は、その電位での目標値 x_\infty(V) へ、時定数 \tau_x(V) で緩和する」—という読み方ができる。時定数も目標値も電位に依存する、というのが HH モデルの心臓部である。
各変数の役割
三つのゲート変数の性格を、はっきりさせておこう。
| 変数 | 対象 | 電位が上がると | 速さ |
|---|---|---|---|
| m | \mathrm{Na}^+ 活性化 | 増える(開く) | 速い(〜0.1 ms) |
| h | \mathrm{Na}^+ 不活性化 | 減る(閉じる) | 遅い(〜1 ms) |
| n | \mathrm{K}^+ 活性化 | 増える(開く) | 遅い(〜1 ms) |
\mathrm{Na}^+ 電流の項に m^3 h という積が入っていることに注目してほしい。m が大きく、かつ h が大きいときだけ電流が流れる。m は速く上がり、h は遅れて下がる。だから \mathrm{Na}^+ 電流は「一瞬だけ流れて、すぐ止まる」。
活動電位が起きる筋書き
以上を組み合わせると、活動電位の物語が読める。一手ずつ追ってほしい。
- 外部電流で V が少し上がる
- m が速く増える → \mathrm{Na}^+ が流れ込む → V がさらに上がる
- これは正のフィードバックである。V が上がるほど m が増え、m が増えるほど V が上がる。爆発的に V が跳ね上がる(立ち上がり)
- 遅れて h が減る → \mathrm{Na}^+ 電流が止まる
- 同じく遅れて n が増える → \mathrm{K}^+ が流れ出る → V が下がる(立ち下がり)
- V が静止電位を下回る(過分極)。h が戻りきらず、n もゆっくり戻るあいだは、次の発火が起きにくい(不応期)。発火直後にまったく発火できない時期があるのは、h が下がったままだからである
ここがポイント
活動電位の本質は、速い正のフィードバック(m)と、遅い負のフィードバック(h と n)の組み合わせである。速い自己増強と遅い抑制—この構造は神経系のあちこちに現れ、第3章では集団レベルの振動を生む同じ構造を見ることになる。
そして「閾値」は、方程式のどこにも書かれていない。正のフィードバックが暴走を始める境目が、結果として閾値のように見えるだけである。LIF モデルが手で入れたものを、HH モデルは導出した—これが「機構モデル」の力である(第1章第2節)。
HH 方程式をオイラー法で解いてみよう。パラメータは原論文の値、刻み幅は \Delta t = 0.01 ms である。コードは コード/02_hh.py にある。
C, gNa, gK, gL = 1.0, 120.0, 36.0, 0.3 # μF/cm², mS/cm²
ENa, EK, EL = 50.0, -77.0, -54.387 # mV
DT = 0.01 # ms
def simulate(I, T=300.0):
V = -65.0
am, bm, ah, bh, an, bn = rates(V)
m = am/(am+bm) # 静止状態から始める
h = ah/(ah+bh)
n = an/(an+bn)
spikes = []
for i in range(int(T/DT)):
am, bm, ah, bh, an, bn = rates(V)
m += DT*(am*(1-m) - bm*m)
h += DT*(ah*(1-h) - bh*h)
n += DT*(an*(1-n) - bn*n)
I_ion = (gNa*m**3*h*(V-ENa) + gK*n**4*(V-EK)
+ gL*(V-EL))
V_new = V + DT*(I - I_ion)/C
if V < 0 <= V_new: # 0 mV の上向き通過=スパイク
spikes.append(i*DT)
V = V_new
return spikes(rates はゲート変数の \alpha, \beta を返す関数。全体は コード/02_hh.py を見てほしい。)
I_{\text{ext}} を階段状に与えると、活動電位が出る。V(t) と一緒に m(t), h(t), n(t) も描いてみてほしい。上の1〜6の筋書きが、グラフの上で読めるはずである。そして I_{\text{ext}} を少しずつ上げて、発火が始まる電流を探す。結果はこうなる。
def firing_rate(I, T=300.0, warmup=100.0):
"""定常状態での発火率 [Hz]。最初の 100 ms は捨てる。"""
s = [t for t in simulate(I, T) if t > warmup]
if len(s) < 2:
return 0.0
return 1000.0*(len(s) - 1)/(s[-1] - s[0])
# 発火が始まる電流を二分法で詰める
lo, hi = 0.0, 10.0
for _ in range(40):
mid = (lo + hi)/2
if firing_rate(mid) > 0: hi = mid
else: lo = mid
print(hi, firing_rate(hi + 0.1))| I [μA/cm²] | 2.0 | 5.0 | 6.2 | 6.3 | 10 | 20 | 40 | 80 |
|---|---|---|---|---|---|---|---|---|
| 発火率 [Hz] | 0 | 0 | 0 | 53 | 69 | 87 | 109 | 0 |
発火率がゼロから 53 Hz へ飛んでいる。いくらでも遅く発火することはできない—これが第5節で扱う Type II の特徴である。ただし、ここで測った発火の境目と、静止状態そのものが不安定になる境目は別である。どの状態から電流を加えるかも効いてくる。
表のいちばん右に注目してほしい。電流をさらに上げると、数えられるスパイクが止まる。全体を描いたのが 図 1 である。
コード/02_hh.py の出力)。発火率はゼロから飛び、電流とともに上がる。だが高電流域ではスパイクの頂点が 0 mV に届かなくなり、この数え方での発火率はゼロになる。
このコードは、膜電位 V が 0 mV を上向きに横切った回数を数えている。I=80 では、その高さに届かなくても膜電位は振動し続ける。さらに電流を上げると振動自体が消え、膜が脱分極したまま固定される脱分極ブロックに至る。
だから表の「0」には二種類ある。弱い入力での静止と、振動はあるが測定の高さを横切らない状態だ。指標がゼロになったことと、現象が消えたことは違う。境界値と計算時間による違いは、第2章第3節「手を動かす」(Web版の追加実験)で確かめられる。
自分で V(t) を描いて確かめてほしい。そのうえで問い—強く脱分極させ続けると振幅が縮むのは、なぜだろうか。h の働きから考えてみよう。
観測時間を変えて境目を確かめる
コード/02_hh.pyで、同じ電流について計算時間を延ばし、0 mV の通過回数と膜電位の最大・最小を別々に記録しよう。コードの初期値と刻み幅は上の例にそろえる。
| 電流密度 I [μA/cm²] | 点検すること |
|---|---|
| 約6.21 | この初期値から発火しても、固定点はまだ安定 |
| 約9.78 | 固定点の固有値の実部が符号を変える |
| 約64 | 頂点が0 mVに届かなくなり、数えた発火率がゼロになる |
| 約154.5 | 高電流側の固定点が安定になり、振動自体が消える |
安定性の判定は第4〜5節を読んでから確かめてほしい。I=6.21 の主固有値(実部が最大の固有値)の実部は約 -0.071 ms^{-1} である。高電流側では短い計算だと約158まで振動が残るように見える。計算を延ばして振幅が縮み続けるなら、それを持続する振動と数えてはいけない。これらは配布コードの設定での境界値であり、実際の細胞一般に共通する数ではない。
4. 相平面解析と次元削減
HH 方程式は4変数(V, m, h, n)の非線形系である。4次元は絵に描けない。2次元に落とせないだろうか。
縮約の論理
第3節の表を見返すと、手がかりがある。
- m は非常に速い。他の変数が動く時間スケールでは、m はすでに目標値に達している。だから m \approx m_\infty(V) と置いてよい(断熱近似)
- h と n は同じくらいの速さで、しかも h が減るときに n が増える。この動作範囲では、数値的に h + n \approx \text{const} が成り立つ。だから一つの変数にまとめられる
こうして4変数が2変数に減る。V と、まとめた回復変数 w である。
FitzHugh–Nagumo モデル
この縮約の考え方—速い興奮変数と遅い回復変数の二つで書く—を、数学的に扱いやすい形に置き直したのが FitzHugh–Nagumo モデルである。
先に断っておく。次の式は、上の二つの近似を HH 方程式に代入して導かれるものではない。w は h と n をまとめた実効的な回復変数であり、三次式 V - V^3/3 は「速い正のフィードバック」を扱いやすく表した便宜的な形である。HH の定量的な縮約ではなく、同じ定性的構造を持つ別のモデルだと思ってほしい。第1章第2節の言葉で言えば、HH が機構モデルなのに対し、こちらは現象の再現をねらう記述モデルに近い。
\begin{aligned} \frac{dV}{dt} &= V - \frac{V^3}{3} - w + I \\ \tau_w \frac{dw}{dt} &= V + a - b w \end{aligned}
ここでは尺度を取り直しているので、V、w、I、t はいずれも物理的な単位を持たないモデル上の量である。V は電位に対応する興奮の程度、w はそれを抑える回復の程度、I は外部入力、t はモデル時間を表す。HH 方程式と同じ記号を使っているが、V をそのまま mV、I を μA/cm² と読むことはできない。定数 a、b、\tau_w も単位を持たず、a は回復変数への駆動をずらし、b はその自己減衰の強さを決め、\tau_w は回復変数の時間尺度を定める。
生物学的な意味は薄れたが、構造は保たれている。V の式の V - V^3/3 が正のフィードバック(小さい V では V が自己増強する)を表し、-w が負のフィードバックである。\tau_w が大きいので w は遅い。
相平面
2変数なら、(V, w) 平面に軌道を描ける。これが相平面である。
平面の上に、二本の特別な曲線を引く。ヌルクライン(nullcline)である。
- V-ヌルクライン — dV/dt = 0 となる曲線。ここでは軌道が垂直に動く
- w-ヌルクライン — dw/dt = 0 となる曲線。ここでは軌道が水平に動く
FitzHugh–Nagumo なら、V-ヌルクラインは w = V - V^3/3 + I という三次曲線(N 字形)、w-ヌルクラインは w = (V+a)/b という直線である。
二本の交点が固定点である。そこでは両方の微分がゼロだから、系は動かない。
FitzHugh–Nagumo モデルの相平面を描いたのが 図 2 である。
固定点の安定性
固定点があっても、それが安定とは限らない。少しずらしたとき、戻ってくるのか離れていくのか。
調べ方は、固定点のまわりで線形化することである。固定点のすぐ近くでは、変化率を「固定点からのずれに比例する部分」で近似できる。これが線形化だ。記号を先に決めておく。\dot{V} は dV/dt の略記である(以後この書き方を使う)。固定点を (V^*, w^*) とし、そこからのずれを \delta V = V - V^*、\delta w = w - w^* と置いて、縦に並べたベクトルを \mathbf{u} = (\delta V, \delta w)^\top \in \mathbb{R}^2 と書く(\top は横の並びを縦にする転置である)。このあと出てくる J \in \mathbb{R}^{2 \times 2} は、ずれが変化率にどう伝わるかを並べた表で、行が上から \dot{V}, \dot{w}、列が左から V, w に対応する。たとえば右上の \partial \dot{V}/\partial w は、V を止めて w を少し動かしたときに \dot{V} がどれだけ変わるかである。四つの偏微分は、すべて固定点 (V^*, w^*) で評価する。すると、
\frac{d}{dt}\begin{pmatrix}\delta V \\ \delta w\end{pmatrix} \approx \underbrace{\begin{pmatrix} \partial \dot{V}/\partial V & \partial \dot{V}/\partial w \\ \partial \dot{w}/\partial V & \partial \dot{w}/\partial w \end{pmatrix}}_{= \, J \ \textsf{(ヤコビ行列)}} \begin{pmatrix}\delta V \\ \delta w\end{pmatrix}
ヤコビ行列 J が現れた。本書で五つの役目を負うことになる J の、最初の登場である(第1章第3節の表を参照)。ここでの役目は固定点の安定性を決めることである。
覚えておいてほしいのは、J の列が「何を動かすか」、行が「何の変化を見るか」を表すことだ。ここでは状態 (V, w) から変化率 (\dot{V}, \dot{w}) への対応を見ている。第9章ではパラメータから出力への対応、第11章では刺激から応答への対応、第15章では生成のための変数変換を見ることになる。大きさも成分も変わるが、微小な変化の伝わり方を並べるという作り方は同じである。
線形の方程式 \dot{\mathbf{u}} = J\mathbf{u} の解は、J の固有値を使って書ける。固有値とは、J\mathbf{v} = \lambda\mathbf{v} を満たす \lambda のことである(\mathbf{v} \ne \mathbf{0} をその固有ベクトルと呼ぶ)。
固有値が実数のときは話が簡単だ。その固有ベクトルの向きに置かれたずれは、向きを変えずに e^{\lambda t} 倍に伸び縮みするだけである。\lambda < 0 なら縮み、\lambda > 0 なら伸びる。
固有値は複素数にもなりうる。J の成分が実数なら、複素固有値は必ず共役な対 \alpha \pm i\beta で現れる(i は虚数単位で i^2 = -1。\alpha を実部、i に掛かる \beta を虚部と呼び、虚部の符号を反転したものを共役と呼ぶ)。この場合、向きを保つ実の方向はもうない。代わりに二次元の平面の上で、ずれが回転しながら伸び縮みする。e^{(\alpha + i\beta)t} = e^{\alpha t}(\cos\beta t + i\sin\beta t) という形から読めるとおり、e^{\alpha t} が伸び縮みを、三角関数が回転を担う。振動の角周波数は |\beta| である(対のどちらを +\beta と呼ぶかは決められないので、大きさだけを見る。もとの座標では軌道は楕円になるので、「一定の角速度で回る」と言えるのは適当な座標に取り替えたときである)。
判定は単純である。
- すべての固有値の実部が負 → ずれは減衰する → 安定
- 一つでも実部が正 → ずれは増大する → 不安定
- 複素固有値 → 実部が振動の減衰・増大、虚部の大きさが振動の速さを与える
実際に作ってみよう。FitzHugh–Nagumo の右辺を偏微分すると、固定点 (V^*, w^*) でのヤコビ行列は
J = \begin{pmatrix} 1 - V^{*2} & -1 \\ 1/\tau_w & -b/\tau_w \end{pmatrix}
となる。四つの成分は、どれも「一方の変数を少し動かしたとき、もう一方の変化率がどれだけ変わるか」という局所的な傾きである。四隅を順に読もう。
- 左上 1 - V^{*2} — V のずれが V 自身の変化率に返る効果。|V^*| < 1 なら正で、ずれが自分自身を増幅する(速い正のフィードバック)。|V^*| > 1 なら負で、ずれは自分で収まる
- 右上 -1 — 回復変数がブレーキとして効く
- 左下 1/\tau_w — V が上がると回復変数が追いかけて上がる
- 右下 -b/\tau_w — 回復変数が自分で減衰する。\tau_w が大きいほどゆっくり
四隅を合わせて読むと、第3節で見た「速い正のフィードバックと遅い負のフィードバック」という描像が、そのまま数字になっていることが分かる。あとは固有値の実部の符号を読めばよい。数を入れて一度通してみよう。a = 1/3、b = 1、\tau_w = 10、I = 0 とすると (V^*, w^*) = (-1, -2/3) が固定点である(両式の右辺が -1 + 1/3 + 2/3 = 0 になる)。V^* = -1 を上の J に入れると成分は 0, -1, 0.1, -0.1 で、固有値は \lambda^2 + 0.1\lambda + 0.1 = 0 を解いて -0.05 \pm 0.312\,i。実部が負で虚部があるから、この固定点は回り込みながら戻る安定な点である。
この段階で身につけてほしいのは「ヤコビ行列を作り、固有値の実部の符号を読む」という手順だけである。一般次元での固有値の理論や、後で出てくる分岐の非退化条件は、いまは直感でよい。
ここがポイント
固定点の安定性は、そこでのヤコビ行列の固有値の実部の符号で決まる。すべての実部が負なら安定、一つでも正なら不安定である(実部がちょうどゼロのときは、この行列だけでは決まらない)。これは本書で最初に現れるヤコビアンの用法である。第3章でネットワークに、第9章で学習ダイナミクスに、第11章で表現の幾何に、同じ行列が別の顔で現れる。
5. 分岐 ── 発火が始まる瞬間の数学
第3節の「手を動かす」で問うたことに答える。入力電流を上げていくと、発火はどう始まるのか。
分岐とは
パラメータ(ここでは I)を連続的に変えていくと、固定点の位置や安定性が変わる。ある値で、解の質的な性質が変わることがある。これを分岐(bifurcation)と呼ぶ。
「質的に変わる」というのは、固定点の数が変わる、安定だったものが不安定になる、周期解が生まれる—といったことである。
二つの型
神経モデルで主に現れるのは、次の二つである。一つはサドルノード分岐で、安定な固定点とある向きからは近づき、別の向きには離れる不安定な固定点(サドル)が近づいてきて、衝突して消滅する。
ただし、固定点が消えただけでは周期軌道が現れる保証はない—別の固定点に落ち着いてしまうこともある。発火が始まるときに起きているのは、この衝突が閉じた軌道の上で起こる場合である。これを不変円上のサドルノード分岐(saddle-node on invariant circle, SNIC)と呼ぶ。二つの固定点は円周上で衝突して消え、あとには円周に沿った流れだけが残る。だから系はその円を回り続ける—それが周期発火になる。このとき、発火率はゼロから連続的に立ち上がる。消滅の直前には、軌道が「元固定点があった場所」の近くで長く足踏みするため、周期が非常に長くなる(=発火率が非常に低い)からだ。任意に低い発火率が実現できる。
もう一つは Hopf 分岐である。固定点は残るが、複素固有値の実部が負から正に変わり、不安定になる。そのまわりに周期軌道(リミットサイクル)が生まれる場合が、Type II の発火を生む代表的な筋書きである。このとき、発火率はゼロでない値から突然始まる。振動の周波数は固有値の虚部で決まっており、それは分岐点で有限の値を持つからだ。
二つの分岐の違いは、f–I 曲線の形の違いとして現れる(図 3)。
Type I と Type II
この違いは、実際のニューロンで観察されている。カニの軸索を電流で駆動して、なめらかに発火率が上がる型と、ある電流で一定の周波数から発火が始まる型を区別したのは Hodgkin (1948) である。
| 分岐の型 | f–I 曲線 | 発火開始 | |
|---|---|---|---|
| Type I | SNIC(不変円上のサドルノード) | 連続的に立ち上がる | 任意に低い発火率が可能 |
| Type II | 代表例:Hopf | 不連続に飛ぶ | 最低発火率が存在する |
Type II のニューロンは、低い発火率を出せない。これは機能的な違いを生む。Type I は入力の強さを発火率で細かく符号化できる(積分器として働く)が、Type II は特定の周波数の入力に共鳴しやすい(共鳴器として働く)。そして興味深いことに、同じニューロンでも、細胞の応答しやすさや伝達の強さを調節する神経修飾物質などの作用で Type I と Type II を行き来することが知られている。回路の計算特性が、その場で切り替わりうるわけである。
FitzHugh–Nagumo モデルで分岐図を描いてみよう。コードは コード/02_fhn.py にある。
DT = 0.005
def run(I, a=0.7, b=0.8, tau=12.5, T=2000.0, warm=1000.0):
V, w = -1.0, -0.5
n, nw = int(T/DT), int(warm/DT)
Vs = np.empty(n - nw); sp = []
for i in range(n):
dV = V - V**3/3 - w + I # 速い変数
dw = (V + a - b*w)/tau # 遅い変数
Vn = V + DT*dV
w += DT*dw
if i >= nw: # 前半は過渡なので捨てる
Vs[i-nw] = Vn
if V < 0 <= Vn: # 上向きのゼロ交差=発火
sp.append(i*DT)
V = Vn
f = 0.0 if len(sp) < 2 else \
1000.0*(len(sp) - 1)/(sp[-1] - sp[0])
return Vs.min(), Vs.max(), f # 分岐図の上下と発火率I を動かし、過渡が収まったあとの V の最大値と最小値を記録する。静止なら一点、振動なら二点に分かれる。モデルの時間は無次元だが、ここでは1単位を 1 ms と読んで Hz に換算する。絶対値より変化の仕方を見よう。
| I | 0.32 | 0.33 | 0.80 | 1.50 |
|---|---|---|---|---|
| V の最小 | -0.98 | -1.99 | -1.93 | 1.03 |
| V の最大 | -0.98 | 1.76 | 1.91 | 1.03 |
| 発火率 [Hz] | 0 | 21 | 27 | 0 |
コード/02_fhn.py の出力)。上は定常状態の V の最大と最小、下は発火率。固定した初期値からは、振動する範囲で一点が二点に分かれる。固定点そのものが不安定になる境目とは区別する。
この初期値 (V,w)=(-1,-0.5) から振動へ向かう境目は I\approx0.324、直上の発火率は約19 Hzで、ゼロから飛ぶ。一方、固定点の安定性は第4節のヤコビ行列で調べられる。二つの固有値の和であるトレース(対角成分の和)は
\operatorname{tr}J=1-V^{*2}-b/\tau_w
だった。これがゼロになる低電流側の境目は I\approx0.3313 で、先ほどとは違う。この点で固有値の虚部は約 \pm0.276 とゼロではなく、有限の周波数で始まる Type II に対応する。
二つの境目のあいだ、たとえば I=0.33 では、安定な固定点とリミットサイクルが共存する。初期値を固定点の近くに変えると、表と同じ電流でも静止に落ち着くはずだ。二つの落ち着き先を自分で確かめてほしい。「この初期値から振動した」と「固定点が不安定になった」は別なのである。
I=1.50 では再び振動が止まる。HH の脱分極ブロックと似た現象だが、FHN に Na チャネルの不活性化の変数はなく、機構まで同じとは言えない。細かな刻みの結果は、第2章第5節「手を動かす」(Web版の追加実験)で確かめられる。問い—a,b,\tau_w を変えて Type I にするには、二本のヌルクラインの交わり方をどう変えればよいだろうか。
初期値を変えて二つの落ち着き先を見る
上のコード/02_fhn.pyでは、a=0.7,b=0.8,\tau_w=12.5、刻み幅0.005で2000モデル時間単位を計算し、最初の1000を捨てる。I=0.33 のまま、本文の初期値と固定点の近くの初期値を比べよう。初期値は run 関数の中で固定されているので、その代入を変えて計算する。
固定点 V^*,w^* は
w^*=(V^*+a)/b,\qquad V^*-V^{*3}/3-w^*+I=0
から求める。この点のヤコビ行列は
J=\begin{pmatrix}1-V^{*2}&-1\\1/\tau_w&-b/\tau_w\end{pmatrix}
となる。トレースがゼロになる電流は約0.3313と1.4187で、行列式はどちらでも正。固定した初期値から振動へ向かう境目の約0.324とは異なることを確認できる。
毎回同じ初期値へ戻す方法と、直前の計算の終状態を次の初期値にする方法も比べてほしい。前者は同じ出発点の行き先、後者は一つの落ち着き先をたどっている。入力を上げる場合と下げる場合を含め、どの手順で境目を測ったかを結果に添えよう。
6. スパイク列の統計
ここまでは決定的な方程式だった。だが実際のニューロンの応答は、同じ刺激を繰り返しても毎回違う。
発火率という記述
スパイクの時刻を t_1, t_2, \dots とする(添字はスパイクが起きた順番である)。この列そのものを扱うのは大変なので、多くの場合、単位時間あたりの平均的なスパイク数である発火率に縮約する。時間とともに変わる発火率を r(t) と書く。前節までの f は周期発火の頻度だったが、r(t) は試行ごとに揺らぐスパイク数の期待値から定める量である。時間を秒で測れば、単位はどちらも Hz になる。
r(t) = \lim_{\Delta t \to 0} \frac{\textsf{区間 } [t, t+\Delta t) \textsf{ のスパイク数の期待値}}{\Delta t}
「期待値」と書いたことに注意してほしい。一回の試行だけからは推定できない。多数の試行にわたる平均か、時間方向の平均(定常であることに加えて、時間平均が試行平均に一致するというエルゴード性を仮定する)が必要である。
ポアソン過程
もっとも単純なモデルは、重ならない時間区間のスパイク数が互いに独立で、十分短い区間ではスパイクが1個起きる確率が r\,\Delta t に比例し、2個以上はそれより高位の微小量になる、というものである。これがポアソン過程である。いまのように率 r が一定のものを斉次、率が時間とともに変わる r(t) のものを非斉次と呼ぶ。刺激に応じて発火率が変わる場面では、後者を使うことになる。
以下の三つの性質は、率 r が一定の斉次の場合について確認する。性質1:時間 T の間のスパイク数 n はポアソン分布に従う。ここでの n は非負の整数で、第3節のゲート変数 n とは別物である。P(n) はスパイクがちょうど n 個入る確率、rT はその区間の平均スパイク数を表す(r を Hz で表すなら T は秒で測る)。
次式の n! は階乗、つまり 1 から n までを掛けた値である(3! = 6、0! = 1 と約束する)。
P(n) = \frac{(rT)^n}{n!}\, e^{-rT}
性質2:平均と分散が等しい。\mathbb{E}[n] = \mathrm{Var}[n] = rT。ここで分散はばらつきの二乗の平均 \mathrm{Var}[n] = \mathbb{E}[(n - \mathbb{E}[n])^2] で、その平方根が標準偏差である。
ばらつきの大きさが、平均の大きさで決まってしまう—独立なコイン投げの重ね合わせだから当然だが、この性質は第6章で決定的に効く。性質3:隣り合うスパイクの時間間隔(ISI)は指数分布に従う。この間隔を \tau \ge 0 と書く(第1節の時定数 \tau とは別の量である)。p(\tau) は間隔の確率密度で、間隔が \tau から \tau + d\tau の狭い範囲に入る確率が p(\tau)\, d\tau で近似できる、という意味だ。\tau を秒で測れば p(\tau) の単位は秒の逆数になる。
p(\tau) = r\, e^{-r\tau}
導出は簡単だ。「次のスパイクまで \tau 以上待つ」確率は、区間 [0,\tau) にスパイクが一つもない確率、つまり P(0) = e^{-r\tau} である。これを \tau で微分して符号を変えれば密度が出る。
指数分布の密度は \tau = 0 で最大である。短い間隔がいちばん出やすい、ということだ。
ここで一つ、取り違えやすい点がある。密度が減っていくことは、「待つほど発火しにくくなる」という意味ではない。まだ発火していないという条件のもとでの瞬間の発火率(ハザード)を計算すると、p(\tau)/\int_\tau^\infty p(u)\,du = r で、\tau によらず一定である。密度が減るのは、そこまで待てる確率が減るからにすぎない。
実際のニューロンはこうなっていない。発火した直後は不応期のせいでハザードが下がる。だから短い間隔は指数分布ほど出ない。
変動性の指標
スパイク列のばらつきを測る指標として、変動係数(coefficient of variation)が使われる。スパイク間隔の平均を \mu_{\text{ISI}}、標準偏差を \sigma_{\text{ISI}} と書き、標準偏差を平均で割った値を C_V とする。分子と分母を同じ時間の単位で測るので、C_V は単位を持たない。平均の間隔に対して間隔がどれだけばらつくかを表す指標である。
C_V = \frac{\sigma_{\text{ISI}}}{\mu_{\text{ISI}}}
指数分布では \sigma = \mu = 1/r だから、C_V = 1 である。ポアソン過程が基準値 1 を与える。
- C_V < 1 — ポアソンより規則的(不応期がある、周期的に発火する)
- C_V > 1 — ポアソンより不規則(短い間に複数のスパイクをまとめて出す、バースト発火をする)
皮質のニューロンを記録すると、C_V はしばしば 1 に近い。皮質は驚くほど不規則に発火している。これは長らく謎とされ、「興奮性入力と抑制性入力がほぼ釣り合っていて、その揺らぎで発火している」という説明(バランス状態)が提案されている。第3章で触れる。
発火率という近似はいつ妥当か
本書の以降の章では、ほとんどの場面で発火率を使う。だがそれが妥当でない場面もあることは、意識しておきたい。
- スパイクの正確なタイミングが情報を担っている場合(音の振動の決まったタイミングで発火する、聴覚系の位相同期など)
- 発火が非常にまばらで、一回の試行に数個しかスパイクがない場合
- 同期が重要な場合(複数の細胞が揃って発火することに意味がある場合)
発火率への縮約は、こうした情報をすべて捨てている。捨ててよいかは、問いによる—第1章第2節の主題が、ここでも効いている。
7. 活性化関数はどこから来たのか
章の締めくくりとして、深層学習へ橋を架ける。
f–I 曲線と活性化関数
第2節で LIF モデルの f–I 曲線を計算した。形は「閾値までゼロ、超えたら急に立ち上がり、やがて緩やかになる」だった。
人工ニューロンの活性化関数は、この曲線の近似である。以下の u は、第1章第1節で定義した活性化関数に入れる量、つまり u = \sum_i w_i x_i + b で、u も出力も単位を持たない実数として扱う。だからシグモイドの上限 1 は、最大発火率が 1 Hz だという意味ではない—電流や発火率に対応させるには、入力と出力の尺度を別に決める必要がある。
シグモイド関数
\varphi(u) = \frac{1}{1 + e^{-u}}
滑らかに立ち上がり、上下に飽和する。飽和は、実際のニューロンが不応期のせいで最大発火率を持つことに対応している。
第2節で残した宿題を、ここで片づけておこう。LIF は電流を上げれば発火率がいくらでも上がるので、そのままでは飽和しない。だが一回の発火に、閾値まで充電する時間 T のほかに、休んでいる時間—不応期 \tau_{\text{ref}}—がかかるとすれば、
f = \frac{1}{T + \tau_{\text{ref}}}
となる。電流を上げると T \to 0 だが、\tau_{\text{ref}} は残る。だから
f \to \frac{1}{\tau_{\text{ref}}}
で頭打ちになる。\tau_{\text{ref}} = 2 ms なら上限は 500 Hz である。シグモイドの上の飽和は、ここから来ている。
ReLU(正規化線形関数)
\varphi(u) = \max(0, u)
閾値以下でゼロ、以上で線形。f–I 曲線の閾値より上の動作範囲を、直線で置き換えたものと読める(閾値のごく近くは、上で見たとおり直線にならない)。飽和を捨てているが、実際のニューロンも通常の動作範囲では飽和まで行かないことが多い。そして ReLU は、勾配が消えないという計算上の利点を持つ(第8章第6節)。生物学的な妥当性と計算上の都合が、たまたま一致した例である。
対応の限界
ただし、この対応を過大評価してはいけない。深層学習で活性化関数が選ばれる理由は、ほとんどが計算上の都合である。ReLU が普及したのは、誤差の影響をさかのぼって伝える途中で微分の値が極端に小さくなる勾配消失を避けられるからであって、生物学的に正しいからではない。GELU や Swish といった近年の関数は、生物学的な動機をまったく持たない。
逆向きの対応も怪しい。実際のニューロンの入出力関係は、入力の時間的な構造や樹状突起での非線形な統合に依存し、単一の静的な関数では書けない。それでも出所を知る価値はある。「活性化関数」という部品が、どこかから降ってきたのではなく、実在するものの単純化として生まれたことを知っていると、モデルを設計するときの判断が変わる。何を捨てているかが見えるからだ。
本書が置いている前提
ここでもう一歩踏み込んでおきたい。活性化関数だけの話ではなく、本書全体が何を前提にしているかである。
本書はこの先ずっと、情報はスパイクとして運ばれ、信号を運ぶ化学物質(神経伝達物質)を放出する接点、化学シナプスで受け渡され、発火率に要約される—という描像のうえに数学を組み立てる。これはニューロン説(neuron doctrine)と呼ばれる考え方の延長にある。神経系の機能単位は個々の細胞であり、細胞は接点を介して信号をやりとりする、という前提だ。これは確定した事実ではなく、記述の選択である。実際の脳には、この枠に収まらない伝達の仕組みがいくつも見つかっている。
- 樹状突起の演算。樹状突起は入力を素直に足し合わせる導線ではない。局所的に活動電位を起こし、枝ごとに違う非線形演算をする(London と Häusser 2005)
- ギャップ結合(gap junction)。細胞膜が直接つながってイオンが行き来する。化学シナプスを経ないので速く、集団の同期に効く
- 体積伝達(volume transmission)。伝達物質が細胞外液に拡散し、接点を持たない多数の細胞に届く。宛先が一つに決まらない伝達である
- 電場介在伝達(ephaptic transmission)。細胞外の電場そのものが、隣の細胞の興奮性を変える
- グリア伝達(gliotransmission)。ニューロンでない細胞が伝達に加わる
さらに、スパイクを使わない伝達もある。眼の奥で光を受ける薄い組織、網膜の光受容細胞(光を電気信号に変える細胞)や双極細胞(その信号を次へ中継する細胞)では、膜電位の緩やかな変化がそのまま伝達物質の放出量に反映される。「情報はスパイクの列である」という前提は、脳のすべての場所で成り立つわけではない。この一覧は Bullock ら (2005) の整理に沿っている。彼らはニューロン説を捨てよと言ったのではない。書き換えが必要だと論じたのである。
もっと大きな水準では、記述の切り方そのものを変える提案もある。皮質を離散した単位の集まりと見るのではなく、連続した媒質と見て、活動を進行波や形状の共鳴モードとして書く立場だ。これは次章の話題に直接関わるので、そこで改めて触れる(第3章第7節)。
では、なぜ本書はそれでもニューロン説の描像を採るのか。
正直に言えば、そこでしか脳とAIを同じ数学に載せられないからである。人工ニューラルネットワークは「ユニットと重み」の言葉で書かれている。脳の側を同じ言葉に翻訳して初めて、両者を並べて比べられる。本書の主題はその比較にあるから、この抽象化を採る。
ただし前提を前提として意識しているかどうかで、結論の読み方は変わる。脳とAIが「同じ数学で書ける」と確かめられたとして、それは脳をこの抽象化で切り取った限りでの話である。切り取り方を変えれば、似ている度合いも変わりうる。第1章第2節で述べた「どの水準で似ていると言っているのか」という問いは、実はもう一段手前から始まっている—どの水準で脳を記述すると決めたのか、である。
本書を読み終えたとき、あなたはこの前提のどこを外してみたいと思うだろうか。
次章へ
本章で扱ったのは一個のニューロンだった。だが脳の計算は、明らかに集団の性質である。
次章では、多数のニューロンの集団を一つの変数(集団発火率)に縮約し、その相互作用を微分方程式で書く。本章で使った道具—微分方程式、相平面、固定点、分岐—が、そのまま集団のレベルで使える。記述の水準が変わっても数学は変わらない、という本書の主題の最初の実例になる。
確認問題
[導出]\tau\,dx/dt = -x + I(I は定数)を解き、t = \tau で初期値と最終値の差が何倍になるかを示せ。(第1節)
[導出]HH 方程式のゲート変数の式 dx/dt = \alpha_x(1-x) - \beta_x x を \tau_x\,dx/dt = -x + x_\infty の形に変形し、\tau_x と x_\infty を \alpha_x, \beta_x で表せ。(第3節)
[確認]\mathrm{Na}^+ 電流の項が m^3h という積を含むことの意味を述べよ。m と h の速さの違いが、活動電位のどの局面に対応するか。(第3節)
[考える]HH 方程式には閾値が明示的に書かれていない。それにもかかわらず閾値のような振る舞いが現れるのはなぜか。(第3節・第5節)
[確認]2変数系の固定点の安定性が、ヤコビ行列の固有値でどう判定されるかを述べよ。(第4節)
[確認]SNIC(不変円上のサドルノード分岐)と Hopf 分岐の違いを、f–I 曲線の形の違いとして説明せよ。一般のサドルノード分岐では周期発火が保証されないのはなぜか。(第5節)
[導出]ポアソン過程において、スパイク間隔が指数分布に従うことを導出せよ。(第6節)
[考える]シグモイド関数と ReLU が、f–I 曲線のどの性質を捉え、どの性質を捨てているかを述べよ。(第7節)
参考文献
- 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 — 本章第3節の原典。読みやすい論文なので一度は当たってほしい
- Gerstner, W., Kistler, W. M., Naud, R., & Paninski, L. (2014). Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. Cambridge University Press. — 単一ニューロンモデルの標準的な教科書。無料で公開されている
- 宮川博義・井上雅司(2013)『ニューロンの生物物理』丸善出版. — 生物物理の側から。膜とイオンチャネルの記述が丁寧[3節・コラム C2]
- Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT Press. — 分岐による発火の分類を徹底的に扱う[4〜5節]
- Bullock, T. H., et al. (2005). The neuron doctrine, redux. Science, 310(5749), 791–793. https://doi.org/10.1126/science.1114394 — 本書が置いている前提を外から眺めるために。ギャップ結合・体積伝達・電場介在伝達など、化学シナプス以外の伝達を一望できる短い論文[7節]
- London, M., & Häusser, M. (2005). Dendritic computation. Annual Review of Neuroscience, 28(1), 503–532. https://doi.org/10.1146/annurev.neuro.28.061604.135703 — 樹状突起を「導線」ではなく演算装置として見る[7節・コラム C2]
- Hodgkin, A. L. (1948). The local electric changes associated with repetitive action in a non-medullated axon. The Journal of Physiology, 107(2), 165–181. https://doi.org/10.1113/jphysiol.1948.sp004260 — Type I/Type II の区別の原典[5節]