神経活動を微分方程式で表してみよう


 脳には様々なタイプの神経細胞が存在しています。タイプによって細胞体や神経突起の形が異なったり、発現しているたんぱく質(マーカー蛋白)が異なったりします。そして何より重要なことに、発火のパターンが異なります。
 神経細胞の挙動を模倣するモデルを設計するとき、できる限り簡単な数式で、かつ、できる限り多彩な挙動を示すモデルが望まれます。2003年にユージン・イジケヴィッチ(Eugene M. Izhikevich)が提唱したモデルは、計算負荷の軽減と挙動の多様性において、理想に近い数理モデルです(欠点は数式に電気生理学的な意味を見出すのが困難(?)なことでしょうか)。ここではそのモデルを取り上げたいと思います。



 ここで使われる数式は4つの変数(v, t, I, u)を含む簡単な常微分方程式(ODE)です。二つの式で成り立っていますので、いわゆる連立線形常微分方程式と呼ばれるタイプになります。

一つ目の数式は

  dv/dt = 0.04 v2 + 5 v + 140 - u + I

です。v の2次の項まで含まれています。v は膜電位を表しています。単位は mV です。t は時間を表します。単位は ms が妥当でしょう。

I は神経細胞に流れ込むシナプス電流で、これは神経細胞を興奮させたり抑制させたりする駆動力になります。

 しかし、膜電位に変化を引き起こす駆動力はそれだけでありません。内因的な膜特性も発火挙動に大きな影響を与えます。それを表す変数が u です。
 上の数式を見てもわかるように、u は膜電位を安定化する方向に働きます(“- u”として負の方向に寄与している)。つまり、K+チャネルの活性化やNa+チャネルの不活性化を表す回復変数になります。u は次の式で表されます。これが二つ目の数式です。

  du/dt = a ( b v - u )

u の時間微分もまた v u によって表現されています。ここに出てきた ab はパラメータです。単位はありません。

a u の時間スケールを決定します。典型的なケースでは a = 0.02 ですが、一般に a が小さいほど膜電位の戻りは小さくなりますので、時間差のあるシナプス入力を加重する傾向が強くなります。

buv への感受性を表します。典型的なケースでは b = 0.2 です。b が大きくなると vu の連関が高まりますので膜電位が振動しやすくなり、LTS細胞などに見られる過分極後のリバウンド発火など、独特の現象が現れます(以下参照)。

a > ba < b の境目は相転移が生じる臨界点なのですが、初心者には難しいのでここでは省略します。



 以上が微分方程式の説明ですが、それだけでは神経細胞の挙動を説明したことにはなりません。

 この数式の場合、膜電位が高くなった(つまり、神経細胞が発火した)あとのリセットを考えなければいけません。
 ここでは v が +30 mV に到達したら、それを“発火した”と解釈して、その直後に v を低い値に引き戻すことにしましょう。つまり、膜電位 vc という値を代入するのです。典型的なケースでは c = - 65 mV です。

 発火したときに起こる現象は、しかし、それだけではありません。
 発火するとA型K+チャネルが開口して膜電位を戻す力が強まります。つまり回復変数 u が一時的に増加するのです。これは、u + d のように値を加算してやることで表現できます。d が大きいほど回復の影響が尾を引き、次の発火までの時間(不応期)が長くなります。典型的なケースでは d = 2 です。

以上の発火のリセットを補式で書くと

  If v > 30 mV, then v ← c, u ← u + d

となります。

 つまり、この数理モデルは、(a, b, c, d)という4つのパラメータ制限のもとに動く4つの変数(v, t, I, u)で成り立っているわけです。

注: 30 mVという閾値が、いわゆる発火閾値でなく、生じた発火に対するリセット点であるところは、他のモデル(たとえばintegrate-and-fire model)と異なる大きな特徴です。



 これで数式の説明はおしまいです。とても簡単な数式なので驚いた人も多いと思います。この程度の数式できちんと神経細胞の発火が表現できるのか疑いたくもなります。ところが、実際に試してみると、様々な神経細胞の発火パターンをスマートに模倣できることがわかります。
 ここでは、実際に体験するためのアプリケーションを用意しましたので、皆さんもダウンロードして試してみましょう。
ここをクリックしてダウンロード (download)

なお、このプログラムはVisual Basicで書かれていますので、起動するためにはVB6ランタイムが必要です。ランタイムはこちら1こちら2などから無料でダウンロードできますので、必ずインストールしてください。


プログラムを起動すると下のような画面が現れます。各説明は以下の通りです。





 早速、実験してみましょう。まず「DC」と書かれたボタンを押してみましょう。これはモデル細胞に定電流を入力するものです。ここでは I = 10 で一定です。



下の窓に入力した電流 I が、上の窓にはこの入力に応じて変化した膜電位 v の様子が表示されています。この例では I がステップした後に7回の発火が確認できます。これは大脳皮質の錐体細胞の発火パターンを模倣したもので、最初の発火間隔が狭く、それ以降は間隔が広くなるのが特徴です(アダプテーションといいます)。


 つぎに、隣の「Pulse」と書かれたボタンを押してみましょう。これは短い矩形波を連続して入力するときの神経細胞の反応になります。



上の画像のように、一つ一つの入力が弱くても、連続で与えられたときは脱分極が加重されて(この場合では3回目の入力で)発火することがわかります。このように「Pulse」ボタンは細胞の入力統合能を調べるために使います。



 さて、それではいよいよパラメータをいじってみましょう。まず手始めに、c d の二つの変数を、c = 50d = 2 のように変えてみましょう。つまり、発火後の膜電位の回復力を弱めてみるわけです。「DC」ボタンを押すと次のように表示されます。



短い時間内に集積した多重発火(バースト)が定期的に観察されます。これはチャッタリング細胞(chattering cell)と呼ばれる神経細胞の発火パターンです。

もう一つ試してみましょう。
数値を元の値に戻すために右下の「Reset」ボタンを押してください。そして今度は b = 0.25d = 2 と変更します。「HyperP」と書かれたボタンを押してください。これは外向きの電流 I、つまり抑制方向の電流を与える実験です。結果はこうなります。



面白いことに、抑制が強いと、その反動(リバウンド)で神経細胞が発火することがわかります。これは抑制性介在ニューロンの一種であるlow threshold spiking cellに典型的な性質です。

 これまで簡単に見てきましたように、神経細胞の発火パターンは、簡単な数式(ここでは連立線形常微分方程式)で表すことができます。皆さんもいろいろなパラメータの組み合わせを試して、数理モデルの潜在性を確認してみましょう。



参考事項

1. このモデルを使って海馬回路をシミュレートする研究を、現在、個人的な趣味で
   進行させています。
2. c = - 55d = 4 でintrinsically bursting cellが模倣されます。
3. a = 0.1d = 2 でfast spiking cellが模倣されます。
4. a = 0.2b = 2c = - 56d = - 16I = - 99 でカオス(chaos)が出現しますが、
   このアプリケーションでは時間分解能が足りず観察できません。


参考文献


Izhikevich EM (2003) Simple Model of Spiking Neurons. IEEE Transactions on Neural Networks, 14:1569-1572.

Izhikevich EM (2004) Which Model to Use for Cortical Spiking Neurons? IEEE Transactions on Neural Networks, 15:1063-1070.

 上のページに戻る