电磁激励下具有混合突触耦合的Chay神经元系统动力学特性研究

Research on the Dynamical Characteristics of Chay Neuron Systems with Hybrid Synaptic Coupling under Electromagnetic Excitation

摘要: 为了研究混合突触对神经元放电模式和同步状态的影响，考虑电磁激励下Chay神经元模型，采用变步长四阶Runge-Kutta算法进行数值计算。利用峰峰间期分岔图、相图、时间历程图以及双参数分岔图研究不同参数对耦合系统分岔的影响，得出耦合神经元系统在不同的参数值下可产生周期放电、逆倍周期放电和混沌簇放电。利用相似函数同步统计量，研究耦合神经系统同步行为，得出增大耦合强度可促进耦合系统的同步行为，可揭示耦合系统神经信息的编码和传递机制，为治疗神经性疾病提供理论基础。

Abstract: To investigate the effects of mixed synaptic on neuronal firing patterns and synchrony states, the Chay neuronal model under electromagnetic excitation is considered, and numerical calculations are performed using a variable-step fourth-order Runge-Kutta algorithm. Using peak-to-peak bifurcation diagrams, phase diagrams, time-history plots, and two-parameter bifurcation diagrams, the effects of different parameters on bifurcation in coupled systems are investigated. It was found that coupled neural systems can generate periodic discharges, inverse-periodic discharges, and chaotic cluster discharges under different parameter values. Synchronizing statistics using similarity functions, the synchronization behavior of coupled neural systems is investigated, and it is concluded that increasing coupling strength promotes synchronization in coupled systems. This reveals the encoding and transmission mechanisms of neural information in coupled systems, providing a theoretical foundation for treating neurological disorders.

1. 引言

神经元作为神经系统的结构与功能基本单元，通过复杂的突触连接及耦合机制实现神经信息的传递。有生物实验[1]证实，神经信息传递呈现显著的协同特性，其过程不仅与大脑的结构功能紧密相关，还深度参与多种疾病的致病机理调控。为了模拟生物神经系统的分岔与同步机制，学者们建立了人工神经元和神经网络模型，如HH模型[2]、ML神经元模型[3]、HR神经元模型[4]、Chay神经元模型[5]等。随着研究的持续深入，当向神经系统施加脉冲[6]、时滞[7]、电磁[8]等不同外部刺激时，会产生丰富的动力学特性。

当神经元通过电突触[9]、化学突触[10]、混合突触[11]的动态耦合时，可实现信息的编码与传递。针对不同耦合构型下神经信息传递机制的研究，目前已成为研究热点。由于神经系统异常放电时会导致神经性疾病，如帕金森疾病、阿尔茨海默症、癫痫等，故探究神经系统放电行为可为这些疾病的治疗提供理论基础。薛睿等[12]利用哈密顿能量控制的方法，对电磁作用下HR神经元模型进行分岔分析和控制研究，揭示在不同反馈增益下神经元系统会有多种放电模式。刘畅等[13]对电磁耦合Chay神经元模型分岔和同步研究，得出增大磁耦合强度能促进耦合系统的同步行为。张洁[14]等人研究磁作用下化学突触耦合环状Morris-Lecar神经网络模型的分岔与同步，揭示化学耦合强度对神经信息传递过程影响的内在机制。高月月等[15]在电磁感应下研究耦合强度对电突触耦合Chay神经元模型和全局耦合神经网络同步状态的影响。文献[16]构建抑制性化学突触耦合HH神经元模型，探讨了耦合系统在不同参数条件下的簇发同步、反相位同步以及鲁棒性同步问题，通过调节突触中某些参数值，能够改变、降低或可抑制引起致病效应的簇发同步。Chay神经元是基于钾离子通道子通道建立的模型，通过量化钠、钾、钙等离子跨膜流动的电导变化，可以精准复现神经元的发放节律与电生理响应特性，其动力学模型在揭示神经元兴奋传递、周期放电及混沌行为等方面具有独特优势。然而，当前研究多集中于单一突触耦合模式或无外部场激励的场景，鲜少关注电磁激励下混合突触耦合Chay神经系统的分岔机制与同步动力学特性。

基于上述分析，本文基于三维Chay神经元模型，施加电磁感应电流，利用混合突触建立电磁激励下具有混合突触耦合的神经元模型，并采用C语言编写变步长四阶Runge-Kutta算法进行数值计算，利用Grapher软件绘制分岔图、相图、时间历程图和双参数分岔图，得出电耦合强度对耦合系统的分岔行为几乎没有影响，而化学耦合强度值能促进新周期放电类型的出现。分别研究了电耦合强度和化学耦合强度对耦合系统的同步行为的影响，得出当两种耦合强度增大到一定范围时，耦合系统都能达到同步状态，这为神经性疾病的治疗提供理论价值。

