Pre-Bötzinger复合体的从簇到峰放电的同步转迁及分岔机制*

杨永霞 李玉叶† 古华光

1) (赤峰学院, 数学与计算机科学学院, 赤峰 024000)

2) (赤峰学院, 应用数学研究所, 赤峰 024000)

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

Pre−Bötzinger复合体是兴奋性耦合的神经元网络, 通过产生复杂的放电节律和节律模式的同步转迁参与调控呼吸节律.本文选用复杂簇和峰放电节律的单神经元数学模型构建复合体模型, 仿真了与生物学实验相关的多类同步节律模式及其复杂转迁历程, 并利用快慢变量分离揭示了相应的分岔机制.当初值相同时,随着兴奋性耦合强度的增加, 复合体模型依次表现出完全同步的“fold/homoclinic”, “subHopf/subHopf”簇放电和周期1峰放电.当初值不同时, 随耦合强度增加, 表现为由“fold/homoclinic”, 到“fold/fold limit cycle”、到“subHopf/subHopf”与“fold/fold limit cycle”的混合簇放电、再到“subHopf/subHopf”簇放电的相位同步转迁, 最后到反相同步周期1峰放电.完全(同相)同步和反相同步的周期1节律表现出了不同分岔机制.反相峰同步行为给出了与强兴奋性耦合容易诱发同相同步这一传统观念不同的新示例.研究结果给出了pre−Bötzinger复合体的从簇到峰放电节律的同步转迁规律及复杂分岔机制, 反常同步行为丰富了非线性动力学的内涵.

1 引 言

神经系统是神经元构成的复杂的非线性网络系统, 表现出复杂的非线性电活动或时空行为, 比如分岔、混沌和同步及同步转迁等, 参与感觉、运动控制和脑高级功能的实现[1−6].其中, 呼吸节律对生命维持具有重要意义, 也是四大生理指标之一 .脑 内 被 称 为 pre−Bötzinger复 合 体 (pre−Bötzinger complex, pre−BötC)的网络与呼吸节律的产生与调控有关[7−12].目前, 实验已经发现pre−BötC内的单个神经元会表现出复杂的放电节律模式特别是簇放电节律模式, pre−BötC会表现出复杂的同步行为[13−18].因此, 利用非线性科学揭示pre−BötC的单神经元的复杂节律和其同步行为是重要的科学问题, 有助于从理论上认识实验发现的复杂多样的节律及其动力学、复杂同步转迁的动力学规律, 对于神经科学和非线性物理都有重要的意义.

生物学实验已经发现, pre−BötC是一个兴奋性突触耦合的神经元网络, 可以产生同步的簇振荡[13,17−21], 以及复杂的同步节律模式和节律模式的同步转迁[13−16,19−21].例如, Smith 等[13,19−21]在pre−BötC离体实验观察到许多不同模式的簇放电.此外, pre−BötC的单神经元会表现出复杂的节律模式及其转迁[13,15,16,22−26].例如, Dunmyre等[25]在pre−BötC的离体实验中发现了一系列的放电模式, 包括静息态、峰放电、方波簇放电、去极化阻滞(depolarization block, DB)以及混合簇放电(方波簇与DB簇); Negro等[26]在离体实验中发现增大细胞外液的钾离子浓度([K+]o), 神经元会经历从静息到簇放电、再到峰放电的节律转迁模式.

针对实验发现大量的非线性动力学现象, 已经有相应的理论研究.描述pre−BötC神经元电活动的离子通道理论模型被建立, 并模拟了不同的放电模式及其转迁[24,26−32].例如, Butera等[24]建立神经元模型并模拟了静息态、簇放电和峰放电.Negro等[26]在Butera模型的基础上, 调节与钾离子通道相关参数, 神经元会经历从静息到簇放电、再到峰放电的节律转迁模式, 也就是会经历复杂的分岔过程.郭豆豆和吕卓生[32]在修改的Butera模型的基础上加入电磁辐射和外强迫电流研究它们对混合簇的影响.进一步, 兴奋性耦合的pre−BötC网络模型的同步行为被研究, 并且已经获得了一些符合常规概念的结果, 如增大兴奋性突触耦合强度引起同步簇放电频率增加, 降低兴奋性突触耦合引起同步簇放电频率降低[27], 兴奋性突触耦合强度变化诱发不同簇放电的现象[33].Dunmyre等[25]构建了包括多个离子流如钙电流、持续钠电流和自耦合的pre−BötC的理论模型, 模拟了静息态、峰放电、方波簇放电、DB簇以及混合簇放电(方波簇与DB簇混合).刘深泉等[34,35]研究了持续钠电流的平衡电位、延迟钾电流的平衡电位以及它们的电导系数对pre−BötC神经元膜电位的影响.因为pre−BötC单神经元就表现出十分复杂的非线性动力学和多时间尺度动力学, 并且神经元间的耦合是化学耦合, 也是强非线性的, 与复杂网络动力学经常研究的电耦合(线性耦合)是不同的, 这就给 pre−BötC动力学的研究带来困难.因此,Best等[29]和段利霞等[36−40]利用快慢变量分离方法来揭示同步簇放电节律模式的复杂非线性动力学.例如, Best等[29]研究了不同突触耦合强度下的对称与非对称的同步簇放电和峰放电模式; 段利霞等[36−40]利用快慢变量分离和双参数分岔分析继续研宄了不同参数(如钠、钾、钙和IP3等)对同步放电模式的影响, 并且区分了不同钠、钾、钙离子电导下的同步簇放电模式, 如混合簇[39]和簇内同向同步和反向同步, 揭示了簇内的峰(Spike)的同相同步和反相同步的动力学与初值的关系: 簇内峰的同相同步对应于相同初值, 簇内峰的反相同步对应不同初值[40].Rubin等[41−42]运用数值模拟和快慢分离变量方法, 描述了神经网络的行为, 同时,分析了产生和控制呼吸节律的状态依赖机制.此外, 对于神经元的异质性、Na+/K+泵、突触性质、网络相互作用对同步的影响也有研究[33,43].但是,与复杂网络的同步动力学的研究相比较, pre−BötC的同步动力学的研究还相对初步, 对同步模式和类型以及同步转迁的过程的关注还较少, 用来刻画同步的指标也相对简单(大多采用的是同步因子), 大多数的研究关注了某个具体参数值下的同步动力学特性, 而随关键参数如耦合强度连续的变化规律研究较少等.

