基于空间相关各向异性规则的仿真小脑模型的神经网络优化研究

Neural Network Optimization of Simulated Cerebellar Model Based on Spatial Correlation Anisotropy Rule

摘要: 仿真小脑模型的颗粒层神经元通过随机投射网络生成时空活动序列，然而网络中的空间相关性和连通性未能充分反映真实生物组织的结构原则。本文针对该问题提出一种基于新的连通规则的优化方案，借助Perlin噪声的纹理特性让神经元的连接具有优先方向偏好和相邻神经元有相似的偏好，从而模拟生物神经元连接在空间上的不对称性；此外，本文还借助高斯分布设定神经元之间的连接强度，使其随着神经元之间距离的增加而呈指数衰减，这与生物学中兴奋性突触后电流随突触距离神经细胞体中心越远而减弱的现象相吻合为了验证改进的有效性，在相同条件下通过对仿真视动反射任务的对比分析，研究了优化前后的小脑模型在记忆学习能力上的差异，并对两种模型下细胞电生理变化进行了深入探讨。该优化的小脑模型不仅对小脑运动学习记忆形成有着理论意义，而且为小脑控制器提供工程上的应用价值。

Abstract: The granular layer neurons of the simulated cerebellar model generate spatiotemporal activity sequences through a random projection network. However, the spatial correlation and connectivity in the network do not fully reflect the structural principles of real biological tissues. Aiming at this problem, this paper proposes an optimization scheme based on new connectivity rules. By using the texture characteristics of Perlin noise, the connection of neurons has the preference of preferred direction and the similar preference of adjacent neurons, so as to simulate the asymmetry of biological neuron connection in space; In addition, we also use the Gaussian distribution to set the connection strength between neurons, so that it decays exponentially with the increase of the distance between neurons, which is consistent with the phenomenon in biology that the excitatory postsynaptic current decreases with the distance from the synapse to the center of the nerve cell body. In order to verify the effectiveness of the improvement, through the comparative analysis of the simulated optokinetic reflex task under the same conditions, the differences in memory and learning ability of the cerebellar model before and after optimization were studied, and the changes of cell electrophysiology under the two models were discussed. The optimized cerebellar model not only has theoretical significance for the formation of cerebellar motor learning and memory, but also provides engineering application value for cerebellar controller.

1. 引言

有序的动作序列对于任何有意义行为的构建至关重要，这意味着大脑区域内的任务特异性神经元在执行任务时必须按精确的时间顺序激发活动 [1] 。人的时间感知受到脑部特定区域，特别是右侧顶叶与壳核的神经活动的深刻影响。这些活动紧密关联着我们的感官输入信号、学习以及记忆过程，同时伴随着多巴胺水平的微妙变化。而大脑对于短时间刻度的测量，则依赖于海马体和内嗅皮层这两个负责记忆与空间定位的关键部位 [2] 。例如，海马体中的时间细胞通过链式信号精确地反映出事件发生的时间顺序 [3] 。多年来，神经元放电的速率被视为衡量神经活动的标准，但准确的时间也很重要，而相位进动这一机制在情景记忆和空间导航中发挥着至关重要的作用。相位进动(phase precession)是进行情景记忆和空间导航的关键机制，有研究人员发现，在学习陌生移动路线的老鼠中存在相位进动的现象 [4] 。处理序列图像也涉及相位进动 [5] 。另外已有多项研究在不同脑区记录中观察到了时间编码的神经活动序列，特别是在小脑协调连贯身体运动中的关键作用上 [6] ，体现出了从几十到几百毫秒的时间尺度动态。然而，在局部随机、稀疏连接的复杂神经网络环境中，如何产生这种有序活动模式是一个核心问题。

一种理论观点认为，神经元或一组神经元可与刺激开始后的特定时间窗口相关联，形成潜在的时间表达机制(pass of time, POT)。Yamazaki的研究揭示，在连续刺激下，颗粒细胞活性状态随机切换，活性序列非循环且各具特定时间点映射 [7] ，并以此建立了小脑尖峰网络模型，其中颗粒层的连接结构采用随机投射(Random project, RP)网络。

为解决随机投射网络中缺乏生物组织原则的问题，科学家们探索新的连通性方案。例如，Spreizer [8] 提出的模型强调通过非对称性和强度差异的空间连接来自发生成有序活动序列，具体表现为：神经元轴突倾向于向某一方向优先投射，且相邻神经元的投射偏好相似，这反映了生物神经元在空间上的不对称分布特性。

