SPINE is a simple Spiking Neuron simulator

version python python

Author: Hiroshi ARAKI @Hiro Lab.
  araki@hirlab.net (@を半角に直してください)
@hirlab_a
HiroshiARAKI/spine: SPINE is a simple Spiking Neuron simulator

目次 Contents

コンセプト

spine


SPINEは一言でいうと,超シンプルなニューロンシミュレータです. Pythonで動作し,必要なパッケージはnumpyとmatplotlibのみ.

もともとは個人使用で,勉強用で作り始めたものなので,クオリティは求めてはいけません

また,日本産で必要最低限の機能しかないので,簡単な実験にぴったりです. 神経科学,スパイキングニューロンなどの研究の入門用として使用ください.

ちなみにキャッチコピーは,

SPINE is a simple Spiking Neuron simulator.


keywords: spiking, neuron, simulator, python, Hodgkin, Huxley, FitzHugh, Nagumo, LIF, スパイキング, ニューロン, シミュレータ, シミュレーション


使い方

Pythonが使える前提での解説になります.ご了承ください.

1. Githubからインストール

https://github.com/HiroshiARAKI/spine からzip形式でダウンロード及び解凍(展開)するか, 以下のコマンドでパッケージを取得しインストールします.

            $ wget https://github.com/HiroshiARAKI/spine/tarball/master
            $ tar xvf master
        

gitが使える場合はcloneしても良いでしょう.

            $ git clone https://github.com/HiroshiARAKI/spine.git
        

もし,パッケージとしてpipで管理したい場合や,いろんなファイルからSPINEを使いたい場合は

            $ python setup.py install
        

でインストール可能です.実際に確認してみます.

            $ pip list
            Package                       Version
            ----------------------------- -----------
            spine                         2.x
        

いい感じです.


2. 使用してみる

インストールが終わればすぐにシミュレーションが可能です.

試しに,.pyファイルを作成して以下のコードを書いてみましょう.

            from spine import HodgkinHuxley
            import numpy as np


            if __name__ == '__main__':

                # Hodgkin-Huxleyニューロンを生成
                neu = HodgkinHuxley(time=100,  # 実験時間 [ms]
                                    dt=0.01)   # タイムステップ [ms]

                # 入力データ (今回はノイズを載せた矩形正弦波)
                input_data = np.sin(0.5 * np.arange(0, neu.time, neu.dt))
                input_data = np.where(input_data > 0, 20, -5) + 10 * np.random.rand(int(neu.time / neu.dt))

                # 入力データを与えて膜電位を計算
                neu.calc_v(input_data)

                # 実験結果を描画
                neu.plot_monitor(save=True, filename='sample.png')
        

Hodgkin-Huxley

簡単にHodgkin-Huxleyニューロンのシミュレーションができました.


3. 実験条件を変更してみる

それでは少し条件を変えてみましょう.

                from spine import HodgkinHuxley
                import numpy as np


                if __name__ == '__main__':

                    neu = HodgkinHuxley(time=100,
                                        dt=0.01,
                                        n=0.4,  # コンダクタンスのパラメータを変えてみる
                                        m=0.05,
                                        h=0.4)

                    # 入力も少し変えてみる
                    input_data = np.sin(0.2 * np.arange(0, neu.time, neu.dt))
                    input_data = np.where(input_data > 0, 10, -5) + 10 * np.random.rand(int(neu.time / neu.dt))

                    neu.calc_v(input_data)

                    neu.plot_monitor(save=True, filename='sample2.png')
            

Hodgkin-Huxley

出力結果が変わりました.このようにニューロン生成時にパラメータを細かく設定することが可能です.


4. その他のモデルやツール

Hodgkin-Huxleyモデル意外にも,LIFモデルやFitzHugh-Nagumoなどのモデルも使用可能です.

