基于一种自适应突触学习的动态相关系数与相位同步分析

严传魁，王如彬

（1.杭州师范大学理学院，浙江 杭州310036；2.华东理工大学理学院认知神经动力学研究所，上海200237）

随着科技的发展，人们发现同步在工程领域应用广泛，在此背景下网络同步研究取得了较大的进展［1-7］，通过对系统状态变量的控制，使系统达到所需要的同步.另外，在大脑神经系统的生理实验过程中，科学家们发现大脑中也存在神经元同步放电的模式［8-10］，这种放电模式可能与大脑的某些功能相关，所以神经系统的同步研究对人们了解大脑意义重大.

神经网络系统的同步研究，也经历了从两个耦合神经元到网络层次的发展，研究中常见的网络拓扑结构一般以规则网络为主，如链式，网格式，环式［1-3，7］.另外，还有更符合大脑神经元解剖结构的小世界网络连接的神经元系统［11-13］.但是这些网络的规则性导致它们与实际的生物网络具有较大的差异，实际上生物网络的突触连接具有更大的随机性.而且目前的大部分研究受工程应用中动力系统的影响，还有个共同的特征就是耦合的对称性，即结点间的耦合作用是相互的，但是实际神经元网络的突触连接是单向的.本文的研究将这两个问题做改进和推广.

既然神经网络的同步放电可能实现大脑的某些功能，那么一个很自然的问题就是这种同步是如何达到的？对现实中动力系统的研究发现，可以通过对系统中的某些状态变量进行外部控制使得系统达到同步［14-16］，但是对于神经系统这显然不是必要的.大脑具有自组织性，神经元之间的突触连接具有可塑性，神经系统通过突触学习和适应达到想要的结果.本文提出了一种自适应的突触学习规则来模拟生物系统中时间点依赖的突触可塑性.

1 模型介绍

1.1 网络模型

进行神经网络的研究带来的一个很基本问题是网络的拓扑结构是怎样的.链式、环式、网格式等规则网络完全无法体现神经系统网络结构的复杂性.因为神经网络中神经元的突触连接具有方向性，并且连接不具有规则图形的对称性.更接近于神经网络的是NW 小世界网络，实际神经系统的突触连接应该介于完全连接和完全不连接之间，小世界网络具有一定的随机性和不规则性，所以神经网络应该是具有小世界网络特性的.

考虑神经系统突触连接的一些特性，用一个经过改造的NW 小世界网络来模拟神经网络.在对NW小世界网络节点进行随机连接时由原来的对称性耦合变成单向性突触连接，网络构造步骤如下：

1）为了保证网络的连通性，让系统中的n个神经元进行环式连接，每个节点只和相邻的两个点对称连接；

2）对于网络中的任意两个不相邻的节点：第i和第j个节点神经元，产生两个［0，1］上均匀分布的随机数pij，pji；

3）然后对这两个节点进行两次概率事件的判定，不妨设小世界的连接概率为p.如果pij＜p表示网络连接中第i个神经元激励第j个神经元，否则表示第i个神经元没有激励第j个神经元；反过来，若pji＜p表示网络连接中第j个神经元激励第i个神经元，否则无突触连接.

1.2 学习模型

考虑一个由n个神经元组成的一般网络，其拓扑连接矩阵A＝（aij）∈Rn×n，其中aij＝1表示第j个神经元激励第i个神经元，否则aij＝0，且aij＝0，i＝j.由于神经网络与一般动力系统的差异性，突触连接是非对称的，所以aij不一定等于aji.C（t）＝（cij（t））∈Rn×n表示突触连接的强度矩阵.

神经元的同步一般指的是在第一个状态变量膜电位上的同步.若单神经元模型为：＝（f1（Xi），f2（Xi），f3（Xi），…，fm（Xi））T，其中表示该神经元的m个状态变量，表示膜电位.突触连接采用电突触模型，则神经网络的数学模型为：

记X＝（X1，X2，…，Xn）∈Rm×n为神经系统的状态.

针对神经元i和j设计如下的突触自适应学习规则：

1）计算突触学习概率

2）在区间［0，1］上产生服从均匀分布的随机数r；

3）比较学习概率与随机数.

若Pij＞r，认为突触发生学习，调整强度；否则，不变.调整规则为：

其中η表示学习速率，Δt表示时间窗长度.