文献[40]中揭示的兴奋性耦合网络中的簇内峰的反相同步的节律与常规的兴奋性作用容易引起同相同步有所不同, 给出了一种不同于传统观念的新示例.近期, 与传统观念不同的新的非线性现象和机制引起了越来越多的关注[44−50].比如, 抑制性作用引起单神经元从静息变为放电、增强放电频率、引起神经元网络的同相同步[44−48], 不同于常规概念的抑制性作用压制放电和引起反相同步; 兴奋性作用可以引起放电簇内峰的个数降低[49]、引起神经元网络从放电变为静息[50].因此, 进一步揭示pre−BötC是否存在兴奋性作用引起的与常规概念不同的新示例, 是一个重要的非线性科学问题.

基于此, 本文将按照研究网络同步的一般方法和流程, 研究pre−BötC的复杂节律的同步转迁及其分岔机制.采用了多个指标刻画同步(除同步因子外, 还采用的峰相位差、簇相位差以及峰峰间期(interspike intervals, ISIs)), 关注于关键生理参数兴奋性耦合强度的连续变化, 并着重于通过快慢变量分离和分岔分析来揭示簇放电的内在非线性动力学, 区分了相同初值与不同初值对同步节律的影响, 关注了与常规概念不同的非线性现象.研究结果给出了随着兴奋性耦合强度增加pre−BötC的从簇到峰放电的同步转迁规律, 并利用快慢变量分离区分了节律转迁的分岔机制; 强兴奋性耦合强度下反相峰同步与常规的兴奋性作用诱发同相同步的概念不同, 并利用快慢变量分离给出了这一反常现象的分岔机制.全面、深入认识了呼吸相关节律的动力学机制, 丰富了非线性动力学的内涵, 为调控pre−BötC的复杂同步节律及转迁提供了参考.

2 神经元模型介绍及同步指标的判定

2.1 神经元模型介绍

2.1.1 单神经元模型

Butera 等[24]以 Hodgkin−Huxley 模型为基础建立了pre−BötC的离子通道单腔室模型, 方程如下:

其中C是膜电容, V是膜电位, h和n分别代表Na+通道和K+通道的打开概率.(1)式右端5项分别代表持续性钠电流、钠电流、钾电流、漏电流和细胞受到兴奋性刺激所产生的耦合电流; gNap, gNa,gK, gL和 gtonic-e分别为各通道的最大电导; ENa, EK,EL和 Etonic-e分别是相应的平衡电位.(2)式中 ε 是调节参数, 以确保变量h是慢变量.x∞(V) 是x的稳态电压依赖激(失)活函数; τx(V) 是电压依赖时间常数; 其表达式形式如下:

其中 gK是分岔参数, 模型中的其他参数值和量纲如表1所示.

2.1.2 耦合神经元模型

Pre−BötC的兴奋性耦合双神经元模型由Butera等提出[24,27], 其是由单腔室pre−BötC模型经化学突触耦合的网络模型, 方程如下:

其中si表示门控变量; 对于i, j = 1, 2且

gsyn-esi(Vi−Esyn-e)表示耦合神经元网络的耦合电流, 记为

其中 gsyn-e为耦合强度, Esyn-e为突触的反转电位.因为 si与Vj有关, 就实现了神经元间的耦合.当Vi是振荡或放电时,

不再可能为常值0, 这是与电突触不同的, 可能会引起复杂的动力学.设 Esyn-e=0mV 确保突触电流是 兴 奋 性 的.参 数 αs=0.2ms−1,σs= −5mV , θs= −10mV , τs=5ms .其他参数值和量纲与单神经元模型相同.

表1 理论模型中的参数值Table 1.Parameter values used in the theoretical model.

2.2 同步指标判定

2.2.1 耦合神经元的峰相位函数及相位差

神经元的峰相位函数常用来刻画同步, 耦合神经元l (l =1,2)的峰相位函数 ϕl(t) 定义如下:

2.2.2 耦合神经元的簇相位函数及相位差

与峰相位函数相一致, 耦合神经元l(l=1,2)的簇相位函数 Φl(t) 定义如下:

簇相位差定义为:∆Φ(t)=|Φ1(t)−Φ2(t)|, 按照 max(∆Φ(t)) 值定义相位同步: (1)max(∆Φ(t))=0, 簇同相同步; (2) max(∆Φ(t))=π , 簇反相同步;(3) 0 2π , 簇的非相位同步.

