利用双慢变量的快慢变量分离分析新脑皮层神经元Wilson 模型的复杂电活动*

梁艳美 陆博† 古华光

1) (河南科技学院数学科学学院,新乡 453003)

2) (同济大学航空航天与力学学院,上海 200092)

具有重要、广谱的生理功能的大脑新皮层神经元表现出规则和快速峰放电、内源性和连续型簇放电,受到两个慢变量 (T 型钙和钙激活的钾离子通道) 的调控 (Wilson 模型).若只选取更慢的慢变量进行单慢变量的快慢变量分离,不能获得Wilson 模型的簇放电的动力学机制,因为此时的快子系统还含有另一个慢变量.因此,本文采用双慢变量的快慢变量分离方法进行分析.首先获得快子系统在参数平面上的分岔及分岔曲线与簇放电相轨线的交点;基于两个慢变量都不足够慢导致的只有部分交点对应的分岔与簇放电相关,本文进一步获得三维空间中分岔曲线与相轨迹的位置关系,识别出与簇放电相关的分岔,并排除无关的分岔.结果发现,内源性簇放电与不变环上的鞍结分岔有关、与鞍结分岔无关,连续型簇放电与不变环上的鞍结分岔和超临界霍普夫分岔有关、与鞍结分岔无关.本研究全面深入认识了新脑皮层神经元簇放电的动力学机制,为调控放电模式奠定了基础.

1 引言

神经元是神经系统的基本结构单位,通过电活动实现各种生理功能.神经元电活动主要包括峰放电 (spiking) 模式、簇放电 (bursting) 模式及静息(resting state) 等[1],在感觉、认知、运动控制和疾病等方面发挥着重要作用[2−5].峰放电模式表现出连续快速放电,而簇放电模式表现为簇 (burst,连续快速的放电峰(spike)) 与休止期的交替.揭示放电模式随参数变化的转迁规律或分岔序列有助于认识神经系统的功能.例如,温度和钙离子诱发神经元电活动从簇放电到峰放电的转迁规律,分别与温度和痛觉相关[6,7];钾离子浓度调控峰放电到簇放电 (扩散性抑制) 的转迁与偏头疼相关[8].除了放电模式随参数的转迁规律,放电模式产生的内在动力学及在外界调控下的复杂动力学也是神经动力学的重要环节.

簇放电模式的产生与慢变量调控有关,快慢变量分离结合分岔已经被广泛用于揭示簇放电的内在机制[9−12].神经元系统可以分解为快子系统和慢变量,如Hindmarsh-Rose 模型和Chay 模型,是由二维快子系统和一个慢变量构成的.这些模型的簇放电的动力学机制已经获得较为全面和深入的理论认识.快子系统会表现出峰放电、静息以及峰放电与静息之间的分岔;慢变量调控下簇放电行为会在快子系统的峰放电和静息之间转迁,簇的起始和结束对应于快子系统的分岔[9].因此,不同的簇放电模式就是由簇的起始和结束的分岔不同所决定的.快慢变量分离已经广泛用于揭示不同调控下的簇放电的复杂动力学机制,如pre-Bötzinger 复合体 (pre-Bötzinger complex) 的簇放电[12−14].此外,具有时滞的突触流与噪声可以作用到快子系统的分岔点附近通过改变运行轨迹来引起新的节律[15−20],而不具有时滞的多类调控可以通过改变快子系统分岔引起复杂的节律[21,22].这些复杂的放电节律还对兴奋性作用引起放电增强或抑制性作用引起放电降低的传统观念提出了挑战,丰富了非线性动力学的内涵[23−26].

虽然快慢变量分离能成功用于揭示只有单个慢变量的簇放电的动力学,但是,当慢变量不足够慢时,簇的起始和结束相位与快子系统分岔点的位置关系会有些偏差[27].更为重要的是,在实际的神经系统中,往往会受到两个及两个以上的慢变量的调控.例如,胰腺β细胞受到两个慢速钾电流的调控[28],偏头疼出现时的海马锥体神经元受到钾和钠离子浓度两个慢变量的调控[29],损伤神经纤维受到细胞内和内质网的钙离子浓度两个慢变量的调控[30].此外,还有描述新脑皮层神经元电活动的理论模型 (本文称为Wilson 模型[31]),受到低阈值T 型钙流和钙激活的钾离子超极化电流的调控等.如何对多个慢变量调控的簇放电进行快慢变量分离获得其内在动力学机制仍是有待解决的重要问题.最近,已有研究关注了两种双慢变量调控的簇放电动力学.一种是以胰腺β细胞[28]的簇放电为代表,两个慢变量中的一个比另一个慢很多,此时,将更慢的慢变量看作常数来简化模型,再进行单慢变量的快慢变量分离.另一种是以海马锥体神经元中的簇放电为代表,两个慢变量 (钾和钠离子浓度) 都特别慢.若只选择其中一个慢变量进行快慢变量分离,快子系统的分岔不能与簇的起始和结束相匹配,说明单慢变量的快慢变量分离对于分析具有双慢变量的簇放电已经失效[29].因此,文献[29]将两个慢变量同时看作分岔参数,对快子系统进行双参数分岔分析,获得了多类余维–1 分岔曲线,簇放电的相轨迹与余维–1 分岔曲线的交点对应着簇的起始和结束.实际上,受两个慢变量调控的簇放电还有很多不同的情形,例如,描述新脑皮层神经元放电的理论模型—Wilson 模型的两个慢变量都不足够慢,与上述两种都不同,可能需要提出新的双慢变量的快慢变量分离的流程.