该算法设定了一个概率判定，只有当Pij＞r时学习事件发生.当突触前后神经元都处于阈下活动时非常小，因此突触学习概率Pij非常小，甚至接近于0，于是P＞r发生的概率非常小，所以突触学习发生的可能性也非常小；当它们先后出现非同步的放电越大，Pij越大，越可能出现突触学习.这两个特性非常符合实际生理实验中出现的突触学习规则：突触学习具有时间点依赖性，前后动作电位顺序发生会显著地增强连接［15］，Hebb原则指出阈下活动不易产生突触的强度变化.

T（t）是一个减函数且，类似于模拟退火算法中的降温函数，所以也可以取成其他形式，在温度较高时更容易发生突触学习.

由学习概率的定义知：

这个结果表明当学习时间充分长之后，只要ε取得足够小，学习概率Pij会小于一个任意给定的正数，这会使得学习事件的发生变成一个小概率事件.这也符合生理现象，因为网络的突触学习不可能是无限下去的，当网络学习到一定的程度会趋于成熟，所以这种设计具有较好的生理意义.

当突触学习的计算步长趋向零时，即连续情形下，此刻学习规则为：

1.3 动态相关系数

计算两个序列同步的一个常用统计量为相关系数，如果两个放电序列达到稳态，计算稳态后的相关系数公式为：

该公式能够计算稳态后两个序列的同步程度，该数值越大说明同步程度越高，越小说明同步情况越差，如果等于1说明它们达到了完全同步.

由于自适应学习算法使得同步化成为一个动态的过程，为了刻画这种动态的同步过程，把上面的相关系数计算公式进行改进和推广.

假设样本时间点T＝｛t1，t2，t3，……｝，第j与第k个神经元的放电序列为，详细的符号意义和计算步骤如下：

1）选取时间窗口.不妨设第i个时间窗口长度为Δti，即窗口为［ti，ti＋Δti］.

2）在每个时间窗口里进行采样，不妨设第i个时间窗口的样本共采样Ni个样本点.

3）定义第i个时间窗口内的两个放电序列的相关系数为时刻ti的动态相关系数.即：CCjk（ti）表示在ti时刻第j与第k个神经元的动态相关系数.用这个相关系数序列就可以刻画它们间突触学习过程中的一种同步动态过程.

1.4 相位同步

在两个耦合不全同的振子系统，两个系统可能可以达到同时振荡，但是幅度可以完全不同，称为相位同步.为了研究各种动力系统的相位同步现象，首先是要定义系统的相位.相位的定义方法较多，如解析信号逼近法，旋转振子的直观方法，基于Poincare映射定义法.

神经系统的运动是不规则的，但是神经系统的动作电位有全有或者全无的性质，因此膜电位在每次放电都会达到一个峰值，从非线性理论的角度可以在这个峰值处构造一个Poincare截面，然后可以记录相轨线从一个固定方向穿过此截面的时间点.每两次穿过截面的时间间隔就是完成一个振荡周期的时间，于是可以定义周期的平均.在实际操作中，可以取定一个足够长的时间窗口τ，根据Poincare映射记录达到峰值这个事件发生的次数为N，则平均周期和平均频率可定义为：

于是可定义相位如下：

其中tn表示相轨线第n次穿过Poincare截面的时间点.

C是一个常数，它是由在计算的时间窗口内，最后达到相同步时相差的动作电位数决定的，等于动作电位差乘以2π.

2 结果分析

网络中单个神经元采用Morris-Lecar（ML）模型.其对应的微分方程为：

其中V表示膜电位，W是一个恢复变量，表示钾离子通道开放概率的演化过程，C是膜电容，φ是表示神经元快慢尺度之间的变化，gCa，gK，gL分别表示钙，钾和漏电流的最大电导，VCa，VK，VL分别表示钙，钾和漏电流通道的反转电压.IDC表示来自环境的刺激和突触后电流总和，m∞（V），W（V）分别表示钙离子和钾离子通道打开的概率稳态值.

为了在ML神经系统产生实验中出现的丰富放电模式，引入一个慢子系统如下：

这里的μ表示一个较小的常数，由于系统原来的两个变量m∞（V），W（V）都是快变量，加入的I（V）是一个慢变量系统，可以组成一个双尺度系统，该系统能够产生丰富的放电模式.将原系统改进为：

