背景介绍

计算神经科学中的一大假设是神经元通过集体的交互处理信息，即使是小型网络的神经元如何交互依然未知。尽管实验中常常能观测到多神经元发射模式，但大脑皮层神经元之间的联系是十分微弱的，这种联系可以被相关性[1]或突触力量[2]度量。两个团队[3,4]的近来工作，已经证明了在神经元之间交互十分微弱的视网膜，可以成功预测网络相关性的90~99%。他们都采用了一个二阶的最大熵模型来预测，模型的输入只需要发射率和成对交互。这种模型只能预测一个时间点的多神经元发射模式，而很多的实验室报告声称体外[5]和体内[6]相关状态横跨多个时间点连续出现。

作者应用二阶最大熵模型于广泛的体外皮质网络，包括人和老鼠组织的急性切片，还有老鼠的器官型和分离型培养，使用预测峰值和本地电势相关状态来衡量模型的准确性，将模型输出的状态组成序列，与实际数据中观察到的连续状态相比较。

方法

神经元\(i\)在一个确定的时间段里的活动用\(\sigma_{i}\)表示，\(\sigma_{i}\)可以取两个值：+1表示出现尖峰（一次或多次），-1表示未出现尖峰。类似地，使用\(\sigma_{i\cdot}\)的值表示节点\(i\)处的\(LFP\)是否超出阈值。\(\sigma_{i}\)和\(\sigma_{i\cdot}\)计算方法可以类比，下面仅给出\(\sigma_{i}\)的计算公式。给定发射速率，\(\sigma_i\)的期望值可以由下式确定：

其中\(\sigma_{i}^{t}\)表示节点\(i\)在时间\(t\)的活动，\(T\)是记录的时间盒子总数。节点\(i\)和节点\(j\)之间的对相互作用的期望值可以表示为：

在一个特定的时间点，具有\(N\)个神经元的网络的状态向量\(V\)是：

对于\(N\)个神经元，在一次记录中可能存在\(2^{N}\)个可能的状态可以被观察到，这些状态不是同频率的被观察到，二阶的最大熵模型的目标就是在给定发射速率\(\left\langle\sigma_{i}\right\rangle\)和对相互作用\(\left\langle\sigma_{i} \sigma_{j}\right\rangle\)，预测观察到每一个状态的概率。为了预测状态的概率分布，神经元的活动被映射到物理学的\(Ising\)模型[6]。这个模型假设每一个神经元表现像一个小条形磁铁，可以和本地的磁场和周围的邻居条形磁体交互。本地节点\(i\)周围的磁场可以表示为\(h_i\)，产物\(\sigma_i h_i\)给出交互相联系的能量。\(\sigma_i h_i\)为正值时，表示节点\(i\)与周围的磁场同谐，这种情况是能量有利的并且更可能的。节点\(i\)和节点\(j\)之间的交互用\(J_{ij}\)表示，产物$J_{ij}\sigma_i\sigma_j \(类似给出这种交互的能量，\)J_{ij}\sigma_i\sigma_j \(为正值表示神经元与它周围的邻居同谐，这种情况是能量有利的并且更可能的。需要注意，作者提出当\)J_{ij}\(为正时，节点\)i\(和节点\)j\(之间的类似活动是有利的；当\)J_{ij}$为负时，不同的活动有利。本地的磁场和交互可以做出下列的近似：

正如\(Ising\)模型中的磁体集合具有一个能量，处于状态\(V\)的\(N\)个神经元的集合的能量\(E\)为：

方程(5)中的第二个总和的前提是\(i \not= j\)。

下一步是将能量映射到概率，为了最少化不必要的假设，选择了最大熵的概率分布。对于一个给定均值的数据集的最大熵分布是一个指数分布[7]，所以观察到状态\(V_j\)的概率为：

其中分母是归一化项（统计力学中的配分函数），在\(2^N\)个状态上求和，这种概率分布中高能量的状态出现概率小，低能量的状态出现概率大。根据这个概率分布，可以计算得到最大熵模型下发射速率和对交互的均值，下标\(m\)表示模型：

其中\(\sigma_{i}\left(V_{j}\right)\)表示网络处于状态\(V_j\)下\(\sigma_i\)的值，根据模型得到的数值和来自数据的真实数值相比较。为了增强模型参数和数据参数之间的一致性，本地磁场\(h_i\)和交互\(J_{ij}\)被一种迭代式的缩放算法[8]调整，调整方法如下：

其中常数\(\alpha < 1\)被用来防止算法变得不稳定，作者发现当\(\alpha = 0.75\)时最快稳定充分。经过一次调整，一组新的能量和概率分布会计算出来，对于每一个网络调整大约\(50,000\)次迭代直到本地磁场和交互与交互线之间的误差小于0.1%。得到新的\(h_i\)和\(J_{ij}\)之后，他们会被回插入到方程(5)中计算得到新的状态能量值，新的状态能量值再回插到方程(6)中来计算观察到的每个状态的概率分布。

模型评估

单个时间点

一阶模型只考虑发射速率，概率分布为\(P_1\)；二阶模型考虑发射速率和成对交互，概率分布为\(P_2\)；\(N\)阶模型考虑发射速率和交互（从二阶直到N阶），可以捕获一个\(N\)个神经元网络的全部交互，它的概率分布\(P_N\)和数据中的概率分布完全一致。熵\(S\)的计算公式如下：