哺乳动物新脑皮层在神经系统高级功能的实现中起着广泛和重要的作用.在新脑皮层神经元观察到4 类基本的放电形式[32−35]: 规则峰放电 (regular spiking,RS)、快速峰放电 (fast spiking,FS)、连续型簇放电 (continuous bursting,CB) 和内源性簇放电 (intrinsic bursting,IB),涉及多类离子电流[36],例如钠流、钾流、T 型钙流 (IT) 和钙激活的钾离子超极化电流(IAHP).RS 神经元主要是第4 层的棘星细胞 (spiny stellate cells) 和第2,3,5 和6 层的锥体神经元,起兴奋性调控作用.FS 神经元是非锥体细胞,例如,γ氨基丁酸能中间神经元 (GABAergic interneurons)[37],与感觉及皮层节律有关.FS 表现出后超极化 (after-hyperpolarization,AHP) 和高频特性,是因为钾流在复极化和AHP 中起重要作用[38].CB 神经元主要是棘星细胞和第2,3,4 层的锥体神经元,产生具有高频峰的簇放电,簇内峰的幅值明显降低,簇通过后去极化 (after-depolarization,ADP) 结束[39].视觉诱发的GAMMA 振荡与CB 神经元的同步活动相关[40].IB 神经元在强直流作用下生成簇放电,然后在电流的作用下表现为伴随ADP 的低频峰放电.IB 神经元对应兴奋性锥体神经元,特别是在第5 层最为丰富.例如,IB 簇放电与癫痫[41]或慢波睡眠[42]有关.因此,了解皮层神经元的不同电活动对于了解大脑皮层的正常功能和病理生理过程至关重要.

理论上,基于主要的离子流INa,IK,IT和IAHP的新皮层神经元模型 (wilson model)[31]被提出,可以模拟FS,RS,IB 和CB 四种放电形式.已有诸多关于该模型的网络动力学和单神经元放电的动力学的研究.例如,利用RS 或IB 神经元网络研究动作电位波形依赖于慢电导刺激[43]、多巴胺影响GAMMA 同步振荡与注意力的关系[44],以及小世界耦合的RS 和IB 神经元网络的同步振荡机制[45].在单神经元中,研究了不同调控下的Wilson 模型的动力学行为.例如,只调控直流激励或同时调控IT和IAHP,诱发FS,RS,Bursting 和混沌等不同放电模式之间的跃迁,表明混沌可以通过加周期或倍周期分岔产生[46].但是,由于该模型为四维模型,含有两个快变量和两个慢变量,比一般的三维 (两个快变量和一个慢变量) 的描述簇放电的动力学模型复杂.因此,到目前还没有对四维Wilson 模型的4 类放电特别是簇放电的内在动力学的认识.该模型的复杂特性对于放电动力学分析,特别是簇放电的分析带来了困难,因此有许多研究对Wilson模型进行简化,简化为两维或三维模型,并用于研究皮层神经元的电活动[47−52].但是,二维Wilson模型需要其他调控才能产生簇放电.到目前为止,对于Wilson 模型的复杂放电特别是簇放电的内在动力学 (与快子系统的分岔的关系) 还没有全面和深入研究,对于外界调控下的放电模式的诸多研究也未与内在动力学建立联系,这是亟待解决的问题.

本文提出了适合于分析Wilson 模型簇放电的双慢变量的快慢变量分离流程,利用单和双慢变量的快慢变量分离来研究四维Wilson 模型峰放电和簇放电产生的动力学机制.首先,仿真了4 类放电及放电模式随参数的变化规律.其次,利用单参数的快慢变量分离揭示了两类峰放电的动力学机制.再次,发现单慢变量的快慢变量分离不能有效分析两类簇放电,簇的开始或结束相位与快速子系统的分岔点不匹配,簇内峰的幅值与极限环的幅值也不匹配.这是因为快子系统含有另一个慢变量不能作为快子系统而导致.为了克服这个问题,提出了新的双慢变量的快慢变量分离的操作流程.将慢变量T和H都作为快子系统的分岔参数,在 (H,T)空间得到余维–2 分岔点 (Bogdanov-Takens (BT)分岔和Cusp 分岔) 和多个余维–1 分岔曲线,包括不变环上的鞍结分岔 (saddle-node on invariant circle bifurcation,SNIC)、亚临界霍普夫分岔 (subcritical Andronov-Hopf bifurcation,SubHopf)、超临界霍普夫分岔 (supercritical Andronov-Hopf bifurcation,SupHopf) 和鞍-同宿轨分岔 (saddle-Homoclinic orbit bifurcation,SHC).然后,在获得簇放电的相轨迹与余维–1 分岔曲线在 (H,T) 平面上的交点的基础上,进一步得到它们在 (H,T,V) 三维空间中的位置关系,识别了交点中与簇的起始和结束的局部动力学相关的分岔,并排除了无关的分岔,获得了簇内峰的幅值与稳定极限环幅值的关系.本研究与两个慢变量都很慢的锥体神经元的结果[29]有两处不同,一是簇放电的轨线与分岔点的位置有些许偏差,这由慢变量不够慢导致;二是要在 (H,T,V) 三维空间中进一步识别簇放电相轨迹和分岔的位置关系,识别出与簇放电相关的分岔并排除无关的分岔,给出了一类适用于分析两个慢变量不是很慢的簇放电的方法.本研究具有两个方面的意义,一是给出一类新的双慢变量的快慢变量分离的操作流程,适用于两个慢变量都不是很慢的情形;二是深入认识新脑皮层神经元的两类峰放电和两类簇放电的复杂内在动力学.该研究可以为揭示不同调控下的放电模式的复杂动力学奠定基础.

2 模型与方法

2.1 Wilson 模型

Wilson 模型是常用来描述新皮层神经元峰放电和簇放电的理论模型[31],其方程如下:

该模型有4 个变量,分别为电压变量V,钾通道恢复变量R,与钙离子电流有关的电导变量T,与钙激活的钾离子电流有关的电导变量H.C代表神经元的电容.τR,τT和τH分别代表IK,IT和IAHP流的时间常数.I表示去极化电流.其中INa,IK,IT和IAHP分别代表钠流、钾流、低阈值T 型钙流和钙激活的钾离子超极化电流,对应的方程分别为:INa=m∞(V-ENa),IK=gRR(V-EK),IT=gT T×(V-ECa)和IAHP=gHH(V-EK).式中gR,gT和gH分别是钾离子、钙离子和钙激活的钾离子的最大电导,ENa,EK和ECa分别是钠离子、钾离子和钙离子的反转电位,函数m∞,R∞和T∞是关于膜电位V的多项式,各项表达式为:m∞=17.8+47.6V+33.8V2,R∞=1.24+3.7V+3.2V2和T∞=8(V+0.725)2.

