中文 Part 2 Izhikevich's simple neuron: Stability Analysis

$$\newcommand{\Vminus}{V_-} \newcommand{\Wminus}{W_-} \newcommand{\Vplus}{V_+} \newcommand{\dvdv}{\frac{\partial\frac{dV}{dt}}{\partial V}} \newcommand{\dvdw}{\frac{\partial\frac{dV}{dt}}{\partial W}} \newcommand{\dwdv}{\frac{\partial\frac{dW}{dt}}{\partial V}} \newcommand{\dwdw}{\frac{\partial\frac{dW}{dt}}{\partial W}}$$

For English Version press here

撰稿:Alex White
翻譯:高暐哲

許多神經元是高度非線性的，Izhikevich神經元模型也不例外。電壓方程中的$V^2$項在神經元模型的動力學中起著重要作用。使用非線性動力學來研究Izhikevich模型是很自然的。奇怪的是，非線性動力學中最重要的概念之一是線性化，線性化是將非線性系統逼近為線性系統的過程，這應該是違反直覺的。包括Izhikevich模型在內，非線性模型在神經科學領域取得的巨大成功的原因是它們可以表現的龐大的現象，其中包括全有或全無動作電位基本上也是一種非線性現象。現在，我聲稱理解這些複雜的非線性現象的秘訣是通過線性化近似非線性項$V^2$，我怎麼得出這個奇怪的結論？

線性化在神經科學中是一個非常強大的工具，因為電壓在靜息態平衡上花費了相當多的時間。線性化告訴我們，當電壓接近任何平衡時，非線性效應很弱以至於可以在數學上忽略。當我們做線性化時，我們可以問一些非常重要的問題：

1.發生神經脈衝之後神經元多快回復至靜息態？

2.回復至靜息態之後神經元會震盪嗎？

3.當神經元接受興奮性刺激後為什會離開靜息態的平衡？

4.神經元多快脫離靜息態？

5.為什麼脫離靜息態後就回不來了

許多這些看似簡單的問題只能通過對系統進行初步線性化來解決。我們將在這篇文章中回答第一個問題，並在續篇中回答第二個問題。後面的一些問題需要分岔理論（bifurcation theory）提供一些解釋（我們將在後面的文章中討論這個理論）。 但是，所有這些問題都可以通過首先將平衡點周圍的區域近似為線性系統來解答，並利用這些近似值來理解非線性動力學給系統增加的額外行為。

在這篇文章中，我想展示於各種細節中使用線性化。這些細節大部分是通過計算機的幾個簡短的指令完成的，但我覺得完全理解整個程序如何運作，並且代入所有複雜的代數是非常重要的。

為了進一步解釋，我們先回答一個看似簡單的問題。我們如何知道當Izhikevich模型在靜息態平衡時是穩定的？這意味著在沒有刺激輸入的情況下，Izhikevich模型將會一直持續維持靜息態。我們將確切地定義穩定的意義，並且我們還將展示平衡點存在的條件，當它們確實存在時系統總是穩定的。

讓我們通過生物近似的Izhikevich模型的線性化。 如果您完全不熟悉特徵值，我建議您從3blue1brown觀看此影片。我建議準備一支鉛筆和紙，這將有助於代數上的理解。

現在我們將忽略神經脈衝機制並關注在微分方程。

$$C_m \frac{dV}{dt} = k (V-V_r)(V-V_{th}) - W + I$$

$$\frac{dW}{dt} = a(b(V-V_r) -  W)$$

如果你想回顧參數代表的意涵，可以參考第1篇文章。

首先我們需要找到平衡點。回想一下，平衡點是導數為零$\frac{dV}{dt}=\frac{dW} {dt}=0$的那些點。我們可以通過繪製零斜率線來將其視覺化。藍線表示$\frac{dV}{dt}=0$的位置，同樣的，黃線代表$\frac{dW}{dt}=0$。它們在兩個地方交叉，最左邊的均衡點$Eq_{-}=(V_{ - },W_{-})$，以及最右邊的平衡點$Eq_{+}=(V_{+},W_{+})$。

圖1：Izhikevich模型的零斜率線。藍線為$\frac{dV}{dt}=0$，黃線為$\frac{dW}{dt}=0$。注意大多數點走向 $Eq_{-}$。這意味著它是一個穩定的靜息態。然而有些點的位置超出畫面的範圍。這對應於鞍點$Eq_{+}$。

