放电频率细究

全文参考《Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems》（MIT Press，2005），作者为Peter Dayan和L.F. Abbott

图 1 猕猴下颞叶皮层神经元对视觉刺激的反应

上图为猕猴观看视频时其下颞叶皮层（inferotemporal cortex，ITC）单个神经元的尖峰放电（下文用“spike”或“放电”指代尖峰放电），横坐标对应的总时间为3s，一个小竖线就代表该神经元的一次放电。

接下来，我们将计算这段神经元放电的频率。

最容易想到的计算方法是用总的放电数目（counts）除以总的时间，这种方法得到的是该神经元在3s内的平均放电频率。

显然，我们不满足于此，我们还想知道每一个瞬间神经元的反应，然而，再短的瞬间也有长度，我们将这些“瞬间”称作时间窗（time window）。

针对3s的总时长，我们先选择一个100ms的时间窗，滑动该时间窗，使其首尾相连，计数每一个时间窗内的神经元放电数目，再将其在时间窗内等分，就得到了其在该“瞬间”的放电频率，将所有这些时间窗的放电频率依次拼接在一起，我们就得到了下图。

图 2 “时间窗法”计算放电频率

若需要提升图2的时间分辨率，只需进一步缩小时间窗即可。理论上来说，时间窗越小，我们对神经元活动的描述就越精确。但实际上，时间窗不能小于一次spike的持续时间（0.3至2毫秒）；不仅如此，进一步缩小的时间窗会让spike变得稀疏，从而极大地降低最后计算得到的放电频率的意义。

过大的时间窗会漏掉细节，过小的时间窗又难凸显放电率的差异，为此，我们常采用两个办法：1.实时调整时间窗大小；2.让时间窗交叠起来。

先说第一个方法，我们完全可以调转思路，从固定时间窗大小，改变放电数目到固定放电数目，改变时间窗大小，借此就能让每个时间窗中都存在足够数目的spike事件，由此计算得到的放电频率会更具可比性。

再说第二个方法，前文中提到的所有计算方法都是将时间窗首尾相连，然而，我们还可以使其交叠，让时间窗滑动的距离小于其自身长度，这样不仅能取到较小的时间窗，还能保证每一个时间窗中都有足够数目的spike。下图就是图1采用滑动时间窗计算得到的结果。

图 3 滑动时间窗计算放电频率

图2和图3中的放电频率图呈现“阶梯”状，而我们有时更想了解神经元活动的变化规律和趋势，也就是说，我们需要一条平滑的曲线来描述放电频率。

此处我们转变思路，不借助上述阶梯图来直接进行平滑操作，而是从神经元的放电特性出发，重新建立公式，绘制放电频率曲线。为了建立这样的公式，有几个基本的概念需要事先说明。

狄拉克函数的特性为：在自变量为0时取无穷大，而在其他所有自变量值时取0。上式中i为第i个spike，n为总的放电个数，t_i为第i个神经元放电的时刻，t为我们关注的时刻。也就是说，只有当t_i累加过程中取值刚好等于t时，函数才有值，在t_i为其它值时函数值都为0。也就是说，如果我们关注的时刻刚好是某个放电时刻，函数累加就有一个值，且该值为1；反之，若关注的为非放电时刻，则函数值为0，该特性使其函数非常好地契合了神经元spike的表现形式。

等式左侧为在t时刻的预估的放电频率r_approx，右侧为借助线性过滤的思想得到的积分形式。

现在我们重点关注这个过滤核，它可以是所有离散变量的加和，也就是说，所有放电频率都为1/Δt的整数倍，也可以是连续变量的积分，离t时刻越近的放电理应对放电频率的贡献更大，为了相对准确地描述距离和贡献的关系，我们选用高斯函数（高斯分布对于短时和长时活动的平滑效果较均衡）：

以高斯分布作为过滤核，代入前文的积分计算并作图，我们就能得到如下图所示的一条平滑的放电频率曲线。

图 4 “高斯过滤核“卷积求取神经元放电频率

让我们再次回到神经元放电的生物学意义上来，前文中所有的计算过程都默认取到了时刻t前后时间窗内的spike事件，但事实上，我们忽视了一个朴素的事实：特定时刻的放电只取决于在它之前的神经元放电事件。也就是说，在t时刻之后的放电不应被用来计算t时刻的放电频率。

用上式中的过滤核做积分运算并画图，我们再一次得到了一个有些不同的神经元放电频率曲线。

图 5 “半波整流过滤核”卷积求取神经元放电频率

从同样一段3s的神经元原始放电数据出发，借助不同的计算方法，我们得到了相似却又有着细微差异的神经元放电频率图。

图 6 不同方法求取神经元放电频率

前文中所有求取放电频率的基础都是我们已经获得了单神经元的放电数据，但还有一些场景，我们无法或难以直接测量得到单神经元放电的原始数据，这时候就需要借助其它易于检测的参数来预估神经元的放电频率了；不仅如此，在已获得了神经元放电原始数据后，我们仍可将预估的放电频率与真实记录得到的放电频率进行比较，从而发现二者的差异，而这些差异往往就能揭示某些被忽略的神经元活动机制。

相比单神经元的放电规律，诱发神经元活动的外界刺激（stimulus）更容易被准确测量得到。声光热电力，我们脑中的许多神经元都会对某个特定的外界刺激稳定地作出回应（response），只要准确测得了这些外界刺激的参数，我们就能预估神经元的活动规律。

上式左侧表示t时刻预估的放电频率，右侧第一项r₀表示在没有任何外界刺激或外界刺激未发生显著变化时的神经元基线放电频率，第二项的积分中，s函数表示在关注时刻t-τ到t时间内的外界刺激变化情况，而D函数则表示外界刺激到神经元放电频率的“换算关系”。（注意上式中预估的神经元放电频率r_est和外界刺激s间的关系为线性关系，后文中会再次提及这一点）

许多神经科学研究的目标就是找到这个D函数，即外界刺激对神经元活动的影响方式，这也是我们下文中讨论的重点。

等式左侧为我们翘首以盼的D函数，右侧则是一系列复杂的表达式：中间那一个分数中的分子为刺激和放电频率的相关系数，它是“时间窗长度τ”的函数；右项的分子中“r”为真实记录得到的神经元放电频率，“C函数”为每次放电前τ时刻的外界刺激均值（spike-triggered average，可参考小邹的另一篇文章Spike-Triggered Average：在神经元放电之前）。将所有的分母项看作常数即可。

图 7 由外界刺激预估单神经元的放电频率

图7中的A图为一段真实的外界图像的移动速度图；B图为某个“图像敏感”神经元在两次“观看图片”的过程中的放电raster图；C图中的实线部分为预估的该神经元的放电频率，虚线则为实际记录到的该神经元的放电频率。

上式被称作Volterra展开（Volterra expansion），该式在功能上类似于泰勒级数展开，随着等式右侧高次项的增多，预估的放电频率也会愈加趋近于实际的放电频率。

深入上式的细节：第一项为基线放电频率r₀，第二项为线性项（二者一起被称作Wiener kernel），第三项则考虑了放电前两个不同时刻外界刺激间的共同作用对神经元活动的影响，其中D₂函数也顺理成章地变成了两个时刻的函数，依此类推，也就是说，神经元响应的是在时间轴上有着互动关系的多个外界变量。

至此，我们了解了由外界刺激出发预估神经元放电频率的基本思路和过程。