Wilson 模型采用的参数值为[31]:τT=14 ms,τH=45 ms,C=1 μF/cm2,gR=26 mS/cm2.ENa=0,EK=–0.95,ECa=1.2.方程计算出的膜电位、反转电位的值乘以100 后的值的量纲对应mV,电流INa,IK,IT,IAHP和I的值除以10 后的值对应的量纲为 nA/cm2.为了描述简便,后续本文的电压和电流值不再给出量纲.

2.2 快慢变量分离

本文研究的Wilson 模型,τR=4.2 ms 或τR=1.5 ms,τT=14 ms,τH=45 ms .τT和τH大于τR,说明T和H是两个慢变量;但τT和τH跟τR的差距不够大,在一个量级,说明两个慢变量不足够慢 (根据经验,时间常数差别在2 个数量级及以上往往是足够慢).传统的快慢变量分离只有一个慢变量.本文研究的簇放电模式受两个慢变量影响,需要发展双慢变量的快慢变量分离方法进行分析.方程 (1a)—(1d) 构成全系统.

2.2.1 单慢变量

对于RS,CB 和IB,因为τH=45 ms 较大,以H为慢变量,方程 (1a)—(1c) 为快子系统.对于FS,因为gH=0,所以以T为慢变量,方程 (1a)和 (1b) 为快子系统.对快子系统,获得以单慢变量为参数的分岔图.将全系统峰(簇)放电的相轨迹叠加到快子系统的分岔图上,研究峰(簇)放电的起始和结束相位与分岔的关系.

2.2.2 双慢变量

选取T和H同时为慢变量,则方程 (1a) 和(1b) 构成快子系统.T和H均作为分岔参数,在(H,T) 参数空间上得到快子系统的双参数分岔.将全系统峰 (簇) 放电的相轨迹叠加到双参数分岔平面上,得到相轨迹与分岔曲线的交点,交点对应的分岔就是可能的与簇放电的起始和结束相位相关的分岔.进一步,在三维空间探寻分岔与相轨线的准确位置关系,获得与簇放电相关的分岔,排除与簇放电无关的分岔: 距离轨线近的就是与簇放电相关的分岔,距离远的就是与簇放电无关的分岔.

2.3 计算方法

运用四阶Runge-Kutta 的方法求解方程,积分步长为0.0001 ms,积分时间为20000 ms 以保证获得稳态放电 (图1 中的时间0 ms 对应积分时间的10000 ms).分岔计算通过XPPAUT[53]软件获得.

图1 Wilson 模型在不同参数下的膜电位 (a) RS,其中 I=10,gT=0.1 mS/cm2,gH=5 mS/cm2,τR=4.2 ms;(b) FS,其中 I=3.5,gT=0.25 mS/cm2,gH=0 mS/cm2,τR=1.5 ms;(c) CB,其中 I=8.5,gT=2.25 mS/cm2,gH=9.5 mS/cm2,τR=4.2 ms;(d) IB,其中 I=3.5,gT=1.2 mS/cm2,gH=3.4 mS/cm2,τR=4.2 ms;图 (c),(d) 中,点A 和点B 分别代表第一个簇(蓝色)的起始与结束相位Fig.1.Membrane potential of Wilson model at different parameter values: (a) RS,I=10,gT=0.1 mS/cm2,gH=5 mS/cm2,τR=4.2 ms;(b) FS,I=3.5,gT=0.25 mS/cm2,gH=0 mS/cm2,τR=1.5 ms;(c) CB,I=8.5,gT=2.25 mS/cm2,gH=9.5 mS/cm2,τR=4.2 ms;(d) IB,I=3.5,gT=1.2 mS/cm2,gH=3.4 mS/cm2,τR=4.2 ms.In Fig.1 (c),(d),the points A and B represent the start and termination phases of first burst (blue),respectively.

3 结果

3.1 四类放电

通过对τR,gT,gH和I赋予不同的参数值(见图1),模拟RS,FS,CB,IB 四种不同的电活动,膜电位的时间历程分别如图1(a)—(d) 所示.

RS 和FS 均表现出紧张性放电,如图1(a),(b)所示,FS 比RS 放电频率更高,主要是因为FS 复极化速度更快.IB 和CB 表现为周期性的簇放电模式,如图1(c),(d)所示.正如生理学上观察到的那样,CB 每个簇内峰的高度依次降低,并且每个簇以去极化后电位 (ADP) 终止.

3.2 峰峰间期分岔图

为进一步描述4 类电活动随参数的变化,不同放电模式的峰峰间期 (inter-spike intervals,ISIs)分别随参数I,gT和gH的变化如图2—图4 所示.

RS,FS,CB 和IB 电活动的ISIs 随I变化的分岔分别如图2(a)—(d) 所示.RS 和FS 的ISIs随I的递增而减小.CB 的ISIs 随着I的变化会产生分岔.类似地,IB 的ISIs 也表现出了分岔,说明两类簇放电在不同参数调控下会有周期数的变化.

4 类电活动的ISIs 随gH改变的分岔如图3(a)—(d) 所示.峰放电的ISIs 都随gH的递增而增长,不同的是FS 的ISIs 随gH增大会产生倍周期分岔.CB和IB 随gH变化表现出加周期分岔.IB 的ISIs 在gH较小时表现为周期1 的峰放电.

参数gT的改变对于神经元4 类电活动的影响如图4(a)—(d) 所示.RS,CB 和IB 的ISIs 随gT变化趋势类似,该模型会经历复杂的分岔,由峰放电转换为不同周期数的簇放电.而FS 的ISIs 随gT递增逐渐减小.

