基于条件互信息的动态贝叶斯法探明生物神经元网络连接结构

任婧雯，董朝轶，朱美佳，白鹏辉，赵肖懿，马爽，贾婷婷

1.内蒙古工业大学电力学院，内蒙古呼和浩特010080；2.内蒙古机电控制重点实验室，内蒙古呼和浩特010051

前言

大脑的连接结构错综复杂，是迄今为止人类已知最为复杂的系统之一，为了进一步探明大脑连接结构并深入研究与大脑相关的生理与病理功能，多种探究脑网络的技术手段，如脑电图（Electroencephalograph,EEG）［1］，脑磁图（Magnetoencephalograph, MEG）［2-3］，功能磁共振成像（functional Magnetic Resonance Imaging,fMRI）［4］，双光子显微镜（Two-Photon Microscope, TPM）［5-6］和多电极阵列（Multielectrode Array, MEA）［7］等得到了快速发展。目前，通过以上技术手段，可获得高通量具有较高时空分辨率的生物神经元网络（Biological Neural Network,BNN）多通道时间序列数据。计算神经生物学者们则利用这些脑网络响应数据，发展出多种BNN 功能性连接结构反向辨识算法，试图建立大脑的功能性图谱。目前，针对BNN 的连接结构辨识方法有动态贝叶斯网络算法［8］、隐马尔可夫模型算法［9］、线性格兰特因果性（Linear Grant Causality,LGC）算法［10］和互信息算法［11］等。虽然这些传统方法对辨识生物调控网络的功能性连接结构均具有一定的有效性，但对于复杂的、具有神经元点火的本质非线性和突触连接塑性的BNN 而言，大部分方法难以获得较高的辨识精度，且计算运行时间较长。而互信息算法和动态贝叶斯网络算法相结合探明BNN连接结构克服了以上算法的弊端，对于以后探究大脑的功能性连接结构和构建人工智能网络具有重要意义。

本研究采用基于条件互信息的动态贝叶斯网络法（Dynamic Bayesian Network Method Based on Conditional Mutual Information, DBNM-BCMI）以 辨识人工构造的BNN连接结构。主要思想为：（1）利用积分点火（Integrate-and-Fire, IF）原理构建具有一定生物真实性的脉冲神经元网络（Spike Neural Network, SNN）。因为其能够以较高精度模拟BNN的IF、传输延时等机制，所以SNN 被称为“第三代人工神经网络”，已被广泛应用于计算神经生物学研究［10］。通过调整建模参数使其接近真实的BNN，在MATLAB环境下运行网络仿真程序获得SNN多通道脉冲响应序列数据。（2）利用提出的DBNM-BCMI 算法探明SNN 的功能性连接结构。（3）对比互信息算法、LGC 算法和DBNM-BCMI算法的辨识效果，结果表明DBNM-BCMI 算法优于互信息算法和LGC 算法，其对4、10、20、50、100 节点SNN 连接结构的辨识准确率均高于互信息算法和LGC 算法，且运行时间较短，验证了所提出算法的有效性。

1 方法

1.1 互信息算法

假设两个随机变量X、Y的联合概率密度函数表示为p(x，y)，变量的边缘概率密度函数为p(x)、p(y)。定义信息熵H，表示变量的不确定度。

其中，H(X,Y)称为X、Y的联合信息熵，H(X|Y)和H(Y|X)被称为条件信息熵。联合信息熵与边缘信息熵和条件信息熵之间具有如下关系［12］：

互信息I(X,Y)用于量化任意两个随机变量X、Y间的相互作用关系，图1描绘了信息熵和互信息间的关系［13］。

图1 信息熵和互信息的关系Fig.1 Relationship between information entropy and mutual information

由式（6）可得到互信息的计算公式，即式（7）、（8）：

互信息算法通过计算BNN在放电活动中各网络节点（神经元）数据携带信息的相关性来辨识网络连接结构。即使神经元在放电活动中没有明显的同步性，仍可利用变量间作用的延时性，辨识出因果性的动态连接关系。因突触前神经元点火的时间一般发生在突触后神经元点火之前，基于互信息的因果性辨识法，可同时辨识出突触连接的作用方向。因此，式（8）可修正为：

表示在动态贝叶斯网络框架下（k为当前时刻，τ为突触前神经元到突触后神经元动态作用的纯滞后，两者均为正数），节点X指向节点Y的因果性互信息。这里约定本研究中所有涉及到的互信息计算均具有方向性的因果性互信息。