因为记录的时间大概是一个小时长，网络的规模相对较小(\(N\leq 10\))，估算熵值的错误应当十分小[9]。应当注意，一阶模型的熵值\(S_1\)永远比所有更高阶的模型熵值大，增加的交互可以减小熵值[10]。互信息量\(I_N\)是一个网络产生的熵值的总量，可以用一阶模型和N阶模型的熵值差异表示[11]：

类似地可以得到二阶模型得到的信息总量为：

二阶最大熵模型的表现可以通过互信息量所占部分来衡量，用\(f\)表示：

这个比例在[0,1]范围内变动，当取1时表示模型取得完美的表现。长期的测量记录，比例还可以表示如下[12]:

其中\(D_1\)表示\(P_1\)和\(P_N\)之间的相对熵，\(D_2\)表示\(P_2\)和\(P_N\)之间的相对熵：

多个时间点

活跃状态是指至少存在一个神经元是活动的，与之相对的不活跃状态指所有的神经元都处于平稳状态，一个序列被定义为连续的活跃状态[13]。考虑序列的两个简单统计指标：序列长度和序列大小，序列长度是连续的活跃状态的数目，序列大小是序列中出现的峰值或\(LFPs\)的次数。尽管[3]设计的模型本意不是预测序列的，但模型在单一时间点预测网络状态的成功促使作者想将其应用于网络状态序列的预测上，如果预测的效果不理想，证明状态序列时间上不相互独立，这同样表明模型需要修改以更好预测序列。

对于长度为\(T\)的记录序列，从\(\frac{T}{2}\)的位置一分为二，选取其中一半可以得到一个序列，对250个神经元网络做相同的操作。为了估计真实数据集中的多样性，作者将每一对序列点对点相减得到差异，将这些差异想加求和，最后的和值可以衡量数据的多样性。对于模型生成的序列，作者从\(P_2\)中按照概率分布取出状态，再将其随机组合成序列长度为\(\frac{T}{2}\)，类似的可以得到模型的差异和值，来衡量模型序列的多样性。模型和真实数据之间的多样性，通过模型序列点对点减去真实序列的方法得到差异和值。为了比较三种不同差异和值的分布，它们被按照升序的方式排列，可以得到累积概率分布。接着使用\(Kolmogorov–Smirnov \quad test\)来探测不同累积分布之间的巨大差异性。所有点中的最大差异性用\(D\)表示。

显著相关性

为了评估一个相关性是否是显著的，作者使用超几何分布[14]来计算神经元\(A\)和神经元\(B\)同时出现峰值的总次数\(M\)：

其中\(nA\)和\(nB\)分别表示神经元\(A\)和\(B\)各自产生的峰值数量，\(nb\)是记录的时间次数。这个分布的标准差:

当峰值和时间记录很大时（正如实验中的情况），超几何分布会趋近于高斯分布，作者就利用高斯分布来估算峰值数量的1%截止值：

基于上面得到的数据，可以计算得到1%截止值出的相关系数：

其中\(\bar{A}=\frac{n A}{n b}, \bar{B}=\frac{n B}{n b}\)，\(N\)表示\(N_{+1 \%}\)或者\(N_{-1 \%}\)。作者计算时间滞后相关性，\(\left\langle\sigma_{i}^{t} \sigma_{j}^{t+1}\right\rangle\)，作为一个可以用来将空间模型扩展成时间空间模型的参数，\(\sigma_{j}^{t+1}\)表示节点\(j\)在时间点\(t+1\)时的活动。

参考文献

[1] Bartho P, Hirase H, Monoconduit L, Zugaro M, Harris KD, BuzsakiG (2004) Characterization of neocortical principal cells and interneurons by network interactions and extracellular features. J Neurophys 92:600–608.

[2] Song S, Sjostrom PJ, Reigl M, Nelson S, Chklovskii DB (2005) Highly nonrandom features of synaptic connectivity in local cortical circuits. PLoS Biol 3:e68.

[3] Schneidman E, Berry MJ, Segev R, BialekW (2006) Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440:1007–1012.

[4] Shlens J, Field GD, Gauthier JL, Grivich MI, Petrusca D, Sher A, Litke A, Chichilnisky EJ (2006) The structure of multi-neuron firing patterns in primate retina. J Neurosci 26:8254–8266.

[5] Beggs JM, Plenz D (2004) Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J Neurosci 24:5216 –5229.

[6] Landau LD, Lifshitz EM (1958) Statistical physics. London: Pergamon.

[7] Jaynes ET (1957) Information theory and statistical mechanics. Phys Rev 106:62–79.

[8] Darroch HN, Ratcliff D (1972) Generalized iterative scaling for log-linear models. Ann Math Stat 43:1470 –1480.

[9] Strong SP, Koberle R, de Ruyter van Steveninck RR, Bialek W (1998) Entropy and information in neural spike trains. Phys Rev Lett 80:197.

[10] Cover TM, Thomas JA (1991) Elements of information theory. New York: Wiley.

[11] Schneidman E, Still S, Berry MJ, BialekW (2003) Network information and connected correlations. Phys Rev Lett 91:238701.

[12] Shlens J, Field GD, Gauthier JL, Grivich MI, Petrusca D, Sher A, Litke A, Chichilnisky EJ (2006) The structure of multi-neuron firing patterns in primate retina. J Neurosci 26:8254–8266.

[13] Beggs JM, Plenz D (2003) Neuronal avalanches in neocortical circuits. J Neurosci 23:11167–11177.

[14] Johnson NL, Kotz S, Balakrishnan N (1997) Discrete multivariate distributions. New York: Wiley.