为了深入探究小脑在运动学习记忆形成中的作用，本文着手优化这类网络，对比改进前后的小脑峰值模型在网络投射结构下的时空活动序列及细胞电生理变化。通过仿真分析和关键指标比较，该研究旨在阐明在不依赖明确学习机制或吸引子动力学的情况下，如何通过新颖的连通规则实现随机神经网络中有序时空活动序列的有效生成，从而深化对小脑运动学习记忆原理的理解及其在实际应用中的价值。

2. 实现方法

2.1. 模型结构

参照Yamazaki小脑神经网络 [9] ，以10,000个颗粒细胞规模的小脑模型作为例子。如图1所示，箭头象征着兴奋信号的传递路径，而圆头则代表了抑制信号的流动方向。模型有三种突触纤维接受刺激，分别为苔藓纤维(mossy fibers, MF)、平行纤维(Parallel fibers, PF)和攀延纤维(climbing fiber, CF)，还有六种细胞组成分别为颗粒细胞(granule cells, GRs)、高尔基细胞(golgi cells, GOs)、蓝状细胞(basket cell, BS)、浦肯野细胞(Purkinje cells, PC)、下橄榄核(Inferior olive, IO)和前庭核(vestibular nuclei, VN)。刺激信号由苔藓纤维产生，之后输入分成两个去向，一个是输入到前庭核，同时输入到颗粒细胞层。颗粒细胞层中分别用10 × 10的矩阵来表示高尔基细胞和颗粒细胞蔟，每个颗粒细胞蔟有100个GR，其中每个高尔基细胞以50%的概率接收以自身为中心，周围7 × 7的颗粒细胞蔟的兴奋性输入；每个拥有4个树突，每个树突以2.5%的概率接收自身周围9 × 9的高尔基细胞的抑制性输入。处理完的信号输入到平行纤维中，其中突触可塑性表达线性层的存在，突触权重经过训练自行学习。最后浦肯野细胞的抑制信号和苔藓纤维传输来的兴奋输入共同作用前庭核，前庭核的神经元尖峰作为模型的最后输出。

指导信号和误差由下橄榄核的神经元活动表达，通过攀延纤维影响高尔基细胞活动，攀延纤维的信号只作为突触权重调节的诱导，实现了长时程增强(long-term potentiation, LTP)和长时程抑制(long-term depression, LTD)这两种双向可塑性变化。改变PF的突触权重达到学习的效果。苔藓纤维信号输入到前庭核，类似于输出的时候加入一个正则项。

Figure 1. Cerebellar model network architecture

图1. 小脑模型网络架构

2.2. 计算机制

本模型神经元的膜电位更新基于Yamazaki模型 [8] ，如式(1)：

$C\frac{\text{d}}{\text{d}t}=\underset{c}{{\displaystyle \sum }}{g}_c\left(E_c-V\left(t\right)\right)$ (1)

其中 $V\left(t\right)$ 和C分别是t时刻的膜电位和电容，E是电势。在给定时间的电流下，用电导和触发电势计算更新后的电位。电导的计算通过卷积模型进行，如式(2)：

$g_c=\overline{g_c}\underset{j}{{\displaystyle \sum }}{w}_j{\displaystyle {\int }_{-\infty }^t\alpha \left(t-s\right){\delta }_j\left(s\right)\text{d}s}$ (2)

其中下标c表示不同的离子通道， $\overline{g_c}$ 表示最大电导， $w_j$ 表示来自突触前神经元j的突触权重。指数函数 $\alpha \left(t-s\right)$ 表示在t时刻的电势衰减， ${\delta }_j\left(s\right)$ 代表的是峰值事件，是否触发了动作电位用0和1来表示，如式(3)：