目前，基于互信息的网络功能性连接辨识算法的详细流程为：通过将两个神经元节点变量间的互信息值与给定的阈值δ进行比较。若该互信息值大于δ，说明两者之间具有相互依赖关系，用数字“1”表示，反之，如果互信息的值小于δ，说明两者间相互独立，用数字“0”表示［9］。这样就可以得到表示网络连接关系的关联矩阵，构造出神经元网络的连接结构。由于该算法计算了突触前神经元数据和延时了突触传递纯滞后时间后的突触后神经元数据，所以该算法能够辨识网络连接的方向。但是该算法存在显著的不足：忽略了在复杂非线性BNN中，其余神经元节点对这两个神经元节点的影响和作用，即未能区分直接连接和间接连接。因此，在实际应用中，存在一定的局限性。

1.2 改进互信息算法

为解决BNN 中第三方神经元节点对两个神经元节点间互信息计算值的影响。引入条件互信息的概念，对已有的互信息算法进行改进，提出DBNMBCMI 算法。该算法采用动态贝叶斯网络作为初始网络构架，利用互信息或条件互信息计算值对网络连接进行修剪，从而得到正确的网络连接图。在动态贝叶斯网络中，存在3 种基本连接结构：依赖结构（A→B→C）、非依赖结构（A←B→C）、导出依赖结构（A→B←C）。构建动态贝叶斯网络结构时，要利用互信息或条件互信息保留依赖结构，删除非依赖结构和导出依赖结构［14-16］。

具体算法描述如下：（1）对于一个含有N个神经元节点SNN，通过网络仿真产生多通道动态响应数据D。（2）计算两两神经元节点X、Y在数据集D条件下的因果互信息值，与给定的阈值δ1进行比较。若I(X，Y|D)＜δ1，说明这两个神经元节点之间无边连接，关联矩阵中相应位置的元素置为“0”；如果I(X,Y|D)＞δ1，说明神经元节点之间有边连接，记录关联矩阵数字为“1”。逐次考察并计算所有两节点间的互信息值，得到关联矩阵T(X,Y)，依据关联矩阵中的连接情况，可获得初始动态贝叶斯网络图O。（3）在O 的基础上增加其它依赖结构的边。如果在O 中发现X、Y之间不存在边，则用Zx(X,Y)表示在X的邻域中，X和Y的通路上的节点集合；Zy(X,Y)表示在Y的邻域中，X和Y的通路上的节点集合。选择Zx(X,Y)和Zy(X,Y)中边数较少的一方，用Z(X,Y)表示。通过与给定的阈值δ2进行比较，如果I(X,Y|Z(X,Y))＞δ2，则在O 中增加边X→Y，记录数字为“1”，反之不变，直到O 中所有符合条件的神经元节点检验过为止。将记录的数据存入新的关联矩阵T1(X,Y)中，并得到新的动态贝叶斯网络图O1，即得到增加依赖结构边后的贝叶斯网络图。（4）在O1的基础上剔除非依赖结构边。如果在O1中发现X、Y之间存在边，则用Z1x(X,Y)表示在X的邻域中，X和Y的通路上的节点集合；Z1y(X,Y)表示在Y的邻域中，X和Y的通路上的节点集合。选择Z1x(X,Y)和Z1y(X,Y)中边数较少的一方，用Z1(X,Y)表示。通过与给定的阈值δ3进行比较，如果I(X,Y|Z1(X,Y))＜δ3，则在O1中剔除边X→Y，记录数字为“0”，反之不变，直到O1 中所有符合条件的神经元节点检验过为止。将记录的数据存入关联矩阵T2(X,Y)，并生成关联矩阵图O2，即得到剔除非依赖结构边后的动态贝叶斯网络图。（5）在O2的基础上剔除导出依赖结构的边。对于任意的神经元节点对X，Y，如果I(X,Y)＜δ4，与X和Y可能形成导出依赖结构的结点可能为Xv1，…，Xvt，需要对每一个可能的导出依赖结构进行碰撞识别［14-16］，假设：

如果X、Y和Xvh形成导出依赖结构，判断是否存在X→Xvh和Y→Xvh，如果存在，则剔除边X→Y［12］，记录关联矩阵相应位置数字为“0”，反之不变，直到将O2中所有符合条件的神经元节点对全部检验过为止。将记录的数据存入新的关联矩阵T3(X,Y)并生成剔除导出依赖结构边后的动态贝叶斯网络图O3。

在每一步计算中，将需要保留边的最小互信息值和需要删除边的最大互信息值平均后，确定为阈值δ1、δ2、δ3、δ4的取值，这样就可保证较高的网络结构辨识率。