2. 构建模型

离子跨膜运动会产生电位差，参照文献[17]中的双曲正切磁控忆阻器模型，来模拟电磁感应效应，其状态方程为

$\{\begin{array}{l}i=\frac{\text{d}q\left(\phi \right)}{\text{d}t}=\frac{\text{d}q\left(\phi \right)}{\text{d}\phi }\frac{\text{d}\phi }{\text{d}t}=W\left(\phi \right)V=-\left(\alpha -2\beta \tanh \left(\phi \right)\right)V\\ \frac{\text{d}\phi }{\text{d}t}=-k\phi +V\end{array}$ (1)

其中，$V$ 、$i$ 分别表示穿过忆阻器的输入电压和输出电流，$\phi$ 表示磁控忆阻器的内部状态变量，$\alpha$ 、$\beta$ 表示磁控忆阻器的内部参数，$W\left(\phi \right)=-\left(\alpha -2\beta \tanh \phi \right)$ 为双曲正切磁控忆导，正实数$k$ 为忆阻反馈系数。

神经元之间是通过化学突触与电突触耦合进行信息编码和传递的，完成人体各器官机能的调节与控制。本文在三维Chay [5]神经元模型的基础上，施加电磁激励并利用电突触与化学突触将两个Chay神经元耦合形成具有混合突触耦合的Chay神经元系统，模型如方程(2)所示：

$\{\begin{array}{l}\frac{\text{d}{V}_{\text{1,2}}}{\text{d}t}=g_I{m}_{\infty }^3{h}_{\infty }\left(V_I-V_{1,2}\right)+g_{KV}\left(V_K-V_{1,2}\right)n_{1,2}^4+g_{KC}\frac{C_{1,2}}{1+C_{1,2}}\left(V_K-V_{1,2}\right)+g_L\left(V_L-V_{1,2}\right)+I_{ext}\\ -k_1\left(\alpha -2\beta \tanh {\phi }_{1,2}\right)V_{1,2}+D\left(V_{2,1}-V_{1,2}\right)+\frac{H_{syn}\left(V_{syn}-V_{\text{1,2}}\right)}{1+{\text{e}}^{-\sigma \left(V_{2,1}-\theta \right)}}\\ \frac{\text{d}{n}_{1,2}}{\text{d}t}=\frac{n_{\infty }-n_{1,2}}{{\tau }_n}\\ \frac{\text{d}{C}_{1,2}}{\text{d}t}=\rho \left(m_{\infty }^3{h}_{\infty }\left(V_C-V_{1,2}\right)-k_C{C}_{1,2}\right)\\ \frac{\text{d}{\phi }_{1,2}}{\text{d}t}=-k_2{\phi }_{1,2}+V_{1,2}\end{array}$ (2)

其中，$m_{\infty }=\frac{{\alpha }_m}{{\alpha }_m+{\beta }_m}$ ，$h_{\infty }=\frac{{\alpha }_h}{{\alpha }_h+{\beta }_h}$ ，$n_{\infty }=\frac{{\alpha }_n}{{\alpha }_n+{\beta }_n}$ ，${\tau }_n=\frac{1}{{\lambda }_n\left({\alpha }_n+{\beta }_n\right)}$ ，${\alpha }_m=\frac{0.1\left(25+V_{1,2}\right)}{1-{\mathrm{e}}^{-0.1V_{1,2}-2.5}}$ ，${\beta }_m=4{\mathrm{e}}^{-\frac{V_{1,2}+50}{18}}$ ，${\alpha }_h=0.07{\mathrm{e}}^{-0.05V_{1,2}-2.5}$ ，${\beta }_h=\frac{1}{1+{\mathrm{e}}^{-0.1V_{1,2}-2}}$ ，${\alpha }_n=\frac{0.01\left(20+V_{1,2}\right)}{1-{\mathrm{e}}^{-0.1V_{1,2}-2}}$ ，${\beta }_n=0.125{\mathrm{e}}^{\frac{-\left(V_{1,2}+30\right)}{80}}$ 。

在本文的研究中，系统参数取值为$V_I=100\text{mV}$ 、$V_K=-75\text{mV}$ 、$V_L=-40\text{mV}$ 、$V_C=145\text{mV}$ 、$g_I=1800\text{mS}/{\text{cm}}^{\text{2}}$ 、$g_{KV}=1700\text{mS}/{\text{cm}}^{\text{2}}$ 、$g_{KC}=10\text{mS}/{\text{cm}}^{\text{2}}$ 、$g_L=7\text{mS}/{\text{cm}}^{\text{2}}$ 、$I_{ext}=6$ 、$k_C=3.3/18$ 、$\rho =0.27$ 、${\lambda }_n=230$ 、$V_{syn}=-15$ 、$\sigma =5$ 、$\theta =-40$ 。$k_1$ 是磁耦合强度，$k_2$ 是忆阻反馈系数，$D$ 为电耦合强度，$H_{syn}$ 为化学耦合强度。

3. 分岔分析

神经系统中信号传递的强弱都是通过耦合强度进行控制的，研究耦合强度对系统动力学特性的影响至关重要。选取参数耦合强度$D$ 、$H_{syn}$ 作为研究参数，初始条件设为$\left(V_1,n_1,C_1,{\phi }_1,V_2,n_2,C_2,{\phi }_2\right)=$ $\left(0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01\right)$ ，即两个神经元的初始条件一样时，采用变步长四阶Runge-Kutta算法进行数值仿真。由于电耦合属于线性耦合，并经过多次数值实验，结果发现改变电耦合强度$D$ 的值对耦合系统的分岔特性和放电模式的影响微小。

由于化学突触是非线性耦合，下面对化学耦合强度$H_{syn}$ 进行数值计算，得到如图1所示的峰峰间期分岔图，图1(a)中所示的p-n表示周期数n放电模式。从图1(a)可以得出，化学耦合强度从零开始增大时，耦合系统处于周期2峰放电。当化学耦合强度增大到1.48时发生倍周期分岔，出现周期4放电，随后再次发生倍周期分岔进入周期8簇放电，耦合系统随着$H_{syn}$ 的递增，最终经过倍周期序列进入混沌放电，这一变化机制也叫吸引子合并激变[18]，即系统从峰放电进入簇发放电状态，如图1(b)所示的局部放大图，并且在图2(c)所示的时间序列图得到验证。当$H_{syn}=2.23$ 时，耦合系统从混沌序列经过边界激变进入周期7簇发放电，峰放电转迁为簇放电。之后随着化学耦合强度$H_{syn}$ 的增大，混沌放电与周期簇放电交替出现，周期数在递增，且混沌窗口变窄，使得当$H_{syn}$ 增大到一定范围时，混沌放电序列消失。因此，增大化学耦合强度的值，系统会出现长时间段的静息状态，即簇发放电，抑制放电节律。

(a) 参数$H_{syn}$ (b) 图(a)的局部放大图

Figure 1. Bifurcation diagram of the peak-to-peak interval of parameter $H_{syn}$

图1. 参数$H_{syn}$ 的峰峰间期分岔图

为了验证图1中，当控制参数变化时耦合系统所产生的放电模式。数值计算出$H_{syn}$ 分别为1、2、3时的相图和时间历程图，如图2所示。先对系统进行去瞬态处理，即取时间段为[100, 130]时，就可得到的耦合系统(2)的数值解。当$H_{syn}=1$ 时，系统处于周期2峰放电；当$H_{syn}=2$ 时，系统经过倍周期分岔进入混沌序列；当$H_{syn}=3$ 时，耦合系统从混沌序列经过边界激变进入周期簇放电，处于周期7簇放电。

Figure 2. Phase diagram and time history diagram

图2. 相图和时间历程图

由于双参数图能够更加清晰的、直观的表现出神经元系统丰富的动力学特性，即参数的取值范围、是否处于周期放电或者混沌放电。所以基于对化学耦合强度$H_{syn}$ 的单参数分岔分析，现在将$H_{syn}$ 作为控制参数，将磁耦合强度$k_1$ 、参数$g_{KV}$ 同时作为研究参数，研究$k_1$ 、$g_{KV}$ 同时变化时对神经元系统动力学特性的影响。经过数值计算得到如图3所示耦合系统在参数平面$\left(k_1,g_{KV}\right)$ 上的放电模式、发生区域和边界特征，右侧颜色块和数字表示对应的周期数，周期数大于等于17时认为是混沌放电，图中的颜色区域表示发生区域。

在参数平面上，周期放电区域都呈现出阿诺德舌[19]形状，方向朝右侧，这些舌状区域代表了系统处于稳定状态。可以得出，不同$H_{syn}$ 值下，周期放电区域的阿诺德舌分布范围存在差异。从图3(a)可以得出，当固定$g_{KV}$ 值，增大$k_1$ 时，耦合系统从周期p-n放电模式经过加周期分岔进入p-(n + 1)放电，再经过倍周期分岔进入短暂的p-2(n + 1)放电模式，后经过鞍结分岔进入伴有周期窗口的混沌区域，在混沌区域产生吸引子合并激变进入周期4放电，最后经历逆倍周期分岔进入周期1尖峰放电。当固定$k_1$ 值，增大$g_{KV}$ 时，耦合系统从p-n周期经历多次逆加周期分岔进入p-(n-i)放电(i表示经历分岔的次数)，后经过鞍结分岔进入伴有周期窗口的混沌区域，同样在混沌区域产生吸引子合并激变进入周期4放电，最后经历逆倍周期分岔进入周期1尖峰放电。若沿着图中黑色直线的方向来看，耦合系统处于伴有混沌放电的加周期分岔模式。在图3(b)，图3(c)中，分岔类型与图3(a)中的类似。整体来看，随着耦合强度$H_{syn}$ 的增大，耦合系统会出现新的周期放电类型，如$H_{syn}=2$ 时，出现p-13放电类型；$H_{syn}=\text{3}$ 时，出现周期p-14、p-15、p-16放电类型。结果表明，增大$H_{syn}$ 能促进新周期的出现，也反映出化学耦合强度对神经元系统动力学稳定状态具有调控的动态作用。

(a) $H_{syn}=1$

(b) $H_{syn}=2$ (c) $H_{syn}=\text{3}$

Figure 3. Bifurcation diagram of the planar $\left(k_1,g_{KV}\right)$ two-parameter system under different chemical coupling strength values

图3. 在不同化学耦合强度值下，平面$\left(k_1,g_{KV}\right)$ 双参数分岔图

现将${\lambda }_n$ 、参数$\alpha$ 同时作为研究参数，仍然将$H_{syn}$ 作为控制参数数值计算出如图4所示的双参数分岔图。同样，在参数平面上也产生了阿诺德舌形状，方向朝上。从图4(a)可以得出，当固定参数$\alpha$ 时，随着增大参数${\lambda }_n$ ，耦合系统从周期1尖峰状态经历倍周期分岔进入周期4放电，后产生吸引子合并激变进入混沌放电区域，从混沌放电区域经历鞍结分岔进入p-n周期放电，最后经历多次加周期分岔进入p-(n + i)放电模式。当固定参数${\lambda }_n$ 时，随着增大参数$\alpha$ ，从下往上看，耦合系统先从p-n经过多次加周期分岔进入p-(n + i)放电模式，再经历鞍结分岔进入混沌区域，后产生吸引子合并激变进入周期4放电，最后经历逆倍周期分岔进入周期1尖峰放电。若沿着图中黑色直线的方向来看，耦合系统处于伴有混沌放电的加周期分岔模式。在图4(b)，图4(c)中，分岔类型与图4(a)中的类似。整体来看，随着$H_{syn}$ 的增大，阿诺德舌形状向左下角运动，在右上角位置处出现周期递增的行为，且递增至混沌区域，也进一步反映出增大化学耦合强度值可使耦合系统呈现出丰富的动力学特性，这为理解神经元信息编码与传递的动力学机制提供更精准的参数依据，也为解决神经系统疾病提供了创新的见解和治疗策略，并完善了神经元系统动力学理论体系。

(a) $H_{syn}=1$

(b) $H_{syn}=2$ (c) $H_{syn}=\text{3}$

Figure 4. Bifurcation diagram of the planar $\left({\lambda }_n,\alpha \right)$ two-parameter system under different chemical coupling strength values

图4. 在不同化学耦合强度值下，平面$\left({\lambda }_n,\alpha \right)$ 双参数分岔图

4. 同步分析

为了研究耦合强度对耦合Chay神经元系统的同步行为的影响，引入同步统计量相似函数$S\left(0\right)$ [20]，其表达式如下所示：

$S\left(0\right)={\left[\frac{\langle {\left(x_1\left(t\right)-x_2\left(t\right)\right)}^2\rangle }{{\left(\langle x_1^2\left(t\right)\rangle \langle x_2^2\left(t\right)\rangle \right)}^{\frac{1}{2}}}\right]}^{\frac{1}{2}}$ (3)

式中，$x_1$ 和$x_2$ 表示的是耦合系统的膜电位，符号$\langle ·\rangle$ 表示的是均值。从(3)式可以得出，当耦合神经元之间的放电规律是相互独立时，都有$S\left(0\right)\ne 0$ ；当耦合神经元之间的放电规律完全同步时，有$S\left(0\right)=0$ ；且$S\left(0\right)$ 值越小，耦合神经元之间的同步化程度越高。

Figure 5. Synchronous graph of the similarity function for chemical coupling strength $H_{syn}$

图5. 化学耦合强度$H_{syn}$ 的相似函数同步图

设定电耦合强度$D=2$ 时，初始条件设为$\left(V_1,n_1,C_1,{\phi }_1,V_2,n_2,C_2,{\phi }_2\right)=(0.01,0.01,0.01,0.01,-0.01,-\text{0}\text{.01}$ $-0.01,-0.01)$ ，即两个神经元的初始条件不相等时，数值计算关于化学耦合强度$H_{syn}$ 的相似函数同步图，如图5所示。随着$H_{syn}$ 的不断增大，相似函数$S\left(0\right)>0$ ，说明耦合系统始终处于完全不同步状态；当$H_{syn}=1.62$ 时，相似函数迅速下降，说明耦合系统在微小的区域内是不完全同步的；之后相似函数逐减小，直到$H_{syn}=6$ 时，$S\left(0\right)\approx 0$ ，耦合系统处于完全同步状态，所以两个神经元之间实现同步状态需要更强的化学耦合强度。

为了验证耦合系统的同步状态，在图5中任意选取两个分岔参数值，通过数值计算得到耦合系统(2)如图6所示的同步验证图。图6(a)是$H_{syn}=1$ 时的相图和时间历程图，从图中可以得出，相轨迹分布杂乱，第一个神经元的膜电位值$V_1$ 与第二个神经元的膜电位$V_2$ 值不相等，即$V_{\text{1}}\ne V_2$ ，耦合系统处于异步状态，与图5对应一致，在时间区间[100, 130]之间，耦合系统处于混沌簇放电模式，且有两个簇；选取$H_{syn}=6$ 时，相轨迹图是分布一条斜率为1对角线上，即$V_1=V_2$ ，耦合系统处于同步状态，在时间区间[100, 130]之间，耦合系统处于混沌簇放电模式，只有一个簇，簇的形状比较宽，说明在该条件下耦合系统持续放电的时间较长。

(a1) $H_{syn}=1$ 相图 (a2) $H_{syn}=1$ 时间历程图

(b1) $H_{syn}=6$ 相图 (b2) $H_{syn}=6$ 时间历程图

Figure 6. Synchronous verification diagram

图6. 同步验证图

以电耦合强度$D$ 作为研究对象，在基准参数值下数值计算出如图7所示的相似函数同步图。可以得出，随着耦合强度值从0开始增大时，相似函数$S\left(0\right)$ 的值一直处于0.23附近，变化微小且耦合系统处于异步状态；直到$D=2.76$ 时，$S\left(0\right)$ 的值大幅度振荡，在这期间耦合系统从异步状态运动到完全同步状态，即$S\left(0\right)$ 的值随时间的变化渐近收敛于0。

Figure 7. Synchronization diagram of similarity function for electrical coupling strength $D$

图7. 电耦合强度$D$ 的相似函数同步图

同理，在图7中任意选取两个分岔参数值$D=1$ 、$D=2$ ，数值计算得到如图8所示的同步验证图。从图8(a)可以得出，当电耦合强度值等于1时，每个神经元处于混沌簇放电状态，两个神经元是异步状态；当电耦合强度等于4时，对比图8(b)，每一个神经元由原来的簇放电转迁为峰放电，且处于完全同步状态，说明在两个神经元初始条件不相等的状态下，耦合系统处于峰放电所需要的电耦合强度值更大。

(a1) $D=1$ 相图 (a2) $D=1$ 时间历程图

(b1) $D=4$ 相图 (b2) $D=4$ 时间历程图

Figure 8. Synchronous verification diagram

图8. 同步验证图

5. 结论

本文建立了在电磁激励下具有混合突触耦合的Chay神经系统。利用数值计算的方法对耦合系统进行单参数和双参数分析，得出不同参数值可以使耦合系统呈现出混沌簇放电模式、倍周期分岔和逆倍周期分岔等丰富的动力学特性。利用相似函数同步统计因子，对耦合神经元系统进行同步研究，得出随着耦合强度增大，耦合系统逐渐趋于完全同步，并利用相图和时间历程图得到验证，这为进一步探究治疗大脑各类疾病提供了理论基础。