神経科学  ·  2018/06/23

心臓震盪とリミットサイクル

（著）山たー

　生理学の講義であった話を検証してみようとしたという内容です。検証できてはいないので悪しからず。あと自分は心臓について講義で聞いた程度の知識しかないので、適当なことを書いていると思います。

 

心臓震盪

心臓震盪というのは胸部に衝撃が加わることで、時に致死的な不整脈が生じて血液を全身に送れなくなり、死に至るという症状です。

 

心臓は心筋が周期的に収縮することで拍動を起こし、血液を送り出すポンプの役割を果たしています。この周期性は、非線形物理の用語を用いれば、心筋の振動状態がリミットサイクル(limit cycle)上にあると言えます。リミットサイクルというのは位相空間上で孤立した閉軌道のことです。もし、心筋がリミットサイクル周辺の状態であっても、時間経過でリミットサイクルに引き込まれます。

 

普通は多少リズムがずれてもリミットサイクルに戻るのですが、衝撃を受ける位置、衝撃の強さ、タイミングによってはリミットサイクルに戻れなくなり、それが心臓震盪の原因ではないか、という話です。

 

以下では簡易的なモデルとしてFitzHugh-南雲モデルを用い、この話を考えてみます。

 

パルス摂動を加えたFitzHugh-南雲モデル

FitzHugh-南雲モデル(以下FHNモデル)は神経細胞の簡易的なモデルです。詳細は「FitzHugh-Nagumoモデルをアニメーションで見る」を参照してください。

 

このFHNモデルにパルス摂動がある瞬間において加えられるというモデルを考えます。パルス摂動(pulse perturbation)は心臓に加わる衝撃のことです。微分方程式は次のようになります。

\begin{align*} \frac{du}{dt} &= c\left(u-\frac{u^3}{3}-v+I\right)+\color{red}{p(t)}\\ \frac{dv}{dt} &= u-bv+a \end{align*} $u$は膜電位、$v$はイオンチャネルの開閉などの隠れ変数です。また、$p(t)$はパルス摂動であり、次のように表されます。 $$ p(t) = \begin{cases} \epsilon & (t_p \lt t\lt t_p+\sigma) \\ 0& (otherwise) \end{cases} $$ ここで$\epsilon$はパルス摂動の強度、$t_p$はパルス摂動の開始時刻、$\sigma$は摂動の持続時間です。

 

まず、摂動が無いときのリミットサイクルは次のようになります。振動子の状態はリミットサイクルを反時計回りします。

周期$T$はおおよそ$4.095s$です。

$t_p=1$のときに摂動($\epsilon=10, \sigma=0.1$)を加えると、振動子の状態は次の青い線上を運動します。

一次的にリミットサイクルから外れますが、その後、すぐにリミットサイクルに戻っていることが分かります。

 

$\epsilon$と$t_p$をランダムにすると次のようになりました。
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import random

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

def FHN_PP(state, t, epsilon, tp, sigma):
    """
    FitzHugh-Nagumo added pulse perturbation model
    u : the membrane potential
    v : a recovery variable
    t  : time array
    epsilon :パルス強度
    tp : pulse perturbation start time
    sigma : パルスの持続時間
    """
    u, v = state
    
    #pulse perturbation パルス摂動
    p = epsilon*(t > tp) - epsilon*(t >tp + sigma)
    dudt = c * (-v + u - pow(u,3)/3 + I) + p
    dvdt = u - b * v + a
    return dudt, dvdt

# the initial state
u0 = 1.710
v0 = 0.374
X0 = [u0, v0]

#pertubation value
epsilon = 10
tp = 1
sigma = 0.1
dt = 0.001
idt = 1/dt
t = np.arange(0, 20, dt)

#limitcycle
Xlc = integrate.odeint(FHN_PP, X0, t, args=(0,0,0))

#perturbated
Xp = integrate.odeint(FHN_PP, X0, t, args=(epsilon,tp,sigma))

