Brian2で神経活動シミュレーション


神経活動をコンピュータ上でシミュレーションする研究は以前から行われています。
シミュレータとしてはイェール大学の[NEURON](https://neuron.yale.edu/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でも出来るようです。

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

次項から[ドキュメントのチュートリアル](https://brian2.readthedocs.io/en/stable/resources/tutorials/index.html)を参考にしながら基本的なところをやってみます。

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

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

```Python
import brian2
```

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

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

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

```Python
brian2.start_scope()
```

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

```Python
eqs = 'dv/dt = (v0-v)/tau : volt (unless refractory)'
```
上記のvoltは変数の物理単位です。SI単位で表します。
ここでは膜電位が変数になります。
(unless refractory)は不応期を使う場合に必要です。

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

```Python
G = brian2.NeuronGroup(1, eqs, threshold='v > -40*mV', reset='v = -60*mV', refractory=5*brian2.ms, method='exact')
```
第一引数は神経細胞の数です。今回は１個で設定しました。
第二引数は膜電位の微分方程式です。
他の引数は
thresholdはスパイクになる閾値電位
resetはスパイク後のリセット電位
refractoryは不応期
methodは数値計算の方法
を設定しています。

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

```Python
G.v = -70 * brian2.mV
```

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

```Python
M = brian2.StateMonitor(G, 'v', record=True)
```
シミュレーションは`run(<時間>)`によって開始されます。

```Python
brian2.run(50*brian2.ms)
```
これでシミュレーションできました。

結果をプロットします。

```Python
brian2.plot(M.t/brian2.ms, M.v[0]/brian2.mV)
brian2.xlabel('Time (ms)')
brian2.ylabel('Ptential (mV)');
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/521288/9d52ce80-47fc-f56f-624b-603d8ce710f8.png)



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

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

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

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

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

```Python
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`を使います。

```Python
M = SpikeMonitor(G)
```
各細胞に異なる入力を定義します。

```Python
I_max = 40*mV
G.I = 'i*I_max/(N-1)'
```
iは細胞のインデックスになります。これで異なる入力を定義できました。

runします。

```Python
duration = 1000*ms
run(duration)
```

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

```Python
plot(M.t/ms, M.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index');
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/521288/a09e72ce-65bf-a4db-0adb-9b2ecb6bfe56.png)

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

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

```Python
plot(G.I/mV, M.count/duration)
xlabel('Input (mV)')
ylabel('Firing rate (sp/s)');
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/521288/c3065cc6-00a4-c685-de21-ec731edafc02.png)
出来ました。
これが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/