耦合演化神经元网络中的自组织临界性研究

Self-Organized Criticality of the Co-Evolutionary Neuronal Network

摘要: 本文借助类Bak-Sneppen (BS)神经元，对稳定的耦合演化网络结构进行了自组织临界性研究。网络中的节点为类BS神经元，节点之间的连边表示神经元之间的突触。在演化过程中，将均匀分布的障碍值统一随机分配给神经元，选择具备最小障碍值即最不稳定的神经元进行“放电”，从而改变自己和最近邻神经元的障碍值，放电神经元随后将处于不应期。本研究主要模拟了耦合神经元网络结构在稳定条件下的雪崩规模分布，首回归时间分布和全回归时间分布，发现其均具有幂律分布特性。另外，在雪崩参数取值在0.45附近时，雪崩规模分布的幂律系数为−1.5，与已有实验结果相吻合。研究结果表明，当耦合演化网络具有稳定结构时，类BS神经元表现出自组织临界性。

Abstract: This study investigates the self-organized criticality (SOC) of stable coupled evolutionary neuronal network by employing Bak-Sneppen (BS)-like neurons. In the network, nodes represent BS-like neurons, while edges denote synapses between neurons. During the evolutionary process, barrier values uniformly distributed in the range are randomly assigned to neurons. The neuron with the smallest barrier value, i.e., the most unstable one, is selected to “fire”, thereby altering its own barrier value as well as those of its nearest neighbors. The firing neuron subsequently enters a refractory period. This research primarily simulates the avalanche size distribution, first return time distribution, and all return time distribution under stable conditions in coupled evolutionary neuronal network, revealing their power-law characteristics. Notably, when the avalanche parameter is approximately 0.45, the power-law exponent of the avalanche size distribution is −1.5, consistent with experimental results. These findings indicate that BS-like neurons exhibit self-organized criticality when co-evolutionary neuronal network possess stable structures.

1. 引言

复杂网络广泛应用于现实生活中，包括人际关系网络[1]，城市道路网[2]，计算机网络[3]等等。复杂网络将这些现实网络抽象为节点和连边，着重研究其内在动力学问题。在复杂网络的发展过程中，提出了多种网络连接方式，包括随机网络[4]，Watts-Strogatz (WS)小世界网络[5]，Newman-Watts (NW)小世界网络[6]，Barabási-Albert (BA)无标度网络[7]等等。近年来，复杂网络在模拟生物神经网络的形成机制方面也取得了进展。Jarman等人[8]在小世界网络的密集区域创建捷径的同时，消除未充分利用的连接，以此对神经网络的形成过程进行模拟。Landmann等人[9]提出了节点状态与连边耦合演化的网络，根据节点的平均活跃度，添加或删除该节点的入边。基于该重连边规则，网络由初始结构向临界状态演化，且临界状态与初始结构无关。

自组织临界性(SOC)的概念最早由Bak等人[10]提出，用来解释产生1/f噪声的基本机制。SOC主要用于描述一类具有多个个体的复杂系统，这类系统仅在个体间相互作用下，可以动态变化到临界状态。在临界状态下，系统中的雪崩行为表现出长程时空关联。较常见的自组织临界模型包括沙堆模型[10]、地震模型[11]、种群演化模型[12]和森林火灾模型[13]。这些系统由于内部组元之间的非线性相互作用，导致一个小的局部扰动可能被放大并传播到整个系统，形成类似于“多米诺骨牌效应”的大规模雪崩，雪崩规模通常遵循幂律分布[14]，表现出自组织临界性。

Bak-Sneppen (BS)模型[15]最早用于研究生态学中不同物种之间的相互作用。Lin等人[16]首次将BS模型应用到研究具有NW小世界网络结构的一维环状神经元网络中，发现神经元节点的雪崩规模具有幂律分布的特性，由此，类BS模型被用于研究神经网络结构中的自组织临界性。随后Lee等人[17]经过研究发现，BA无标度网络在类BS模型中同样具有神经元雪崩规模呈现幂律分布的特性。Zeng等人[18]在WS小世界网络的基础上提出的二维WS小世界网络结构同样具有这样的特性。如今，类BS神经元模型已经成为研究网络结构自组织临界性的一般性方法，该模型的主要优点是计算复杂度低，能够以最短的运行时间体现出神经元网络的不同物理特征。

在人类的大脑发育过程中，尤其从第20周到两岁这一阶段，白质纤维束分数各向异性与年龄增长有关，且增长梯度逐渐变缓[19]。为此，本文模拟了结构趋于稳定的耦合演化网络，并结合类BS神经元模型，研究该稳定网络结构中的自组织临界性。

2. 网络结构