我們的目標是證明當存在最左邊的平衡點$(V_{-},W_{-})$時，它總是穩定的。我們首先根據所有參數找出$V_{-},W_{-}$的值。因此，我們設$C_m \frac{dV}{dt}= 0$和$\frac{dW}{dt} = 0$，我們得到

$$0 = k (V-V_r)(V-V_{th}) - W + I$$

$$0 = a(b(V-V_r) -  W)$$

$$W = b(V- V_r)$$

$$0 = k (V-V_r)(V-V_{th}) - b(V- V_r) + I$$

請注意，這個方程在$V$項中是二次方的，並且它獨立於$W$。我們現在可以將其視為一個二次方程，並使用這個二次方程來找到

$$V=\frac{b+ k(V_r+V_{th})\pm \sqrt{(b+k(V_r+V_{th}))^2- 4 k(I-b V_r-k V_r V_{th})}}{2k}$$

在這裡，我們希望取負根$\Vminus$，所以我們取根號項為負

$$\Vminus=\frac{b+ k(V_r+V_{th})- \sqrt{(b+k(V_r+V_{th}))^2- 4 k(I-b V_r-k V_r V_{th})}}{2k}$$

我們先在這裡停一下，看看這個方程式代表的含意。因為我們想在這裡顯示兩件事。 $\Vminus \leq \Vplus$，以及當根消失時。$V_o \leq plus$

讓我們先建立我們的幾何直覺。看看GIF圖。從$I=0$開始。隨著$I$上移，$\Vplus$和$\Vminus$靠近在一起並最後於$V_o$相撞。 我們想證明，如果$b \geq 0$那$\Vminus \leq V_o$。這件事稍後會很有用。

圖2:隨著$I$增加，兩個平衡點越來越接近，直到它們最後於$V_o$相撞。進一步推動$I$消滅平衡點，神經元維持在激發態。

下面我們將需要用代數的方法操作根號內容。

$$\Vminus= \frac{b+ k(V_r+V_{th}) - \sqrt{[b+k(V_r+V_{th})]^2- 4 k(I-b V_r-k V_r V_{th})}}{2k}$$

現在我們單獨著重在根號的內容。我們希望簡化它，方便稍後在一些不等式中使用。為了做到這一點，我們將首先展開根號內部以便消掉部份的項。我們首先單獨看根號內

$$- \sqrt{[b+k(V_r+V_{th})]^2- 4 k(I-b V_r-k V_r V_{th})}$$

我們展開以獲得

$$- \sqrt{[b^2 + 2 b k V_r +k^2 V_r^2 +2 b k V_{th} + 2 k^2 V_r V_{th}+V_{th}^2] +(-4b k V_r- 4 b k V_r V_{th} -4 I k)}$$

然後收集類似的項

$$- \sqrt{[b^2 - 2 b k V_r +k^2 V_r^2 +2 b k V_{th} - 2 k^2 V_r V_{th}+Vth^2] - (4 I k)}$$

現在注意到除了符號不同外，$[]$括號中的所有內容都是相同的嗎？這意味著我們可以將括號內重新排列成這種形式

$$- \sqrt{[b+k(-V_r+V_{th})]^2- 4 k I}$$

將以下代入回$\Vminus$的表達式中得到

$$\Vminus= \frac{b+ k(V_r+V_{th}) - \sqrt{[b+k(-V_r+V_{th})]^2- 4 k I}}{2k}$$

這個方程比前面的$\Vminus$表達式簡單得多。 您現在會注意到在這個式子裡更清楚的解釋了輸入電流$I$對系統的影響。隨著$I$的增加，根號內部的值越來越小。這意味著兩個平衡點$\Vplus$和$\Vminus$越來越接近。為了更清楚地表明這一點，我們先設置$I=0$。我們得到下面的式子：

$$\Vminus= \frac{b+ k(V_r+V_{th}) - \sqrt{[b+k(-V_r+V_{th})]^2}}{2k}$$

$$\Vminus= \frac{b+ k(V_r+V_{th}) - {[b+k(-V_r+V_{th})]}}{2k}$$

$$\Vminus =V_r$$