2.2.3 耦合神经元膜电位的相关系数

两个耦合神经元的膜电位的相关系数也常用于刻画同步, 其定义如下:

相关系数ρ越大, 两神经元膜电位的关联度越高, 即同步程度越高.ρ =1 表示完全同步,ρ=0表示不同步; | ρ |<0.6 , 0 .6 ≤ |ρ|<0.8 和0.8 ≤ |ρ|<1分别代表同步程度的弱、中和强.

2.3 神经元模型的快慢变量分离

2.3.1 单神经元模型的快慢变量分离

对单神经元模型(1)—(3)进行快慢变量分离.(1)和(2)式作为快子系统, (3)式描述的h是慢子系统, 选h为快子系统的分岔参数, 研究单神经元的快子系统的平衡点和极限环分岔类型.将单神经元全系统的簇放电的放电轨迹镶嵌到快子系统的分岔图上, 通过轨线与分岔的位置关系来识别不同的簇放电模式.

2.3.2 耦合神经元模型的快慢变量分离

耦合神经元模型的快慢变量分离与单神经元类似.对于(4)—(7)式描述的耦合神经元模型,(5)式描述的 hi(i = 1, 2)变量为慢子系统.若研究神经元i的簇放电的快慢变量分离, 将h=h1=h2选为参数作为慢变量, 剩余的6维方程为快子系统.首先, 将h选为参数, 研究快子系统的平衡点和极限环的分岔.然后, 将神经元1(或2)的簇放电的放电轨迹镶嵌到快子系统的分岔图上, 通过轨线与分岔的位置关系来识别不同的簇放电模式.

在本文中, 采用的是四阶龙格−库塔(Runge−Kutta)数值积分法, 积分步长为0.001 ms.分岔计算采用软件是XPPAUT[51]进行分岔的数值计算.

3 仿真结果

3.1 单神经元模型放电行为及ISI分岔图

3.1.1 单神经元放电行为

在 pre−BötC 单神经元模型, 当 gK分别取为7.1, 7.8, 10.0和25.0 nS时, 单神经元系统表现出DB簇放电、周期18方波簇放电、周期12方波簇放电和周期3方波簇放电, 分别如图1(a), 图1(b),图1(c)和图1(d)所示.

图1 不同 gK 下单神经元放电在(h, V)相平面的轨迹 (a) gK=7.1nS ; (b) gK=7.8nS ; (c) gK=10.0nS ; (d)gK=25.0nSFig.1.The (h, V) trajectory of the single neuron at different g K values: (a) gK=7.1nS ; (b) gK=7.8nS ; (c) gK=10.0nS ;(d) gK=25.0nS .

图2 单神经元模型的随 g K 的分岔 (a) ISIs分岔序列; (b)图(a)左下角方框的局部放大Fig.2.Bifurcation of the single neuron model with increasing g K : (a) Bifurcations of ISIs; (b) the enlargement of ISIs within the square at the down−left corner of fig (a).

3.1.2 单神经元模型的复杂ISI分岔序列

在 pre−BötC 单神经元模型, 以 gK作为分岔参数, 发现随着 gK的增大, 放电经历复杂的转迁规律,由静息进入DB簇放电、再进入高周期数的方波簇放电, 经过逆加周期分岔进入周期1方波簇放电.对应的峰峰间期(ISIs)的分岔图如图2(a)所示,图2(b)是图2(a)的局部放大图, 与Negro等[26]在实验中发现的增大细胞外的钾离子浓度引起的pre−BötC神经元会经历从静息到簇放电、再到峰放电的节律转迁历程相近似.

3.1.3 单神经元的快慢变量分离及簇放电模式分类

在不同分岔参数 gK下, 对单神经元系统(1)—(3)进行快慢变量分离.以 gK=7.1nS 为例(当 gK=7.8,10.0,25.0nS 时分岔结构类似, 关键点参数值参见表2), 随着慢变量h的增加, 快子系统的平衡点在(h, V)平面上展示了“S”型曲线, 其中下支和上支的蓝色实线分别为稳定结点和稳定焦点, 上支和中支的黑色虚线分别是不稳定的焦点和鞍点, “S”型曲线上的拐点 F1 (h ≈ 0.4928)代表平衡点的鞍结分岔, F2 位于h ≈ —1.6780, 如图3所示.

表2 不同 gK 下快子系统中关键点的慢变量h的值Table 2.The values of slow variable h of the bifurcation or key points at different g K values.

图3 当 gK=7.1nS 时, 单神经元模型的快子系统随着慢变量h变化的分岔Fig.3.Bifurcations of the fast−subsystem of the single neur−on with respect to h when gK=7.1nS .

在“S”型曲线的上支, 随着慢变量h增加, 当h ≈ 0.2128时, 不稳定焦点变为稳定焦点, 同时出现了不稳定极限环, 发生了亚临界霍普夫分岔(即subh点); 当 h ≈ 0.3265时, 发生了鞍−同宿轨分岔(即HC), 随着h继续增大, 有稳定的极限环出现; 当h ≈ 0.4308时, 稳定的极限环与不稳定的极限环相碰发生了极限环的鞍结分岔(即LPC), 如图3所示.快子系统的放电区域(即稳定极限环存在)是从HC点到LPC点, 于是h ∈ (0.3265, 0.4308)是快子系统的稳定的结点与极限环(放电)共存的区域, 红色实线和虚线分别代表稳定及不稳定极限环的最大、最小值.