詳しくはDocumentationを参照

                # インポートできるもの
                
                from spine import (
                    HodgkinHuxley,
                    FitzHughNagumo,
                    LIF,
                    IF,
                    PoissonSpike
                )
                # kernel['single] や kernel['double']のように呼び出せる
                from spine.tools import kernel
                from spine.tools.plotting import plot_spike_scatter

            

5. 各クラスのメンバ変数について

各ニューロンクラスには,monitorと言う名の辞書型変数が存在します.

この変数には膜電位を計算した後に,その実験結果の情報が格納されています.

もし,自分独自でグラフの描画を細かく設定したい場合はこのmonitorを使用してください.

                from spine import HodgkinHuxley
                import numpy as np
                import matplotlib.pyplot as plt


                if __name__ == '__main__':

                    neu = HodgkinHuxley(time=100, dt=0.01)

                    input_data = np.sin(0.2 * np.arange(0, neu.time, neu.dt))
                    input_data = np.where(input_data > 0, 10, -5) + 10 * np.random.rand(int(neu.time / neu.dt))

                    neu.calc_v(input_data)

                    # monitorから膜電位だけ取得
                    v = neu.monitor['v']

                    plt.plot(v)
                    plt.show()
            

Hodgkin-Huxley

このように,膜電位だけ取得し扱うことができます.
monitor['v'] は膜電位 (FitzHugh-Nagumoモデルの場合はrecovery voltage),
monitor['n'], monitor['m'], monitor['h'] はコンダクタンスパラメータ(HHモデル限定),
monitor['s'] は出力スパイク列 (LIFモデル限定),
monitor['u'] はFitzHugh-Nagumoモデルの場合の膜電位
です.

その他の機能についてはコードを直接みていただいた方が良いかもしれません.

もっと機能が増えたら,マニュアルを作成します.


6. 超簡単に実験したいとき

コードもなにも書きたくない場合,すでにお試し実験コードが用意されています.

例えばFitzHugh-Nagumoモデルであれば,

                $ python fhn_sample.py
            

FitzHugh-Nagumo

でOKです.fhn.pyの他にhh_sample.py, lif_sample.py, dlif_sample.pyにも同じようにサンプルコードが書かれています.


スポンサードリンク


簡易ドキュメント

Neuron (parent class)

※これ単体で使うことはない.

                class Neuron(object):
                    def __init__(self, time: int, dt: float):
                        """
                        Neuron class constructor
                        :param time:
                        :param dt:
                        """
                        self.time = time
                        self.dt = dt
            

以下共通のメンバ変数

メンバ変数名概要
time観測時間 [ms]
dtタイムステップ [ms]

LIF

よくあるニューロンモデル.一番使い勝手が良いかもしれない.

本LIFモデルは,入出力をスパイク列にしているためエンジニアリングの観点からも扱いやすい.

スパイクから電流への変換ではカーネルを選べる.

$$\tau_{m}\frac{dV(t)}{dt}=(-V(t) + E_{rest}) + I(t)$$ $$I(t) = \sum_{i \in N} w_i (S_i(t) * K(t))$$ $$ \begin{eqnarray} K(t) = \begin{cases} \exp(-t/\tau) & (\mathrm{single \ exponential \ filter})\\ \exp(-t/\tau_1) - \exp(-t/\tau_2) & (\mathrm{double \ exponential \ filter}) \end{cases} \end{eqnarray} $$

                class LIF(Neuron):
                    """
                    LIF: leaky integrate-and-fire model
                    """

                    def __init__(self,
                                 time: int,
                                 dt: float = 1.0,
                                 rest=-65,
                                 th=-40,
                                 ref=3,
                                 tc_decay=100,
                                 k='single',
                                 tau: tuple = (20, ),
                                 **kwargs):
                        """
                        Initialize Neuron parameters
                        :param time:      experimental time
                        :param dt:        time step
                        :param rest:      resting potential
                        :param th:        threshold
                        :param ref:       refractory period
                        :param tc_decay:  time constance
                        :param k:         kernel {'single', 'double'}
                        :param tau:       exponential decays as tuple(tau_1 ,tau_2) or float
                        """
            
