一类四维神经元模型的放电模式研究

Study on Firing Patterns for a Four-Dimensional Neuron Model

摘要: 本文研究一类描述坐骨神经慢性压迫性损伤的四维神经元模型。利用快慢动力学方法，并结合数值模拟结果来分析此模型的各种放电模式。通过作出快子系统随慢变量变化的单参数分岔图，以及对相应于放电状态的稳定极限环的讨论，来得到不同放电模式出现的动力学机制。

Abstract: This paper studies a sciatic nerve chronic constriction injury model. This model consists of four coupled nonlinear differential equations. By using the fast-slow dynamical methods and numerical simulations, we analyze the firing patterns of the model. We give the single parameter bifurcation diagrams of the fast subsystem, and discuss the stable limit cycles which are corresponding to the firing state of the neuron. Dynamical mechanisms of different firing patterns are obtained.

1. 引言

对应不同的放电模式，神经元有着不同的编码和信息处理方式。利用动力系统的分岔理论和快慢动 力学方法来研究神经元模型的不同类型的放电模式，是近几十年被大量采用和行之有效的方法，可参见文献 [1] - [6] 等。如最早由Rinzel [1] 在研究神经元簇放电模式时提出的“切片法”(dissection)，将神经元模型看成快慢动力系统来研究，即将模型中的变量分为快变量和慢变量，快变量负责动作电位的产生，而慢变量起着调节快变量的作用，其动力学行为表现为静息状态与反复放电状态的相互转迁，从而产生出簇放电。

为了与试验数据相吻合，文 [7] 将描述坐骨神经慢性压迫性损伤的三维Chay模型 [8] [9] 进行了改进，变为如下四维神经元模型