1.3 动态贝叶斯网络算法

在DBNM-BCMI 算法进行网络辨识的过程中，对网络初始结构的假设来源于动态贝叶斯网络结构形式。动态贝叶斯网络是静态贝叶斯网络的扩展，在一阶马尔科夫过程的假设下，所有网络变量的联合概率密度可表示为在父节点作用下子节点发生的条件概率密度乘积［17］，其中父子节点的数据间存在一个纯滞后的时间推移。

假设BNN有N个神经元节点，节点变量集合为X={X1,X2,X3,…,Xn}，用Pa(Xi)表示节点Xi的父结点集合。静态贝叶斯网络在变量集X的联合概率密度表示为：

在贝叶斯网络基础上加入时间推移，演变为动态贝叶斯网络，式（11）变为：

图2和图3分别表示了具有相同网络结构的静态贝叶斯网络和动态贝叶斯网络的示意图。在动态贝叶斯网络中，每个变量之间在前后相邻时间片间的因果性概率依赖关系描述了SNN 的网络拓扑结构，可采用定量化推理方法对动态贝叶斯网络进行动态建模、分析、预测和控制，这一网络结构推理方法已被广泛应用于语音识别、系统控制、演化分析中［18］。采用动态贝叶斯网络算法辨识BNN 的结构不仅可以减少模型结构的复杂度，而且可以增加学习推理的效率。

图2 静态贝叶斯网络的先验图Fig.2 Predictive diagram of static Bayesian network

图3 动态贝叶斯网络图Fig.3 Dynamic Bayesian network diagram

2 IF模型

为验证DBNM-BCMI 算法的有效性，本研究采用SNN 对BNN 进行仿真模拟。常见的SNN 建模和仿真模型有IF 模型、Hodgkin-Huxley 模型、二维神经元模型等［19］。其中，Hodgkin-Huxley 模型和二维神经元模型基于多变量非线性微分方程组，求解过程复杂，不适用于大规模网络仿真；而IF 模型在考虑所有突触前神经元脉冲电压对突触后神经元膜电位的累积效应的同时，仅涉及到一元积分，计算量小，所以更适合于大规模BNN 仿真，相对于其他两种模型更适合本研究。

IF 模型由Louis Lapicque 在1907年提出，主要用于BNN 网络脉冲编码、记忆和网络动态研究。本研究利用IF 模型构建SNN，从而实现对大规模BNN 进行仿真、模拟以及网络结构反向辨识算法验证［20］。

IF 模型将神经元网络视为开关网络，对突触后神经元i的动作电位进行建模，膜电位u(t)受钠、钾离子通道和氯离子通道影响，表示为：

其中，τm为神经元膜电位动态时间常数，R和C分别为膜电阻和膜电容。输入电流I(t)由突触前神经元电活动产生，每一突触前尖峰都可产生一个突触后电流脉冲。若突触前神经元j在时刻点火，神经脉冲信号经轴突传播后，突触后神经元i将获得输入电流此时，神经元i的总输入电流为所有输入电流之和：

多个突触前神经元j对突触后神经元i膜电位的贡献度不同，这种不同效应用权值wij表示。当膜电位高于阈值电压θ时：

突触后神经元i被点火。点火后，膜电位迅速降落到静息电位水平：

在点火后的绝对不应期内，即使膜电位再次高于θ，突触后神经元也不会被点火。IF 模型具有结构简单、计算量小等特点，相比二维或多维微分方程模型，更适用于大规模BNN 仿真和建模。本研究采用IF模型对大规模BNN 进行仿真，最终将所提出的DBNM-BCMI 算法应用于网络动态响应数据，验证算法的正确性。

3 仿真与分析

3.1 10个神经元节点SNN结构辨识

以10 个神经元节点SNN 为例，分别验证互信息算法、LGC 算法和DBNM-BCMI 算法对SNN 结构的辨识精度。另外矩阵B和矩阵T3分别表示真实SNN的关联矩阵和辨识出的网络关联矩阵，维数均为N×N。辨识精度被定义为：满足条件B(i,j)=1 且T3(i,j)=1 或B(i,j)=0 且T3(i,j)=0（i,j=1,2,3,…,N）的所有有向支路数，除以N节点网络的全连接支路数N(N-1)。根据BNN 的连接特点，网络连接度通常在10%～30%之间［19］，本研究所构造SNN 的随机连接比为0.2，即每一神经元同网络中20%的其他神经元随机相连。假设每个神经元均受到当地噪声电流Iin(t)作用，Iin(t)均为高斯白噪声。仿真参数为：突触传递动态时间常数τs= 3.5 ms，胞膜动态时间常数τm= 8 ms，突触传递纯时延τ= 50 ms，神经元的绝对不应期Trefractory=2 ms，膜电位的点火阈值θ= -40 mV，突触后神经元静息电位ur= -65 mV。

