神経科学 · 2018/06/24

FitzHugh-Nagumoモデルの位相場を求める (Python)

(著)山たー

FitzHugh-Nagumoモデルの位相場(Phase Field)を中尾裕也氏の『非線形振動子の位相縮約理論とその応用』に書いていることを元に実装してみました。ただし、FHNモデル中のvとuの定義を入れ替えています。用いたモデルは「FitzHugh-Nagumoモデルをアニメーションで見る」と同じです。

計算方法

リミットサイクル上の点の位相

リミットサイクルの周期を$T$, 角速度$\dot{\theta}=\omega$とします。つまり$T=2\pi\omega$が成り立っています。リミットサイクル上の点$\boldsymbol{X_0}$が位相$\theta=\theta(\boldsymbol{X_0})$を持つとし、$\boldsymbol{X_0}$から出発した軌道が位相の原点 $\theta=0$ を初めに通過するまで時間 $\tau\ (0\lt \tau \leq T)$ かかるとします。このとき$\theta=\theta(\boldsymbol{X_0})$は $$ \theta=2\pi-\omega\tau=2\pi\left(1-\frac{\tau}{T}\right) $$ で計算できます（下図）。

リミットサイクル外の点の位相

リミットサイクル外の点$\boldsymbol{X}$が位相$\theta=\theta(\boldsymbol{X})$を持つとします。$\boldsymbol{X}$から軌道をスタートさせ、 $T$ の整数倍 $nT$ だけ軌道を計算します。このとき位相の原点に最も近づいた時刻を $\tau$ とし、 $\tau$ から $mT\ (m\in\mathbb{N})$ を除いた時間を $\tau'$ とします。つまり \begin{align*} \tau&=\tau'+mT\\ \therefore \tau'&=\tau\%T \end{align*} が成り立ちます。$\%$は割り算の余りを計算することを意味します。このとき$\theta=\theta(\boldsymbol{X})$は $$ \theta=2\pi\left(1-\frac{\tau'}{T}\right) $$ で計算できます（下図）。

まとめ

実はリミットサイクル上でも外の点の場合と同じ計算式で計算できます。いずれにせよリミットサイクル上の点を厳密に決めるのは難しいです（数値的に揺らぐため）。

コード

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.integrate as integrate

I = 0.34 #external stimulus
a = 0.7
b = 0.8
c = 10

# 原点
u0 = 1.710
v0 = 0.374
X0 = [u0, v0]

#time
dt = 0.001
idt = 1/dt
t = np.arange(0, 40, dt)
T = 4095 #msec

def FHN(state, t):
    """
    FitzHugh-Nagumo added pulse perturbation model
    u : the membrane potential
    v : a recovery variable
    """
    u, v = state
    dudt = c * (-v + u - pow(u,3)/3 + I)
    dvdt = u - b * v + a
    return dudt, dvdt

def PhaseField(u,v):
    #軌道を計算
    Xp = integrate.odeint(FHN, [u,v], t)
    
    #軌道中の点と原点の距離を計算
    d = Xp-X0
    L2 = d[:,0]**2 + d[:,1]**2
    
    #軌道中の点で最も原点に近い最初の点のindex(=time)を取得
    tau_0 = np.argmin(L2)
    
    #周期で割って余りを出す
    tau = tau_0 % T
    
    #位相を2piで割った値
    theta_per_2pi = 1 - tau / T
    
    #位相
    #theta = theta_per_2pi * 2*np.pi
    return theta_per_2pi

#print(PhaseField(u0,v0)) #1
ds = 0.01
U, V = np.meshgrid(np.arange(-3, 3, ds),
                   np.arange(2, -1, -ds))
ulen, vlen = U.shape
arr_PhaseField = np.zeros((ulen,vlen))

#網羅的に位相を計算
for i in range (ulen):
    for j in range(vlen):
        u = U[i,j]
        v = V[i,j]
        arr_PhaseField[i,j] = PhaseField(u,v)
        print("save : "+str(i*vlen+j))

# ヒートマップを出力
plt.figure(figsize=(5,5))
sns.heatmap(arr_PhaseField, cmap = "hsv",
            xticklabels = False, yticklabels = False,
            cbar = False)
#plt.savefig('FHN_PhaseField.svg')
plt.show()

描画にはかなり時間がかかります。