這正是我們所期望的$\Vminus$值。 我們一樣對$\Vplus$作同樣的作法

$$\Vplus= \frac{b+ k(V_r+V_{th}) + {[b+k(-V_r+V_{th})]}}{2k}$$

$$\Vplus=\frac{2b+ 2k(V_{th})}{2k}$$

請注意$\Vminus = V_r$，但是$\Vplus \neq V_{th}$。 因此，$b$越來越正，$\Vplus$越來越正。 然而正如我們所看到的，隨著$I$增加$\Vminus$和$\Vplus$一起移動，並且最終根部將相撞於某個值$V_o$，然後消失。 事實上，當 $$\sqrt{[b+k(-Vr+Vth)]^2- 4 k I}=0$$ 時，根會發生碰撞。

解$I$這個方程式，我們看到這在什麼時候發生

$$I = \frac{[b+k(-Vr+Vth)]^2}{4 k}$$

這意味著我們得到

$$V_o = \frac{b+ k(Vr+Vth)}{2k}$$

小記$V_o$被稱為鞍點分岔，我們將在以後的文章中討論它。但現在你可以直接想成，當$I \geq  \frac{[b+k(-Vr+Vth)]^2}{- 4 k}$時，神經元總是在產生脈衝，因為沒有靜息態的平衡點。

現在我們可以開始進行線性化。回想一下，我們的目標是找到圍繞這個休息均衡定點的近似線性系統。接近平衡點意味著什麼？我附上了兩張圖來說明這一點。一個是標明$Eq_-$和$Eq_-$的整個相平面。我還附上了一張圖，放大了$Eq_-$附近的區域。 你會注意到，在這個尺度上，藍色的二次零斜率線看起來像一條線。這就意味著它幾乎是線性的。放大時系統近似線性。

我們可以使用這個近似線性來確認系統是否平衡穩定。若且唯若線性系統的所有特徵值具有負實部時，平衡點才是穩定的。

圖3：上圖為完整的系統。再次，彩虹點是系統在相空間移動的各種解。

下圖為放大於$Eq_{-}$附近的平衡點。注意，黃色和藍色的曲線看起來接近直線。因此我們可以將它線性化。紅色和黑色線是特徵向量。請注意，解首先轉向黑色特徵值，然後轉向靜息態。

為了釐清當我們放大觀察時這個線性系統的行為，我們需要找到一個雅可比矩陣。 雅可比矩陣正是這些直的藍色和黃色零斜率線的數學描述。在數學上，當$\vec{x}\approx <\Vminus,\Wminus>$. 時，雅可比矩陣 $J$ 滿足$\frac{d\vec{x}}{dt} = J\vec{x}$。為了找到 $J$ ，我們用位於$\Vminus$的平衡點切線來逼近零斜率線。為此，我們需要計算相對於$V$零斜率線和$W$零斜率線的$V$和$W$的偏導數。我們由此可以得到方程式，矩陣外觀像為：

$$J = \left( \begin{array}{cc} \dvdv & \dvdw \\ \dwdv & \dwdw \end{array} \right)$$

其中$\frac{\partial}{\partial V}$ and $\frac{\partial}{\partial W}$就是每個函數相對於$V$和$W$的導數。實際運算代數我們得到以下結果：

$$\dvdv = \frac{\partial}{\partial V}\frac{k (V-V_r)(V-V_{th}) - W + I}{Cm} = \frac{2 V-k(V_r+V_{th})}{C_m}$$

$$\dvdw = \frac{\partial}{\partial W}\frac{k (V-V_r)(V-V_{th}) - W + I}{Cm} = -\frac{1}{C_m}$$

$$\dwdv = \frac{\partial}{\partial V} a(b(V-V_r) - W) = ab$$

$$\dwdw = \frac{\partial}{\partial W} a(b(V-V_r) - W) = -a$$

設$V= \Vminus$，我們現在可以寫下我們的雅可比矩陣

$$J = \left( \begin{array}{cc}  \frac{2 \Vminus-k(V_r+V_{th})}{C_m} & -\frac{1}{C_m}\\ ab & -a \end{array} \right)$$