图2—图4 展示了Wilson 模型在不同参数下展示出不同的放电模式,本文将以图1 展示的4 类放电模式为代表进行快慢动力学分析.

图2 4 类放电的ISIs 关于I 的分岔 (a) RS;(b) FS;(c) CB;(d) IB.(a)—(d) 的竖虚线对应图1(a)—(d)Fig.2.Bifurcation of ISIs with respect to I of four types of firing patterns: (a) RS;(b) FS;(c) CB;(d) IB.The vertical dotted line in panels (a)–(d) corresponds to Fig.1(a)–(d).

图3 ISIs 关于 gH 的分岔 (a) RS;(b) FS;(c) CB;(d) IB.(a)—(d) 的竖虚线对应图 1(a)—(d)Fig.3.Bifurcations of ISIs with respect to gH : (a) RS;(b) FS;(c) CB;(d) IB.The vertical dotted line in panels (a)–(d) corresponds to Fig.1(a)–(d).

图4 ISIs 关于 gT 的分岔 (a) RS;(b) FS;(c) CB;(d) IB.(a)—(d) 的竖虚线对应图 1(a)—(d)Fig.4.Bifurcations of ISIs with respect to gT : (a) RS;(b) FS;(c) CB;(d) IB.The vertical dotted line in panels (a)—(d) corresponds to Fig.1(a)–(d).

3.3 峰放电和簇放电的快慢变量分离: 单慢变量

图1(a) 的RS 的单慢变量的快慢变量分离的结果如图5(a) 所示.快子系统的平衡点曲线呈“Z”形,由稳定平衡点 (红色实线) 和不稳定平衡点(黑色虚线) 组成.红色实线为稳定结点 (右下分支),具有较低膜电位,对应于静息状态;黑色虚线代表鞍点 (中间分支) 和不稳定焦点 (上分支).在H≈0.568757 处存在不变环上的鞍结分岔 (SNIC),稳定极限环 (最大振幅Vmax和最小振幅Vmin分别由上绿线和下绿线表示) 消失.稳定的极限环对应于放电.此外,“Z”形曲线的上拐点 (SN),出现在鞍点 (中间分支) 和不稳定焦点 (上分支) 的交点处.由图5(a) 可以看出,快子系统极限环幅值的最大值与全系统峰放电轨线 (蓝色轨迹) 幅值的极大值对应,全系统轨线都在快子系统的放电区域内,不在静息的区域内.

图5 峰放电的快子系统随H 或T 变化的分岔和峰放电的相轨迹 (蓝色实线) (a) RS,SNIC 代表H≈0.568757 处不变环上的鞍结分岔;(b) FS,SNIC 和SupHopf 分别代表H≈–0.365159 处不变环上的鞍结分岔和H≈7.012722 处超临界霍普夫分岔.红色实线、黑色虚线分别代表稳定和不稳定平衡点.上 (Vmax) 和下 (Vmin) 绿线代表稳定极限环的最大和最小幅值Fig.5.Bifurcation of the fast subsystem of the spiking with respect to H or T and the phase trajectory of spiking (solid blue curves): (a) RS,SNIC represents saddle-node on invariant circle bifurcation at H≈0.568757;(b) FS,SNIC and SupHopf represent saddle-node on invariant circle bifurcation at H≈-0.365159 and the supercritical Andronov-Hopf bifurcation at H≈7.012722,respectively.The red solid curves and the black dashed curves represent the stable and unstable equilibrium points,respectively.Upper (Vmax) and lower (Vmin) green curves represent the maximum and minimum amplitudes of the stable limit cycle.

对于图1(b) 的FS,gH=0,因此,将慢变量T作为分岔参数.FS 的单慢变量的快慢变量分离的结果如图5(b) 所示.快速子系统的平衡点曲线呈“S”形,由稳定平衡点 (红色实线) 和不稳定平衡点 (黑色虚线) 组成.左红色实线代表稳定结点,具有较低膜电位,对应于静息状态;右红色实线代表稳定焦点,具有较高膜电位,对应于去极化阻滞.黑色虚线代表鞍点 (中间分支) 和不稳定焦点 (上分支).稳定极限环 (绿色) 出现在SNIC (T≈–0.365159) 和 SupHopf (T≈7.012722) 之间,即放电区域.可以看出,FS 的相轨迹 (蓝色轨迹) 的极值与极限环极值对应,位于快子系统的稳定极限环区间内.

图1(c) 的CB 以慢变量H为分岔参数的快慢变量分离的结果如图6(a),(b) 所示.快子系统的分岔如图6(a) 所示,平衡点曲线呈“Z”形,上支和下支红色实线代表稳定结点,中支黑色短划线代表不稳定的鞍点.稳定结点与鞍点相碰形成两个鞍结分岔点,SN1和SN2,分别位于H≈0.287483 和H≈0.718103 处.将CB 的相轨迹叠加到图6(a)上,得到图6(b).将CB 的相轨迹叠加到图6(a),得到图6(b).浅蓝色和灰色实线代表CB 的簇和休止期.可以看出,簇的起始相位点A与快子系统的SN1有一定关系,休止期与快子系统的静息有关;但是,簇的结束相位点B与快子系统的分岔无关,簇内峰的幅值与稳定极限环的幅值无关.以上结果表明,单慢变量的快慢变量分离不能解释CB 的动力学机制,这是因为此时快子系统含有慢变量T.

