两耦合HR混沌神经元同步研究

吴望生　(长江大学物理与光电工程学院，湖北 荆州 434023)

唐国宁　(广西师范大学物理科学与技术学院，广西 桂林 541004)



两耦合HR混沌神经元同步研究

吴望生(长江大学物理与光电工程学院，湖北 荆州 434023)

唐国宁(广西师范大学物理科学与技术学院，广西 桂林 541004)

[摘要]基于Hindmarsh-Rose(HR)神经元动力学模型，通过数值计算耦合神经元系统的同步差和峰峰间期(ISI)随耦合强度的变化研究了2个耦合混沌HR神经元在近邻反馈耦合下的同步过程。结果表明，初态混沌的HR神经元在耦合中存在着复杂的同步过程；伴随耦合强度的增大，神经元会交替出现混沌发放、周期发放等放电现象；在达到精确的完全同步之前，会交替出现簇同步、近似同步和完全同步等同步现象；在弱耦合条件下，神经元会出现反同步和去同步。

[关键词]Hindmarsh-Rose神经元；同步；反同步

自17世纪Huygens发现同步现象[1]后的很长一段时间里，有关同步的研究始终是建立在周期运动的基础之上。20世纪90年代，Pecora和Carro[2,3]首次在实验室的电子线路上实现混沌同步，从此混沌同步研究在许多学科领域深入展开。特别是在神经科学领域，同步成为当前研究的热点问题。生理实验证实，大脑中存在着神经元系统的局部同步模式。1989年Gray等[4]发现，当给予同一个物体的轮廓作为共同刺激时，分布在相互分离模块的视觉皮层细胞常常会产生同步化放电；1997年Riehle等[5]在猴的运动皮质区发现不同神经元之间出现了峰同步，证实了大脑局部同步区的存在。结合生物神经系统中发现的许多同步现象，Pouget等[6]提出神经系统的大量信息是由神经元系统共同处理。神经元系统是一个由快慢变量组成的多尺度系统，慢变量是神经元的簇行为产生的机制[7]，导致神经元会呈现丰富的发放行为，形成各种周期和混沌模式，因而耦合神经元系统的同步过程也非常复杂。研究表明，神经元在达到稳定的完全同步之前，存在着不同类型的同步过程。胡岗等[8]研究了广义同步和相位同步之间的同步过程，指出相位同步可能优先也可能落后于广义同步；Dhamala等[9]发现耦合神经元达到完全同步前，会出现神经元簇的同步和峰的同步。目前来看，对不同类型的同步过程仍有待进一步研究，分析耦合神经元系统的复杂同步过程对理解现实的神经元同步行为具有重要的指导意义。为此，笔者选用近邻反馈耦合，研究了2个初态混沌的Hindmarsh-Rose(HR)神经元的同步过程。

1动力学模型

2个初态混沌的耦合HR神经元的动力学方程[10]如下：

(1)

(2)

(3)

(4)

(5)

(6)

其中，角标1(或2)分别代表HR神经元1(或2)，统一简写为HR1和HR2； u代表神经细胞的膜电位；v是神经元内部的快速恢复电流；w是慢变调节电流(u和v是快变量，w是慢变量，共同构成一个双尺度系统)；a、b、c、d、γ、s及χ都是系统参数；Iext表示外界直流激励； k是耦合强度。各参数取值如下：a=1，b=3，c=1，d=5，γ=0.006，s=4，χ=-1.6和Iext=3.0，去暂态后单个神经元进入混沌初态。

在数值模拟中，采用四阶龙格库塔法求解耦合方程组，积分步长Δt=0.001。神经元实现精确的完全同步需要时间依赖耦合强度，为保证2个神经元基本达到同步的渐进态，总积分时间取为10000时间单位。同时为了保证数值计算的准确性，需通过消除一定时间单位的暂态过程而得到稳态结果。笔者通过对2个神经元赋不同的随机初值(-0.5,0.5)，去除500时间单位的暂态，让各神经元单独演化进入混沌初态。图1展示了去暂态时间段内HR神经元的动力学行为，图1(a)显示暂态末期2个神经元先后进入混沌发放模式，图1(b)可看出2个神经元在(ui,wi)相平面上的轨道先后进入混沌吸引子。

图1　去暂态时间段内HR神经元的动力学行为

为了研究2个神经元的同步程度，将同步差ei(t)及平均同步差δi定义为：