经网络仿真后，以产生的10个神经元节点中前5神经元节点的神经脉冲响应数据（经10 ms 等周期采样，记录平均点火率）为例，如图4所示。所构造的随机连接的10 个神经元节点SNN 的关联矩阵如表1所示（突触连接的方向遵循“列作用于行”的原则），其基于IF模型构建的SNN结构如图5所示。

图4 10个神经元节点中前5个神经元节点SNN仿真所产生的脉冲序列Fig.4 Pulse sequences generated by spike neural network (SNN)simulation of the first 5 neuron nodes of 10 neuron nodes

表1 在IF模型构建的SNN结构B矩阵Tab.1 SNN structure B matrix constructed in integrate-andfire(IF)model

图5 基于IF模型构建的SNN结构Fig.5 SNN structure constructed based on IF model

采用互信息算法、LGC 算法和DBNM-BCMI 算法对该SNN 结构进行辨识。互信息算法的阈值δ为1.5；LGC 算法的阈值为0.2；DBNM-BCMI 算法的阈值分别取δ1=0.3、δ2=0.1、δ3=0.1、δ4=0.3。辨识后的结果如图6所示。其中，粗虚线表示真实网络结构中边未连，而辨识网络结构中边相连；虚点线表示真实网络结构中边相连，而辨识网络结构中边未连。

图6 3 种不同方法辨识出的SNN 结构Fig.6 SNN structures identified by 3 different methods

DBNM-BCMI 算法辨识的准确率为0.98，LGC算法辨识的准确率为0.93，MI 算法辨识的准确率0.92。由此可见，DBNM-BCMI 算法对该10 个神经元节点SNN 结构的辨识准确率明显高于另外两种方法。

3.2 不同规模的SNN结构辨识

分别对4、10、20、50、100节点的SNN结构进行辨识。对每一种规模网络，进行10轮，每轮100次测试。在每一次测试中，以0.1、0.2、0.3网络连接比各随机产生一个SNN。最后，统计在每一网络规模和每一网络连接比下的平均辨识准确率，结果如表2～4所示。

表2 连接比为0.1时不同方法辨识准确率的对比Tab.2 Comparison of identification accuracy among different methods when the connection ratio is 0.1

表3 连接比为0.2时不同方法辨识准确率的对比Tab.3 Comparison of identification accuracy among different methods when the connection ratio is 0.2

表4 连接比为0.3时不同方法辨识准确率的对比Tab.4 Comparison of identification accuracy among different methods when the connection ratio is 0.3

结果表明：（1）无论对于传统方法还是改进方法，在同一连接比下，随着网络规模的增大，辨识准确率降低，且消耗计算时间显著增加；同时，在同一网络规模下，网络连接比越大，辨识准确率显著降低，且消耗计算时间增加。（2）虽然随着网络规模增加，辨识准确率降低，但50个神经元节点网络和100个神经元节点网络相比，辨识准确率变化不大，提示对于大规模网络而言，辨识准确度趋于平台期，下降并不显著。（3）DBNM-BCMI 算法的辨识准确率显著高于互信息算法和LGC 算法，且对于大规模网络，其辨识准确率均保持在80%以上。

4 结论

辨识复杂BNN 的功能性连接结构，从而探明各种生理、病理功能的网络结构特点。本研究对传统的互信息网络结构辨识法进行改进，提出DBNMBCMI 算法，对人工构造的具有一定生物真实性的SNN 结构进行辨识。DBNM-BCMI 算法通过计算条件互信息，删除动态贝叶斯网络中非依赖结构和导出依赖结构，最终辨识出动态贝叶斯网络的连接结构。结果表明相比较传统的网络结构辨识法，DBNMBCMI 算法对于BNN 具有较高的辨识正确率，验证了该算法的有效性。

除此之外，本研究提出的算法源于条件互信息计算和动态贝叶斯网络结构，仅需要网络多通道响应的时间序列数据，不依赖于网络结构和模型阶数的先验知识，且能较好地适应生物神经元网络的本质非线性动态，是一种较有效的非线性概率方法。相比较传统的线性方法，尤其适用于BNN的网络结构辨识任务。