图6 H 为慢变量时,CB 的快慢变量分离 (a) 快子系统的平衡点和极限环的分岔,在 H≈0.287483 和H≈0.718103 处的SN1 和SN2 代表鞍结分岔点,红色实线和黑色虚线分别代表稳定和不稳定平衡点;(b) CB 相轨迹与图 (a) 的叠加,浅蓝色和灰色实线分别代表簇和休止期的相轨迹,点A 和点B 分别代表簇的起始和结束Fig.6.The fast-slow variable dissection to CB with H taken as the slow variable: (a) The bifurcations of the equilibrium point and the limit cycle of the fast subsystem.SN1 at H≈0.287483 and SN2 at H≈0.718103 represent saddle-node bifurcation points;respectively,the red solid curves and the black dashed curves represent stable and unstable equilibrium points,respectively;(b) phase trajectory of CB plotted with Fig.6(a),the light blue and gray solid curves represent the phase trajectory of the burst and the quiescent state,respectively,and points A and B represent the start and termination phases of burst,respectively.

图1(d) 的IB 以H为慢变量的快慢变量分离如图7(a),(b)所示.快子系统的动力学行为和分岔如图7(a),包括静息 (右下支红色实线),放电(上、下绿色实线) 和去极化阻滞 (左上支红色实线),多个平衡点的分岔 (SupHopf (H≈0.331783)和SubHopf (H≈0.919389));鞍结分岔 (SN1和SN2),极限环的鞍结分岔 (FLC1和FLC2).浅蓝色和灰色实线代表簇和休止期的相轨迹,点A和点B分别代表簇的开始和结束.虽然休止期与快子系统的静息有关,而且簇的起始相位 (点A) 比较靠近SN1分岔,结束相位 (点B) 比较靠近FLC1,但簇内峰的幅值与稳定极限环的最大、最小值 (上、下绿色实线) 并不匹配.不匹配的原因是此时系统还有慢变量T.

图7 H 为慢变量时,IB 的快慢变量分离 (a) 快子系统的平衡点和极限环的分岔;(b) IB 相轨迹与图 7 (a) 的叠加,上、下支红色实线和中支黑色短划线、点线分别代表稳定和不稳定平衡点.上、下绿色 (深蓝色) 实线分别代表稳定 (不稳定) 极限环,浅蓝色和灰色实线代表簇和休止期的相轨迹,点A 和B 分别代表簇的起始与结束.SN,SubHopf,SupHopf 和FLC 分别代表鞍结分岔、亚临界霍普夫分岔、超临界霍普夫分岔和极限环的鞍结分岔Fig.7.The fast-slow variable dissection to IB with H taken as the slow variable: (a) Bifurcation of the equilibrium point and limit cycle of fast subsystem;(b) phase trajectory of IB plotted with Fig.7(a).The upper and lower red solid curves represent stable equilibrium points,and the middle black dashed and dotted curves represent unstable equilibrium points,respectively.The upper and lower green (dark blue) solid curves represent the stable (unstable) limit cycles.The light blue and gray solid curves represent the phase trajectory of the burst and the quiescent state,and points A and B represent the start and termination phases of burst,respectively.SN,SubHopf,SupHopf,and FLC represent the saddle-node bifurcation,the subcritical Andronov-Hopf bifurcation,the supercritical Andronov-Hopf bifurcation,and the fold limit cycle bifurcation,respectively.

3.4 簇放电的快慢变量分离: 双慢变量

对于CB 和IB,H和T同时作为慢变量,由变量V和R构成的方程 (1a)和(1b) 作为快子系统,快子系统中的H和T作为分岔参数.在 (H,T) 空间上获得的快子系统的双参数分岔结构;然后,获得簇放电相轨迹与快子系统的分岔曲线关系;依据位置关系和分岔来揭示簇放电的动力学机理.

3.4.1 CB

对于CB,快子系统的双参分岔平面 (H,T)上的分岔曲线 (粗线) 和簇放电的相轨迹 (细线)如图8(a) 所示.红色实线是SN 分岔曲线,绿色实线是SNIC 分岔曲线,黑色实线是SupHopf 分岔曲线.CB 的相轨迹 (细线) 与分岔曲线共有A,B和C三个交点,分别对应SNIC,SupHopf 和SN,需要对A,B和C三点处的分岔分别进行分析来获得簇放电的复杂动力学.其中,在A和C点上方的参数空间,快子系统表现出复杂的分岔,如图8(b) 所示,与簇放电的相轨线距离较远,不影响簇放电的动力学行为.存在余维–2 分岔点,包括退化分岔点P2[54,55],BT 和Cusp 分岔点.左、右鞍结分岔曲线相碰产生Cusp 分岔 (H≈0.322795 和T≈0.0592512).右SN 曲线 (右红色) 从Cusp 发出,向左下运行至BT 分岔点 (H≈0.3182227 和T≈0.0552198),同时出现SupHopf 分岔曲线 (黑色) 和SHC 分岔曲线 (蓝色),因此BT 点是导致SupHopf 和SHC 分岔的原因;SHC 与SNIC 曲线相交于点P2,SHC 在点P2处结束,同时SNIC 分岔曲线变为SN 分岔曲线.

图8 CB 的快子系统在 (H,T ) 平面上的双参分岔和相轨迹 (a) SNIC (绿色实线)、SN (红色实线) 和SupHopf (黑色实线) 分别代表不变环上的鞍结分岔曲线、鞍结分岔曲线和超临界霍普夫分岔曲线,BT 代表余维–2 分岔点;浅蓝色和灰色实线代表簇和休止期的相轨迹;(b) (a) 中BT 分岔点附近的放大.P2 是退化分岔点,SHC (蓝色实线) 代表鞍-同宿轨分岔Fig.8.The bifurcations of fast subsystem of CB in two-parameter plane (H,T ) and the phase trajectory: (a) SNIC (green solid line),SN (red solid line) and SupHopf (black solid line) represent the saddle-node on invariant circle bifurcation curve,the saddlenode bifurcation curve and the supercritical Andronov-Hopf bifurcation curve respectively,BT represents the codimension-2 bifurcation point;the light blue and gray solid curves represent the phase trajectory of the burst and the quiescent state;(b) the enlargement near BT bifurcation point in Fig.(a).P2 is a degenerate bifurcation point,and SHC (blue solid line) represents the Saddle-Homoclinic orbit bifurcation.