接着把ML神经元模型应用到网络上去，考虑一个由n个神经元组成的一般网络，其拓扑连接矩阵A＝（aij）∈Rn×n，C（t）＝（cij（t））∈Rn×n表示突触连接的强度.每一个神经元用一个ML模型表示，加上随机的突触连接.

神经网络动力学模型如下：

其中Vi表示膜电位，Wi表示钾离子通道的活化概率，A＝（aij）n×n表示连接矩阵，C（t）＝（cij（t））n×n表示突触连接强度变化矩阵.

以VCa作为控制参数，不失一般性，取神经元数目n＝20，突触连接概率p＝0.2，其余参数取值如下：

为了模拟神经元个体的非全同性，取网络中所有参数为符合正态分布的随机变量.特别的，随机参数VCa～N（μca，σca）的不同取值可以使得网络中神经元处于各种放电节律下（图1（a）），这显然更符合生物神经系统的集群放电特征.

观察非全同网络中神经元的放电历程和主要的计算结果（图1）.前500ms为弱耦合网络在没有经过学习的时空放电模式，在这段时间，由于网络的非全同性，神经元在未学习前，全部处于各自的放电节律，有的处于周期1放电，同时有的神经元却处于2周期，3周期甚至混沌放电（图1（a））.这个结果有别于全同神经网络的结果，全同非对称神经网络在学习前全部神经元处于共同的放电周期，区别在于相位不同.从横坐标原点自适应学习开始，随着学习的进行，发现网络经过几次密集的簇爆发（图1（b））后迅速进入同步（图1（c））.在某个时间点同步地出现第一个动作电位，然后进入一个共同的周期2峰放电节律，达到完全同步，非全同系统中的神经元同步转迁过程呈现复杂的吸引子过渡（图1（d））.

图2给出了非全同神经网络在突触学习过程中神经元动态相关系数的变化，其中每个时间窗口采样点为1000个.结果表明在学习之前，网络神经元之间的相关性较低.当学习开始之后，动态相关系数迅速走高，神经元间的联系逐渐加强，最后趋向于1，达到同步.从动态相关系数形态上看，非全同神经网络同步过程相对平缓.

相位随着放电次数的增加逐渐变大，每当神经元释放一个动作电位，其状态点在相空间里正好完成一条闭轨线，然后在相位上增加2π.而在阈下活动时相位明显增加比较缓慢，因为动作电位发生时状态发生瞬间的急剧变化，所以有一个明显的增加.

为了考察网络的相位同步情况，先考虑两个神经元.根据定义1知道两个神经元之间的相位同步应该满足条件

网络的相同步是指神经网络中所有的神经元，也就是任意两个神经元都必须满足相位同步的条件，于是可以按如下方式来定义.

定义2 一个有N个神经元的网络系统，对于神经网络中神经元i和j定义相位差Δφij（t）＝若有神经系统的平均相位差满足，则称该网络相位同步.

其中，平均相位差可按下式进行计算：

网络相位差计算结果如图3，结果显示任意两个神经元之间的相位差随着时间变化趋于常数（图3（a）），即网络中任意两个神经元出现相同步，网络平均相位差计算如图3（b），神经网络平均相位差趋于常数，由定义知神经网络出现全局的相位同步.相同步是比完全同步稍弱的同步现象，因为出现相同步时由于非全同的原因可能会出现神经元相位相同但是幅值不同的情形.

图1 （a）网络中非全同神经元的各种放电过程；（b）神经网络放电时空模式，μca＝0.875，σca＝0.03；（c）V1-V9 平面上的相图；（d）吸引子Fig.1 （a）Discharge of neurons in a non-identical neural network；（b）Spatiotemporal discharge pattern of non-identical neural network，μca＝0.875，σca＝0.03；（c）Phase diagram on V1-V9；（d）Attractor

图2 （a）第1与第9个神经元的动态相关系数；（b）非全同网络的平均动态相关系数Fig.2 （a）The dynamic correlation coefficient between No.1and No.9neurons in non-identical neural network；（b）Average dynamic correlation coefficient of non-identical neural network

图3 （a）网络中任取两个神经元的相位差变化；（b）神经网络整体上的平均相位差Fig.3 （a）Difference of phase between two neurons；（b）Average difference of phase in network