plt.figure(figsize=(5,5))
plt.plot(Xp[:,0], Xp[:,1])
plt.plot(Xlc[:,0], Xlc[:,1])
plt.plot(u0, v0,"ko")
plt.plot(Xlc[int(tp*idt),0],Xlc[int(tp*idt),1],"ro")
plt.annotate("t=0",xy=(u0+0.05,v0),size=10)
plt.annotate("tp="+str(tp),xy=(Xlc[int(tp*idt),0]+0.1,Xlc[int(tp*idt),1]),size=10)
plt.xlabel("u: Membrane Potential")
plt.ylabel("v")
plt.title("Perturbated FitzHugh-Nagumo model ")
plt.grid()
plt.show()

#ランダムに生成
plt.figure(figsize=(5,5))
for i in range(10):
    tp = random.uniform(0, 4)
    epsilon = random.uniform(0,20)
    Xp = integrate.odeint(FHN_PP, X0, t, args=(epsilon,tp,sigma))
    plt.plot(Xp[:,0], Xp[:,1], "--")
    plt.plot(Xlc[int(tp*idt),0],Xlc[int(tp*idt),1],"ro")
    plt.annotate("tp="+"{0:.1f}".format(tp),xy=(Xlc[int(tp*idt),0]+0.1,Xlc[int(tp*idt),1]),size=10)

plt.plot(Xlc[:,0], Xlc[:,1],color="#ff7f0e", linewidth = 2.0)
plt.annotate("t=0",xy=(u0+0.05,v0),size=10)
plt.plot(u0, v0,"ko")
plt.xlabel("u: Membrane Potential")
plt.ylabel("v")
plt.title("Perturbated FitzHugh-Nagumo model ")
plt.grid()
plt.show()

 

なんやかんやで全てリミットサイクルに戻っていることが分かります。摂動が起こることで位相が変化していますが、元の位相と新しい位相を求めるにはアイソクロン(isochron)という考えを導入しなければなりません。詳細は「FitzHugh-Nagumoモデルの位相場を求める (Python) 」に譲るとして、位相場は次のようになります。

位相が確定しない、位相空間上の特異点（＝固定点）に状態が入ってしまったとき、心臓震盪が起こるのではという話だそうです。

 

元の位相(=摂動を受けたタイミングでの位相)と新位相、及び摂動の強度により、位相リセット曲線(Phase Resetting Curve; PRC)が描けます（次図）。 

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
arr_PhaseField = np.load("arr_PhaseField.npy")

def FHN_PP(state, t, epsilon, tp, sigma):
    """
    FitzHugh-Nagumo added pulse perturbation model
    u : the membrane potential
    v : a recovery variable
    t  : time array
    epsilon :パルス強度
    tp : pulse perturbation start time
    sigma : パルスの持続時間
    """
    u, v = state
    
    #pulse perturbation パルス摂動
    p = epsilon*(t > tp) - epsilon*(t >tp + sigma)
    dudt = c * (-v + u - pow(u,3)/3 + I) + p
    dvdt = u - b * v + a
    return dudt, dvdt

def readPhase(u,v):
    """
    u,v：座標　⇒ i,j
    v=2 j=0
    v=1.99 j=1
    v=-0.99 j=299
    """
    i = int(u*100 + 300)
    j = int(v*(-100) + 200)
    #print(i)
    #print(j)
    return arr_PhaseField[j,i]
   
# the initial state
u0 = 1.710
v0 = 0.374
X0 = [u0, v0]

#pertubation value
epsilon = 10
tp = 1
sigma = 0.1
dt = 0.001
idt = 1/dt
T = 4.095 #msec
t = np.arange(0, 20, dt)
tp = np.arange(0, T, dt)

phase = np.zeros((len(tp),2)) #old and new phase

#limitcycle
Xlc = integrate.odeint(FHN_PP, X0, tp, args=(0,0,0))

pertubated_pos = np.zeros((len(tp),2))

for i in range(len(tp)):
    Xp = integrate.odeint(FHN_PP, X0, t, args=(epsilon,tp[i],sigma))
    pertubated_pos[i] =  Xp[int(tp[i]*1000 + 100),:]
    print(i)
    
#位相を意味込むために変換
pertubated_pos = np.round(pertubated_pos,2)

for i in range(len(tp)):
    phase[i,0] = readPhase(Xlc[i,0],Xlc[i,1])
    phase[i,1] = readPhase(pertubated_pos[i,0],pertubated_pos[i,1])
    print(i)
    