eu(t)=|u1(t)-u2(t)|

(7)

ew(t)=|w1(t)-w2(t)|

(8)

(9)

其中，eu(t)描述快变量的空间同步程度；ew(t)描述慢变量的空间同步程度。即eu(t)描述了t时刻2个神经元膜电位u的同步差，ew(t)描述了t时刻2个神经元慢变调节电流w的同步差。同步差越小，表明同步程度越高。通常同步差是振荡式下降，根据其大小，一般不容易判断不同参数下神经元的同步程度，所以笔者用δi描述去掉同步暂态后的一段时间内神经元之间的平均同步情况。不失一般性，选择时间范围[9000,10000]。当δi=0时，2个神经元实现精确的完全同步，然而有限的时间内δi一般是逐步趋于零，并且快变量δu和慢变量δw的同步速率不同。笔者设定，当快变量和慢变量均小于等于0.001时，即δi≤10-3时，则2个神经元达到完全同步；当快变量和慢变量其中之一小于等于0.001时，即δu≤10-3或δw≤10-3时，2个神经元达到近似同步。

2数值模拟

采用近邻反馈耦合，对耦合混沌的神经元的峰峰间期(ISI)和同步差进行数值计算。图2给出了HR神经元1的峰峰间期(ISI)随耦合强度k变化的分岔图，图3给出了耦合混沌神经元快变量的平均同步差δi随耦合强度k的变化关系。由图2和图3可以看出，随着耦合强度k的逐步增大，HR神经元呈现丰富的放电活动。随着耦合强度的增加，图2被交替分割成多个带状密集区和线型轨迹区，表明混沌耦合的神经元在混沌发放和周期发放2个状态反复切换，直至进入混沌发放达到完全同步。

图2　HR1的ISI随k变化的分岔图　　　　　　　　　　　　　　图3　平均同步差δi随k的变化

图2的第1段线型轨迹区对应耦合强度区间0.008≤k≤0.028，表明神经元HR1由开始的混沌发放模式变成了规则的周期发放，神经元HR2-ISI分岔图在此区间也有同类现象，这表明耦合神经元发生了簇同步。不失一般性，取耦合强度k=0.024来分析神经元的簇同步状态，图4展示了2个耦合神经元膜电位的时间历程和在(u1,w1)和(u1,u2)平面上所对应的相图。结合图4(a)和(c)可以发现，2个耦合神经元的周期簇同时发生，出现了簇同步；而图4(b)显示峰错位和图4(d)中相位定位在二四象限，都表明在神经元簇里的峰是不相关的，未发生峰同步。

图4　k=0.024时，耦合神经元簇同步状态的分析图

从图2中可以清楚的看到，当耦合较弱且参数范围较小时，耦合强度可以抑制混沌而变为周期的发放行为。但结合图3和图4，虽然神经元由于耦合强度增大出现倍周期轨道，并出现簇同步，但2个神经元之间并未达到峰同步，直到k大于0.42后，2耦合神经元才出现峰同步现象。这说明混沌发放被抑制变为周期发放，出现簇同步，只是峰同步前的一个中间过程，还需要继续增大耦合强度，神经元之间才有可能出现峰同步，并最终达成完全同步。

随着耦合强度逐步增强，周期发放的神经元又进入混沌状态，且不同步，但当0.151≤k≤0.180时，2个耦合神经元随耦合强度的变化交替出现倍周期轨道、混沌轨道，不同步和反同步等复杂的过渡模式。具体来说，当0.151≤k≤0.157时，神经元会出现周期8的倍周期轨道，且出现反相的周期性同步状态。下面，取k=0.151来分析倍周期和反同步现象(见图5)。

图5　k=0.151时，耦合神经元反同步状态的分析图

从图5(a)中可以看到，耦合神经元HR1出现了周期8的倍周期轨道，比图4(c)中周期数增加了，虽然从混沌轨道变成了周期轨道，但2个耦合神经元并未出现如图4(a)和(b)所示的簇同步，而是出现了反同步状态；从图5(b～d)中可以清楚看到神经元3个快慢变量均处于反相的周期性同步状态。经过分析，反同步所对应的平均同步差δi均出现较大的变化，如膜电位同步差δu>1就会出现反同步，这正好解释了图3在0.151≤k≤0.180区间内出现的突然增大的同步差峰包，其δu值均大于1，说明在峰包对应的耦合强度范围内，神经元出现了反相的周期性同步。