3 结 论

本文提出了一种神经网络的自适应突触学习规则，这种学习法具有较强的生理意义，不但符合突触学习的生理实验结果，而且对于各种神经网络拓扑结构都具有较好的适应性.

定义了一种动态相关系数计算方法，用移动时间窗口的方式给出了在神经系统同步过程中的同步程度变化，学习之前系统中的神经元都处于相同的周期轨道上，但是相位不同，经过学习完全同步到具有相同相位的周期轨道.针对网络中神经元的非全同性，对参数做随机化处理.假设所有的参数符合正态分布，特别是控制参数的随机化会导致神经系统放电的多样性，这显然更符合实际生物网络的特性.从神经元时空放电模式中发现在学习过程中出现几次密集的簇爆发然后进入同步.在这种突触学习算法下，模拟结果表明即使处于完全不同周期轨道上的神经元组成的集群依然可以同步到任意的轨道.在非全同网络拓扑结构下，处于不同放电节律的神经元可以同步到同一个放电节律下，同一种放电节律的神经元在不同拓扑特性的集群中能够同步到其他的放电节律，最终的集群放电节律由网络拓扑特性决定，如控制参数正态分布的均值和方差.为了刻画网络在整体上的相同步提出了一种网络平均相位差计算方法，在该方法下模拟了网络中任意两个神经元的相位同步和网络的整体相同步，计算结果表明这种方法是有效的.

［1］Belykh I，Enno de Lange，Hasler M.Synchronization of bursting neurons：what matters in the network topology［J］.Phys Rev Lett，2005，94：188101.

［2］Bazhenov M，Huerta R，Rabinovich MI，etal.Cooperative behavior of a chain of synaptically coupled chaotic neurons［J］.Physica D，1998，116（3-4）：392-400.

［3］Yoshioka M.Cluster synchronization in an ensemble of neurons interacting through chemical synapses［J］.Phys Rev E，2005，71：061914.

［4］Wang Q Y，Lu Q S，Chen G R，etal.Chaos synchronization of coupled neurons with gap junctions［J］.Physics Letters A，2006，356（1）：17-25.

［5］Li C G，Chen G R.Synchronization in general complex dynamical networks with coupling delays［J］.Physica A，2004，343：263-278.

［6］Dhamala M，Jirsa V K，Ding M.Enhancement of neural synchrony by time delay［J］.Phys Rev Lett，2004，92：074104.

［7］Drogoi V，Grosu I.Synchronization of locally coupled neural oscillators［J］.Neural Proc Lett，1998，7（3）：199-210.

［8］Neiman A B，Russell D F，Yakusheva T A，etal.Response clustering in transient stochastic synchronization and desynchronization of coupled neuronal bursters［J］.Phys Rev E，2007，76：021908.

［9］Selverston A I，Rabinovich M I，Abarbanel H D I，etal.Reliable circuits from irregular neurons：a dynamical approach to understanding central pattern generators［J］.Journal of Physiology，2002，94：357-374.

［10］Regina M G.D-glucose sensitive neurosecretory cells of the crab cancer borealis and negative feedback regulation of blood glucose level［J］.The Journal of Experimental Biology，1997，200：1421-1431.

［11］Zheng Y H，Lu Q S.Spatiotemporal patterns and chaotic burst synchronization in a small-world neuronal network［J］.Physica A，2008，387：3719-3728.

［12］Wang Q Y，Duan Z S，Perc M，etal.Synchronization transitions on small-world neuronal networks：effects of information transmission delay and rewiring probability［J］.Europhysics Letters，2008，83：50008.

［13］Perc M.Stochastic resonance on excitable small-world networks via a pacemaker［J］.Phys Rev E，2007，76：066203.

［14］Fradkov A L，Andrievsky B，Evans R J.Adaptive observer based synchronization of chaotic systems with first order coder in the presence of information constraints［J］.IEEE Trans on Circuit Syst I：Fundamental Theory and Applications，2008，55（6）：1685-1694.

［15］Li C P，Sun W G，Kurths J.Synchronization between tow coupled complex network［J］.Phys Rev E，2007，76：046204.

［16］Li R，Duan Z S，Chen G R.Global synchronization of drive-response dynamical networks subject to input nonlinearity［J］.J Phys A：Math Theor，2008，41：385103.