在参数平面 (H,T) 上,对应A,B和C三点的T分别取值为0.00522,0.535852 和0.004636.对二维快子系统,将T分别取为A,B和C三点的值,获得对参数H的单参数分岔,将CB 的相轨迹叠加到单参数分岔图上,如图9(a)—(c) 所示,能分别解释簇放电的相轨迹在A,B和C点附近的局部动力学.从图9(a) 可以看出,簇的起始点A距离SNIC 分岔较近;图9(b) 中,簇的结束点B距离SupHopf 分岔较近;但是,簇放电的幅值与极限环的幅值并不对应,说明只能揭示簇的起始点和结束点附近的局部动力学.对于交点C(T≈0.004636),快慢变量分离如图9(c) 所示,点C距离SN 分岔点比较远,说明C点不是真正的分岔曲线和簇放电的相轨迹的交点;因此,C点对应的SN 分岔与簇放电无关.簇放电的相轨迹和余维–1 分岔曲线是否相交还要在更高维的空间内才能得到准确的认识,比如 (H,T,V) 空间内.

图9 在对应A,B 和C 三个T 值处的CB 的快慢变量分离 (a) 对应CB 的起始点 A(T≈0.00522),A 点对应SNIC 分岔点;(b) 对应CB 的结束点 B(T≈0.535852),B 点对应SupHopf 分岔点;(c) 对应点 C(T≈0.004636);点C 与SN 分岔点距离较远.浅蓝色和灰色实线代表簇和休止期的相轨迹.红色实线、黑色虚线和绿色实线分别代表稳定平衡点、不稳定平衡点和稳定极限环.SN,SNIC 和SupHopf 分别代表鞍结分岔、不变环上的鞍结分岔和超临界霍普夫分岔Fig.9.The fast-slow variable dissection to CB at three T values corresponding to points A,B and C: (a) Corresponding to initial point A of CB (T≈0.00522).Point A corresponds to the SNIC bifurcation point;(b) corresponding to termination point B of CB(T≈0.535852).Point B corresponds to the SupHopf bifurcation point;(c) corresponding to point C (T≈0.004636).Point C is far away from the SN bifurcation point.The light blue and gray solid curves represent the phase trajectory of burst and quiescent state,respectively.The red solid curves,black dashed curves,and green solid curves represent the stable equilibrium points,the unstable equilibrium points,and the stable limit cycle,respectively.SN,SNIC,and SupHopf represent the saddle-node bifurcation,the saddle-node on invariant circle bifurcation,and the supercritical Andronov-Hopf bifurcation,respectively.

为全面和准确揭示簇放电的相轨迹和快子系统的分岔的关系,在三维空间 (H,T,V) 中,将CB的相轨迹与快子系统的分岔叠加,如图10(a)—(d)所示.如图10(a) 所示,蓝色、黄色和粉色曲面构成平衡点曲面,对应不稳定焦点、鞍点和稳定结点,上、下绿面代表稳定极限环的极大值 (Vmax) 和极小值 (Vmin),绿、红和黑粗线分别代表SNIC,SN和SupHopf 分岔曲线 (局部放大见图10(b)—(d)),黑色细线代表CB 的相轨迹.可以看出,簇内峰的幅值与极限环的幅值 (绿面) 也是基本符合的,说明簇放电确实与快子系统的放电有关;簇放电的休止期与粉色面 (静息) 很近,说明簇放电确实与快子系统的静息有关.

图10 (H,T,V) 空间中CB 的双参分岔和相轨迹.蓝色、黄色和粉色面分别由不稳定焦点、鞍点和稳定结点构成,红色、绿色、黑色粗线分别代表鞍结 (SN) 分岔、不变环上的鞍结 (SNIC) 分岔和超临界霍普夫 (SupHopf) 分岔,上、下绿面代表稳定极限环的极大值 (Vmax) 和极小值 (Vmin),黑色细线代表簇的相轨迹 (a) 全图;(b) 图(a)中B 点附近的放大图,B 点位于SupHopf 分岔曲线正下方,对应簇的结束相位;(c) 图(a)中BT 分岔点附近的放大;(d) 图(a)中A 点和C 点附近的放大,A 点在SNIC 分岔曲线正下方,对应CB 轨迹的起始相位,C 点距离SN 分岔曲线较远Fig.10.Two-parameter bifurcation of fast system and phase trajectory of CB in the (H,T,V ) space.The blue,yellow,and pink surfaces are composed of the unstable focus,the saddle and the stable node,respectively.The red,green,and black thick curves represent the saddle-node (SN) bifurcation,the saddle-node on invariant circle (SNIC) bifurcation,and the supercritical Andronov-Hopf (SupHopf) bifurcation.The upper and lower green surfaces represent the maximum value (Vmax) and the minimum value(Vmin) of the stable limit cycle,respectively,and the black thin curve represents the phase trajectory of the bursting: (a) Global view;(b) the enlargement near point B in Fig.(a),point B is located just below the SupHopf bifurcation curve,corresponding to the termination phase of the burst;(c) the enlargement near BT bifurcation point in Fig.(a);(d) the enlargement near points A and C in Fig.(a).The point A is just below the SNIC bifurcation curve,corresponding to the start phase of burst,and point C is far from the SN bifurcation curve.

图10(b) 为簇的结束点B附近的放大,可以发现,簇的结束B靠近黑粗线SupHopf 分岔,但略有偏差,说明簇放电结束于SupHopf 分岔,但因慢变量不够慢使得结束点与分岔点有些许偏差.簇的起始点A和簇放电的相轨迹的关系见图10(c),(d)(图10(a)的局部放大).余维–2 分岔BT,以及余维–1 分岔SNIC,SupHopf,SN 和 SHC 等的空间位置关系清晰可见,如图10(c) 所示.图10(d) 中,簇的起始点A与SNIC 分岔曲线距离较近,说明簇起始于SNIC 分岔.可以明显看出,簇放电的C点与SNIC 分岔曲线距离很远,并无关系.