当耦合强度离开这个反同步区，神经元的反同步被破坏，神经元重新由周期轨道进入混沌轨道，直到耦合强度k≥0.422时，2个混沌耦合神经元又出现近似同步。随着耦合强度的的进一步增加，当耦合强度位于0.422≤k≤0.434时，耦合神经元在这个区间存在着复杂的过渡模式。

图6　对数化后平均同步差δi随k的变化

具体的表现行为通过图6可以分析清楚。将同步差δi对数化后可清晰看出，随着耦合强度的增加，同步差震荡式下降直至为0。蓝色虚线表示δi=10-3，为完全同步分界线，线下为完全同步区。在δi与蓝色虚线交界处，即耦合强度位于0.422≤k≤0.434时，有时δu<10-3而δw>10-3或δu>10-3而δw<10-3，说明出现了近似同步；有时快慢变量的同步平均差均小于0.001，即δi<10-3，发生了完全同步。上述情况表明，耦合强度较大时，神经元在完全同步之前，即可能出现峰同步(即近似同步)，也可能出现簇和峰均同步(即完全同步)，而且上述同步现象随耦合强度的变化而交替出现，表明耦合混沌神经元同步之前的过渡状态中包含着复杂的放电信息。

为了进一步分析这些复杂的同步过程，考察完全同步和近似同步的不同之处，图7在不同的耦合强度下研究神经元快慢变量的各种相图，其中黑色轨迹对应k=0.423，红色轨迹对应k=0.433。图7(a)和(b)显示，黑色轨迹的相位与一三象限角平分线完全重合，表明k=0.423时，神经元出现了完全同步。红色轨迹的相位则定位在一三象限角平分线附近的小邻域内，这表明k=0.433时，神经元产生了近似同步且处于混沌状态。需要说明的是，在耦合强度0.422≤k≤0.434的这个区间内，没有发现之前弱耦合状态下，神经元的峰不同步而簇同步的周期放电现象，神经元不论同步与否，均处于混沌放电状态。

图7　神经元在不同耦合强度下的平面相图

从图3和图6中2条曲线完全归零的位置可查出，当耦合强度k≥0.474时，平均同步差δi=0，这时2个耦合神经元达到了完全的峰同步和簇同步。图8显示，当k=0.474时，2个耦合神经元的混沌放电时间历程图8(a)和混沌吸引子图8(b)完全重合，表明神经元达成完全同步后仍然保持了稳定的混沌状态，结合图1说明初态随机的2个混沌神经元通过近邻反馈耦合达成了精确的完全同步。

图8　k=0.474时，耦合神经元完全同步状态的分析图

3同步过程分析

HR神经元是一个具有快慢变量的双尺度系统，其中快变量产生了神经元的峰，慢变量导致了神经元的簇的形成，在不同的参数条件下，呈现出丰富的动力学行为，如峰峰间期分岔、混沌吸引子等。下面，笔者来比较不同耦合条件下神经元在(u,w)相平面上的吸引子变化情况。

图1(b)显示了未耦合状态下，神经元的混沌吸引子；图4(c)展示耦合强度k=0.024时，神经元HR1的周期5轨道；图5展现了k=0.151时，神经元HR1的周期8轨道；图8则表明神经元达到精确同步后仍然处于混沌吸引子中。随着耦合强度的逐步增大，耦合神经元吸引子轨道由混沌先后过渡到周期5、周期8，然后又进入混沌状态，直至以混沌状态进入精确的完全同步。

对于这种吸引子变化的原因可以用文献[11]中耦合的Lorenz系统的动力变化的解释：在弱耦合的情况下，混沌的耦合破坏了原来Lorenz系统双圈吸引子来回切换所遵循的一致性，结果造成原来吸引子的破坏而产生了新的更复杂的吸引子。上述解释对于Lorenz系统是成立的，对同样具有多圈吸引子的其他耦合系统也是成立的，因此在笔者所研究的HR神经元模型中，神经元在给定的初始参数下呈现了如图1(b)所示的多圈混沌吸引子，随着耦合强度的增加，混沌轨道被耦合破坏而进入倍周期轨道(见图4(c)和图5(a))，继续增大耦合强度，神经元又脱离周期轨道而进入混沌，直至最终达成精确的完全同步仍然保持混沌状态。