${\delta }_j\left(s\right)=\{\begin{array}{c}1\text{if}{V}_c\left(t\right)\ge h_c\\ 0\text{if}{V}_c\left(t\right)<h_c\end{array}$ (3)

其中 $V_c\left(t\right)$ 是t时刻的电位， $h_c$ 是膜电位阈值。对于发生突触权重更新的浦肯野细胞，其电位在模型中前向计算的公式如式(4)：

$V_{PC}\left(t\right)=\underset{i}{{\displaystyle \sum }}{w}_i^{PC\leftarrow GR}{V}_{MF}^i\left(t\right)$ (4)

其中 $V_{PC}\left(t\right)$ 是苔藓纤维的刺激信号途经随机投射网络、平行纤维和浦肯野细胞后的状态， $w_i^{PC\leftarrow GR}$ 是第i条苔藓纤维对应连接的突触权重， $V_{MF}^i\left(t\right)$ 是在t时刻第i条苔藓纤维的刺激信号。 $w_i^{PC\leftarrow GR}$ 的计算公式如式(5)：

$w_i^{PC\leftarrow GR}={\delta }_{IO}\left(t\right)w_i^{PC\leftarrow GR}{C}_{LTD}{}^{WinG{R}_i}+{\delta }_{G{R}_i}\left(t\right)\left(w_i^{PC\leftarrow GR}+C_{LTP}\left(1-w_i^{PC\leftarrow GR}\right)\right)$ (5)

其中 ${\delta }_{IO}\left(t\right)$ 表示下橄榄体在t时刻有无脉冲， $C_{LTD}$ 是减小突触权重的速率， $WinG{R}_i$ 是单个颗粒细胞在给定时间窗口下产生脉冲的次数，窗口长度设为50 ms， ${\delta }_{G{R}_i}\left(t\right)$ 表示第i个颗粒细胞在t时刻有无脉冲， $C_{LTP}$ 是增大突触权重的速率， $\left(1-w_i^{PC\leftarrow GR}\right)$ 是防止权重在接近1时增长过快。

模型最终的输出结果受苔藓纤维的刺激和浦肯野细胞的抑制，如式(6)：

$V_{VN}\left(t\right)=\underset{i}{{\displaystyle \sum }}{w}_i^{VN\leftarrow MF}{V}_{MF}^i\left(t\right)-\underset{i}{{\displaystyle \sum }}{w}_i^{VN\leftarrow PC}{V}_{PC}\left(t\right)$ (6)

前庭核的电位由苔藓纤维的刺激和浦肯野细胞的抑制决定，因此第二项取负号。模型的最终输出对应前庭核的输出脉冲，如式(7)：

${\delta }_{VN}\left(t\right)=\{\begin{array}{c}1\text{if}{V}_{VN}\left(t\right)\ge h_{VN}\\ 0\text{if}{V}_{VN}\left(t\right)<h_{VN}\end{array}$ (7)

根据任务的不同，我们可以对输出作不同的处理来提取其信息。在视动反射实验中，针对训练后单个周期的前庭核脉冲序列分析，本文能够通过拟合其脉冲频率至正弦波形来揭示运动学习的效果。具体而言，所得到正弦波形的振幅强度直接反映了运动学习过程中产生的增益效应，即系统对视觉刺激引起的运动反应增强或减弱的程度。通过这种方式，得以量化并直观地展示出在模型的学习适应性变化。

2.3. 空间相关各向异性规则

空间相关各向异性(Spatially correlated anisotropic network, SCA)网络规则规定神经元在特定方向投射轴突，且相邻神经元有相似投影方向。这种各向异性但相关的连通性，在纯抑制网络(inhibition network, I-network)和兴奋性与抑制性神经元网络(excitation-inhibition network, EI-network)中都能产生时空活动序列。小脑尖峰模型的颗粒层作为EI网络，针对原来随机投射网络的结构与功能优化，引入了柏林分布以指导神经元对刺激的选择性响应机制。高尔基细胞接收颗粒细胞刺激时，其方向选择遵循柏林分布，刺激概率与距离呈高斯分布，即中心附近刺激概率更高。同样，颗粒细胞处理高尔基细胞输入时，也遵循此原理。这种机制模拟了神经元间刺激的随机传递，同时考虑了空间方向性和距离依赖性的选择性接收。这与生物神经元的方向偏好轴突、树突以及相近神经元突触连接强度大的特征相符 [10] 。因此，颗粒层不仅模拟了神经元刺激的随机性，还引入了空间方向性和距离依赖性的选择性机制。

2.4. 评估标准

参考yamazaki定义的指标，运动幅度或速度的空间信息称为“增益”，这里以输出的正弦波振幅为评估增益的标准，单位是hz。另外对颗粒细胞簇的活动模式随着时间的推移而演变进行评估。对于颗粒细胞簇i，在t时刻的种群平均活性如式(6)：

$z_i\left(t\right)=\frac{1}{{\tau }_{PC}}\underset{s=0}{\overset{t}{{\displaystyle \sum }}}\exp \left(-\left(t-s\right)/{\tau }_{PC}\right)\left(\frac{1}{N_c}\underset{j=1}{\overset{N_c}{{\displaystyle \sum }}}{\delta }_{ij}\left(s\right)\right)$ (8)

其中 ${\tau }_{PC}$ 是8.3 ms， $N_c$ 是100，表示GR层中一个蔟的GR细胞数量。 ${\delta }_{ij}$ 表示第i个蔟的第j个GR细胞有无产生电位尖峰，若有则为1，反之为0。之后进一步计算在 $t$ 和 $t+\Delta t$ 时刻GR蔟活动模式的自相关，如式(7)：

$C\left(t,t+\Delta t\right)=\frac{{{\displaystyle \sum }}_i{z}_i\left(t\right)z_i\left(t+\Delta t\right)}{\sqrt{{{\displaystyle \sum }}_i{z}_i^2\left(t\right)}\sqrt{{{\displaystyle \sum }}_i{z}_i^2\left(t+\Delta t\right)}}$ (9)

其中分子表示GR细胞簇在时间和t时的种群向量的内积，分母将向量长度归一化。因为只取正值，所以相关性取一个介于0和1之间的值。如果时刻和t的种群向量相同，相关性为1，如果向量正交，相关性为0，表明活动种群没有重叠。之后计算相似度指数，如式(8)：

$S\left(\Delta t\right)=\frac{1}{T}\underset{t=0}{\overset{T}{{\displaystyle \sum }}}C\left(t,t+\Delta t\right)$ (10)

其中，T表示MF输入的振荡频率的倒数即周期2 s。如果GR蔟在时间推移中出现重复的活动模式，那么随着 $\Delta t$ 增加， $S\left(\Delta t\right)$ 的曲线呈现震荡，反之 $S\left(\Delta t\right)$ 单调下降，表示GR蔟的时间序列是非周期性的，活跃细胞蔟随着时间的推移而演变为不相关的活动模式。

另外一个指标是GR细胞蔟在不同的MF信号振荡周期中活性模式的重现性，首先计算了两个连续周期在t时刻下活动模式的相关性，如式(9)：

$C^{\left(k\right)}\left(t\right)=\frac{{{\displaystyle \sum }}_i{z}_i^{\left(k\right)}\left(t\right)z_i^{\left(k+1\right)}\left(t\right)}{\sqrt{{{\displaystyle \sum }}_i{z}_i^{\left(k\right)2}\left(t\right)}\sqrt{{{\displaystyle \sum }}_i{z}_i^{\left(k+1\right)2}\left(t\right)}}$ (11)

其中 $z_i^{\left(k\right)}\left(t\right)$ 表示颗粒细胞簇i在第k个周期下t时刻的种群平均活性，以此计算再现性指数 $R\left(t\right)$ ，它表明了两种活动模式是如何在时间上被区分的，作为与连续K周期的相关性的平均值如式(10)：

$R\left(t\right)=\frac{2}{KT}{\displaystyle {\sum }_t\left(C^{\left(1\right)}\left(t\right)+C^{\left(3\right)}\left(t\right)+\cdots +C^{\left(k-1\right)}\left(t\right)\right)}$ (12)

其作者建议用10对连续的周期，即取K = 20，来计算R。

3. 对比与分析

3.1. 实验方案

以视动反射(optokinetic eye movement response, OKR)实验作为例子，用该模型进行仿真，观察视动反射的运动学习记忆形成过程。刺激信号如图2所示，一次训练周期有2000 ms，刺激信号和误差信号的脉冲分布符合正弦信号的频率，分别是30 hz和3 hz，最后两者均在2000 ms时结束，这样的流程对应一次训练，共训练200周期，观察各种细胞的膜电位变化。

Figure 2. Stimulus signal of optokinetic eye movement response

图2. 视动反射刺激信号

在构建颗粒层网络结构时，本文对随机投射网络和空间相关各向异性网络采用了不同的处理策略。对于随机投射网络部分，每个颗粒细胞接收周围9 × 9个高尔基细胞的刺激信号，在方向和距离特性上遵循完全随机分布原则，即各个方向以及不同距离的高尔基细胞传递来的刺激具有同等可能性。

然而，在设计空间相关各向异性网络时，本文引入了一种新颖且生物启发式的连接规则：每一神经元仅将其轴突的小部分(约2%至5%)优先指向特定的方向，并且这一方向选择在相邻神经元之间呈现出显著的相关性。也就是说，邻近的神经元更倾向于共享相似的优先刺激接收方向，而距离较远的神经元则表现出独立的方向选择偏好。

因此，在该空间相关各向异性网络中，针对颗粒细胞优先接收信号的方向特征，本文运用了Perlin分布以模拟自然界的平滑噪声与连续变化趋势。同时，在考虑距离因素时，则采用高斯分布模型，使得距离颗粒细胞越近的高尔基细胞向其传递信号的概率越高。

3.2. 试验分析

图3展示了100 × 100个神经元矩阵的优先投射方向分布，左图是RP网络，右图是SCA网络，颜色代表投射方向。可以看出随机投射网络的投射方向没有特别的偏好，而SCA网络揭示了这样的规则：神经元优先在给定方向投射其一小部分连接，且相邻神经元的优先投射方向相似。

图4所示的图表描绘了在模型训练流程中，配备有SCA网络时PC和VN两种神经元的尖峰放电活动分布。PC和VN的脉冲活动则表现出明显的动态变化特性，尤其是在IO激活频率较高的阶段(1000 ms附近)。第一次训练时，PC神经元的激活率处于较高水平。随着训练过程的推进，其激活率呈现出逐渐且平缓的下降态势。到第30次训练时，已经能够观察到明显的学习迹象。到第200次时，PC出现了尖峰放电暂停。这些现象与采用RP网络的模型的仿真结果相吻合，进一步验证了SCA网络生成的时空活动模式对于小脑运动学习过程具有实际效用。

Figure 3. Distribution of neuron projection direction

图3. 神经元投射方向分布

Figure 4. Potential changes during training

图4. 训练过程的电位变化

图5是两种网络下PC和VN神经元放电频率随训练次数的演变情况。对比分析上半部分，在相同训练阶段下，由于误差信号在数量和时间位置上随机，训练刚开始时其频率的分布还不明显，trail1两种网络的PC放电频率差异较大，但都维持在较高水平。训练后可以看到SCA网络模型中PC的最小放电频率比随机投射网络模型更快地下降。

对比下半部分关于VN的放电频率，训练刚开始时两种模型的VN的放电频率差异较大。当训练进行到trial100时，两者的最大放电频率相近；但进一步延续至200次训练迭代(即trial200)，基于SCA网络颗粒层调控的VN展现出了更高的最大放电频率。

综上所述，在不同的颗粒层配置中，采用SCA网络结构的模型在整个训练过程中体现出更为快速的学习速度和动态调整能力，尤其在后期训练阶段，其对VN神经元放电行为的增益效果相较于随机投射网络模型更为显著。

Figure 5. Frequency change of learning process

图5. 学习过程的频率变化

在图6左图中，本文对比了RP网络与SCA网络两种颗粒层的再现性指数。结果显示，这两种网络在再现性指数曲线上的表现相似，这表明选用不同的网络架构对颗粒细胞簇内部活动分布的再现性能并无显著影响。转至右图，本文进一步分析了两种网络模型相似性指数。研究发现，无论是随机投射网络还是SCA网络，在随着∆t增大时，各自内部两个不同时刻GR细胞簇的放电空间分布相似度均呈现逐渐下降的趋势。值得注意的是，相较于随机投射网络，SCA网络的相似性指数随着∆t的增长表现出更快的衰减速度。综上所述，在用不同活动模式的GR细胞蔟去对应不同时间间隔时，SCA网络能够更有效地捕捉并展示细胞簇的活动模式随时间演变的差异，让其作为时空活动序列在表达POT上有更好的效果。

Figure 6. Index comparison of two kinds of granular layer networks

图6. 两种颗粒层网络的指标对比

4. 讨论

针对小脑神经元模型中颗粒层使用随机投射网络产生时空活动序列表达POT的问题，本文通过使用空间相关各向异性网络作对比分析。通过试验仿真模型能学习视动反射，并仿真各种细胞的脉冲状态变化。学习过程离不开突触可塑性的形成，平行纤维是调节突触权重来学习视动反射的关键位点，受颗粒层生成的时空活动序列调控。最后对比使用空间相关各向异性网络拥有更快的学习速度，这个效果离不开GR细胞蔟的时空活动序列在POT表达上的优化。

本模型相对于Yamazaki的模型而言，颗粒层的GRGO连接结构使用空间相关各向异性网络，具体来说，通过Perlin分布生成颗粒层神经元的优先投射方向，使相邻神经元有相似的连接偏好，从而模拟轴突和树突在空间上的不对称；同时神经元映射周围细胞的距离与概率满足高斯分布，距离越近概率越高，模拟相距较短的两个神经元之间，其突触连接通常具有更大的生理效应或传递强度，从而在模型中实现对真实生物系统中空间相关性的精确模拟和刻画及其对信号传递效率的影响。

5. 结语

本论文主要研究空间相关各向异性网络的特性以优化仿真小脑模型，在空间特性上更接近生理神经元突触分布，借此提高由GR蔟的活动模式生成的时空活动序列表达POT的能力。同时通过视动反射模型的仿真验证对比两种网络对模型学习能力的效果。现在模型在平行纤维突触权重仍不能很好地从生理方面解释权重更新，在采用可塑性反应级联后空间相关各向异性网络是否仍然可以提升POT的表达能力是未来研究的方向。