视觉皮层具有层级组织结构，然而，广泛存在的循环通路和平行通路使得解析神经元群体间的信号传递方式变得颇具挑战。在此，我们追踪了从小鼠视觉层级中六个相互连接层级记录到的锋电位活动流。通过分析同步记录神经元对之间超前与滞后的锋电位时序关系，我们构建了一个细胞尺度的有向网络图。利用模块检测算法，基于共享的连接模式对神经元进行聚类，我们发现了两个分布在整个层级结构中的多区域通信模块。这些模块之间及模块内部的信号流动方向、层和区域分布的差异，以及不同的时间动力学特征，表明其中一个模块传递前馈感觉信号，而另一个模块则整合输入以进行循环处理。这些结果表明，多区域功能模块可能是皮层区域之外组织结构的一个基本特征，有助于信号在层级循环网络中的传播。
一、引言
新皮层中的信息处理涉及处理层级或模块之间的信号表征、转换和传递。感觉系统按解剖学层级组织。不同层级区域的定位被认为对应于由前馈和反馈投射连接的不同处理模块。尽管刺激表征（如单神经元调谐偏好和群体编码）已被广泛研究，但基于锋电位的信号在处理层级间的神经元群体中传递的原理，人们了解得少得多。
追踪分布式网络中的信号传播，需要在行为动物体内同时测量大量区域内及区域间相互作用的神经元。先前在毫秒级时间尺度上检查区间神经元通讯的研究，揭示了这些相互作用的刺激依赖性和功能特异性。其他研究使用局部场电位作为群体活动总量的指标，发现了不同频带与前馈（伽马）和反馈（贝塔）信号之间的关联。然而，局部场电位测量无法分辨单个神经元，限制了这些实验监测细胞尺度信号流的能力。宽场钙成像技术可同时获取多个皮层区域的信息，但这些信号动态变化慢得多，无法捕捉毫秒级时间尺度上的信号传输。因此，目前缺乏跨越多层级感觉处理的网络通讯的综合视图，因为具有单神经元和单锋电位精度的多区域记录直到最近才因紧凑型高密度植入式电极（如Neuropixels）的引入而变得易于操作。
在我们先前的工作中，我们构建了一个使用多根Neuropixels探针的记录平台，能够同时从小鼠视觉层级的多达八个层级中测量数百个发放的神经元（图1A）。通过分析发放单元之间的相关性，我们发现自下而上的感觉驱动过程中信号流的平均方向遵循解剖学定义的层级。与此一致，我们观察到在更高层级有更长的反应潜伏期。然而，充分的解剖学证据表明，皮层层级之间存在平行通路和循环连接。在小鼠中，初级视觉皮层（V1）直接投射到所有更高级视觉区，甚至单个细胞也能通过分支轴突并行投射到多个区域。此外，局部子网络表现出聚类连接性，这可能支持神经元群组（或集合）的共激活。
图1. 多区域锋电位记录中功能模块的识别
（A）上图：实验装置示意图。下图：来自六个视觉皮层区域的Neuropixels记录。视网膜拓扑符号图叠加在血管图像上以引导区域定位。右图：根据解剖学层级得分重新绘制的图。
（B）在漂移光栅刺激期间536个同时记录神经元的栅格图（显示了六个连续试次，每个试次包含2秒光栅呈现，随后为1秒灰度期）。
（C）六个视觉皮层区域对光栅刺激响应的归一化刺激时间直方图（PSTH）。
（D）一个V1神经元与一个LM神经元之间经抖动校正的互相关图（CCG）示例，具有尖锐峰值。连接权重为0时滞前后各13毫秒窗口内CCG幅值平均值的差值，这反映了CCG中的不对称性（方向性）。
（E）一只示例动物的有向细胞尺度连接矩阵。神经元按区域和深度排序。
（F）同一数据的聚类过程及结果。
（G）同一只小鼠的三个功能簇的连接轮廓，源神经元按区域和深度排列。
（H）按簇组织的邻接矩阵，显示模块化结构。
（I）同一示例小鼠中，簇1神经元到各簇神经元的经抖动校正的CCG。
（J）簇2神经元到各簇的经抖动校正的CCG。
（K）簇3神经元到各簇的经抖动校正的CCG。顶行：经抖动校正的CCG按连接权重从高到低自上而下排列。底行：平均CCG，误差条表示均值的标准误差。
鉴于局部和区间连接的复杂性，识别相关的信号传输模块颇具挑战性。经典的解释是每个皮层区域对应一个处理阶段，执行局部计算并将输出信号依次从一个区域传到下一个区域。然而，这种标准观点可能仅是一阶描述。由于密集的侧向和循环连接，主要的通讯枢纽可能实际上跨越多个解剖区域。因此，一个关键的基础步骤是在多区域细胞层面上将系统细分为功能相关的处理模块。
深圳市富临神通科技有限公司是Open Ephys中国代理商，我们在神经科学领域拥有丰富经验。我们为客户提供：Headstages、脑机接口芯片、Miniscope、记录系统、刺激/记录系统、SPI线。
在此，我们采用这样一种方法，即基于功能连接将皮层网络分解为模块。利用来自六个相互连接、层级组织的皮层区域的多区域记录，我们基于神经元间超前-滞后锋电位时序的统计量来推断信号流。将每个神经元视为一个节点，我们构建了一个带方向权重的邻接矩阵，然后使用无监督聚类算法，根据神经元共享的功能连接模式来揭示多区域神经元集合。该方法揭示了两个跨越皮层层级的主要视觉参与模块。多条证据表明，一个模块主要参与感觉信息的前馈分发，而另一个模块则更多地参与信息的循环整合。
二、结果
我们部署了多达六根Neuropixels探针，从小鼠视觉皮层的六个区域记录神经元。每个区域都有各自的视觉空间映射图，并处于由解剖学和生理学共同确定的不同层级位置。初级视觉皮层（V1）位于层级底部，其后依次为嘴外侧区（RL）、外侧内侧区（LM）、前外侧区（AL）、后内侧区（PM）和前内侧区（AM）（图1A）。总体而言，本研究中的记录批次每项实验获得632 ± 18个同时记录的神经元（即分选的“单元”），分布于各皮层层次和区域中（n = 19只小鼠，均值±标准误）（图1B）。我们使用全场漂移光栅刺激来提供强烈的自下而上的感觉输入，以诱发单位时间内大量锋电位的发放。与已知的小鼠视觉层级一致，各区域的平均反应潜伏期遵循顺序递进（图1C）。然而，所有区域在视觉反应的大部分时间内都处于共同激活状态，从而为循环交互提供了机会。为便于功能连接分析，我们依据最低发放率（>2个锋电位/秒）和感受野位置（感受野中心距显示器边缘至少10度）对数据集中的神经元进行了筛选。筛选后剩余n = 3,487个单元，占所有小鼠记录单元总数的29%（详见STAR方法）。
01.多区域功能模块的识别
为表征与信号传输相关的快速时间尺度交互作用，我们使用经抖动校正的互相关图（CCG）分析量化了每只小鼠中所有神经元对之间的锋电位时序互相关性；该方法捕获抖动窗口（25毫秒）内两个神经元之间的相对锋电位时序，但去除刺激锁定信号及长于抖动窗口的相关性。对于每一对神经元，我们通过计算0时滞前后各13毫秒窗口（约抖动窗口的一半）内CCG曲线下积分的差值来确定其连接权重（图1D），该指标反映了短时滞下CCG的不对称性。因此，该权重的符号描述神经元对之间的“信号流方向”（时序上超前或滞后）。请注意，我们未根据CCG的形状（如峰值尖锐度）进行筛选，因此该测量可能包含所选时间窗口内的单突触和多突触连接。通过计算所有神经元对的连接权重，我们生成了一个描述所有同时记录神经元有向功能交互的连接矩阵（图1E）。这些功能连接矩阵显示出非随机结构（图S1）。观察表明存在具有相似功能连接模式的可分离神经元群组。接下来，我们试图通过对其连接轮廓进行聚类来算法性地揭示这些群组。
为系统识别具有相似连接的源神经元集合，我们将交互矩阵按每个源神经元到记录网络中所有神经元的连接作为特征进行聚类（图1E和1F；图S2）。该过程产生了三个稳健的源神经元簇（图S3）：簇1大多为弱连接；簇2中的神经元以强正连接权重为主，表明它们倾向于引领（或驱动）网络活动；而簇3中的神经元以强负连接权重为主，表明它们倾向于跟随（或被驱动于）网络活动（图1G和1H）。可视化各不同簇中所有神经元之间的经抖动校正的CCG（图1I–1K）显示：（1）与簇1单元的CCG通常非常弱；（2）簇2与簇3之间的CCG不对称，簇2单元领先于簇3单元的活动；（3）簇2或簇3内部的CCG既有领先对也有滞后对，平均而言更为平衡。
这三个簇在我们检测的每只小鼠中均被一致识别，表明它们是视觉皮层的一个核心组织特征（图2A）。鉴于各簇正负连接权重的偏向性（及其方向性的功能含义），我们将簇2源神经元称为“驱动模块”，簇3源神经元称为“被驱动模块”（图2B）。支持这些簇稳健性的是，我们使用谱聚类和双聚类算法也观察到了类似的网络模块（图S3）。驱动和被驱动单元的相对比例较为相似（图2C）。这些模块在高对比度漂移光栅刺激期间显而易见（图2A和2B），但在自发活动期间（灰屏）则不太明显（图2D和2E）。然而，我们确实观察到在自发条件下存在一种与漂移光栅条件相似的功能连接模式，当我们扩大y轴尺度以揭示这些弱连接时（图2E，右图）。这表明，导致驱动和被驱动模块分离的网络结构很可能在没有视觉刺激的情况下也存在，但在自发条件下由于发放率低一个数量级而较不明显（自发活动期间平均发放率为0.51 ± 0.0062个锋电位/秒，而漂移光栅期间为4.7 ± 0.059个锋电位/秒）。
图2. 多只小鼠中观察到的一致有向功能模块
（A）三只示例小鼠在光栅刺激期间的邻接矩阵（单元按簇ID分组）。
（B）漂移光栅刺激期间各簇的平均功能连接权重，跨小鼠取平均。簇2和簇3的值显著不同于0（双尾单样本t检验与0比较，簇1：T = 1.1，p = 0.29；簇2：T = 10.9，p = 2.3e−9；簇3：T = −7.6，p = 4.7e−7）。
（C）每只小鼠中各簇单元所占比例。统计采用学生配对t检验（经Bonferroni校正）。–，不显著；∗，p < 0.05；∗∗∗，p < 0.001。
（D）同样三只小鼠在自发活动期间的邻接矩阵（单元排序与A相同）。
（E）自发期间（灰屏）的平均权重。双尾单样本t检验与0比较，簇1：T = 3.4，p = 0.0036；簇2：T = 6.8，p = 2.5e−6；簇3：T = −2.5，p = 0.024。右图为放大视图。灰点表示各只小鼠。数据表示为跨小鼠（n = 19只）的均值±标准误。
这两个模块中的神经元位于哪些皮层区域和层次？有趣的是，驱动模块和被驱动模块中的神经元在空间上分布于皮层层级的所有层级，而非局限于特定区域（图3A）。尽管如此，两个模块中神经元的比例显示出区域偏向性。总体而言，驱动神经元沿层级递减（图3A；与各区域解剖学层级得分的斯皮尔曼相关：r = −0.94，p = 0.005），而被驱动神经元则递增（被驱动模块的斯皮尔曼相关 = 0.94，p = 0.005）（图S4A和S4B）。驱动和被驱动神经元均存在于所有皮层层次，但也显示出层次偏向性。驱动神经元富集于中层和浅层，而被驱动神经元在较深的皮层层次中更为常见（图3B）。先前的解剖学示踪实验表明，介导前馈投射的神经元往往起源于浅层，且浅层中前馈投射神经元的比例沿层级递减。在我们的分析中，驱动模块中的神经元遵循相同模式，表明该模块可能参与前馈处理，而被驱动模块可能更多地参与循环处理。
图3. 不同模块的区域分布、发散度和收敛度存在差异
（A）驱动和被驱动模块在视觉皮层层级中的分布（见STAR方法）。
（B）所有小鼠中驱动和被驱动模块的皮层层次偏向性。（A）和（B）中的误差条表示自助法标准误（见STAR方法）。
（C）驱动模块（红点）和被驱动模块（蓝点）中神经元相对比例（以点大小表示）在视觉皮层层级中的示意图。箭头示意驱动和被驱动模块在汇聚和发散模式上的不同。
（D）驱动和被驱动模块跨区域的发散度分布箱线图。箱体边界为第一至第三四分位数，中间线表示分布的中位数。须线对应5%和95%百分位数。
（E）驱动和被驱动模块的收敛度分布箱线图。
02.不同模块的信号整合与分发
节点收敛度（来自其他神经元到一个神经元的“输入”）和发散度（从一个神经元到其他神经元的“输出”）是两种网络属性，分别差异化支持信号的整合与分发。驱动神经元的连接权重以正为主，而被驱动神经元则以负为主，这表明驱动神经元平均而言应具有更高的发散度（相对于输入，有更多的输出投射）。为量化这一点，我们计算了每个神经元的发散度，即正连接（从该神经元出发的输出投射）数量相对于网络规模的比例。同样，神经元的收敛度计算为负连接所占比例。在随机连接的全对全网络中，发散度和收敛度的默认值为0.5。驱动模块神经元比被驱动神经元具有更高的发散度（图3D；跨区域双因素方差分析，模块间F = 1029.5，p = 5.4e−184；区域间F = 35.4，p = 7.9e−35；交互作用F = 3.6，p = 0.003）。相比之下，被驱动模块神经元具有更高的收敛度（图3E；模块间F = 1029.5，p = 5.3e−184；区域间F = 35.3，p = 8.2e−35；交互作用F = 3.6，p = 0.003）。因此，从网络角度来看，驱动模块中的神经元更适合分发信息，而被驱动模块中的神经元则更适合信号整合。我们还发现发散度和收敛度在视觉层级上存在显著相关性（图S4）。例如，驱动神经元的发散度沿视觉层级逐渐降低（斯皮尔曼r = −0.94，p = 0.0048；皮尔逊r = −0.84，p = 0.038）（图S4C）。这表明早期视觉区域中的驱动单元比更高视觉区域中的对应单元具有更强的信号分发能力。这种层级相关性并非我们聚类的简单结果，因为区域标签打乱的控制分析未显示出与层级的显著相关性（区域间方差分析F = 0.60，p = 0.70；斯皮尔曼相关r = −0.09，p = 0.87），但驱动和被驱动神经元之间在收敛度上仍保留了显著差异（发散度模块间方差分析：F = 1553.7，p = 1.6e−254；收敛度模块间方差分析：F = 311.7，p = 3.9e−65）。
被驱动神经元具有更高的收敛度表明它们在视觉处理过程中整合信号。在视觉系统中，受漂移光栅相位调制的多个类简单神经元的汇聚可产生对不同光栅相位具有容忍度的类复杂神经元。类简单和类复杂神经元可通过调制指数（MI）来量化，该指数描述神经元在其偏好时间频率下的调制程度；类简单细胞显示较大的MI，而类复杂细胞的MI较小（极限情况下接近于零）。为检验被驱动模块是否包含更多与它们相对驱动模块占主导的汇聚输入相一致的类复杂神经元，我们计算了单个神经元的MI并比较了两个模块（图4A和4B）。总体而言，被驱动群体比驱动群体更偏复杂（p = 3.3e−24；曼-惠特尼U检验）。由于V1中存在更多类简单神经元（反映在图4C的双峰分布中），我们还仅在比较高级视觉区域时检验了驱动和被驱动模块之间的MI。我们发现即使在移除V1后，驱动神经元的MI仍显著高于被驱动神经元（p = 0.00011；曼-惠特尼U检验）。该结果与驱动神经元具有更高发散度一致，并表明它们在视觉处理通路中处于功能上更深的位置。
图4. 被驱动模块中类复杂神经元多于驱动模块
（A）漂移光栅刺激示意图及来自驱动和被驱动模块的两个示例单元的栅格图。
（B）驱动（红色）和被驱动（蓝色）神经元的相位调制指数（MI）分布。
（C）不同脑区中属于不同模块的神经元的平均MI。误差条表示标准误。显著性采用曼-惠特尼U检验。∗，p < 0.05；∗∗∗，p < 0.001。
03.模块间的信号传递
在定义了两个跨越视觉层级的稳健网络模块后，我们接下来试图进一步探索这些可分离网络集群内部及之间的信号流。为此，我们对先前定义的有向权重矩阵的子网络进行了可视化和分析（图1和图2）。我们首先通过检查从驱动到被驱动神经元（图5A，上图）以及从被驱动到驱动神经元（图5A，下图）的连接所定义的子网络，聚焦于模块间的信号传输。正如预期，这些模块之间的连接在很大程度上是单向的；即，来自各源区域（包括视觉层级顶部的区域，如PM和AM）的驱动神经元向被驱动神经元发出输出连接，而各源区域中的被驱动神经元则接收来自驱动神经元的输入连接（图5B）。这些连接的不对称性表明信号从驱动模块到被驱动模块大体呈单向流动。
图5. 不同模块之间的单向信号流
（A）一只示例小鼠的模块间子网络有向连接矩阵（上图：驱动［n = 74］→ 被驱动［n = 81］；下图：被驱动 → 驱动）。
（B）（A）中模块间子网络源区域（顶部标记）与所有其他区域之间强连接（为可视化目的取权重值 > 10⁻⁶）的图表示。图中每个节点代表皮层中的一个记录位点。视觉区域身份以节点颜色表示，并按解剖位置排列。各区域内的节点按浅层至深层顺时针排列。线条颜色表示连接的方向性（红色，输出投射；蓝色，输入投射）。
（C）模块间子网络跨小鼠（n = 19只）平均的区域出入指数。
（D）源区域V1、LM和AM的模块间子网络（驱动到被驱动和被驱动到驱动）叠加的图表示。左图：一只示例小鼠源区域V1的模块间子网络有向图表示。强连接（为可视化目的取权重值 > 2∗10⁻⁶）的图表示。线条颜色表示连接的方向性。红色表示V1中驱动神经元的输出 → 其他区域中的被驱动神经元；蓝色表示其他区域中驱动神经元的输入 → V1中的被驱动神经元。源区域V1的输入和输出有明显分离，且这种分离依赖于皮层深度。中图：源区域LM。右图：源区域AM。
（E）两个模块归一化PSTH的群体平均值。误差条表示标准误，因太小而不可见。
在确定功能模块时，我们的聚类算法仅考虑了连接模式，而忽略了它们的区域标签。然而，皮层区域的层级组织与信号传播相关是被普遍接受的。因此，为从各子网络中单个皮层区域的角度量化这种信号流，我们定义了一个称为“区域出入指数”的度量，该指数描述给定子网络中来自源区域的输入连接与输出连接的相对比例。出入指数为1表示所有连接均为对该源区域的输入，指数为−1表示所有连接均为输出（见STAR方法）。对于驱动到被驱动子网络的所有区域，出入指数接近−1（图5C，上图），表明几乎所有连接都是输出。相比之下，跨小鼠的被驱动到驱动子网络中每个区域的出入指数接近+1（图5C，下图）。这些数值是信号从驱动模块向被驱动模块单向流动的结果，并且与对权重矩阵的观察预期一致。然而，对模块间图（图5B）的更仔细评估表明，连接在各区域的层次分布并不均等。为比较这些连接的解剖布局，我们将各源区域的两个子网络（驱动到被驱动和被驱动到驱动）叠加；结果显示在皮层深度上输入和输出投射有明显的分离（图5D），输出投射集中在源区域的中/浅层，而输入投射集中在深层。这些观察反映了我们先前发现的偏向性，即驱动神经元在浅层丰富，而被驱动神经元则偏向深层（图3B）。
从驱动到被驱动模块的方向性通讯表明，这两个模块在视觉刺激期间可能被顺序激活。为验证这一点，我们量化了每个模块的刺激反应潜伏期，发现驱动模块中的锋电位早于被驱动模块（驱动模块达到峰值的时间为60.5 ± 0.3毫秒，而被驱动模块为80.0 ± 0.6毫秒；秩和检验统计量 = −29，p = 7.1e−186）（图5E）。我们进行了仿真模拟，表明这些模块中神经元之间存在短暂时间尺度的相关性并不必然意味着模块间在刺激触发平均响应中存在时间偏移，反之亦然（图S5）。因此，模块间子网络所暗示的方向性通讯得到了模块间信号时间进程的支持。
04.模块内部的信号传递
虽然模块间的连接在很大程度上是单向的，但每个模块内部的连接包含两个方向，表明存在更多的循环处理可能性（图6A和6B）。这得到了模块内网络平均CCG的支持，其在0时滞处显示出一个峰值（图1J和1K）。为进一步探索模块内通讯结构的区域特异性，我们分析了模块内子网络（图6C）。有趣的是，对于驱动模块，区域出入指数沿层级递增：V1具有负的出入指数，主要向其他区域中的驱动神经元发出输出连接；相比之下，AM中的驱动神经元接收到的输入多于输出，如正的出入指数所示（图6C，上图；与层级的斯皮尔曼相关为0.83，p = 0.04）。在被驱动模块内部，连接更为平衡：各区域的出入指数接近0，且与解剖学层级得分无显著相关（斯皮尔曼相关p = 0.21；图6C，下图）。因此，驱动神经元之间的信号流（而非被驱动神经元之间的信号流）与视觉层级相关；这一结果不能用驱动神经元跨区域的偏向性分布来解释，因为被驱动神经元的比例也跨区域存在偏向性（图S4A和S4B）。
图6. 各模块内部的信号流表明前馈与循环过程的分离
（A）与图5同一示例小鼠的模块内子网络有向连接矩阵（上图：驱动 → 驱动 ［n = 74］；下图：被驱动 → 被驱动 ［n = 81］）。
（B）（A）中模块内子网络源区域（顶部标记）与所有其他区域之间强连接（权重值 > 10⁻⁶）的图表示。
（C）所有19只小鼠模块内子网络的平均出入指数。
（D）所有19只小鼠不同模块跨区域的反应潜伏期（首次锋电位时间）。
（E）数据揭示的信息流三个时间阶段示意图。颜色深浅表示激活水平（越深表示越强）。
（C）和（D）中的误差条表示标准误。
这些模块内的连接模式表明，驱动模块神经元沿层级将前馈信号中继给其他区域中的驱动神经元，而被驱动模块内的双向连接则定位于介导区域间的循环交互。与此观察一致的是，驱动模块内的平均视觉诱发反应潜伏期沿区域层级系统性增加（图6D；与解剖学层级得分的相关：斯皮尔曼r = 0.83，p = 0.04）。相比之下，被驱动模块中神经元的视觉潜伏期相对于驱动模块有所延迟，但未表现出跨层级的组织性递进（斯皮尔曼r = 0.14，p = 0.79）。
综合起来，这些结果表明一个工作模型，其中这两个模块支持信号传播的不同阶段：一个模块沿层级传递关于外部刺激的前馈信号；另一个模块在接收驱动模块输入的基础上整合并处理循环信号（总结于图6E）。
05.不同模块间群体锋电位的时域离散度存在差异
前馈网络的建模研究探讨了锋电位信号如何跨顺序处理模块传播。根据包括连接密度和突触权重在内的各种网络参数，连续阶段可以传播同步活动或发放率的异步波动。培养网络的体外实验报告指出，随着信号跨顺序处理模块中继，模块内锋电位活动的同步性降低。然而，来自行为动物体内网络的证据很少。我们在此描述的分布式功能模块可能代表信号处理的连续阶段。在此背景下，我们逐试次比较了驱动或被驱动模块中神经元群体对首次刺激诱发锋电位的时域离散度（图7A–7C）。驱动群体对刺激起始的响应时域离散度小于被驱动模块（图7C；驱动模块的起始时域离散度为14.0 ± 3.9毫秒，被驱动模块为17.1 ± 3.4毫秒，n = 6,000个试次跨小鼠合并，学生双尾t检验统计量T = 43.6，p < 0.001）。有趣的是，单试次群体响应的起始离散度在驱动和被驱动模块之间跨试次呈弱相关（皮尔逊r = 0.26，p < 0.001，n = 6,000），这表明群体抖动的总体水平可能随试次变化或随时间波动（例如基于脑状态）。我们还计算了逐试次群体响应的时间展宽（图7D–F），发现驱动模块更为紧凑（图7E；学生双尾t检验统计量T = −31.4，p = 3.3e−211，n = 6,000个试次跨19只小鼠）且发生更早（图7F；T = −17.4，p = 1.4e−67）于被驱动模块。这些结果表明，驱动与被驱动模块之间的传输与锋电位时域离散度增加及刺激起始同步性降低相关联。这与在处理链更深层循环交互增多的概念相一致。
图7. 不同模块间群体锋电位的时域离散度存在差异
（A）一只示例小鼠单个试次中驱动和被驱动模块同时记录神经元的栅格图。每行是一个神经元的锋电位时间。红色叉号标记首次刺激诱发锋电位。
（B）（A）所示试次中各模块的首次锋电位潜伏期分布。展宽量化了模块起始潜伏期的离散度。
（C）比较驱动和被驱动模块跨试次（n = 6,000个试次，来自19只小鼠）的模块起始时域离散度联合图。
（D）左图：一只示例小鼠驱动神经元（n = 91）和被驱动神经元（n = 94）的归一化PSTH。右图：同一只小鼠对应的归一化群体响应（驱动：峰值在66毫秒，展宽为15.24毫秒；被驱动：峰值在66毫秒，展宽为15.24毫秒；n = 75个试次）。
（E）群体响应峰值展宽跨所有试次的联合散点图。
（F）群体响应峰值潜伏期跨所有试次的联合散点图。白点表示单只小鼠的均值，带跨试次标准差。
三、讨论
01.分发与整合信号的分离神经元群体
我们的结果提供了小鼠视觉皮层信号流的多区域视角。与每个层级水平沿区域链逐点顺序传输信号的简单前馈模型相反，我们提供了基于模块的信号传播的证据，该传播同时涉及层级中多个水平的神经元。我们称之为驱动模块的神经元集群在感觉处理流中更早参与，而被驱动神经元则较晚被募集；这表明这两个模块可能代表不同的处理阶段：一个模块分发前馈感觉信号，而另一个负责循环整合。这一结论得到了这两个模块之间多种差异的支持，包括它们的功能交互、层次和区域分布、网络收敛与发散度，以及模块内的时间协调。
我们在实验和分析过程中做出了若干生理学和计算假设。首先，我们选择了在漂移光栅刺激驱动下平均发放率大于2个锋电位/秒阈值的神经元；这使我们的分析限定于受刺激激活的单元。其次，由于视觉皮层区域间的解剖连接具有视网膜拓扑性，为最大化区域间可测量的功能连接数量，我们选择了感受野中心距显示器边缘至少10度的单元。该选择标准确保我们所研究的神经元具有受视觉刺激激发的视觉感受野位置，并有利于功能连接估计中感受野的重叠。为验证，我们在应用任何发放率或感受野标准之前分析了所有记录单元的完整功能网络，发现未通过我们选择标准的单元连接权重非常低，因此落入以弱连接为主的簇1；因此，这些选择标准不太可能影响我们关于驱动（簇2）和被驱动（簇3）模块的结论。第三，在根据经抖动校正的CCG估计的功能连接构建邻接矩阵时，我们未基于CCG的形状（如选择假定单突触连接的尖锐峰值轮廓）进行分析。相反，我们在±13毫秒的时间窗口内使用了有向功能权重的宽容定义（CCG在零时移周围的不对称性），该窗口可能包含多突触效应。第四，我们仅使用漂移光栅刺激（以不同方向呈现并多次重复）和灰屏时期（自发）研究了功能网络；这些视觉输入不具备更生态化的视觉输入中存在的丰富时空结构。最后，我们的记录是在被动观看视觉刺激且未参与任务的小鼠中进行的。这可能使功能网络偏向反映自下而上的感觉驱动信号，而非来自高级皮层区域的自上而下信号。确实，在未向小鼠呈现刺激的自发条件下，模块之间的连接权重显著降低（图2），这与我们先前的发现一致。我们的结果和发现应在考虑这些注意事项和限制的条件下进行解读。
尽管如此，我们的研究揭示了小鼠视觉皮层中一种新的组织原则，即分布在层级各水平的神经元子集共享功能连接并参与多区域信号传输模块。从类简单神经元的组成到解剖分布和时间动力学的独立测量均与这些模块的不同处理阶段一致，表明这些多区域模块是小鼠视觉感觉处理的基本组织方式。我们的通用方法可应用于其他脑区和感觉系统，以检验这一组织原则的普遍性。
02.识别多区域功能模块的方法
社区检测方法已用于人脑成像研究，以在毫米尺度上识别分布式功能模块，但这些技术很少应用于皮层的细胞分辨率网络。我们在本研究中开发的聚类算法为识别锋电位网络中的功能模块提供了一种新方法，并可能有助于剖析各种系统和情境中的网络子结构。我们识别功能模块的方法依赖于首先量化记录网络中快速时间尺度的功能连接。我们定义了所有同时记录神经元之间交互的有向而非无向邻接矩阵。然后，我们的无监督聚类方法基于推断的输入和输出连接的共享模式来识别网络中不同的神经元群体。通过允许我们的算法同时对来自所有记录皮层区域的神经元进行聚类，我们揭示了并不直接映射到单个皮层区域或层次的模块，尽管存在区域和层次偏向性。因此，我们的方法有可能识别不受标准解剖分区方案约束的功能相关且相互作用的神经元亚群。这一点至关重要，因为许多行为和认知操作很可能由这种分布式神经元集合介导。
由于兴奋性神经元主要负责区间投射，我们进行了一项控制分析，从功能网络中移除假定的抑制性细胞，并重新计算驱动和被驱动模块的各种指标。基于锋电位波形持续时间，我们识别出约19%的单元为快速发放（FS）（假定的抑制性神经元），其余为规则发放（RS）（假定的兴奋性细胞）。FS单元同时存在于驱动和被驱动模块中。我们发现仅含RS网络和解剖分布及属性与包含RS和FS单元的完整网络无差异（图S6）。
03.与先前工作的比较
先前的工作基于层次特异性前馈和反馈连接提出了皮层层次与区域之间通讯的经典模型。具体而言，在单个柱内，输入接收于中层（初级感觉区）或浅层（高级皮层区域），然后从浅层传播到深层。在皮层区域之间，前馈连接主要起源于浅层，而反馈连接主要来自深层。我们发现驱动模块神经元集中在中层和浅层，其相对比例沿视觉层级逐渐减少，这与前馈投射神经元的解剖分布一致。被驱动神经元更可能出现在深层，接收来自同一区域和其他区域驱动神经元的功能输入，并与其他区域的被驱动神经元形成循环连接。这种分布类似于解剖学定义的反馈神经元。然而，重要的是要记住，解剖连接提供了约束功能网络的框架，但功能网络依赖于刺激和脑状态，可能随时间变化。因此，尽管我们观察到功能定义的驱动和被驱动模块的解剖分布与参与前馈和反馈过程的神经元分布之间存在相似性，但这可能是功能网络的感觉驱动状态的结果，在不同的行为和认知背景下可能并非总是成立。
几种皮层处理理论，包括预测编码和反流假说，假设了分离的前馈和反馈回路。我们描述的功能模块可能对这些理论有所启示。超出前馈和反馈回路之外，涉及分发与整合信号的分离神经元群体的概念在感觉处理层级中尚未得到充分探索，该领域的主导框架是顺序前馈处理层。我们的发现可能对皮层网络架构具有有趣的计算意义。这些模型不是构建具有映射到区域的顺序阶段的层级网络，而是可以包含参与具有不同计算特性的不同子网络的穿插神经元群体。
四、未来方向
我们在此识别的功能模块的解剖学基质尚待确定。在小鼠视觉皮层中，单个V1神经元常常产生分支投射到多个高级视觉区域。如果这些投射是模块特异性的，则这种解剖连接可能构成我们观察到的功能网络结构的基础。有证据表明前馈和反馈回路由分离的神经元群体介导，并且V1神经元之间局部连接的规则取决于它们投射到的区域；这些解剖区室化的观察可能有助于解释驱动和被驱动模块之间的划分。最终，因果扰动实验将是确认这两个模块功能相关性并阐明其信号传输确切机制所必需的。由于我们描述的模块高度分布于不同皮层区域和层次，对这些模块中神经元的精确标记和操作在技术上将具有挑战性。细胞类型和层次特异性遗传工具如Cre系目前不足以研究这些模块。允许随机访问激活分布式和混杂细胞的光刺激方法可能对此目的有用。该方法需要在进行光刺激之前识别驱动和被驱动模块。原则上，这可以通过大规模双光子钙成像来实现，以成像跨视觉皮层区域分布式神经元群体的活动。
我们描述的模块必然受到解剖学的约束；然而，脑区内和脑区之间存在精细的解剖投射，这些投射可能以任务和状态依赖的方式参与信息的灵活动态路由。未来的研究可以使用更复杂的刺激来探究模块网络结构如何依赖于不同类型的输入，特别是自然视觉刺激如时空电影和视觉运动反馈（例如虚拟现实）。此外，行为状态和认知任务条件也可能重构小鼠视觉皮层的功能网络。自上而下的信号是否利用了我们在此观察到的自下而上的功能结构（可能通过调节驱动和被驱动模块之间的传输增益）仍有待确定。自上而下的信号是否或如何重新配置回路以参与不同于我们在被动自下而上刺激期间观察到的驱动和被驱动模块的分离神经元子集也是未知的。最后，解剖层级远不止于视觉皮层。需要对这一网络进行更完整的采样，以确定驱动和被驱动模块在小鼠皮层和丘脑系统中延伸的程度，以及更广泛地，其他物种的大脑中是否存在类似的功能模块。
五、方法详情
01.小鼠
小鼠饲养于艾伦研究所动物设施，并按照艾伦研究所机构动物护理与使用委员会批准的方案进行使用。共使用了四种小鼠基因型：野生型C57BL/6J（杰克逊实验室）（n = 11），或Pvalb-IRES-Cre（n = 1）、Vip-IRES-Cre（n = 2）和Sst-IRES-Cre（n = 5）小鼠（均为内部繁殖，并与Ai32光遗传学通道视紫红质报告基因系杂交）。手术后，所有小鼠单笼饲养，并维持在12小时反相光照周期下。所有实验均在黑暗周期内进行。
02.数据采集
实验数据采集遵循先前描述的操作流程。下文提供了这些方法的概要。
03.手术方法
此处使用的所有手术方法与先前研究相同。简言之，每只动物均植入一个钛制头部固定架，该支架允许以可重复的构型进行头部固定。为植入头架，小鼠首先用5%异氟烷麻醉（1-3分钟），并置于立体定位仪上。手术期间异氟烷维持在1.5%–2.5%。体温保持在37.5°C。使用卡洛芬进行疼痛管理（5-10 mg/kg，皮下注射），并使用阿托品抑制支气管分泌和调节心律（0.02-0.05 mg/kg，皮下注射）。将头架置于颅骨上，并用白色粘固剂固定。待粘固剂干燥后，将小鼠置于定制夹具中，使颅骨旋转至20°角，以便于在视觉皮层上方制作颅窗。移除直径为5毫米的圆形颅骨片，并进行硬脑膜切开术。脑表面覆盖直径为5毫米的圆形玻璃盖玻片，其1毫米边缘延伸至完整颅骨之上。盖玻片底部涂覆一层硅胶以减少与脑表面的粘附。在手术结束、麻醉恢复前，将小鼠转移至拍照记录站，以获取颅窗的空间配准图像。
Neuropixels记录在初次手术至少四周后进行。移除颅盖玻片，并替换为含有与六个视觉皮层区域对齐的孔洞的插入窗。首先，用异氟烷麻醉小鼠（诱导3%–5%，维持1.5%，100%氧气），并用眼部润滑剂保护眼睛。体温维持在37.5°C。轻轻移除颅窗以暴露脑组织。然后将带有探针穿透孔的插入窗置于头架槽内并用粘固剂密封。在窗下注入琼脂糖混合物。该混合物由高EEO琼脂糖、认证低熔点琼脂糖和人工脑脊液组成。该混合物经优化，既足够坚固以稳定脑组织并尽量减少探针漂移，又足够柔韧以允许探针穿过而不弯曲。在插入窗的孔洞上方添加一层硅油以防止琼脂糖干燥。将3D打印的塑料盖旋入头架槽内以防笼内碎屑进入。该程序结束后，小鼠返回原笼中休息1-2小时，然后进行Neuropixels记录。
04.内源信号成像
内源信号成像（ISI）在初次手术后约15天、实验前25天进行。ISI用于生成视野（即刺激显示器上的坐标位置）到各皮层区域的视网膜拓扑映射图（图1A）。这些映射图用于划定功能定义的视觉区域边界，以将Neuropixels探针引导至初级和高级视觉皮层区域中经视网膜拓扑定位的位置。
05.适应训练
小鼠在隔音训练箱中进行为期两周的适应训练，训练箱内配有头架固定器、跑轮和刺激显示器。在两周期间，每只小鼠由同一操作人员处理。第一周，操作人员引导小鼠接触跑轮，并每天逐渐延长头部固定的时间。第二周，小鼠每天暴露于视觉刺激10至50分钟。在记录前两天，小鼠在记录装置上进行75分钟和100分钟的适应训练，期间观看实验中所用刺激集的缩短版本。
06.电生理实验
神经记录使用Neuropixels探针。每根探针有960个记录位点，其中一部分可配置用于任意给定时间的记录。使用最靠近探针尖端的电极位点。这些位点在70微米宽、10毫米长的探针柄上呈棋盘格图案排列。每个记录位点的信号被分成锋电位带（30 kHz采样率，500 Hz高通滤波）和局部场电位带（2.5 kHz采样率，1000 Hz低通滤波）。
实验装置设计用于将六根Neuropixels探针大致垂直于视觉皮层表面插入大脑。每根探针安装在三轴微操作器上，微操作器再安装在固体铝板（探针盒）上。将小鼠置于跑轮上并固定于头架夹持器。每根探针的尖端对准各区域中目标视网膜拓扑位置。在探针完全缩回、到达脑表面以及完全插入后，均拍摄明场照片记录。在右眼前方放置红外二向色镜，使眼动追踪相机能够在不受视觉刺激干扰的情况下工作。随后在装置前方降下黑色幕布，使小鼠处于完全黑暗环境中，仅视觉刺激显示器提供光源。
Neuropixels数据使用Open Ephys图形用户界面以30 kHz（锋电位带）和2.5 kHz（局部场电位带）采集。锋电位带和局部场电位带的增益设置分别为500倍和250倍。每根探针要么连接到专用的通过以太网传输数据的FPGA，要么连接到机箱内的PXIe卡。原始神经数据以压缩格式流式传输存档，分析前提取。
07.皮层区域定位
为确认视觉皮层区域的身份，将实验期间拍摄的探针图像与内源信号成像期间拍摄的脑表面血管图像进行比对。利用血管模式将视觉区域映射图叠加到带有插入探针的脑表面图像上（图1A）。为最大化跨区域可测量的功能连接，我们在所有区域中靶向凝视中心（RL除外，因几何形状原因靶向其质心），并依据视网膜拓扑图引导重叠的感受野。通过使用在屏幕上不同位置呈现的小Gabor斑块来映射所有分选单元的感受野，对靶向定位进行验证（见下文）。所有分析仅限于屏幕边界内具有清晰感受野的神经元。
08.视觉刺激
视觉刺激使用基于PsychoPy的自定义脚本生成，并使用LCD显示器显示，分辨率为1920 × 1200像素（21.93英寸宽，60 Hz刷新率）。刺激为单眼呈现，显示器距离小鼠右眼15厘米，在刺激畸变前覆盖120° × 95°的视野空间。每台显示器均经伽马校正，平均亮度为50 cd/m²。考虑到小鼠的近距离视角，对所有刺激应用球面畸变校正，以确保从小鼠视角观察时，整个显示器上的表观尺寸、速度和空间频率保持恒定。
09.感受野（RF）视觉刺激
感受野位置使用小Gabor斑块映射。该映射刺激由2 Hz、0.04周期/度的漂移光栅（3个方向：0°、45°、90°）组成，带有20°圆形掩模。Gabor斑块随机出现在屏幕上81个位置之一（9 × 9网格，间距10°），每次呈现250毫秒，无空白间隔。
10.用于电流源密度（CSD）的视觉刺激
用于层次估计的电流源密度使用全视野闪光刺激（一系列暗或亮全视野图像，亮度 = 100 cd/m²），每次持续250毫秒，试次间隔为1.75秒。
11.用于功能连接的视觉刺激
刺激驱动条件下的功能连接使用漂移光栅刺激测量，以4个方向（0°、45°、90°、135°）呈现，时间频率为2周/秒，对比度为0.8。每个试次中，光栅呈现2秒，随后为1秒灰屏。每种条件呈现75-100个试次。
12.锋电位分选
在锋电位分选之前，锋电位带数据经过4个步骤：直流偏移去除、中值减法、滤波和白化。首先，减去每个通道的中值以使信号围绕零居中。接下来，减去通道间中值以去除共模噪声。中值减去后的数据文件输入到Kilosort2 MATLAB包中，该包应用150 Hz高通滤波，然后以32通道为一块进行白化。滤波和白化后的数据保存到单独文件中用于锋电位分选步骤。
使用Kilosort2识别锋电位时间并将锋电位分配给单个单元。Kilosort2试图将完整数据集建模为锋电位“模板”的总和。每个模板的形状和位置被迭代细化，直到数据可以从N个模板在M个锋电位时间点以每个模板按幅度a缩放的方式准确重建。Kilosort2的一个关键特性是允许模板随时间改变形状，以解释实验过程中神经元相对于探针的运动。使用琼脂糖填充的塑料窗稳定脑组织几乎消除了与动物跑动相关的探针运动，但在约3小时的实验中仍观察到探针的缓慢漂移。Kilosort2能够随着单元沿探针轴移动而准确追踪，消除了原始版本Kilosort所需的手动合并步骤。锋电位分选步骤在使用NVIDIA GTX 1070 GPU的双处理器Intel 4核2.6 GHz工作站上以接近实时（每批次约3小时）运行。我们使用Kilosort2中的默认参数，初始阈值为12，最终通滤波阈值为8。
Kilosort2算法偶尔会将模板拟合到从原始数据中减去另一个模板后剩余的残差上，导致锋电位被重复计数。这可能造成一个单元的ISI违规数人为偏高，或附近单元之间零时滞同步性人为偏高。为消除这种人为同步污染数据分析的可能性，对Kilosort2的输出进行后处理，移除峰值时间在5个样本（0.16毫秒）内且峰值波形在5个通道（约50微米）内的锋电位。
Kilosort2生成固定长度（2毫秒）的模板，与细胞外检测到的锋电位波形的时间过程匹配。然而，对模板形状没有约束，这意味着该算法常常将模板拟合到具有不可能由动作电位相关电流流动产生的特征的电压波动上。与这些模板相关的单元被视为“噪声”，并根据3个标准自动滤除：扩展（单通道或>25通道）、形状（基于小波分解无峰值和谷值）或多个空间峰值（波形沿探针轴非局域化）。
锋电位分选步骤之后，每个批次的数据上传至实验室信息管理系统。每个数据集通过一套项目特定工作流程运行相同的系列处理步骤，以生成用于进一步分析的NWB文件。
六、量化与统计分析
01.数据集
总计来自19只小鼠的单元被纳入功能连接分析。锋电位分选、质量控制和预处理步骤遵循与先前研究相同的程序。其中19个数据集中的13个已通过AllenSDK在艾伦研究所网站上先前发布。来自19只小鼠的所有原始数据均以NWB文件形式在Zenodo上提供。每只小鼠平均同时记录632 ± 18个分选的皮层单元。我们设定发放率阈值以选择用于功能连接分析的单元。发放率（FR）定义为漂移光栅刺激开始后50毫秒至500毫秒窗口内的平均锋电位数。仅平均FR > 2个锋电位/秒的单元用于成对互相关图（CCG）计算，每只小鼠平均为356 ± 7个单元（总计n = 6773个单元）。由于功能连接随感受野位置变化，我们进一步将数据集限制为感受野中心距显示器边缘至少10度的单元（见下文视觉感受野）。经过感受野位置筛选后，每只小鼠最终有184 ± 8个单元用于最终的聚类程序（总计n = 3487个单元）。在对每只小鼠中受发放率和感受野位置约束的功能连接矩阵应用聚类后，各簇单元总数分别为：簇1 n = 1386，簇2 n = 1131，簇3 n = 970。
02.量化与统计分析
所有分析均在Python中执行。本文使用的主要分析包包括Scipy、scikit-learn、statsmodels和networkx。除非另有说明，误差条以均值的标准误计算。比较两个独立变量之间的差异时，若其分布呈类高斯分布（正态性检验），使用学生t检验；若分布非高斯，使用秩和检验。检验分布是否显著不同于0时，使用单样本t检验。比较皮层区域间各模块之间的变量时，使用双因素方差分析（ANOVA）评估模块间的主效应以及跨区域是否存在交互作用。与先前建立的小鼠解剖学视觉层级比较相似性时，计算测量变量（如首次锋电位潜伏期）与先前计算的解剖学层级得分之间的相关性，使用斯皮尔曼相关估计等级次序显著性。统计详情和p值可在结果部分或图例中找到。
03.视觉感受野
感受野使用Gabor斑块映射（每个20°；3种不同方向（0、45、90），时间频率 = 2周/秒，空间频率 = 0.04周/度），随机出现在120° × 95°显示器（1920 × 1200像素，21.93英寸宽，60 Hz刷新率）上81个不同位置（9 × 9网格，像素中心间距10°），背景为灰色。一个单元的感受野图定义为81个位置处锋电位计数的平均二维直方图，每个像素覆盖10° × 10°正方形。然后将感受野在最大响应的20%处进行阈值处理以去除潜在噪声像素。然后，拟合二维高斯函数。
将拟合阈值处理后的视觉感受野图，以估计感受野位置的中心。
04.刺激时间直方图（PSTH）
为可视化神经元群体的时间动态（图1、图3、图4和图7），将每个神经元的活动按1毫秒分箱，跨试次平均（n = 75），用标准差为3毫秒的高斯滤波器平滑，减去基线（基线期为刺激开始后0至0.03秒），并除以刺激开始后0至1.5秒内响应的最大值进行归一化。单个神经元的归一化PSTH在神经元群体内取平均；误差条表示神经元间均值的标准误。
05.功能连接
我们通过计算锋电位序列互相关图（CCG）来分析同时记录神经元对之间的功能交互。对于锋电位序列为s₁和s₂的一对神经元，CCG定义为：
其中 M 为试次数量，N 为每个试次中的时间箱数量，x₁ⁱ 和 x₂ⁱ 分别为两个单元在第 i 个试次中的锋电位序列，τ 为相对于参考锋电位的时间滞后，λ₁ 和 λ₂ 为两个单元的平均发放率。CCG 本质上为两个锋电位序列之间的滑动点积。T(τ) 为三角窗函数，用于校正由滑动窗口引起的重叠时间箱。为校正发放率依赖性，我们将 CCG 按几何平均锋电位率进行归一化。对每个漂移光栅方向分别计算归一化的 CCG，并跨方向取平均，以获得每对单元的 CCG。
经抖动校正的 CCG 通过减去原始数据集重采样版本（其中锋电位时间在抖动窗口内随机扰动）产生的 CCG 的期望值来构建。给定两个离散锋电位序列 X 和 Y，维度为 N × T，其中 X_{n,t} 为来自试次 n、时间箱 t（宽度 = 1 毫秒）中锋电位序列 X 的锋电位数，Y_{n,t} 为来自试次 n、时间箱 t（宽度 = 1 毫秒）中锋电位序列 Y 的锋电位数。X 的刺激时间直方图表示为：
通过对每个试次 n 的每个抖动窗口 L 内的锋电位数量求和来对锋电位序列进行分箱后，我们获得了一个新的矩阵 X_bin，其维度为 T' = T/L 个时间箱和 N 个试次。该 X_bin 矩阵通过将每个元素除以 X_bin 的行均值进行归一化，得到 X_bin_normalized，这是一个 T' × N 的矩阵。然后将 X_bin_normalized 中的每个时间箱复制 L 次，将 X_bin_normalized 从 T' × N 转换回 T × N 矩阵。随后我们计算原始锋电位序列 X 的刺激时间直方图抖动均值如下：
该 P(X_n|L) 表示重采样锋电位序列的概率分布，是一个 T × N 的矩阵。该矩阵直接用于计算抖动后的互相关图（jittered CCG）：
校正项 是真实的期望值，反映了原始数据集所有可能重采样的平均值。该校正项在从原始CCG中减去之前，先通过几何平均发放率进行归一化。用于生成重采样锋电位概率分布的分析公式由先前研究人员提出。该方法在保持每个抖动窗口内锋电位数量以及跨试次平均PSTH形状的同时，破坏了抖动窗口内的时间相关性。
在我们的测量中，基于先前研究选择了25毫秒的抖动窗口。该抖动校正方法既去除了反应中刺激锁定的成分，也去除了大于抖动窗口的缓慢波动。剩余的快时间尺度相关性更可能与两个神经元之间的信号传播有关。因此，经抖动校正的CCG反映了抖动窗口（25毫秒）内一对神经元之间的时间相关性。
然后，我们通过从经抖动校正CCG的（0至13）毫秒总和中减去（−13至0）毫秒的总和来计算有向连接权重（图1D）。13毫秒窗口定义为我们所用25毫秒抖动窗口的一半，同时也因为小鼠神经元之间实际的功能延迟发生在毫秒至数十毫秒的时间尺度上。所得数值表示功能连接的强度，其符号表示一对神经元之间功能连接的方向性。对所有神经元对进行计算后，为每只小鼠生成了一个有向的细胞分辨率连接矩阵（图1E）。
七、聚类
01.非随机性
我们首先通过计算图谱（基于谱图理论）来检验测量的连接矩阵中是否存在模块化结构（非随机性）。图的特征值定义为其邻接矩阵的特征值。图的特征值集合构成图谱。通过将原始连接矩阵的图谱与其打乱后的连接矩阵（x轴和y轴独立打乱）以及相同大小的随机生成连接矩阵进行比较，来量化矩阵的随机性。我们发现原始矩阵的图谱在顶部特征值上的解释方差显著高于打乱矩阵和随机矩阵，表明测量的连接矩阵具有非随机结构（图S1）。
02.定义聚类数量
使用几种互补方法确定聚类数量（图S3A）：
1. 肘部法则估计给定聚类数k的解释方差百分比。当曲线转为平缓时估计聚类数k̂。以下度量表示给定聚类Cj（包含nj个点）中所有点之间的簇内距离（成对距离）之和：
将归一化的簇内平方和相加，可得到聚类紧凑性的度量，即围绕簇均值的合并簇内平方和：
W_k随聚类数k的增加而单调递增。聚类数选择在边际增益下降（或斜率变化最剧烈）的“肘部”点处。
2. 间隙统计量（Gap statistics）试图将log(W_k)与数据的零参考分布（即无明显聚类的分布）进行比较进行标准化。最优聚类数k的估计值为log(W_k)低于该参考曲线最远的取值。该信息包含在间隙统计量的以下公式中：
其中E*表示来自参考分布、样本量为n的期望值。在考虑抽样分布后，估计值k̂将是使Gap(k)最大化的取值。
3. 聚类密度估计通过计算密度函数F(k)来估计给定k下的数据分布密度。F(k)的值为实际畸变与估计畸变之比。当数据均匀分布时，F(k)的值为1。当数据分布中存在集中区域时，F(k)的值减小。因此，k̂聚类数通过寻找F(k)的最小值来确定。
综合使用上述三种方法对k̂的估计，我们确定最优聚类数为3。
03.聚类方法
为寻找与网络其余部分具有相关连接模式的神经元，我们将每个源神经元到所有目标神经元的权重作为特征，对有向连接矩阵进行聚类（图1F和图S2）。为降低噪声，我们使用主成分分析（PCA）将连接特征投影到低维空间，仅保留解释总方差80%的前几个主成分。随后我们应用基于k-means的一致性聚类方法，以获得不受随机初始条件偏差影响的稳健聚类。首先，我们通过在不同初始条件下运行k-means 100次来构建共聚类关联矩阵（达到稳定的共聚类）。矩阵中的每个条目表示两个单元属于同一簇的概率。然后，我们使用层次聚类对关联矩阵进行聚类以确定簇标签。聚类数使用前一节所述方法确定。
04.不同聚类方法的比较
我们的一致性聚类基于k-means聚类方法，该方法根据降维PCA空间中的特征测量点的紧凑性（见上文）。我们将该聚类方法与另外两种检测邻接矩阵中模块化结构的聚类方法进行比较：谱聚类方法和双聚类方法。
谱聚类根据数据点的连通性确定聚类：相连或紧邻的点被归入同一簇。在谱聚类中，数据点被视为图的节点，聚类被视为图划分问题。然后将节点映射到低维空间，以便容易地分离形成簇。谱聚类分三步进行：1. 计算相似度图（k近邻）。2. 将数据投影到低维空间（计算图拉普拉斯矩阵及其特征值和特征向量）。3. 创建簇（基于对应于第二个特征值的特征向量为每个节点赋值，然后对给定聚类数的节点用k-means进行分割）。
双聚类（或块聚类）是一种同时对矩阵的行和列进行聚类的方法。对于m（样本）× n（特征）矩阵，该算法生成双簇，即在一组列上表现出相似连接模式的行子集。
功能连接矩阵的一致性聚类、谱聚类和双聚类结果如图S3C所示。三种方法在检测属于以不同权重模式为主的三个簇的单元方面显示出相对一致的聚类结果（图S3D）。因此，我们的聚类发现具有普遍性，不依赖于我们所使用的具体聚类方法。
05.聚类质量
我们使用了先前用于评估锋电位分选聚类质量的两种方法，来量化给定不同聚类数下神经元群体聚类的质量（图S3B）。d-prime（d′）使用费希尔线性判别分析计算，以在PC空间中找到最大分离线。d′表示目标簇与所有其他簇之间的无偏可分离性。该值越高，簇之间越可区分。命中率使用最近邻方法（n_neighbors = 3）计算，这是对每个簇中样本污染的非参数估计。对于属于目标簇的每个单元，识别主成分空间中最近的三个单元。“命中率”定义为这些单元中属于目标簇的比例。该指标基于“隔离度”指标。该值越高，每个簇中的污染越少。
06.调制指数
为估计每个单元的锋电位活动如何受漂移光栅刺激时间频率的调制，我们计算了刺激调制指数（MI），其定义为：
PS表示刺激时间直方图（PSTH）的功率谱密度。功率谱使用韦尔奇方法在10毫秒分箱的PSTH上计算。⟨PS⟩f表示所有频率上的平均功率；PS(f_pref)为该单元偏好时间频率（f_pref）下的功率。MI反映了单元在其偏好刺激频率下的视觉诱发响应功率与平均功率谱之间的差异。MI > 3对应锋电位在刺激频率下的强调制（指示类简单细胞样响应），而较小的MI值表示受刺激时间频率的调制较弱（指示类复杂细胞样响应）。
07.模块分布
为估计各模块单元跨区域的分布，我们首先合并所有小鼠的单元以增加各区域的总单元数并最小化记录偏差。然后，由于观察到跨区域的采样偏差，我们在各区域随机欠采样单元以匹配各区域的单元数。之后，通过将各区域驱动单元数除以欠采样后驱动单元总数来估计驱动神经元的分布。同样的计算用于估计被驱动单元的分布。一个模块跨区域的比例之和为1。最终结果为自助法均值（有放回抽样以匹配各区域单元数；n_boot = 100）。误差条表示各样本的自助法标准差。结果仅显示“驱动”和“被驱动”模块。簇1（非显著连接簇）单元未观察到系统性区域偏差（结果未显示）。
每个神经元模块跨层次的分布量化如下：首先根据CSD估计的第4层位置将单元分为浅层、中层和深层。为估计各模块单元跨深度的分布，我们首先合并所有小鼠的单元以增加单元数并最小化跨小鼠的采样偏差。然后，由于单元数存在深度依赖性，我们在每个深度随机欠采样单元以匹配浅层、中层和深层的单元数。之后，通过将各深度驱动单元数除以驱动单元总数来估计驱动单元跨深度的分布。被驱动单元同样处理。一个模块跨深度的比例之和为1。误差条表示各样本的自助法标准差（n_boot = 100）。
08.图生成
为生成图可视化（图5B、图5D和图6B），我们首先通过将具有相同电极上峰值通道的单元合并，将单单元连接矩阵压缩为单记录位点连接矩阵。然后，将每个位点视为图中的一个节点。为了直观表示，属于同一皮层区域的节点相邻，并按浅层到深层顺时针排列。各区域的位置由左半球视觉区域物理位置的俯视图确定（图1A）。图的边表示位点之间的连接，红线表示来自源单元的输出投射（正权重），蓝线表示返回到源单元的输入投射（负权重）。
09.发散度和收敛度
发散度在概念上类似于图的出度。其定义为从源神经元到网络其余部分（N个神经元）的显著正连接（权重 > 10⁻⁶）的比例。k_{i,+}表示从神经元i到网络的正连接数。
收敛度在概念上类似于图的入度。其定义为从网络其余部分到源神经元 i 的显著负连接（权重 < −10⁻⁶）的比例。
八、时间动态分析
01.反应潜伏期
使用两种不同测量方法来估计反应潜伏期。首先，峰值反应潜伏期定义为神经元在刺激开始后反应达到首次峰值的时间。其次，通过在每个试次中查找刺激开始后30毫秒内首次锋电位的时间来估计首次锋电位时间。如果在刺激开始后250毫秒内未检测到锋电位，则该试次不纳入分析。每个单元的总体潜伏期定义为跨试次首次锋电位的平均时间。
02.群体起始时域离散度
我们使用每只小鼠单个试次中一个模块内所有神经元的首次锋电位时间展宽（时域离散度）作为群体反应起始同步性的指标。在获得每只小鼠给定试次中群体内所有神经元首次锋电位时间的分布后，我们计算分布的中位数（t_median）。展宽通过在t_median ± 40毫秒窗口内对该首次锋电位时间分布拟合高斯函数来计算：
其中x为相对于刺激开始的锋电位时间，μ为平均首次锋电位时间的估计值，σ为该试次首次锋电位时间分布展宽的估计值。
03.首次峰值的群体响应展宽
为量化神经元集群的锋电位响应展宽，我们估计了每个试次群体PSTH的宽度。PSTH以2毫秒分箱计算，并用宽度为5毫秒的高斯核进行卷积。首次峰值的属性使用scipy.signal.find_peaks估计。峰值宽度表示峰值半高全宽，峰值高度表示峰值的最大值。神经元群体的锋电位展宽是量化信号如何通过前馈网络传输的重要参数。
04.区域出入指数
为量化从源区域输出的投射相对于返回到源区域的输入的比例，我们将区域出入指数定义为：
其中N_out为给定网络中从源区域到其他区域的连接数，N_in为从其他区域到源区域的连接数。该指数反映源区域出入度的不对称性。当该值接近−1时，源区域以输出投射（正权重）为主。当该值接近1时，源区域以输入投射（负权重）为主。当该值为0时，源区域的输出和输入连接保持平衡。
05.层次定义
我们通过首先使用同时记录的局部场电位（LFP）计算电流源密度（CSD）来估计皮层中层的深度。CSD使用刺激开始后250毫秒内的LFP，按照先前研究人员的方法进行计算。首先，我们计算每个记录位点处的平均诱发（刺激锁定）局部场电位。接下来，我们复制最上方和最下方的场电位轨迹，并跨位点对这些信号进行平滑处理。
其中V为场电位，z为垂直于层次的坐标，h为空间采样间隔（本文中为40微米）。然后，我们计算二阶空间导数。
在生成的CSD图中，电流汇以向下偏转表示，电流源以向上偏转表示。为便于可视化，我们使用二维高斯核（σ_t = 1；σ_d = 2）对CSD进行平滑处理。为找到中层，我们通过在CSD图上搜索局部最大值（第一个汇），即刺激开始后100毫秒内的第一个汇（随后为源），将其定义为输入层（中心通道）。
我们使用中层估计来计算“驱动”和“被驱动”模块的层次分布偏向性。我们将皮层层次划分为三层：中层（中心通道 ± 8个通道，即 ± 40微米）、浅层（中层以上的通道）和深层（中层以下、白质以上的通道）。
06.模拟以检验PSTH形状与CCG尖锐峰值之间的数学关系
由于我们观察到功能连接定义的“驱动”模块响应早于“被驱动”模块（图5E），我们想知道“驱动领先被驱动”的短时间尺度关系是否是平均PSTH所反映的一般潜伏期的结果。尽管我们的抖动校正方法应已去除刺激锁定成分，且观察到的方向性应仅反映短时间尺度的信号传递，但我们仍希望排除CCG中观察到的不对称性仅仅是试次平均PSTH潜伏期反映的可能性。
我们使用简单模拟进行阳性和阴性对照（图S5）。阴性对照检验两个具有相关但时间偏移的PSTH轨迹的神经元是否必然在其经抖动校正的CCG（25毫秒抖动窗口）中显示显著峰值。阳性对照检验如果我们人为引入毫秒级时间尺度相关性，两个不相关PSTH轨迹的神经元是否能在其CCG中产生显著峰值。测试的数学表达式如下：
给定两条PSTH轨迹：r₁(t)和r₂(t)，我们模拟泊松锋电位序列：
在时间 t 上重复100次，其中两条模拟锋电位序列（s₁和s₂）的PSTH分别匹配r₁(t)和r₂(t)的形状。仅在阳性对照中向s₁(t)和s₂(t)引入同步锋电位。计算s₁和s₂之间抖动校正前后的CCG。我们发现两个神经元之间的短时间尺度相关性（以CCG中显著峰值标识）不依赖于它们PSTH的形状和相对时间，而是——如预期——仅反映它们精细时间尺度的时间关系。