耦合的多尺度系统的混沌同步与快慢子系统的动力行为是紧密相关的，耦合神经元达到完全同步前，会出现神经元簇的同步和峰的同步[9]。笔者研究发现，耦合混沌的HR神经元随着耦合强度的增加，会在弱耦合区间出现簇同步，即慢变量的同步(见图4(a)和(b))，慢变量的同步随耦合强度继续增加而被破坏，直至耦合强度大到能使神经元出现峰的同步和簇的同步，即完全同步，意味着快慢变量都达到同步。因此，2个耦合混沌的HR神经元采用近邻反馈耦合，随着耦合强度的增加，会先后出现簇同步、峰同步直至完全同步。耦合强度的增加除了可以改变神经元的吸引子，还可以影响神经元的放电节律形成反同步。如图5所示，当在某些耦合强度区间内出现周期性吸引子时，2个神经元的快慢变量如膜电位和慢电流会出现反相的周期性同步，事实上在这个区间附近，2个耦合神经元随耦合强度的变化交替出现倍周期轨道、混沌轨道，不同步和反同步等复杂的过渡模式。

4结论

研究了采用近邻反馈耦合2个混沌HR神经元的同步问题。在数值模拟中，采用四阶龙格库塔法求解耦合方程组，通过计算耦合神经元系统的同步差和峰峰间期随耦合强度的变化，发现在耦合混沌的HR神经元模型中，存在着复杂的同步过程。

1)随着耦合强度的增大，神经元的吸引子会被改变，会交替出现混沌发放、周期发放等放电现象，表明混沌轨道被耦合破坏而进入倍周期轨道，继续增大耦合强度，神经元又脱离周期轨道而进入混沌，直至最终达成精确的完全同步仍然保持混沌状态。

2)在达到精确的完全同步之前，随着耦合强度的增大，神经元会交替出现混沌发放、周期发放等放电现象，神经元系统会依次出现神经元簇的同步、峰的同步，最终达成完全同步。

3)耦合强度的增加除了可以改变神经元的吸引子，影响神经元的放电节律，甚至在某些弱耦合区间，会导致神经元系统出现反相的周期性同步，抑制同步的发生。

[参考文献]

[1]Huygens C.Horologium oscillatorium [M]. Parisiis, France, 1673.

[2] Pecora L M, Carroll T L. Synchronization in chaotic systems [J]. Physical Review Letters, 1990, 64(8): 821～824.

[3]Pecora L M, Carroll T L. Driving systems with chaotic signals[J]. Physical Review A, 1991, 44(4): 2374～2383.

[4]Gray C M, König P, Engel A K, et al. Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties[J]. Nature, 1989, 338(6213): 334～337.

[5]Riehle A, Grün S, Diesmann M, et al. Spike synchronization and rate modulation differentially involved in motor cortical function[J]. Science (New York, N.Y.), 1997, 278(5345): 1950～1953.

[6]Pouget A, Dayan P, Zemel R. Information processing with population codes[J]. Nature Reviews Neuroscience, 2000, 1(2): 125～132.

[7]王青云, 石霞, 陆启韶. 神经元耦合系统的同步动力学[M]. 北京: 科学出版社, 2008: 38～50.

[8] Zheng Z G, Hu G. Generalized synchronization versus phase synchronization[J]. Physical Review E, 2000, 62(6): 7882～7885.

[9] Dhamala M, Jirsa V K, Ding M Z. Transitions to synchrony in coupled bursting neurons[J]. Physical Review Letters, 2004, 92(2): 02810.

[10]Hindmarsh J L, Rose R M. A model of neuronal bursting using three coupled first order differential equations[J]. Proceedings of the Royal Society of London Series B, 1984, 221(1222): 87～102.

[11]Liu Z H, Lai Y C, Matias M A. Universal scaling of Lyapunov exponents in coupled chaotic oscillators. [J]. Physical Review E, 2003, 67: 045203(R).

[编辑]张涛

[文献标志码]A

[文章编号]1673-1409(2016)07-0022-06

[中图分类号]O415.5

[作者简介]吴望生(1979-)，男，硕士，讲师，现主要从事统计物理方面的教学与研究工作；E-mail:wshwu@yangtzeu.edu.cn。

[基金项目]国家自然科学基金项目(11165004)。

[收稿日期]2015-11-27

[引著格式]吴望生，唐国宁.两耦合HR混沌神经元同步研究[J].长江大学学报(自科版)，2016,13(7)：22～27.