現在我們的下一步是找到這個矩陣的特徵值。如上所述，特徵值告訴我們系統恢復到靜息態的速度有多快。一個大的負特徵值對應於快速的恢復速度。一個二維系統有兩個特徵值（$\lambda_1$ and $\lambda_2$），並因此有兩個方向返回至靜息態，稱為特徵向量$\vec{\xi_1}$ 和 $\vec{\xi_2}$。在上面的圖片中，我繪製了具有淺灰色矢量的矢量場的完整Izhikevich相平面。我也附加了放大版本。這次圖表上有紅色和黑色的線條。黑線和紅線是線性系統的特徵向量。在幾何上，如果一個線性解開始於這些特徵向量之一，那麼它就停留在這些特徵向量上永遠不會脫離。但是大多數解都是從特徵向量開始的。系統沿著紅色特徵向量比黑色特徵向量更快地恢復至靜息態。這意味著紅色特徵值比黑色特徵值更負。

圖4A：放大的$Eq_{-}$。這個系統的兩個特徵值都是負的。然而，與紅色特徵向量的特徵值$\lambda_1$衰減得比與黑色特徵向量的特徵值$\lambda_2$更快。這表示紅$\lambda_1$比$\lambda_2$更負。我們得到$\lambda_1 \leq \lambda_2$。

從圖形上來說，這意味著解很快沿著紅色特徵向量衰減，並沿黑色特徵向量緩慢衰減。

圖4B：放大的$Eq_+$。請注意，與紅色特徵向量的特徵值是負值，另一個黑色特徵向量的特徵值是正的。這意味著系統有一個正的和一個負的特徵值。這種平衡點是一個鞍點。注意，因為黑色特徵值是正的，所以解會遠離平衡點。這意味著這個特徵向量是不穩定的。

另外，這是公式的幾何解釋：

$$\begin{align} V(t) = C_1 \xi_{1_V} e^{{\lambda_1}t} +C_2 \xi_{2_W} e^{{\lambda_2}t}\\ W(t) = C_1 \xi_{1_W} e^{{\lambda_1}t} +C_2 \xi_{2_W} e^{{\lambda_2}t} \end{align}$$

請注意$\vec{\xi_1}=<\xi_{1_V},\xi_{1_W}>$。 我們可以用這個公式來估計我們的神經元在接近靜息態時的行為。通過這個公式我們可以看出特徵值的重要性：定義穩定性。如果特徵值$\lambda$為正，則$t\rightarrow \infty$，然後$e^{\lambda t}\rightarrow \infty$。但是，如果$\lambda$為負數，當$t\rightarrow \infty$，那麼$e^{\lambda t}\rightarrow 0$。這告訴我們解沿著它的特徵向量移動的方向。因此我們知道如果一個系統的特徵值都有負實部，那麼它是穩定的。此外，由這些特徵值可以得知它恢復到靜息態的速度。請注意$\Vplus$有一個正特徵向量。這產生了分界面，也是神經元能夠產生脈衝的原因。

我們現在將為$\Vminus$找到這兩個特徵值$\lambda_1$和$\lambda_2$的值，並證明它們總是為負。我們為這些特徵值找到的公式將用於創造各種有趣的神經元。

為了找到雅可比矩陣的特徵值，我們需要找到這個系統的行列式

$$\det(J-\lambda \boldsymbol{I}) = \left| \begin{array}{cc}  \frac{2 \Vminus-k(V_r+V_{th})}{C_m}-\lambda & -\frac{1}{C_m}\\ ab & -a-\lambda \end{array} \right|= 0$$

其中$$\boldsymbol{I}$是單位矩陣。如果你不熟悉行列式，請看這裡。這個矩陣的行列式產生下面給出的特徵多項式：

$$

0=\left(\frac{2 \Vminus-k(V_r+V_{th})}{C_m}-\lambda\right)\left(-a-\lambda\right) - (-\frac{1}{C_m})(b)

$$

收集所有$\lambda$我們得到

$$

\lambda^2 +\left(a -\frac{k(V_r+V_{th}-2\Vminus)}{C_M}\right)\lambda -\left( \frac{ab - a k(V_r+V_{th}-2\Vminus)}{C_m} \right)

$$

解$\lambda$產生另一個複雜的二次方程式。

$$

\lambda =  -a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus) \pm  \sqrt{(a -\frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}+a\frac{k}{C_m}(V_r+V_{th}-2\Vminus)\right]}

$$