接着, 将 gK=7.1nS 的全系统的簇放电(如图1(a)所示)在(h, V)平面上的轨迹叠加到其快子系统的分岔图上, 如图4(a)所示.从图4(a)可以看出:随着h的增加, 全系统放电轨线从下支的静息态经鞍结分岔F1(h ≈ 0.4923)转迁到上支, 再随着h的减小, 放电轨线围绕稳定焦点作阻尼振荡, 当经过亚临界Hopf分岔点subh时, 振幅逐渐减小到零, 在subh点左侧不稳定焦点的作用下, 系统振荡幅度逐渐增大, 但是没有达到稳定极限环的幅值, 放电超出了快子系统的放电区域最后转变为静息态.根据Izhikevich簇放电分类的介绍[52], 此放电簇称为经“fold/homoclinic”滞后环的“subHopf/subHopf”型簇放电.

当 gK=7.8,10.0,25.0nS 时, 其相应的快子系统在(h, V)平面的平衡点和全系统的簇放电(如图1(b), 图1(c), 图1(d)所示)在(h, V)平面上的轨迹叠加, 分别如图4(b), 图4(c)和图4(d)所示.由图4(b), 图4(c)和图4(d)可知, 它们簇的类型是相同的, 因为它们有相类似的快子系统分岔图且全系统的放电都是处在相应的快子系统放电和稳定结点共存的区域内, 如表2和图4所示.随着h的增加, 全系统放电轨线从下支的静息态经鞍结分岔(F1)转迁到上支重复放电区域, 再随着h的减小, 重复放电通过鞍−同宿(HC)分岔转为下支的静息态, 此类型的簇放电被称为“fold/homoclinic”型簇放电.

3.2 兴奋性化学突触耦合神经元模型的同步转迁

图4 单神经元在不同的 gK 下簇放电模式的快慢变量分离 (a) gK=7.1nS ; (b) gK=7.8nS ; (c) gK=10.0nS ; (d) gK= 25.0 nSFig.4.The fast−slow variable dissection of bursting of single neuron at different g K values: (a) gK=7.1nS ; (b) gK=7.8nS ;(c) gK=10.0nS ; (d) gK=25.0nS .

为了与以前的多个研究的条件相一致, 本文研究的如(4)—(7)式所示的两个pre−BötC神经元经兴奋性化学突触耦合形成的网络也是同质的, 就是两个神经元的参数是相同的.分岔参数 gK取为7.8 nS, 其它参数值与单神经元模型完全一致.考虑两种初值情况: 1)两耦合神经元的初值完全相同; 2)两耦合神经元的初值不同, 对应初值向量中任一分量相差超过10—10.对于相同初值和不同初值, 我们都研究了多组初值.相同初值选取为(V1,h1, n1, s1) = (1.74551, 0.49343, 0.7561, 1.53 ×10—4)和 (V2, h2, n2, s2) = (1.74551, 0.49343,0.7561, 1.53 × 10—4), 不同初值选取为 (V1, h1, n1,s1) = (1.74551, 0.49343, 0.7561, 1.53 × 10—4)和(V2, h2, n2, s2) = (—52.1421, 0.45472, 0.00306,2.81 × 10—4).在这两种情况下改变耦合强度 gsyn-e,研究 gsyn-e对耦合神经元的同步节律以及同步转迁的影响.

3.2.1 相同初值条件

因为两个耦合神经元是同质的, 当两个神经元的初值相同时, 它们接受到的耦合电流=gsyn-esi(Vi−Esyn-e)是相同的, 网络一定是完全同步的.但是, 神经元接受到的耦合电流是非0 的, 随着 gsyn-e增加也会增加.因此, 随着耦合强度 gsyn-e的增加, 网络内的单个神经元的节律是会变化的, 变化结果如图5(a1)—(a5)所示.