因此,与图9 的平面图相比,图10 给出三维空间中CB 相轨迹与快子系统的分岔图可以更全面认识簇放电的内在动力学: 簇与快子系统的放电有关,休止期与静息有关;簇的起始与SNIC 分岔有关,结束与SupHopf 有关;簇放电与SN 分岔无关,虽然平面图 (图9) 提示SN 分岔曲线与簇放电相轨迹相交,但实际上在三维空间并不相交.因此,在三维参数空间 (H,T,V) 中,簇放电的动力学和快子系统的分岔能更为全面和准确的识别簇放电动力学.

3.4.2 IB

对于IB 平面上获得快子系统的分岔类似于CB,如图11(a),(b) 所示.红色、绿色和黑色实线分别对应SN,SNIC 和SupHopf 分岔曲线.IB 的簇 (蓝色线) 和休止期 (灰色线) 的相轨迹位于双参分岔图的左下角,即BT 分岔点的左下方.图11(b)为图11(a) 的放大,并且将放电区域标记为蓝色,空白色区域对应静息态.IB 的相轨迹顺时针穿过SNIC 和SN 曲线,共有A,B,C和D四个交点.点A和B对应SNIC 分岔,点C和D对应SN 分岔.

同样,为探究IB 的簇的起始与结束,需要对4 个交点处的分岔单独进行分析.在 (H,T) 平面上,图11(b) 中A,B,C和D四个交点处的T值分别为0.003489,0.178687,0.138212 和0.025683.在A,B,C和D四个点,将IB 的轨线叠加到二维快子系统随H变化的分岔图上,如图12(a)—(d)所示.图12(a) 中,簇的起始点A与SNIC 相匹配,簇的休止期与静息有关,并且首个峰的幅度可以很好地用快子系统的分岔来解释,但是无法解释簇的结束.簇的结束相位B点与SNIC 分岔点相近但略有差距,是因为慢变量不够慢导致的.从图12(c),(d) 可以看出,点C和点D距离SN 分岔点都较远,说明簇放电与SN 无关.该结果仅能解释IB 起始和结束相位附近的局部动力学.

图12 在对应A,B,C 和D 四个T 值处IB 的快慢变量分离,红色实线、黑色虚线和绿色实线分别代表稳定平衡点、不稳定平衡点和稳定极限环,SN,SNIC 和SupHopf 分别代表鞍结分岔、不变环上的鞍结分岔和超临界霍普夫分岔.浅蓝色和灰色实线代表IB 的相轨迹 (a) 对应点 A(T≈0.003489,簇的起始相位);(b) 对应点 B(T≈0.178687,簇的结束相位);(c) 对应点 C(T≈0.138212,距离SN 分岔曲线较远);(d) 对应点 D(T≈0.025683,距离SN 分岔曲线较远)Fig.12.The fast-slow variable dissection to IB bursting at four T values corresponding to points A,B,C and D.The red solid curve,the black dashed curves,and the green solid curves represent the stable equilibrium point,the unstable equilibrium point,and the stable limit cycle,respectively.SN,SNIC,and SupHopf represent the saddle-node bifurcation,the saddle-node on invariant circle bifurcation,and the supercritical Andronov-Hopf bifurcation,respectively.The light blue,and gray solid curves represent the phase trajectory of IB: (a) Corresponding to point A(T ≈ 0.003489,the start phase of burst);(b) corresponding to point B(T≈0.178687,the termination phase of burst);(c) corresponding to point C(T≈0.138212,far from the SN bifurcation curve);(d) corresponding to point D(T≈0.025683,far from the SN bifurcation curve).

在三维空间 (H,T,V) 中,IB 的相轨迹和快子系统的分岔之间的关系如图13(a),(b) 所示,图13(b) 为图13(a) 的局部放大.平衡点曲面为“Z”型,右下方粉色区域对应静息状态,中支黄色区域对应鞍点,上支蓝色区域对应不稳定的焦点;上、下两个绿面对应稳定极限环的最大值 (Vmax)和最小值 (Vmin);中支与下支平面的交线 (绿色)对应SNIC 分岔曲线,中支与上支平面的交线 (红色) 对应SN 分岔曲线.黑色细线代表IB 的相轨迹.点A,B,C和D对应图11(b) 中的4 个交点.如图13(a) 所示,簇放电的峰的幅值与快子系统的极限环的幅值相对应,簇放电的休止期与快子系统的静息 (粉色曲面) 相对应.

图11 (H,T) 平面上的IB 的相轨迹与快子系统的分岔,浅蓝色和灰色实线代表簇和休止期的相轨迹 (a) SNIC (绿色实线)、SN (红色实线) 和SupHopf (黑色实线) 分别代表不变环上的鞍结分岔曲线、鞍结分岔曲线和超临界霍普夫分岔曲线,BT 代表余维–2 分岔点 (Bogdanov-Takens);(b)图11(a) 中IB 相轨迹附近的放大图,蓝色区域对应放电区域,空白区域对应静息态Fig.11.The phase trajectory of IB bursting and bifurcations of the fast subsystem in phase plane (H,T ).The light blue and gray solid curves represent the phase trajectory of burst and quiescent state,respectively: (a) SNIC (green solid curves),SN (red solid curves),and SupHopf (black solid curves) represent the saddle-node on invariant circle bifurcation curve,the saddle-node bifurcation curve,and the supercritical Andronov-Hopf bifurcation curve,respectively,and BT represents the codimension-2 bifurcation point (Bogdanov-Takens);(b) the enlargement around the phase trajectory of IB in Fig.11(a),the blue area corresponds to the firing behavior,and the blank area corresponds to the resting state.

如图13(b) 所示,簇的起始点A和结束点B与SNIC 分岔曲线相近但略有距离,是慢变量不够慢引起的;C点和D点与SN 分岔曲线较远,说明簇放电与SN 分岔无关.