現在這個等式有很多項，但這是我們之前所做的所有看似無意義的工作在這裡都會得到回報。讓我們先看根號以外的部分。我們想要得到 $$-a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus)\leq 0$$ 。我們注意到$a \geq 0$，而且我們知道$2\Vminus \leq 0$。 因此，我們必須證明 $$2\Vminus \geq V_r +V_{th}$$ 。 當這個不等式成立時，則得到 $$-a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus)\leq 0$$ 。

回想一下，我們有$V_o \geq \Vminus$ and that $V_o = \frac{b+ k(V_r+V_{th})}{2k}$。 現在我們可以說明

$$\frac{b}{2k}+(Vr+Vth) \geq 2\Vminus \geq V_r +V_{th}$$

$$\frac{b}{2k}+(Vr+Vth) \geq V_r +V_{th}$$

$$\frac{b}{2k} \geq 0$$

並且因為$b \geq 0$和$k \geq 0$，此不等式成立。 因此，根號外總是負的。現在這意味著至少其中一個特徵值總是負的。註腳

我們的下一個目標是要表明，當我們取正根時，根號項總是小於根號之外的部分。我們得到公式

$$-a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus) >  \sqrt{(a -\frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}+a\frac{k}{C_m}(V_r+V_{th}-2\Vminus)\right]}$$

現在你可能已經可以看到我們要做什麼了。上面的步驟我們展開並整理了根號內的項。

$$\sqrt{(a -\frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}+a\frac{k}{C_m}(V_r+V_{th}-2\Vminus)\right]}$$

$$\sqrt{(a^2 -2 a \frac{k}{C_m}(V_r+V_{th}-2\Vminus)+ \frac{k^2}{C_m^2}((V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}+a\frac{k}{C_m}(V_r+V_{th}-2\Vminus)\right]}$$

$$\sqrt{(a

+ \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}\right]}$$

代回到我們的lambda方程中

$$

\lambda =  -a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus) +  \sqrt{(a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}\right]}

$$

現在讓我們回過頭來看這代表的意義。記得$b \geq 0$，所以我們可以在$b=0$時得到一些特別的啟發。我們注意到越來越大的$b$，根號項越來越小。這給了我們

$$

\lambda =  -a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus) +  \sqrt{(a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2}

$$

$$

\lambda =  -a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus) + (a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus)

$$

$$\lambda = 0$$

當特徵值為0時，這意味著什麼？ 這意味著它不會沿著這個方向改變。

我們還記得

$$

\frac{dW}{dt} = a(b V- W)

$$

因此，現在$W$不依賴於$V$，系統在某種意義上被破壞了。 然而$b\neq 0$和因為更正的$b使得根號項更小。公式上為：

$$+\sqrt{(a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2} \geq+\sqrt{(a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}\right]}$$

因此我們可以得出結論，我們的系統總是有兩個負的特徵值。

這些是相當初步的數學。但希望它給了我們如何線性化Izhikevich模型一些概念。但我想用這些程序來實現更進階的功能。注意，更大的$b$使根號項越小。事實上，如果

$$(a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2 \leq 4\left[a\frac{b}{C_m}\right]$$
Then
$$(a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}\right] \leq 0$$

這意味著我們正在取一個負數的平方根！我們知道負數的平方根是一個虛數，代表我們的特徵值是複數。我們也知道實部總是為負。 因此我們的兩個特徵值將取$\lambda = \alpha \pm i \beta$ where $\alpha = -a - \frac{k}{C_m} (V_r+V_{th} -2\Vminus)$ and $\beta=|\sqrt{(a + \frac{k}{C_m}(V_r+V_{th}-2\Vminus))^2-4\left[a\frac{b}{C_m}\right]}|$

這在圖形上是什意思？這在圖形上就引發了旋轉行為。由於特徵值存在虛部，解將圍繞平衡點旋轉。 由於$\alpha$的負值，這個螺旋向內緩慢地移向平衡點。在下一篇文章中，我們將展示如何由此來構建一種非常特殊的神經元，它只對$\beta$的特定頻率的輸入有反應。

圖5: 特徵值含虛部的系統。注意$Eq_+$在畫面右側。注意到$Eq_-$的特徵值皆具有負實部。 因此，即使平衡點向靜息態震盪，平衡點也是穩定的。