FitzHugh-Nagumo modelで遊びたかった

　これを読んでいたら後半でモデリングの話が出現。その中でもFitzHugh-Nagumo modelは扱う式が少ないので実装できる気がしたのでやってみた。
　ググって出てきたこちらを参考に、まずは振動子が1つの場合をば。

---------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt

a=0.7
b=0.8
c=10
I=0.33
T=10
dt=0.01

def d_u(u, v):
    return c*(-v+u-u**3/3+I)

def d_v(u, v):
    return u-b*v+a

u=[1,]
v=[1,]

t=0
t+=dt
while t<T:
    du=d_u(u[-1], v[-1])
    dv=d_v(u[-1], v[-1])
    u.append(u[-1]+dt*du)
    v.append(v[-1]+dt*dv)
    t+=dt

plt.plot(np.arange(0, len(u), 1)*dt, u, c='b')
plt.plot(np.arange(0, len(v), 1)*dt, v, c='r')
plt.show()

plt.plot(u, v, c='g')
plt.show()
---------------------------------------------------------

リミットサイクル的な図に。

。。。えっ？
　前述の参考にしたサイトには「I を約 0.342 より 大きくすると (subcritical な) Hopf 分岐が起こり、周期的に変動する解が得られます。」と記載されている。上ではI=0.33としているにも関わらず周期的に変動している。合わない。

　数値シミュレーションの関係と推測し、dt=0.01からdt=0.005に小さくして再度実行。

数値シミュレーションの都合だった模様です。

同様にして次はI=0.35、dt=0.005として再度実行。

こちらはまぁさっき見たのと同じですな。

　で、今度はこの振動子を増やす。まずは振動子の数は2つ。コネクションは1、delayも1とテキトーである。その他a, b, cとかタウとかガンマとか、パラメータのバランスをいい感じでとらないと発散しそうだなと思いつつ、まずは実装する。
---------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt

a=0.7
b=0.8
c=10

N=2
C=np.ones((N, N))
D=np.ones((N, N))
for i in range(N):
    C[i, i]=0
    D[i, i]=0

def floor(vector):
    shape=vector.shape
    vector=vector.reshape(-1)
    for i in range(vector.size):
        vector[i]=np.int(vector[i])
 
    vector.reshape(shape)
    return vector

def i_input(u, i):
    con=C[:, i]
    delay=D[:, i]
    temp=np.zeros(N)+u.shape[1]-1
    temp=floor(temp-delay)
    con[temp<0]=0
    temp[temp<0]=0
    pulse=np.zeros(N)
    for j in range(N):
        pulse[j]=u[j, temp[j]]
 
    iinput=np.dot(con, pulse)
    return iinput

def d_u(u, v):
    base=g(u[:, -1], v[:, -1])
    for i in range(u.shape[0]):
        base[i]+=i_input(u, i)
 
    return base

tau=1
gamma=1
def g(u, v):
    return tau*(v+gamma*u-u**3/3)

def d_v(u, v):
    return h(u[:, -1], v[:, -1])

def h(u, v):
    return -(u-a-b*v)/tau

t=0
T=10
dt=0.005

u=np.array([1,1]).reshape(-1, 1)
v=np.array([1,1]).reshape(-1, 1)

t+=dt
while t<T:
    du=d_u(u, v)
    dv=d_v(u, v)
    temp_u=u[:, -1]+du*dt
    temp_v=v[:, -1]+dv*dt
    u=np.c_[u, temp_u]
    v=np.c_[v, temp_v]
    t+=dt

plt.plot(np.arange(0, T, dt), u[0, :], c='g')
plt.plot(np.arange(0, T, dt), v[0, :], c='r')
plt.show()
---------------------------------------------------------

後半のvの伸びがすごい。
嫌な予感を抱きながらT=15にしてみる。

ダメそう。
実際T=20にしたら悲しみのoverflow。
現時点ではさっぱりパラメータの意義を理解していないので、その辺りを式眺めながら考えていい感じのパラメータ組み合わせを探すのが今後か。お手上げなら絨毯爆撃で。

それにしても前半の1振動子の時の図を見ると、波形から""生き物の鼓動""が感じられて素敵(※一大学生の感想)。蔵本モデルは単純な分波形が無機物的なので、その差を感じられたのが今年のクリスマスプレゼントでしたとさ。


ところで元のレビューのFitzHugh-Nagumoモデルの式ちょいちょい間違ってるよね？


(12/25 17:35追記)
vを小さくするということで、まずはbを半分にしてみたら早速振動してくれた。

ただしその意味するところは不明である。
それとリミットサイクルの回転方向が前半のものと逆というのも不自然。membrane potentialがギュインと上がった直後にrecovery potentialがギュインと上がってもとに戻る、つまりu-v平面でいえば反時計回りの方が自然である。
さらにいえばuもvも負の値に突入してるのもおかしい気がする。ここで負になってることが反時計回りとかとも関係してそうな気がする。u**3とかあるし。

つまりやっぱりわからん。。。
今後2よりも多くの振動子に増やしたいので、あんまりテキトーにその場しのぎしても後で結局戻ってくることになりそうですな。