$\{\begin{array}{l}\frac{\text{d}V}{\text{d}t}=g_i{m}_{\infty }^3{h}_{\infty }\left(v_i-V\right)+g_{kv}\left(v_k-V\right)n^4+g_{kc}\frac{C_{in}}{1+C_{in}}\left(v_k-V\right)+g_l\left(v_l-V\right),\hfill \\ \frac{\text{d}n}{\text{d}t}=\frac{n_{\infty }-n}{{\tau }_n},\hfill \\ \frac{\text{d}{C}_{in}}{\text{d}t}=m_{\infty }^3{h}_{\infty }\left(v_c-V\right)-k_c{C}_{in}-k_{erp}{C}_{in}+\frac{k_{rel}{C}_{in}}{k_{crc}+C_{in}}\left(C_{lum}-C_{in}\right),\hfill \\ \frac{\text{d}{C}_{lum}}{\text{d}t}=k_{erp}{C}_{in}-\frac{k_{rel}{C}_{in}}{k_{crc}+C_{in}}\left(C_{lum}-C_{in}\right),\hfill \end{array}$ (1)

其中V为膜电压，n为钾离子通道激活的概率， $C_{in}$ 为细胞内钙离子浓度， $C_{lum}$ 对应内质网上腔钙浓度。参数 $v_i,v_k,v_l,v_c$ 分别表示混合Na⁺-Ca²⁺通道，电压依赖K⁺通道(或钙依赖K⁺通道)，泄漏离子通道和Ca²⁺通道的反转电位， $g_i,g_{kv},g_{kc},g_l$ 分别代表各通道的最大电导除以膜电容， $m_{\infty },h_{\infty }$ 分别为混合Na⁺-Ca²⁺ 通道激活和失活的概率， $n_{\infty }$ 为n的稳定态值， ${\tau }_n$ 是松弛时间， $k_c$ 表示细胞内钙离子流出的比例常数， $k_{erp}$ 表示钙库中Ca²⁺-ATPase的泵激活率， $k_{rel}$ 和 $k_{crc}$ 分别表示 $C_{lum}$ 的释放率和 $C_{in}$ 的离解常数。用y表示h、

m或n，则由 $y_{\infty }=\frac{{\alpha }_y}{{\alpha }_y+{\beta }_y}$ 可得 $h_{\infty },m_{\infty }$ 和 $n_{\infty }$ ，其中 ${\alpha }_h,{\beta }_h,{\alpha }_m,{\beta }_m,{\alpha }_n,{\beta }_n,{\tau }_n$ 的表达式如下

${\alpha }_m=\frac{0.1\left(25+V\right)}{1-{\text{e}}^{-0.1V-2.5}}$ , ${\beta }_m=4{\text{e}}^{-\frac{V+50}{18}}$ , ${\alpha }_h=0.07{\text{e}}^{-0.05V-2.5}$ , ${\beta }_h=\frac{1}{1+{\text{e}}^{-0.1V-2}}$ ,

${\alpha }_n=\frac{0.01\left(20+V\right)}{1-{\text{e}}^{-0.1V-2}}$ , ${\beta }_n=0.125{\text{e}}^{-\frac{V+30}{80}}$ , ${\tau }_n=\frac{1}{{\lambda }_n\left({\alpha }_n+{\beta }_n\right)}$ .

在本文中，参数值取为

$v_l=-40\text{mV}$ , $v_k=-75\text{mV}$ , $v_i=160\text{mV}$ , $g_i=1080\text{pS}$ , $g_{kv}=1500\text{pS}$ ,

$g_{kc}=16\text{pS}$ , $g_l=\text{7pS}$ , ${\lambda }_n=320$ , $v_c=60$ , $k_c=\frac{3.3}{18}$ , $k_{erp}=15$ , $k_{rel}=0.3$ , $k_{crc}=0.5$ .

在文 [7] 中，作者通过数值模拟，对模型(1)进行了分析，讨论了加周期分岔、混沌在峰放电和簇放电中的变化过程等复杂动力学现象，说明了四维模型(1)与三维Chay模型相比，其放电模式与实验数据更加吻合。在本文中，我们将利用分岔和快慢动力学方法 [1] [10] ，结合分岔软件MATCONT的使用，对产生模型(1)的几种放电模式的动力学机制进行研究。

2. 分岔和快慢动力学分析

由于变量 $C_{lum}$ 改变的速率远慢于其他的变量 $V,n,C_{in}$ ，因此我们将 $V,n,C_{in}$ 作为快变量，将 $C_{lum}$ 作为慢变量，即将模型(1)看成由快子系统

$\{\begin{array}{l}\frac{\text{d}V}{\text{d}t}=g_i{m}_{\infty }^3{h}_{\infty }\left(v_i-V\right)+g_{kv}\left(v_k-V\right)n^4+g_{kc}\frac{C_{in}}{1+C_{in}}\left(v_k-V\right)+g_l\left(v_l-V\right),\hfill \\ \frac{\text{d}n}{\text{d}t}=\frac{n_{\infty }-n}{{\tau }_n}\text{,}\hfill \\ \frac{\text{d}{C}_{in}}{\text{d}t}=m_{\infty }^3{h}_{\infty }\left(v_c-V\right)-k_c{C}_{in}-k_{erp}{C}_{in}+\frac{k_{rel}{C}_{in}}{k_{crc}+C_{in}}\left(C_{lum}-C_{in}\right)\hfill \end{array}$ (2)

和慢子系统

$\frac{\text{d}{C}_{lum}}{\text{d}t}=k_{erp}{C}_{in}-\frac{k_{rel}{C}_{in}}{k_{crc}+C_{in}}\left(C_{lum}-C_{in}\right)$ (3)

构成，将慢变量 $C_{lum}$ 作为快变子系统(2)的分岔参数。

令模型(1)的四个方程的右端为零，得到平衡点满足的方程组。现取模型(1)的后两个方程的右端为零，然后相加，得到 $m_{\infty }^3{h}_{\infty }\left(v_c-V\right)-k_c{C}_{in}=0$ ，由此解出 $C_{in}$ ，代入(3)的右端为零的式子，可求得

$C_{lum}=\frac{1}{k_c{k}_{rel}}\left[m_{\infty }^3{h}_{\infty }\left(v_c-V\right)\left(k_{erp}+k_{rel}\right)+k_{erp}{k}_{crc}{k}_c\right]$ 。 (4)

于是，由(2)的三个方程右端为零和(4)式构成的方程组的解即为(1)的平衡点坐标。因此由系统(2)在 $C_{lum}-V$ 平面作出的随分岔参数 $C_{lum}$ 变化的平衡点曲线与方程(4)确定的曲线的交点即为(1)的平衡点。

下面通过选取几个不同的 $v_k$ 值，利用快慢动力学方法来分析模型(1)的不同的放电模式的动力学机制。本文使用的方法，与文 [11] 中研究三维Chay模型的不同的簇放电模式的方法类似。我们分以下5种情形来讨论。

2.1. v_k = −70

图1(a)为v_k = −70时神经元模型(1)的膜电位V的时间历程图，其放电模式为簇放电模式。为了分析产生此放电模式的动力学机制，我们作出了含有单参数分岔的几何分析图(图1(b))。

在图1(b)中，黑色实线和虚线构成的Z型曲线为快子系统(2)的平衡点曲线，上支和下支的实线部分对应稳定平衡点曲线，中支的虚线部分对应不稳定平衡点曲线，点画线为方程(4)表示的曲线。 $S{N}_1$ 与 $S{N}_2$ 点分别为系统(2)的鞍结分岔点，P点为点画线与平衡点曲线的不稳定部分的交点，对应模型(1)的不稳定平衡点(图1(b)中小空心圆点)。图中附加的蓝色曲线为模型(1)的轨线在 $C_{lum}-V$ 平面的投影曲线，与上支和下支的稳定平衡点曲线有部分重合。

当分岔参数 $C_{lum}$ 减少时，位于Z型平衡点曲线下支附近出发的(1)的轨线，沿着下支到达鞍结分岔点

Figure 1. v_k = −70. (a) Time course of membrane potential; (b) geometric analysis diagram

图1. v_k = −70。(a) 膜电位的时间历程图；(b) 几何分析图

$S{N}_1$ ，然后向上转迁到Z型平衡点曲线上支的稳定平衡点附近，此时分岔参数 $C_{lum}$ 增加，(1)的轨线沿着上支到达鞍结分岔点 $S{N}_2$ ，然后又转迁到Z型平衡点曲线的下支，随着 $C_{lum}$ 减少又重复上述过程，从而产生出图1(a)的簇放电模式。按照文 [3] 的分类，此簇放电模式为“fold/fold”点-点滞后环型簇放电。

以下文中出现的图中各曲线、点的意义与图1中相同，不再重复。

2.2. v_k = −76

图2(a)为v_k = −76时神经元模型(1)的膜电位v的时间历程图，图2(b)为相应的几何分析图。在图2(b)中，Z型平衡点曲线的下支和中支仍然分别为(2)的稳定平衡点曲线和不稳定平衡点曲线，上支中 $H_1$ 和 $H_2$ 皆为超临界的Hopf分岔点(靠近 $H_2$ 右方还有一个次临界的Hopf分岔点未标出)，在 $H_1$ 左边和 $H_2$ 到 $S{N}_2$ 部分对应稳定平衡点曲线， $H_1$ 到 $H_2$ 部分对应不稳定平衡点曲线。红色曲线 $V_{\max },V_{\min }$ 分别表示快子系统(2)随参数 $C_{lum}$ 变化时所产生的一族稳定极限环上的膜电位 $V\left(t\right)$ 的最大值和最小值曲线(在后面的几何分析图中 $V_{\max },V_{\min }$ 表示意义一样)。点画线与平衡点曲线的不稳定部分的交点P，对应模型(1)的不稳定平衡点。

当分岔参数 $C_{lum}$ 减少时，位于Z型平衡点曲线下支附近出发的(1)的轨线，沿着下支到达鞍结分岔点 $S{N}_1$ ，然后向上转迁到Z型平衡点曲线上支的稳定平衡点附近，此时分岔参数 $C_{lum}$ 增加，(1)的轨线沿着上支到达 $H_1$ ，然后环绕 $V_{\max },V_{\min }$ 曲线向右延伸到 $H_2$ , 经过 $H_2$ 后沿着平衡点曲线附近到达 $S{N}_2$ 点，然后转迁到Z型平衡点曲线的下支，随着 $C_{lum}$ 减少又重复上述过程，从而产生出图2(a)的簇放电模式。按照文 [3] 的分类，此簇放电模式为经由“fold/fold”滞后环的“Hopf/Hopf”型簇放电。

2.3. v_k = −80

图3(a)为v_k = −80时神经元模型(1)的膜电位v的时间历程图，图3(b)为相应的几何分析图。在图3(b)中，z型平衡点曲线的下支和中支分别为(2)的稳定平衡点曲线和不稳定平衡点曲线，上支中 $H_1$ 和 $H_2$ 皆为超临界的Hopf分岔点(靠近 $H_2$ 右方还有一个次临界的Hopf分岔点未标出)，在 $H_1$ 左边和 $H_2$ 到 $S{N}_2$ 部分对应稳定平衡点曲线， $H_1$ 到 $H_2$ 部分对应不稳定平衡点曲线。点画线与平衡点曲线的不稳定部分的交点P，对应模型(1)的不稳定平衡点。

当分岔参数 $C_{lum}$ 减少时，位于z型平衡点曲线下支附近出发的(1)的轨线，沿着下支到达鞍结分岔点

Figure 2. v_k = −76. (a) Time course of membrane potential; (b) geometric analysis diagram

图2. v_k = −76。(a) 膜电位的时间历程图；(b) 几何分析图

Figure 3. v_k = −80. (a) Time course of membrane potential; (b) geometric analysis diagram

图3. v_k = −80。(a) 膜电位的时间历程图；(b) 几何分析图

$S{N}_1$ ，然后向上转迁到z型平衡点曲线上支并环绕 $V_{\max },V_{\min }$ 曲线向右延伸到 $H_2$ , 经过 $H_2$ 后沿着平衡点曲线附近到达 $S{N}_2$ 点，然后转迁到z型平衡点曲线的下支，随着 $C_{lum}$ 减少又重复上述过程，从而产生出图3(a)的簇放电模式。由文 [3] 知，此簇放电模式为经由“fold/fold”滞后环的“fold/Hopf”型簇放电。

2.4. v_k = −87

图4(a)为v_k = −87时神经元模型(1)的膜电位V的时间历程图，图4(b)为相应的几何分析图。在图4(b)中，z型平衡点曲线的下支和中支分别为(2)的稳定平衡点曲线和不稳定平衡点曲线，上支中H为超临界的Hopf分岔点(靠近H点右方还有一个次临界的Hopf分岔点未标出)，在H左边部分对应不稳定平衡点曲线，H到 $S{N}_2$ 部分对应稳定平衡点曲线。点画线与平衡点曲线的不稳定部分的交点P，对应模型(1)的不稳定平衡点。

此时根据几何分析图图4(b)进行的快慢动力学分析与v_k = −80时的分析完全一致，得到的簇放电模式也为经由“fold/fold”滞后环的“fold/Hopf”型簇放电。

2.5. v_k = −103

图5(a)为v_k = −103时神经元模型(1)的膜电位V的时间历程图，图5(b)为相应的几何分析图。在图5(b)中，z型平衡点曲线的下支和中支分别为(2)的稳定平衡点曲线和不稳定平衡点曲线，上支中H点为超临界的Hopf分岔点(靠近H右方还有一个次临界的Hopf分岔点未标出)，且在H点左边部分对应不稳定平衡点曲线，在H到 $S{N}_2$ 部分对应稳定平衡点曲线。点画线与平衡点曲线的稳定部分的交点P，对应模型(1)的稳定平衡点(在图5(b)中为实心圆点)。

Figure 4. v_k = −87. (a) Time course of membrane potential; (b) geometric analysis diagram

图4. v_k = −87。(a) 膜电位的时间历程图；(b) 几何分析图

Figure 5. v_k = −103. (a) Time course of membrane potential; (b) geometric analysis diagram

图5. v_k = −103。(a) 膜电位的时间历程图；(b) 几何分析图

我们选取模型(1)的初始点坐标 $\left(V,n,C_{in},C_{lum}\right)=\left(-35,0.3,0.1,25\right)$ ，随着分岔参数 $C_{lum}$ 的增加，(1)的轨线环绕 $V_{\max },V_{\min }$ 曲线向右延伸到H，经过H后沿着平衡点曲线到达 $S{N}_2$ 点，然后转迁到Z型平衡点曲线的下支，随着 $C_{lum}$ 减少最后趋向于稳定平衡点P，神经元最终到达静息状态。在文 [4] 中称此放电状态为簇兴奋模式。

3. 结论

文 [7] 建立了一类坐骨神经慢性压迫性损伤的四维神经元模型，并通过参数的变化，讨论了该模型复杂的放电模式。本文对此四维神经元模型，利用分岔和快慢动力学的方法，选取与文 [7] 不同的参数变化，即通过 $v_k$ 的不同取值，来讨论模型的几种放电模式，即“fold/fold”点-点滞后环型簇放电，经由“fold/fold”滞后环的“Hopf/Hopf”型簇放电和“fold/Hopf”型簇放电，以及簇兴奋模式，得到了产生相应的放电模式的快慢动力学机制。通过对产生神经元不同的放电模式的动力学研究，使我们对神经元及神经网络的信息处理方式有更深入的理解。

基金项目

国家自然科学基金(11572127)。