メンバ変数名概要
rest静止膜電位 $E_{rest}$ [mV]
th発火閾値 [mV]
ref不応期 [ms]
tc_decay減衰時定数 $\tau_{m}$ [ms]
kスパイク応答カーネル ({'single', 'double'}から指定.どちらでもない場合'single'が選択される)
tauカーネルの時定数 [ms] ('single'ならtuple(tau, )'double'ならtuple(tau_1, tau_2):デフォルト値はそれぞれ(20, )(5, 2.5))
monitor記録用辞書変数{'s': [{0,1} spike train], 'v': [membrane voltage], 'f': [firing times]}

膜電位計算

LIF.calc_v(data: Tuple(spikes[N][int(time/dt)], weights[N]))


LIF (double exponential filter)

漏れ(リーク)がない,ただの積分発火モデル.

需要があるのだろうか.LIFより後に作成したので,なぜかLIFが親クラスとなっているちょっとよくないコーディング.

動作にはなんら問題がないので,別にいいけどいつか直す.

$$V(t) = \sum_{i\in Pre} w_i (S_i(t) * K(t)) + E_{rest}$$ $$ \begin{eqnarray} K(t) = \begin{cases} \exp(-t/\tau) & (\mathrm{single \ exponential \ filter})\\ \exp(-t/\tau_1) - \exp(-t/\tau_2) & (\mathrm{double \ exponential \ filter}) \end{cases} \end{eqnarray} $$

                class IF(LIF):
                    """
                    IF: integrate-and-fire model
                    """
                    def __init__(self,
                                 time: int,
                                 dt: float = 1.0,
                                 rest=-65,
                                 th=-40,
                                 ref=3,
                                 k='single',
                                 tau: tuple = (20, ),
                                 **kwargs):
                        """
                        Initialize Neuron parameters
                        :param time:
                        :param dt:
                        :param rest:
                        :param th:
                        :param ref:
                        :param k:
                        :param tau:
                        :param kwargs:
                        """
            
メンバ変数名概要
LIFと同じLIFと同じ

膜電位計算

IF.calc_v(data: Tuple(spikes[N][int(time/dt)], weights[N]))


Hodgkin-Huxley model

一番ニューロンのモデル化を真面目にやったもの.

エンジニアリングの観点から見たら扱いづらいモデル.入力は電流としている点に注意.

$$ 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)$$ $$ \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)$$

※ $(t)$ は省略しています

            class HodgkinHuxley(Neuron):
                """
                Hodgkin-Huxley neuron model
                """

                def __init__(self, time: int,
                             dt: float = 0.01,
                             rest=-65.,
                             Cm=1.0,
                             gNa=120.,
                             gK=36.,
                             gl=0.3,
                             ENa=50.,
                             EK=-77.,
                             El=-54.387,
                             n=0.32,
                             m=0.05,
                             h=0.6,
                             **kwargs):
                    """
                    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
                    :param n:    Conductance param
                    :param m:    Conductance param
                    :param h:    Conductance param
                    """
            
メンバ変数名概要
rest静止膜電位 $E_{rest}$ [mV]
Cm 膜容量 $C_m$ [μF/cm2]
gNa $Na^{+}$チャネルコンダクタンス $g_{Na}$ [mS/cm2]
gK $K^{+}$チャネルコンダクタンス $g_{K}$ [mS/cm2]/td>
gl リークコンダクタンス $g_{l}$ [mS/cm2]
ENa ナトリウム平衡電位 $E_{Na}$ [mV]
EK カリウム平衡電位 $E_{K}$ [mV]
El リーク平衡電位 $E_{l}$ [mV]
n カリウム活性化変数 ($0 \leq n \leq 1$)
m ナトリウム活性化変数 ($0 \leq n \leq 1$)
h ナトリウム不活性化変数 ($0 \leq n \leq 1$)
monitor記録用辞書変数{'v': [membrane voltage], 'n': [n], 'm': [m], 'h': [h]}