x = np.arange(0, 1, 1/4095)
plt.figure(figsize=(6,6))
sns.set_style("ticks")
plt.plot(x,x,"--",alpha=0.5)
plt.plot(x,phase[:,1])
plt.title("Phase Resetting Curve")
plt.xlim([0, 1.0])
plt.ylim([0, 1.0])
plt.xlabel("Old Phase")
plt.ylabel("New Phase")
plt.grid()
plt.tight_layout()
plt.savefig("PRC.png")

 

PRCをまとめると位相リセット地図(Phase Resetting Map)が描けます（次図）。これは摂動の強さを-20~40の間で変化させたものです。摂動を加えたのはuのみですから、x軸の方向に左右に摂動を加えたということです。

  

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

"""
位相を読み込み　
→　https://omedstu.jimdo.com/2018/06/24/fitzhugh-nagumoモデルの位相場を求める-python/
"""
arr_PhaseField = np.load("arr_PhaseField.npy")

def FHN_PP(state, t, epsilon, tp, sigma):
    """
    FitzHugh-Nagumo added pulse perturbation model
    u : the membrane potential
    v : a recovery variable
    t  : time array
    epsilon :パルス強度
    tp : pulse perturbation start time
    sigma : パルスの持続時間
    """
    u, v = state
    
    #pulse perturbation パルス摂動
    p = epsilon*(t > tp) - epsilon*(t >tp + sigma)
    dudt = c * (-v + u - pow(u,3)/3 + I) + p
    dvdt = u - b * v + a
    return dudt, dvdt

def readPhase(X):
    """
    u,v：座標　⇒ i,j
    u:-3→3
    v:2→-1
    """
    i = np.clip(int(X[0]*100 + 300),0,599)
    j = np.clip(int(X[1]*(-100) + 200),0,299)
    return arr_PhaseField[j,i]
   
# the initial state
u0 = 1.710
v0 = 0.374
X0 = [u0, v0]

#pertubation value
sigma = 0.1
dt = 0.001
idt = 1/dt
T = 4.095 #msec
t = np.arange(0, 10, dt)
tp = np.arange(0, T, dt)
epsilon = np.arange(-20, 40, 0.1)

pertubated_phase = np.zeros((len(tp),len(epsilon))) #old and new phase

#limitcycle
Xlc = integrate.odeint(FHN_PP, X0, tp, args=(0,0,0))

for j in range(len(epsilon)):    
    for i in range(len(tp)):
        #print(str(epsilon[j]),str(tp[i]))
        Xp = integrate.odeint(FHN_PP, X0, t, 
                              args=(epsilon[j],tp[i],sigma))
        pertubated_phase[i,j] = readPhase(Xp[int(tp[i]*1000 + 100),:])
        print(j,i)
        
plt.figure(figsize=(5,7))
sns.heatmap(pertubated_phase.T, cmap = "hsv",
            xticklabels = False, yticklabels = False,
            cbar = True, 
            cbar_kws={"orientation": "horizontal"})
plt.gca().invert_yaxis()
plt.xlabel("Old Phase/2π")
plt.ylabel("Pulse Amplitude")
plt.savefig('FHN_PRM.png')
plt.show()

 

このように、特異点に状態が入ってしまうのは衝撃の強度、タイミングが特定の狭い範囲内に入り込んでしまったときです。強い衝撃はもちろん危険ですが、小さな衝撃であってもタイミングが悪ければ心臓震盪が起こるかもしれないと言えます。

 

まとめ

今回はFitzHugh-南雲モデルを用いましたが、「これが心臓震盪のメカニズム」というには、より精度の高いモデルを作り、それを用いてシミュレーションする必要があるということです。再三になりますが、自分は講義内容を検証しただけなので、ちゃんとした心臓震盪のメカニズムは分かっておりません。

 

参考文献

・Steven H. Strogatz, 非線形ダイナミクスとカオス, 2015, 丸善出版

・高松敦子, 生命システム論, 2017 (pdfファイル)

・中尾 裕也, 非線形振動子の位相縮約理論とその応用, 2009 (pdfファイル)

・野村泰伸, 不整脈の力学系モデル (pdfファイル)

・Phase resetting in neurons -Wikipedia