本文从N个孤立神经元开始，通过神经元状态和连边之间的耦合演化，模拟大脑在20周至两岁之间的发育过程，并在其形成的稳定结构上，利用类BS神经元模型验证其自组织临界性。在大脑发育阶段，每个神经元i可以处于激活或静息状态，即${\sigma }_i=\left\{0,1\right\}$ 。神经元之间的连边可以处于激发、抑制或无连边状态，即$C_{ij}=\left\{-1,0,1\right\}$ 。在耦合演化过程中，神经元的状态决定了该神经元入边的增加或减少，而连边的增加或减少同样对下一时刻神经元的状态造成影响。同时，网络连边是有向的，即$C_{ij}$ 并不一定等于$C_{ji}$ 。

耦合演化过程中，神经元之间通过如下过程生成连边：

1. 初始时刻，神经元之间无连边。

2. 耦合演化过程中，神经元状态变化的概率如公式(1)、(2)所示。$t+1$ 时刻，神经元i的激活概率取决于前一时刻与其通过入边相连的所有神经元状态，其中$\beta$ 为温度的倒数，可以表示系统所包含的噪声，

本文在数值模拟过程中，令$\beta =10$ 。定义活跃度参数$A_i\left(t\right)$ ，表示神经元i在t及t时刻之前的时间段$E_T$ 内的平均状态，即平均活跃度$A_i={\langle {\sigma }_i\rangle }_{E_T}$ 。

$Prob\left[{\sigma }_i\left(t+1\right)=1\right]=\frac{1}{1+\exp \left[-2\beta \left(f_i\left(t\right)-0.5\right)\right]}$ (1)

$f_i\left(t\right)={\displaystyle \sum _{j=1}^N{C}_{ij}{\sigma }_j\left(t\right)}$ (2)

3. 根据规则2，选择$E_T=1000$ ，在每个演化时间步(每个时间步等于一个$E_T$ 时间段)中，随机选择神经元i，计算其平均活跃度$A_i$ ，并依据其可能的三种情况选择相应的连边方式，其概念图如图1所示：

1) 如果$A_i\left(t\right)\in \left[0,\xi \right]$ ，表示在之前的$E_T$ 时间段内，神经元i受到的刺激较少，需要添加一个激发态的入边，即$C_{ij}=1$ ，其中j是t时刻与神经元i之间无连边的某个随机选择的神经元；

2) 如果$A_i\left(t\right)\in \left[1-\xi ,1\right]$ ，表示神经元i在$E_T$ 时间段内受到了过多的刺激，需要增加一个抑制性的入边，即$C_{ij}=-1$ ；

3) 如果$A_i\left(t\right)\in \left(\xi ,1-\xi \right)$ ，则表示神经元i处于细致平衡状态。为保证其动态演化过程，随机选择一条入边进行删除操作，其中$\xi \ll 1$ (本文采取$\xi =0.001$ )。

4. 演化一定的时间步之后，网络结构趋于稳定。

Figure 1. The conceptual structure of co-evolutionary neuronal network

图1. 耦合演化神经元网络概念图

3. 类BS神经元模型及其动态演化

网络结构稳定后，在其中测试类BS神经元模型[15]的动态演化，具体过程如下：

1. 为网络中的每个神经元i随机分配一个障碍值$B_i$ ($B_i\in \left(0,1\right)$ 之间均匀分布的随机数)。

2. 在每个时刻t，选择具有最小障碍值$B_{\min }$ 的神经元，进行放电操作，从而改变神经元i及其邻居j的障碍值。放电后，神经元i处于不应期$T_r$ 。

3. 接下来的$T_r$ 时刻内，处于不应期的神经元不会再次放电，但其最小障碍值可以随着其他神经元的放电而发生变化。

4. 演化过程结束后，统计体系中相关物理量的分布，并判断其是否具有自组织临界性。

4. 结果分析

4.1. 耦合演化神经元网络结构

耦合演化神经元网络结构形成的具体过程如第2节所述。图2给出了网络从无连接状态，以耦合演

化的方式过渡到具有稳定结构的过程。横坐标表示演化步数，纵坐标代表平均连边数$\langle k_{+}\rangle$ 、$\langle k_{-}\rangle$ 和参数$\langle \lambda \rangle$ ，其中$\langle k_{+}\rangle$ 和$\langle k_{-}\rangle$ 分别表示激发性入边和抑制性入边占所有连接的比例，$\langle \lambda \rangle$ 为某节点的状态翻转所引起的状态发生变化的平均节点数量。$\langle \lambda \rangle \approx \text{1}$ 表示网络结构处于临界状态；$\langle \lambda \rangle >1$ 表示在若干个时间步长之后，网络中所有神经元的状态可能发生变化；$\langle \lambda \rangle <1$ 则表示单个神经元状态翻转引起的变化可能不会在网络中传播。因此，$\langle \lambda \rangle$ 是衡量神经元网络演化结果的一个重要参量，可称之为分支参数。

