基于模拟的简单神经网络结构推断

《Stochastic Processes and their Applications》：Simulation Based Inference of a Simple Neural Network Structure

　　

# 摘要

- 提出了一种新的用于神经网络胞外记录的统计分析框架：聚焦于实际网络被视为其实现的分布的参数，而不是子图或功能网络重建。
- 在胞外记录中发生的神经网络群体严重欠采样问题，至少在简单情况下，可以通过全网络模拟来解决。
- 少数被记录神经元（约1000个中的5到10个）的放电频率，足以估计一个有向Erd?s-Rényi随机模型中由Galves-L?cherbach模型描述的动态的连接概率，前提是此类设定适用于所研究的网络。

## 引言

在分析神经元动作电位序列时，有一个"房间里的大象"。虽然神经生理学家如今能够从大量胞外电极进行记录，并从原始数据中提取出由多个神经元产生的动作电位或放电序列[1]、[2]；但这些"众多神经元"仅代表构成所研究网络的神经元群体中极小的一部分。利用关联统计量，如交叉相关函数的估计[3]、[4]或随机强度模型[5]、[6]，神经生理学家试图推断所记录神经元组成的网络结构。但这一推断受到神经元群体严重欠采样的阻碍[7]。这产生了一个通常称为功能网络的"网络图像"[8]，其特征强烈依赖于记录条件（如刺激的存在与否），且其与实际网络结构的关系远非显而易见。

鉴于目前对大脑布线发展的认知[9, 第7章]和[8, 第4章]，我们认为合理假设：当神经生理学家在同一物种的多个个体上重复实验，聚焦于某一特定网络时——例如蝗虫的第一嗅觉中继站、大鼠的海马、猴子的纹状皮层——不同个体的网络可被视为同一底层随机图的实现。这是我们工作的核心假设。因此，与其重建所记录神经元组成的网络（这是一个明确的不适定问题），不如聚焦于该图的生成概率分布，因为这才是跨实验可重复的部分。在本文中，我们旨在验证一个概念证明，并将考虑一个实际脑网络的玩具模型：一个有向Erd?s-Rényi模型，涉及两个参数：神经元数量（顶点数）和任意神经元/顶点对之间的连接概率。显然，后续工作需要考虑更复杂的模型，如随机块模型[10]以及包含基序（即非独立连接）的模型[8]。

现在，在分析神经元放电序列时改变所提问题，显然至少不能解决我们最初面临的严重欠采样问题。为解决后者，我们提出使用整个网络的数值模型，从中采样少量放电序列，模拟实际实验中的操作。关键在于，我们对许多实际网络中不同神经元的形态学了解甚多[8]、[11]——这可用于约束随机图模型——我们也了解单个神经元的生理学以及其突触连接的生理学；这正是神经生物学教科书[9]如此厚重的原因！假设能够满意地执行此类全网络模拟，问题便变为从采样放电序列中选择正确的统计量进行计算，以便对所假设随机图模型的参数进行令人满意的推断，其含义将在后文进一步明确。通过这种方式，我们通过尽可能忠实地数值重现欠采样问题来解决该问题。这使我们能够估计统计量的采样分布，并观察后者是否严重依赖于随机图参数。这是我们提出的基于模拟的推断方法的核心，在我们的情况下，它直接受到[12, 第6节, 第208-210页]中最后一个例子和[13]的启发。需要注意的是，"常规"放电序列分析方法完全忽视了数十年实验研究积累起来的解剖学和生理学知识——在我们看来，这是对资源的严重浪费，也是实验者与统计学家之间进行有益对话的重大障碍——而我们的方法能够由设计而纳入此类知识。

我们的数值模型遵循Antonio Galves的"现实主义者"信条，而非涉及大量参数的详细生物物理模型的"自然主义者"方法[14]；这种方法对应于归于爱因斯坦的名言："尽可能简单，但不要太简单"。因此，我们不模拟实际神经元的三维扩展，而使用"点"神经元代替；我们图中的边被视为对实际突触网络复杂性的良好（足够）表征。我们认为神经元是"随机单元"，并采用特定的Galves-L?cherbach模型[15]来描述其动态。我们邀请读者阅读[16, 第1章和附录A]中对这些选择的批判性讨论和实证论证。

本文组织如下：第2节描述随机图模型（第2.1节），然后是神经动态（第2.2-2.4节），随后指定已知参数及其取值以及未知参数及其需从数据中估计的取值（第2.5节）。第3节描述本文中使用的通用基于模拟的推断方法。第一部分（第3.1-3.3节）解释为何经典的最大似然推断无法实现，从而导致使用一个"精心选择"的汇总统计量的采样分布。第3.4节明确说明该方法的推断部分。第4节提供了对所作关于汇总统计量的简化假设的实证论证。第5节使用模拟数据表明，我们的推断方法表现良好，尤其是与更传统的方法相比（第5.1节），并且能为估计参数提供有意义的置信区间（第5.3节）。第6节展示了结论及本工作的展望。附录A解释了为何一个比我们最终使用的汇总统计量更复杂的统计量不能提供关于我们想要估计的参数的更多信息。附录B给出了推断部分实现的所有细节。附录C给出了文中陈述的少数引理的明确证明。附录E简要概述了我们的软件实现和所采用的数值方法，并附有源代码链接。

## 片段：数据生成模型

在此，我们同时呈现神经网络的随机图建模以及此类网络的放电活动模拟作为其随时间的随机演化。该模拟设定产生一组放电序列，其中每一个都是代表一个神经元放电时间顺序的点过程的实现。此模拟设定用于两个目的：

- 首先，模拟大量数据以估计似然，见第3.4节；
- 其次，生成额外数据。

## 片段：蒙特卡洛方法参数估计

本文考虑的典型设定是一个包含$??=1000$个神经元的网络，从中"随机选择"10个神经元，记录$??=10^6$个时间步（如果我们认为离散化时间窗口为1 ms，该步数对应于大约15分钟的实时时间；与通常持续数小时的实际记录相比并不算长）。我们用于后续应用的实际网络是蝗虫的触角叶（第一嗅觉中继站）。

## 片段：估计 $?_p\{φ(Y)=φ(y)\}$

在本节中，我们首先实证验证第3.3节中所作的假设，然后提供证据表明我们可以估计 $\prod_{????=1} ?_p\{φ(Y_j)=φ(y_j)\}$。

## 片段：参数 $p$ 的推断

在此，我们展示我们方法的结果。我们使用了前述两个不同的 $φ$ 函数，并在附录A中介绍，分别是放电频率，标记为spikefreq，以及ISI分布的形状参数，标记为alpha。对 $?_p\{φ(Y)=φ(y)\}$ 的估计可以是参数化的（高斯），标记为gaussian_method，也可以是非参数化的，通过计算直方图，标记为unparam_method。第一个结果（图8）以这四种变体呈现，但在此之后……

## 结论与展望

本文所展示的方法源于两个动机：一是面对神经网络胞外记录的基本局限性——只有实际网络中极少数神经元被记录——的需要；二是希望解决一个在实验重复中答案应可重复的问题：实际观测网络被视为其实现的随机图模型的特征是什么？在讨论图结构问题时，P. Charitat | S. Geffray | C. Pouzat