神経活動をコンピュータ上でシミュレーションする研究は以前から行われています。
シミュレータとしてはイェール大学のNEURONが有名ですが、他にも様々なソフトがあります。

今回はその中でも、去年(2019年)にeLifeで論文発表されたBrian2(バージョン1とは異なる)を使ってみます。
＊間違えやすいですが、名前は**Brain（脳）ではなくてBrian（ブライアン）**ですね。

ホームページ：https://briansimulator.org/
GitHub：https://github.com/brian-team/brian2
論文：https://elifesciences.org/articles/47314
開発元はフランスのソルボンヌ大学などです。

インストール

環境は
Ubuntu 16.04
Anaconda
Python 3.6
です。

インストールはcondaでもpipでも出来るようです。

pip install brian2

これでインストールできました。
バージョンはbrian2-2.3.0.2でした。

次項からドキュメントのチュートリアルを参考にしながら基本的なところをやってみます。

動かしてみる

まずは単純な膜電位シミュレーションをやってみます。

インポートします。
(この項ではどれがモジュールかを分かりやすくするために、あえてマニュアル通りのfrom brian2 import *ではなく普通にインポートします。)

import brian2

Brian2には物理単位がモジュールになっていて、Hzやmsなどをそのまま定義することができます。
これを使うことで、パラメータや定数を簡潔に記述できます。
また、計算時に物理単位が揃わないとエラーになるようになっているので設定ミスを防ぐことができます。

v0 =  -10 * brian2.mV
tau = 10 * brian2.ms

シミュレーションを始めるときにはまずstart_scope()を入れます。

brian2.start_scope()

微分方程式（膜電位の変化を定義する）はString型で定義するようです。

eqs = 'dv/dt = (v0-v)/tau : volt (unless refractory)'

上記のvoltは変数の物理単位です。SI単位で表します。
ここでは膜電位が変数になります。
(unless refractory)は不応期を使う場合に必要です。

神経細胞はNeuronGroupモジュールで定義します。

G = brian2.NeuronGroup(1, eqs, threshold='v > -40*mV', reset='v = -60*mV', refractory=5*brian2.ms, method='exact')

第一引数は神経細胞の数です。今回は１個で設定しました。
第二引数は膜電位の微分方程式です。
他の引数は
thresholdはスパイクになる閾値電位
resetはスパイク後のリセット電位
refractoryは不応期
methodは数値計算の方法
を設定しています。

膜電位の初期値を設定します。

G.v = -70 * brian2.mV

StateMonitorによってモニタリングができます。

M = brian2.StateMonitor(G, 'v', record=True)

シミュレーションはrun(<時間>)によって開始されます。

brian2.run(50*brian2.ms)

これでシミュレーションできました。

結果をプロットします。

brian2.plot(M.t/brian2.ms, M.v[0]/brian2.mV)
brian2.xlabel('Time (ms)')
brian2.ylabel('Ptential (mV)');

この結果はLIFモデルのニューロンのスパイクをシミュレーションした結果になります。(ピーク電位がありませんが。)

これでシミュレーションの流れが分かりました。

発展として、xi(標準正規分布の確率変数)を使うことで、ノイズを含んだ確率微分方程式を扱うこともできるようです。

複数の神経細胞のスパイクシミュレーション

次に複数の神経細胞のスパイクをシミュレーションします。
今回は各細胞に異なる入力電位（I）を与えてみます。

インポートからやっていきます。

from brian2 import *

start_scope()

N = 100 #細胞の数
tau = 10*ms #膜時定数
v0 = -50*mV #静止膜電位

eqs = '''
dv/dt = (v0-v+I)/tau : volt (unless refractory)
I : volt
'''

G = NeuronGroup(N, eqs, threshold='v>-40*mV', reset='v=-60*mV', refractory=1*ms, method='exact')
G.v = -70*mV

膜電位変化ではなくスパイクのみを見たいときはSpikeMonitorを使います。

M = SpikeMonitor(G)

各細胞に異なる入力を定義します。

I_max = 40*mV
G.I = 'i*I_max/(N-1)'

iは細胞のインデックスになります。これで異なる入力を定義できました。

runします。

duration = 1000*ms
run(duration)

ラスタープロットを作成します。

plot(M.t/ms, M.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index');

細胞のインデックスが上がるごとに入力電位が上がるので、スパイクが多くなりました。

よく見るF-Iカーブも作ってみます。

plot(G.I/mV, M.count/duration)
xlabel('Input (mV)')
ylabel('Firing rate (sp/s)');

出来ました。
これがReLU関数に似ていると言われますね。

神経回路をシミュレーションするには

詳細は割愛しますが、神経回路のシミュレーションも簡単に出来るようです。

シナプスはSynapsesモジュールで定義できます。
connectメソッドで細胞の結合が行われます。（どの細胞からどの細胞へつなげるか）

シナプスの強度に微分方程式を定義することもできて、STDPモデルも作れるようです。

他にもNetworkモジュールでネットワークの定義ができたりするようです。

感想

シミュレーションのための機能が充実していて、使いやすいと思いました。

ドキュメントには他にも様々な例が載っています。
https://brian2.readthedocs.io/en/stable/examples/index.html
コンパートメントモデルの例もありました。
これらも試してみたいですね。

＊追記①
今年の1月にBrian2GeNNというものが発表されていました。
https://www.nature.com/articles/s41598-019-54957-7
これでGPUも活用できるようになったみたいです。

＊追記②
またbrian2toolsというライブラリもあり、これを使うと可視化が簡単に出来るようになって便利です。
https://brian2tools.readthedocs.io/en/stable/