スパイキングニューロンのシミュレータを触ってみる記事の２回目。
以前の投稿の続きとして、サンプルにはないものを書いてみた。
izhikevichモデルニューロンの集団発火
サンプルではLIFモデルだったが、izhikevichモデルで同じものを作ってみる。
プログラムはこのサイトのものをほぼパクらせて参考にさせていただいた。この場を借りて感謝したい。Thank you！！！
Google グループ
Google グループ
実際のコード
from brian2 import *
def izhikevich_model():
def draw_normals(n,start,stop):
mu,sigma,numbers = start+(stop-start)/2, (stop-start)/6, zeros(n)
for i in range(n):
s = -1
while (s<start) or (s>stop) :
s = numpy.random.normal(mu,sigma,1)
numbers[i]=s
return numbers
n = 2000
R = 0.8
eqs = Equations('''
dv/dt = (0.04/ms/mV)*v**2 + (5/ms) * v + 140*mV/ms - u + I_syn/ms + I_in/ms : volt
du/dt = a*((b*v) - u) : volt/second
dx/dt = -x/(1*ms) : 1
I_in = ceil(x)*(x>(1/exp(1)))*amplitude : volt
dI_syn/dt = - I_syn/tau : volt
a : 1/second
b : 1/second
c : volt
d : volt/second
amplitude : volt
tau : second
''')
reset = '''
v = c
u += d
'''
P = NeuronGroup(n, model=eqs, threshold='v>30*mvolt', reset=reset, method='euler')
Pe = P[:int(n*R)]
Pi = P[int(n*R):]
Pe.a = 0.02/msecond
Pe.b = 0.2/msecond
Pe.c = (15*draw_normals(int(n*R),float(0),1) - 65) * mvolt
Pe.d = (-6*draw_normals(int(n*R),float(0),1) + 8) * mvolt/msecond
Pe.tau = draw_normals(int(n*R),float(3),15) * msecond
Pi.a = (0.08*draw_normals(n-int(n*R),float(0),1) + 0.02) * 1/msecond
Pi.b = (-0.05*draw_normals(n-int(n*R),float(0),1) + 0.25) * 1/msecond
Pi.c = -65 * mvolt
Pi.d = 2 * mvolt/msecond
Pi.tau = draw_normals(n-int(n*R),float(3),15) * msecond
P.x = 0
P.v = draw_normals(n,float(-50),float(-25)) * mvolt
P.amplitude = draw_normals(n,0,8) * mvolt
Ce = Synapses(Pe, P, on_pre='I_syn+=1.5*mV')
Ce.connect(p=0.05)
Ci = Synapses(Pi, P, on_pre='I_syn-=8*mV')
Ci.connect(p=0.05)
M = SpikeMonitor(P)
V = StateMonitor(P,'v',record=True)
run(500*ms)
subplot(211)
plot(M.t/ms, M.i, '.')
subplot(212)
plot(V.t[:200]/ms, V.v[0][:200]/mV, 'r')
plot(V.t[:200]/ms, V.v[100][:200]/mV, 'g')
plot(V.t[:200]/ms, V.v[200][:200]/mV, 'b')
show()
if __name__=="__main__":
izhikevich_model()
実行結果
以前の記事のLIFモデルのように発火が続かないがパラメータによるところが大きいから気にしなくてもいいと思う。
下の図は適当にとってきたニューロンの膜電位。ラスタープロットだとめっちゃ発火してるように見えるけど、実際はそんな発火してないよねってことを確認したかった。
'''
dx/dt = -x/(1*ms) : 1
I_in = ceil(x)*(x>(1/exp(1)))*amplitude : volt
'''
の部分は外部刺激に対応する項だけど今回は使わなかった。
I_synとI_inをなんでこんな形で書き分けるのか謎。全部I_synで良くない？誰か教えて欲しい。
感想
brianは書き方が独特すぎてしんどい。（自由度が高いとしょうがないよね。）
サンプルを見ながら試行錯誤するか、最悪ソースを直接見に行くかしないと解決しないことがあってつらいね。
変数に単位があるっていうのは嫌いじゃないんだけどもっと良い方法はなかったのかね。1*second って…
それと、サンプルがだいたい
from brian2 import *
の形で書かれてて、いつか変数名で問題起こしそうでひやひやする…
関連記事：
pythonでスパイキングニューロンのシミュレーション (Brianについて） - 技術メモ
に対して、
(2)互いに素な自然数
に対して、
に対して
について、
とする。 自然数
が
ならば、任意の自然数
に対して、
が成り立つ。
と互いに素でない時、明らかである
とも互いに素であり、
かつ
を示せば十分である。 まず、
(kは非負整数)と表せる。 したがって、
フェルマーの小定理から以下が成り立つ
よって