图2表明，在经历约20,000个时间步的演化后，网络中的平均连边数$\langle k_{+}\rangle$ 、$\langle k_{-}\rangle$ 和参数$\langle \lambda \rangle$ 均达到

了动态平衡。图3给出了节点数为4900，演化到第70,000步时的网络度分布，发现在半对数坐标系下，该分布接近于指数分布，表明此时的网络中节点的连接方式处于类随机网络结构。除此之外，耦合演化达到的稳态结构所具有的集聚系数$c<1/1000$ ，而平均路径长度$L=\infty$ ，表明网络中存在孤立神经元。下述关于类BS神经元的自组织临界性模拟研究所采用的网络均为该稳定结构。

Figure 2. Temporal graph of co-evolutionary neuronal network

图2. 耦合演化神经元网络的时序图

Figure 3. Degree distribution graph of co-evolutionary neuronal network

图3. 耦合演化神经元网络的度分布图

4.2. 最小障碍分布

在演化结束后的时间段T ($T=300000$ )中的每个时刻，稳定网络结构上具有最小障碍值$B_{\text{min}}$ 的类BS神经元，将被选中并进行放电，影响其最近邻神经元的障碍值。当演化过程结束后，统计整个过程中的最小障碍值分布$p\left(B_{\text{min}}\right)$ ，如图4所示。结果表明，$p\left(B_{\text{min}}\right)$ 的值随着$B_{\text{min}}$ 的增大而减小，即神经元的障碍值最大，被选中进行放电的概率就越小。由于障碍值为在0~1之间均匀分布的随机数，因此最小障碍值较小的神经元放电的概率较大。放电结束后，放电神经元节点及其最近邻神经元的障碍值都会被重新赋值，因此在某些时刻，体系中神经元的最小障碍值会以较小的概率增加。图4还表明，整个演化过程中，被选中进行放电操作的神经元所具有的$B_{\text{min}}$ 值小于截止障碍值$B_{\text{0}}$ ，$B_{\text{0}}=0.4$ 。

Figure 4. The minimum barrier distribution of co-evolutionary neuronal network

图4. 耦合演化神经元网络的最小障碍值分布图

4.3. 雪崩行为

雪崩行为是研究体系是否具备自组织临界性的重要指标。实验结果表明，在小鼠神经细胞中，雪崩规模分布呈现幂律分布，且幂律指数$\alpha =-1.5$ [20]。对于类BS神经元模型，定义雪崩规模如下：选择雪崩临界值$B_c$ ，演化过程中某时刻某神经元的$B_{\text{min}}$ 大于$B_c$ ，则雪崩开始，在雪崩的持续过程中，当系统中某个神经元的$B_{\text{min}}$ 值再次大于$B_c$ 时，一次雪崩结束。雪崩过程所持续的时间步称为本次雪崩的规模。根据雪崩规模的定义，不难看出，$B_c$ 值直接决定了雪崩规模分布。为检验$B_c$ 值对雪崩规模分布的影响，本文借由图4的$B_{\text{min}}$ 分布图，引入雪崩参数$\psi =P\left(B_{\min }<B_c\right)$ ，当$\psi =1$ 时，$B_c$ 取$B_{\text{min}}$ 分布中的最大值，即截断障碍值$B_{\text{0}}$ 。

图5展示了雪崩参数$\psi$ 取不同参数值时的雪崩规模分布，在双对数图中，雪崩规模分布的幂律趋势随着$\psi$ 的增加而愈发明显：$p\left(S\right)~S^{-\alpha }$ ，幂律指数$\alpha$ 逐渐减小，如图6所示。值得一提的是，当雪崩参数$\psi$ 的取值在0.45附近时，雪崩规模的幂律系数$\alpha$ 趋于1.5，与实验相符[20]。在雪崩参数$\psi$ 较小时，相应的$B_c$ 也会变小，导致系统中障碍值$B>B_c$ 的神经元数量增加，从而呈现出较小的$\alpha$ 值。

Figure 5. The power law distribution graph of avalanche size

图5. 雪崩规模分布图

Figure 6. The power law coefficients of avalanche size change with avalanche coefficients

图6. 雪崩规模幂律系数随雪崩系数变化的趋势图

4.4. 首回归时间