图13 (H,T,V) 空间中快子系统的双参分岔和IB 的相轨迹,H 和T 为慢变量,蓝色、黄色和粉色面分别由不稳定焦点、鞍点和稳定结点构成,红色和绿色粗线分别代表鞍结 (SN) 分岔和不变环上的鞍结分岔 (SNIC),上、下绿面代表稳定极限环的极大值和极小值,黑色细线代表IB 的相轨迹 (a)全图;(b) IB 相轨迹附近的放大,点A 和点B 对应簇的起始相位和结束相位,点C 和点D 距离SN 分岔曲线较远Fig.13.Two-parameter bifurcations of fast subsystem in the (H,T,V) space and phase trajectory of IB with H and T regarded as slow variables.The blue,yellow,and pink surfaces are composed of unstable focus,saddle,and stable node,respectively.The red and green thick curves represent the saddle-node (SN) bifurcation and the saddle-node on invariant circle (SNIC) bifurcation,respectively.The upper and lower green surfaces represent the maximum and minimum values of the stable limit cycle,and the thin black line represents the phase trajectory of IB: (a) Global view;(b) the enlargement around the phase trajectory of IB in Fig.13(a).Points A and B correspond to initial and termination phases of the burst,and points C and D are far from the SN bifurcation curve.

4 结论

新脑皮层神经元具有广谱、重要的生理功能,表现出复杂多样的峰和簇放电模式,依赖于神经元的类型和生理参数.揭示新脑皮层神经元放电的复杂动力学是认识其生理功能的重要环节[9].本文提出新的双慢变量的快慢变量分离的流程,结合传统的单慢变量的快慢变量分离方法,全面、深入地分析了新脑皮层神经元的两类峰放电和两类簇放电的复杂动力学,创新性体现在以下两个方面.

一是提出了新的双慢变量的快慢变量分离的流程.单慢变量的快慢变量分离在分析只受单个慢变量调控的簇放电或峰放电时获得了巨大成功[9],但是,不能用于分析两个慢变量调控的簇放电的动力学机制[28,29],因为此时快子系统含有慢变量.因此,需要根据双慢变量的不同特征开发新的快慢变量分离方法.例如,β细胞的簇放电受到很慢和较慢的慢变量调控[28],固定很慢的慢变量为常数、然后利用较慢的慢变量为分岔参数进行快慢动力学分析.在海马锥体神经元[29],因为两个慢变量都很慢,将两个慢变量都作为快子系统的分岔参数进行双参数分岔分析;因为慢变量足够慢,双参数平面的分岔曲线和簇放电相轨迹的交点与簇放电全都相关.但是,对于本文研究的新脑皮层神经元的簇放电,两个慢变量都不足够慢,这是不同于β细胞和锥体神经元的一种新情况.因此,本文首先参照文献[29]的流程,获得了快子系统的双参数平面的分岔曲线与簇放电的相轨迹的交点.与文献[29]结果不同,部分分岔点与簇放电相关,部分不相关.进一步,本文获得了三维空间的簇放电轨迹与快子系统的分岔曲线、平衡点曲面以及极限环的极值曲面的位置关系,依此识别了与簇放电相关的分岔,排除了与簇放电无关的分岔;与文献[29]中分岔曲线和簇放电相轨线的交点与簇放电都相关不同.

二是全面、深入地认识了新脑皮层神经元峰放电和簇放电的动力学机制.对于峰放电,只受或主要只受一个慢变量的调控,可以利用传统的单快慢变量分离的方法来认识峰放电 (RS 和FS) 的动力学行为,峰放电对应的轨迹只位于快子系统的放电区域内,与静息和分岔无关.而对于簇放电,受到两个慢变量的调控,单慢变量的快慢变量分离不能解释CB 和IB 簇放电的动力学行为,即簇的起始和结束相位与快子系统的分岔点、簇内峰与极限环的幅值或者簇的休止期与快子系统的静息之间不匹配.因此,将两个慢变量H和T作为快子系统的分岔参数,获得双参数平面 (H,T) 上的余维–1 分岔曲线 (SNIC 和SupHopf 分岔曲线) 和余维–2 分岔点 (BT 和Cusp 分岔点),用于分析簇的起始和结束相位.CB 簇放电的相轨迹与分岔曲线相交于A,B和C三点;通过三个交点处的单参分岔分析,可知CB 簇放电轨迹的起始和结束相位仅与SNIC(A点) 和SupHopf (B点) 相关.IB 的相轨迹与分岔曲线相交于A,B,C和D四点,而簇的起始和结束仅与SNIC (A和B点) 相关.这就获得了簇放电起始点和结束点附近的局部动力学,但是,该平面方法也有局限性: 不能认识簇内峰的幅值与稳定极限环的幅值是否匹配,以及簇放电的相轨迹与快子系统的分岔曲线的具体位置关系.因此,本文进一步分析了三维空间 (H,T,V) 中的簇放电轨迹与快子系统平衡点曲面、极限环极值曲面和分岔曲线的位置关系,进一步解释了平面图上SN 分岔曲线与簇放电相轨迹相交是三维空间投影导致的,排除了与簇放电无关的分岔.此外,还展示了簇与快子系统的放电有关,休止期与静息相关.这说明,针对两个慢变量都不足够慢的情况,本文给出了分析其簇放电动力学行为一种有效的双慢变量的快慢变量分离流程.

近期研究提示,针对不同的双慢变量的快慢特征,可以发展不同的快慢变量分离法来分析簇放电的复杂动力学.本文提出了新的流程分析了两个不足够慢的慢变量 (T 型钙和钙激活钾通道) 相关的簇放电的动力学.将利用本文的方法研究Wilson模型在不同参数下的不同放电的复杂动力学,将研究不同调控下 (如不同电流、突触流、自反馈等) 的簇放电的复杂动力学,以及如何调控快子系统的分岔点来影响簇放电的动力学行为.进一步,还可以基于双慢变量的快慢变量分离法研究Wilson 模型在网络中的复杂动力学.