膜電位計算

HodgkinHuxley.calc_v(data: list(input current))


FitzHugh-Nagumo model

正直使いどころがわからないが,シンプルなHodgkin-Huxleyモデルと言える.

とにかく定数や変数が少ないのが魅力的だが,これは膜電位の挙動を見ているだけで,膜電位そのものとは言えないと思う.

これもエンジニアリングの視点からでは正直扱いづらい上に,定数設定がシビア.

$v$ではなく$u$が膜電位にあたるので注意

$$\frac{du(t)}{dt}=c(u(t)-\frac{u(t)^3}{3}-v(t)+I(t))$$ $$\frac{dv(t)}{dt}=u(t)-bv(t)+a$$

                class FitzHughNagumo(Neuron):
                    """
                    FitzHugh-Nagumo neuron model
                    """

                    def __init__(self,
                                 time: int,
                                 dt: float = 0.1,
                                 a=0.7,
                                 b=0.8,
                                 c=10,
                                 **kwargs):
                        """
                        Initialize Neuron parameters
                        :param time:
                        :param dt:
                        :param a:
                        :param b:
                        :param c:
                        """
            
メンバ変数名概要
a定数その1
b定数その2
c定数その3
monitor記録用辞書変数{'u': [membrane voltage], 'v': [recovery variable]}

膜電位計算

FitzHughNagumo.calc_v(data: list(input current))


Poisson Spike Generator

一般的にニューロンの自発的発火周期はポアソン分布に従うと言われている.エンジニアリングでもよく,画像->スパイク列などの前処理として使われる.

入力は画像などの[0,1]実数値データ.これを周波数に変換してスパイク列に変える.

詳しくは,Pythonで実装しながら理解するPoisson Spike (ポアソンスパイク) - Qiita

                class PoissonSpike:
                    def __init__(self, data, time=500, dt=0.1, max_freq=128):
                        """
                        Generate Poisson spike trains
                        :param data:
                        :param time:
                        :param dt:
                        :param max_freq:
                        """
            
メンバ変数名概要
max_freq正規化に使う最大周波数 [Hz] ([0, max_freq]で正規化を行う)
freq_data正規化後の周波数データ [Hz]
norm_data正規化後の周期データ [ms]
spikes出来上がった{0,1}スパイク列 (基本的にはPoissonSpike(data).spikesで事足りる)
fires出来上がった発火時刻リスト
monitor記録用辞書変数 {'s': [{0,1} spike trains], 'f': [firing times]}

Plottong

プロットに関する関数も随時整備していきます.現在は{0,1}スパイク列から,ラスター図(散布図: scatter)でのスパイク描画関数が使えます.

                from spine import PoissonSpike
                from spine.tools.plotting import plot_spike_scatter

                import numpy as np
                import matplotlib.pyplot as plt


                spikes = PoissonSpike(np.random.random((10, 10)), time=300, dt=0.1).spikes
                plot_spike_scatter(spikes, time=300, dt=0.1)
                plt.title('Spike firing timing')
                plt.show()
            


さいごに

このライブラリは正直言って大したものではないので,実際にこんなに専用ページまで作るほどのものではありません.

ただ,作成者である私の自己満で趣味な部分が大きいので,使っていて不自由に感じても文句は言わないでください.笑

どこかで使用・再利用・改変する場合には,特に断りなく使ってください.自由です.

もし,要望等あれば,連絡いただければ対応するかもしれません.

とにかく気楽にお使いください.





Hiroshi ARAKI

とある国立大学院生.
情報工学を専攻し,現在はスパイキングニューラルネットワークについての研究を行っている.
もちろん普通のニューラルネットワークや機械学習にも精通している.というか神経科学よりそっちが専門.
趣味は映画見たり,ドライブしたり,美味しいもの食べたりすること.

    araki@hirlab.net (@を半角に直してください)





SPINE is a simple Spiking Neuron simulator