在现实的神经元实验中，通常会监测参量——峰值时间段[16]，与本文定义的首回归时间T本质上是一致的。在时间段T内，神经元i在某个时间点a放电，并在时间点a后的时间点b再次放电，则神经元i的首回归时间为$t_b-t_a$ 。因此，每个神经元可以有多个首回归时间。在演化过程结束后，对各神经元的首回归时间进行统计，其分布$F\left(T\right)$ 如图7副图所示。在不应期参量$T_r$ 取不同值时，网络中神经元的首回归时间分布稍有不同。由演化过程可知，如果将神经元的不应期$T_r$ 为1个时间步时，则最小首回归时间为2，即最小首回归时间$T=T_r+1$ 。鉴于不同$T_r=\left\{2,3,4,5,6\right\}$ 的取值，首回归时间分布$F\left(T\right)$ 具有一定的相似性，本文对其进行了重标度处理，对应于图7的主图，即将首回归时间T及其分布$F\left(T\right)$ 分别乘以$q^{\eta }$ 和$q^{{\eta }_1}$ ，其中不应期，重标度参数$q=T_r/n_0$ ，$n_0$ 为最小不应期2，$\eta =0.4$ ，${\eta }_1=0.2$ 。因此，在不同的不应期参数条件下，首回归时间及其分布之间存在标度关系：$q^{\eta }F\left(T\right)~{\left(q^{{\eta }_1}T\right)}^{-\tau }$ ，且标度系数为$\tau =1.\text{8}$ 。

Figure 7. The first return time distribution graph of co-evolutionary neuronal network

图7. 耦合演化神经元网络的首回归时间分布图

4.5. 全回归时间

全回归时间$T_a$ 是神经元i在演化时间段T内，任意两个放电时刻之间的时间间隔，因此全回归时间是首回归时间的数学扩展。与首回归时间类似，在整个演化过程中，每个神经元可能有多个全回归时间。如果神经元i在三个时间点上放电，$t_a$ 、$t_b$ 和$t_c$ ，那么该神经元具有三个全回归时间：$t_b-t_a$ 、$t_c-t_b$ 和$t_c-t_a$ ，类似于统计学中的不重复抽样问题，因此如果神经元i在n个时间点上触发，则神经元i有$C_n^2$ 个全回归时间。由此，对耦合稳定网络中所有神经元的全回归时间进行统计，其分布$G\left(T_a\right)$ 如图8副图所示，其中神经元的不应期取值条件下与首回归时间分布一致。但是与图7不同的是，全回归时间分布未做归一化处理。图8的主图对首回归时间分布做了相应的重标度处理。在重标度过程中，将全回归时间

$T_a$ 及其分布$G\left(T_a\right)$ 分别乘以对应的${\theta }^{\delta }$ 和${\theta }^{{\delta }_1}$ ，重标度参数$\theta =T_r/n_0$ ，$n_0$ 为最小不应期1，$\delta =0.5$ ，${\delta }_1=0.1$ ，不应期的取值为$T_r=\left\{1,2,3,4,5\right\}$ 。图8说明，在较小的全回归时间分布中，存在如下的标度关系：${\theta }^{\delta }G\left(T_a\right)~{\left({\theta }^{{\delta }_1}{T}_a\right)}^{-\tau }$ ，且标度系数$\tau =1.4$ 。

Figure 8. The all return time distribution graph of co-evolutionary neuronal network

图8. 耦合演化神经元网络的全回归时间分布图

5. 结论

本文讨论了稳定耦合演化神经元网络的自组织临界性。已有研究表明，大脑在发育过程中，白质纤维束分数的各向异性随年龄增大而逐渐趋于稳定[19]，为此，本文研究了网络结构和神经元状态在耦合演化条件下所实现的稳定结构的特征，并由此结合类BS神经元模型，探究其自组织临界性。

研究发现，类BS神经元模型下，网络结构的各项指标，如雪崩规模S在不同雪崩临界值$\psi$ 下的分布，首回归时间分布和全回归时间分布都表现出幂律行为，表明系统的动态行为具有非线性特征。当网络结构的雪崩参数$\psi$ 接近临界值0.45时，雪崩规模的幂律系数$\alpha$ 趋于1.5，与小鼠神经细胞的实验结果相符[20]，同时，雪崩规模的幂律分布呈现出长尾效应，即大规模雪崩的发生概率较低，同时衰减速度缓慢，这种分布反映了耦合演化神经元网络中极端动态事件的偶然性与复杂性。

本文丰富了现有使用类BS神经元模型对一维NW小世界网络[16]、二维WS小世界网络[18]和BA无标度网络[17]的自组织临界性的研究，并将进一步推进有向网络和多层网络中自组织临界性的研究。

致 谢

感谢Dr Baumgarten在耦合演化网络建模方面提供的宝贵建议和支持。

基金项目

国家自然科学基金(11705097)、南京邮电大学自然科学基金(221101, 222134)。

NOTES

^*通讯作者。