第一部分 (0 nS≤ gsyn-e< 2.45nS 、第二部分(2 .45nS ≤ gsyn-e< 4.45nS)、 第 三 部 分 (4.45nS≤gsyn-e<12.70nS)和第四部分 (12.70nS ≤gsyn-e≤18.0nS)分别为完全同步方波簇放电, 完全同步的DB簇放电, 完全同步的DB簇放电(与第二部分的DB簇存在区别, 此DB簇的放电尖峰的中间没有小的振荡, 几乎是一个固定值), 和完全同步周期1峰放电.四个部分的代表性的放电节律和相应的耦合电流如图6(a) (gsyn-e=0.35nS), 图6(b)(gsyn-e=2.5nS), 图6(c) (gsyn-e=5.0nS)和图6(d)(gsyn-e=18.0nS)的上部实线和下部虚线所示.其中方波簇、DB簇和周期1峰放电与Dunmyre等[25]的实验发现相似.

图5 随着耦合强度 g syn-e 增大, 耦合神经元模型的同步转迁过程.相同初值 (a1)耦合电流平均值 ; (a2)峰相位差的最大值max(∆(ϕ(t))); (a3)簇相位差的最大值 max(∆(Φ(t))) ; (a4)相关系数ρ; (a5)神经元1的ISIs序列.不同初值: (b1)耦合电流平均值; (b2)峰相位差的最大值 max(∆(ϕ(t))) ; (b3)簇相位差的最大值 max(∆(Φ(t))) ; (b4)相关系数ρ; (b5)神经元1的ISIs序列Fig.5.Transitions with respect to g syn-e of coupled neurons model.The same initial values: (a1) The mean values of coupling cur−rent ; (a2) maximum spike phase difference max(∆(ϕ(t))) ; (a3) maximum burst phase difference max(∆(Φ(t))) ; (a4) coefficient ρ; (a5) ISIs of neuron 1.Different initial values: (b1) The mean values of coupling current ; (b2) maximum spike phase difference max(∆(ϕ(t))); (b3) maximum burst phase difference max(∆(Φ(t))) ; (b4) coefficient ρ; (b5) ISIs of neuron 1.

3.2.2 不同初值条件下

第一部分(0 ≤gsyn-e<0.45nS): 两耦合神经元都表现为方波簇放电, 峰相位差的最大值103.3 ＜max(∆ϕ(t))＜ 103.6, 未达到峰相位同步; 簇相位差最大值 max(∆(Φ(t))) ≈ 3.14 ＜ 2π , 达到簇反相同 步; 相 关 系 数 | ρ |<0.2 , 说 明 同 步 程 度 弱.以gsyn-e=0.35nS为例, 此时两耦合神经元都表现为周期18方波簇放电序列, 簇相位差最大值max(∆Φ(t))= 3.14 ＜ 2π , 峰相位差的最大值max(∆ϕ(t))=103.3, 相关系数 ρ = —0.02.神经元1(红)和2(蓝)放电节律(上)及其对应的耦合电流(下)如图7(a)所示.

第 二 部 分(0 .45nS ≤ gsyn-e< 2.0nS): 两 耦 合神经元都表现了方波簇放电模式, 峰相位差的最大值在大参数范围内为 2π ＜ max(∆ϕ(t)) ＜ 50.2,相比第一部分有降低, 簇相位差最大值max(∆(Φ(t)))<0.07＜ 2π , 说明未达峰相位同步、达到了簇相位同步; 相关系数 0.6 ＜ ρ ＜ 0.8, 说明同步程度达到中等.以耦合强度 gsyn-e= 1.5 nS为值 max(∆ϕ(t)) = 910.9, 簇 相 位 差 最 大 值max(∆Φ(t))= 0.06 ＜ 2π , 相关系数 ρ =0.65 .神经元1(红)和2(蓝)放电节律(上)及其对应的耦合电流(下)如图7(c)所示.

图6 初值相同时, 不同耦合强度下神经元1(红)和2(蓝)的膜电位 V (上)及耦合电流 (下), 插图是局部放大 (a)gsyn-e=0.35nS; (b) gsyn-e=2.5nS ; (c) gsyn-e=5.0nS ; (d)gsyn-e=18.0nSFig.6.Membrane potential V (top) and coupling current (low) of neurons 1 (red) and 2 (blue) with the same initial values at different gsyn-e values (Insert figure: the enlargement of bursting): (a) gsyn-e=0.35nS ; (b) gsyn-e=2.5nS ; (c) gsyn-e=5.0nS ;(d) gsyn-e=18.0nS .

第四部分(3 .55nS ≤ gsyn-e< 10.5nS): 耦合神经元表现出DB簇放电模式, 峰相位差的最大值2π ＜ max(∆ϕ(t)) ＜ 653.4, 簇相位差最大值max(∆(Φ(t)))＜ 0.02 ＜ 2π , 相关系数 0.97 ＜ ρ ＜1, 未达到峰相位同步, 而达到簇相位同步, 同步程度强.以 gsyn-e= 5.0 nS为例, 此时簇相位差最大值 max(∆Φ(t)) = 0.99 ＜ 2 π , 峰相位差的最大值max(∆ϕ(t))= 157.9, 同步因子 ρ =0.99 .神经元1(红)和2(蓝)放电节律(上)及其对应的耦合电流(下)随时间t的变化如图7(d)所示.

第五部分(1 0.55nS ≤ gsyn-e≤ 1 8.0nS): 两耦合例, 两神经元都表现为周期23方波簇放电序列,簇相位差最大值 max(∆Φ(t))=0.02 ＜ 2π , 峰相位差的最大值 max(∆ϕ(t)) = 4.1, 相关系数 ρ =0.64 .放电时间序列及耦合电流 Isyn-e随时间t的变化如图7(b)所示.此时两耦合神经元的簇内的峰是反相同步的, 如图7(b)的插图所示, 与Best等[29]的研究类似.图7(b)的插图所示的是单个簇的放大图.

第 三 部 分(2 .0nS ≤ gsyn-e< 3.55nS): 每 个 神经元的电活动随时间变化都表现为方波簇与DB簇的交替的混合簇, 峰相位差的最大值 2π ＜max(∆ϕ(t))＜ 2013.1, 簇 相 位 差 最 大 值max(∆(Φ(t)))<0.08＜ 2π , 相关系数 0.9 ＜ ρ ＜ 1,未达到峰相位同步, 而达到簇相位同步, 同步程度较强.例如, 当 gsyn-e= 2.5 nS时, 峰相位差的最大神经元表现为周期1峰放电, 峰相位差最大值max(∆(ϕ(t)))≈3.14, 达到峰反相同步; 相关系数−0.88＜ ρ ≤ − 0.6 .由于此时不是簇放电, 就没有定义簇相位差.以 gsyn-e= 18.0 nS为例, 峰相位差最大值 max(∆ϕ(t)) = 3 .14 , 同步因子 ρ = —0.88,达到峰反相同步.放电时间序列图及耦合电流Isyn-e随时间t的变化如图7(e)所示.

图7 初值不同时, 不同耦合强度下神经元1(红)和2(蓝)的膜电位V(上)及耦合电流 I syn-e (下), 插图是局部放大 (a) gsyn-e =0.35 nS; (b) g syn-e = 1.5 nS (c) g syn-e = 2.5 nS; (d) g syn-e = 5.0 nS; (e) g syn-e = 18.0 nSFig.7.Membrane potential V (top) and coupling current I syn-e (low) of neurons 1 (red) and 2 (blue) with different initial values at different gsyn-e (Insert figure: the enlargement of bursting): (a) gsyn-e = 0.35 nS; (b) gsyn-e = 1.5 nS; (c) gsyn-e = 2.5 nS;(d) gsyn-e = 5.0 nS; (e) g syn-e = 18.0 nS.

其中方波簇、方波簇与DB簇的交替的混合簇和周期1峰放电与Dunmyre等的实验发现相似[25].

3.3 耦合神经元的快慢变量分离

在不同耦合强度下, 对两耦合神经元系统(4)—(7)进行快慢变量分离.本文以神经元1为例展示结果.

3.3.1 快子系统的分岔的示例

以 gsyn-e= 1.5 nS为例(其它耦合强度下的快子系统的分岔的关键点h值, 如表3所列), 随着慢变量h的增加, 快子系统的平衡点在(h, V1)平面上展示了“S”型曲线, 其中下支和上支的蓝色实线分别为稳定结点和稳定焦点, 黑色虚线是不稳定平衡点, “S”型曲线上的分岔点 F1 都代表平衡点的鞍结分岔, 如图8(a)所示.

在“S”型曲线的上支, 随着慢变量h增加, 不稳定焦点变为稳定焦点, 出现了两个不稳定极限环, 发生了两个亚临界Hopf分岔即subh1点(h≈0.2682)和 subh2点 (h ≈ 0.2857).subh2点的不稳定极限环随着h增大消失, 而subh1点的不稳定极限环随着h的增大, 在h ≈ 0.4549发生了极限环的鞍结分岔(即LPC1), 有稳定的极限环出现,再随着h减小稳定的极限环在h ≈ 0.3321发生了极限环的鞍结分岔(即LPC2), 同时产生了不稳定的极限环, 接着随着h的增大不稳定的极限环消失, 如图8(b)所示, 红色实线和虚线分别代表稳定及不稳定极限环的最大值和最小值.

3.3.2 相同初值条件下的簇放电模式

当耦合强度 gsyn-e= 0.35 nS时, 神经元1处于周期19的方波簇放电状态(如图6(a)所示), 将其簇放电在(h, V1)平面上的轨迹叠加到其快子系统的分岔图上, 如图9(a)所示.随着h的增加, 全系统的轨线从下支的静息态经鞍结分岔(F1)转迁到上支的重复放电, 再随着h的减小, 重复放电经鞍−同宿(HC)分岔转为下支的静息态, 完成了一次周期振荡.根据文献[52]中的簇放电的命名规则, 此类簇放电被称为“fold/homoclinic”型簇放电, 如图9(a)所示.

当 gsyn-e= 2.5 nS时, 神经元1处于DB簇放电状态(如图6(b)所示), 如图9(b)所示.簇放电的相轨线从下支静息态经鞍结分岔(F1)转迁到上支, 但没有围绕着稳定的极限环发生振荡, 放电围绕稳定焦点作衰减振荡, 当经过亚临界Hopf分岔点subh1时, 振幅逐渐减小但没减少到零, 但在subh1点左侧不稳定焦点的作用下, 系统振荡幅度逐渐增大, 越过极限环的鞍结分岔点, 最后上支的放电态转变为下支的静息态, 所以称此放电簇为“fold/fold limit cycle”滞后的“subHopf/subHopf”型簇放电.

表3 不同 g syn-e 下快子系统中关键点的慢变量h的值Table 3.The slow variable h values of the bifurcation or key points at different g syn-e values.

图8 当 gsyn-e = 1.5 nS时, 两耦合神经元的快子系统的分岔, 插图是局部放大 (a)平衡点分岔; (b)平衡点分岔和极限环的分岔Fig.8.Bifurcations of the fast−subsystem of the two coupled neurons with respect to h when gsyn-e = 1.5 nS (Insert figure: the en−largement): (a) Equilibrium points; (b) equilibrium points and limit cycle.

当 gsyn-e= 5.0 nS时, 神经元1仍然处于DB簇放电状态(如图6(c)所示), 如图9(c)所示.簇放电的相轨线从下支的静息态经鞍结分岔(F1)转迁到上支, 绕着稳定的极限环发生振荡, 而放电轨线围绕稳定焦点作阻尼振荡, 当经过亚临界Hopf分岔点subh1时, 振幅减小到零, 但在subh1点左侧不稳定焦点的作用下, 系统振荡幅度逐渐增大, 越过极限环的鞍结分岔点, 最后上支的放电转变为下支的静息, 所以称此放电簇为经“fold/fold limit cycle”滞后环的“subHopf/subHopf”型簇放电.

当 gsyn-e= 18.0 nS时, 神经元1都处于周期1峰放电状态(如图6(d)所示), 如图9(d)所示.相轨线未经“S”型曲线的下支, 这就是峰放电的特征.

3.3.3 不同初值条件下的簇放电模式

不同初值条件下, 神经元1和神经元2的时间历程不同, 但是两个神经元的放电簇模式还是相同的.限于篇幅和避免重复, 本文展示神经元1的结果.

gsyn-e= 0.35 nS, 神经元1处于周期18方波簇放电状态(如图7(a)所示), 图10(a)所示.虽然与初值相同时该耦合强度下耦合系统的分岔类似(如图9(a)所示), 此类型的簇放电被称为“fold/homoclinic”型簇放电, 但是此时的方波簇的峰个数少, 如图10(a)所示.

gsyn-e= 1.5 nS, 神经元1处于周期23方波簇放电状态(如图7(b)所示), 如图10(b)所示.簇放电轨线从静息态经鞍结分岔(F1)转迁到放电态,放电态经极限环的鞍结(LPC2)分岔进入静息态,形成“fold/fold limit cycle”型簇放电, 如图10(b)所示.型簇放电, 如图10(c)所示; 而对于方波簇放电时段, 与初值不同时耦合强度 gsyn-e= 1.5 nS下耦合系统的分岔 (如图10(b))类似, 形成 “fold/fold limit cycle”型簇放电, 如图10(d)所示.

图9 初值相同时, 神经元1在不同耦合强度下簇放电模式的快慢变量分离, 插图是局部放大 (a) gsyn-e = 0.35 nS; (b) gsyn-e =2.5 nS; (c) gsyn-e = 5.0 nS; (d) gsyn-e = 18.0 nSFig.9.The fast−slow variable dissection of neuron 1 for different initial values at different gsyn-e values (Insert figure: the enlarge−ment): (a) g syn-e = 0.35 nS; (b) gsyn-e = 2.5 nS; (c) gsyn-e = 5.0 nS; (d) gsyn-e = 18.0 nS.

图10 初值不同时, 神经元1在不同耦合强度下簇放电模式的快慢变量分离, 插图是局部放大 (a) gsyn-e = 0.35 nS; (b) gsyn-e =1.5 nS; (c)和(d) gsyn-e = 2.5 nS; (e) gsyn-e = 5.0 nS; (f) gsyn-e = 18.0 nSFig.10.The fast−slow variable dissection of neuron 1 for different initial values at different g syn-e values (Insert figure: the enlarge−ment): (a) gsyn-e = 0.35 nS; (b) gsyn-e = 1.5 nS; (c) and (d) gsyn-e = 2.5 nS; (e) gsyn-e = 5.0 nS; (f) gsyn-e = 18.0 nS.

gsyn-e= 5.0 nS, 神经元1处于DB簇放电状态(如图7(d)所示), 如图10(e)所示.快子系统下支的稳定状态在平衡点的分岔与初值相同时该耦

gsyn-e= 2.5 nS, 神经元1处于DB簇与方波簇的混合簇放电状态(如图7(c)所示), 如图10(c)和图10(d)所示.其中的DB簇的放电时段与方波簇放电时段的快慢变量分离的结果分别如图10(c)和图10(d)所示.DB簇放电时段的结果与初值相同时该耦合强度下耦合系统的分岔类似, 形成“fold/fold limit cycle”滞后环的“subHopf/subHopf”合强度下耦合系统的分岔类似, 此放电簇为经“fold/fold limit cycle”滞后环的“subHopf/subHopf”型簇放电.

gsyn-e=18.0nS, 神经元1处于周期1峰放电状态(如图7(e)所示), 如图10(f)所示.簇放电轨线未经“S”型曲线的下支, 快子系统分岔与初值相同时该耦合强度下耦合系统的分岔相同, 称为峰放电, 但是两个峰放电被诱导的原因不同.

3.4 共存的反相同步峰和同相同步峰的不同

图11 反相同步(紫色)和同相同步(绿色)周期1峰放电节律 (a) (h, V1)相平面上的相轨迹图; (b)耦合电流随时间t的变化Fig.11.The anti−phase (purple) and in−phase (green) period−1 spiking: (a) The V-h trajectory; (b) coupling current.

图12 (a)快子系统的平衡点和极限环的分岔; (b)图(a)中极限环分岔处的放大; (c)反相同步(紫色)和同相同步(绿色)周期1峰放电的快慢变量分离; (d)图(c)中反向同步(紫色)和同向(绿色)同步周期1峰放电的放大Fig.12.(a) Bifurcations of equilibrium points and limit cycle of the fast−subsystem; (b) enlargement of (a); (c) fast−slow variable dissection of anti−phase (purple) and in−phase (green) period−1 spiking; (d) enlargement of anti−phase (purple) and in−phase (green)period−1 spiking in Fig.(c).

当耦合强度 gsyn-e= 18.0 nS时, 相同初值时,耦合神经元表现为完全同步周期1峰放电节律(如图6(d)所示), 也是同相峰同步; 初值不同时, 耦合神经元表现出反相同步的周期1峰放电节律(如图7(e)所示).两种放电节律的峰峰间期大小及极限环的幅值不同; 相比于不同初值, 相同初值下的峰峰间期和极限环的幅值都大, 如图11(a)所示.且耦合神经元的耦合流大小不同, 相同初值下的耦合流大于不同初值下的耦合流, 如图11(b)所示.这说明: 当 − 0.099≤h≤0.089 时, 耦合神经元存在共存的周期1峰吸引子.当初值相同时, 耦合神经元表现为同相同步周期1峰放电, 当初值不同时,耦合神经元表现为反相同步周期1峰放电.

利用快慢表变量分离可以进一步区分反相同步和同相同步的周期1放电节律的不同, 如图12所示.可以看出两类峰放电产生的机制不同.快子系统会产生一个大的稳定极限环(粗红线)和一个小的稳定极限环(红线), 如图12(a)所示.图12(b)是极限环分岔的局部放大.现将两共存的周期1放电在(h, V1)平面的相轨线叠加到相应的快子系统的分岔图上, 同相同步的周期1峰放电(绿色)是由快子系统的大的稳定的极限环诱导产生, 而反相同步的周期1峰放电(紫色)是由快子系统的比较小的稳定的极限环诱导产生, 如图12(c)所示.图12(d)是两周期1的极限环处放大.因此同相同步的峰峰间期和极限环的幅值都大于反相同步的峰峰间期和极限环的幅值.

4 结 论

本文研究了pre−BötC单神经元模型的放电规律以及耦合神经元模型由簇到峰的复杂同步节律转迁和分岔机制.在单神经元模型中模拟了由DB簇、到方波簇、到周期1簇放电的节律转迁过程, 与Negro等[26]在实验中发现的pre−BötC神经元会经历从静息到簇放电、再到峰放电的节律转迁历程相近似.对于由pre−BötC的单神经元经兴奋的化学突触耦合建立的两神经元网络, 随着耦合强度增加, 会表现出不同同步节律的转迁.初值相同时, 表现出完全同步的方波簇放电, 到DB簇放电,到周期1峰放电的复杂节律转迁.当两耦合神经元初值不同时, 表现出由簇相位反相同步的方波簇放电, 到簇相位同步但簇内峰反相同步的方波簇放电, 到方波簇与DB簇的混合簇放电, 最后到反相同步的周期1峰放电的复杂节律、同步转迁.其中方波簇, DB簇, 方波簇与DB簇交替的混合簇和周期1峰放电与Dunmyre等[25]的实验发现相似.利用快慢变量分离, 揭示出相同初值下的同步行为依次为“fold/homoclinic”, “subHopf/subHopf”簇放电和周期1峰放电.不同初值时的簇相位同步转迁表现为由“fold/homoclinic”, 到“fold/fold limit cycle”, 到“subHopf/subHopf”与“fold/fold limit cycle”的混合簇放电、再到“subHopf/subHopf”簇放电、最后到反相同步周期1峰放电.这些复杂的节律模式的同步转迁与实验结果相类似, 为全面深入认识实验现象提供了帮助.与已有的关于pre−BötC复合体的同步的数值仿真结果不同, 本文不仅结合了生物学实验结果, 还采用了多个指标刻画了多种同步, 并且重点研究了耦合神经元同步模式随着参数(耦合强度)连续变化的节律转迁规律,不再是只关注某个参数值下的同步模式.

传统观念中, 兴奋性的耦合神经元网络, 兴奋性强度的增加可以促进网络中神经元的同步, 而抑制性神经元网络, 抑制性强度的增加可以抑制网络中神经元同步[53].而本文两耦合神经元网络, 随着兴奋性的增加, 在强耦合强度下不仅发现了两耦合神经元的周期1峰放电的完全同步, 而且还发现了两耦合神经元的周期1峰放电的反向同步, 给出了一种不同于传统观念的新示例.进一步利用快慢变量分离及分岔机制揭示出了同向同步周期1峰放电及反向同步周期1峰放电的不同: 同相同步的周期1峰放电是由快子系统的大的稳定的极限环诱导产生, 而反相同步的周期1峰放电是由快子系统的比较小的稳定的极限环诱导产生.与传统观念不同的新的非线性现象和机制引起了越来越多的关注[44−50], 包括, 抑制性作用引起单神经元从静息变为放电、增强放电频率、引起神经元网络的同相同步; 兴奋性作用可以引起放电频率降低、引起神经元网络从放电变为静息.而本文发现的强兴奋性的强度下的反相同步是一个新的反常现象, 并用分岔给出了理论解释, 丰富了非线性动力学的内涵.

本文及以前的关于pre−BötC的复杂节律及其同步动力学的理论研究[28−40], 有助于从理论上全面和深入研究复杂的、相对离散的实验现象.特别是, pre−BötC的每个神经元的节律就是涉及多时间尺度的、强非线性的复杂动力学, 再加上强非线性的兴奋性化学耦合的影响, 不仅造成了实验现象复杂、难以识别, 还造成了理论研究上的困难, 这就使得pre−BötC的同步动力学的研究比一般的复杂网络同步动力学的研究要更为困难.因此, 快慢变量分离和分岔分析的介入对于这些复杂节律的认识是极为必要的[54−56].今后, 应该继续发挥快慢变量分离和分岔分析的优势, 从多个方面(例如不同离子、考虑自耦合等)开展对pre−BötC的同步动力学的研究, 服务于神经科学和非线性科学的发展.