大脑自发活动在fMRI中表现出丰富的时空特征，但是大脑在不同状态之间整合和重构的原理仍不清楚。本文通过对fMRI进行逐帧聚类将大脑自发活动状态划分成几个代表性的模式，并且证明这些模式在个体水平和群体水平均有较好的稳定性。文章进一步发现这些模式处于全脑信号波动的不同阶段，其波动可以和慢波动力学网络模型良好耦合。最后，文章探究了自闭症相关基因缺陷小鼠非典型的共激活模式和次慢波动力学的改变。本文提出了一套新的大脑静息状态评估方法，并且探究了其在脑疾病中的变化。本文发表在Current Biology杂志。
研究背景：
2012 PNAS文章提出通过静息态时间序列中少数具有高信号强度的时间帧就可以模拟出功能连接网络的结果（2012, PNAS, Xiao Liu & Jeff H. Duyn）。文中，作者先分别提取了后扣带回（PCC）和内侧前额叶（mPFC）的时间序列，找到时间序列中信号强度较高的时间点，对这些时间点的fMRI图像取均值得到的结果与基于相同种子点的功能连接模式相似。本文在此基础上通过逐帧聚类的方法将静息态时间序列划分成六种循环出现的共激活模式，并证明这种新的计算方法对大脑静息态数据的具有较好的代表性和解释力。
A：后扣带回（PCC）和内侧前额叶（mPFC）具有不同的时间序列，其高信号帧分别为Frame 1和Frame 2；
B：对时间序列信号强度取不同阈值和相关性图像进行比较，阈值为99%时相关性约0.5，阈值为85%时相关性约0.8，阈值在20%-85%时相关性基本持平；
C：左图为片层信号平均图，右图为基于种子点的相关性图（2012, PNAS, Xiao Liu & Jeff H. Duyn）
文章思路：
- 重复PNAS文章结果，证明共激活模式与功能连接网络之间的关系
- 对静息态数据集进行逐帧聚类，基于聚类结果提出新的分类模式
- 探究该分类模式在个体和群体中的稳定性
- 探究该分类模式动力学特征
- 基于该分类模式和动力学特征分析自闭症小鼠功能网络变化
实验方法：
文章共使用三个数据集，其中数据集1为本次实验获取，数据集2、数据集3均从已发表文献获取，分别用于证明该分析方法的可重复性及在疾病模型中的应用。数据集1：雄性C57BI6小鼠（12-18周龄），n = 40。数据集2：雄性C57BI6小鼠（12-18周龄），n = 41。数据集3：C57Bl6/J源，15-18周龄 Chd8+/+（雄性9只，雌性14只），Chd8+/-（雄性6只，雌性13只）。
数据集1：
异氟醚麻醉小鼠（5%诱导），插管和人工通气（2%，手术）。左股动脉插管进行持续血压监测和终末动脉采血。手术结束时，停用异氟醚，使用0.75%氟烷（halothane）代替。停止用异氟醚45分钟后开始采集功能像数据。整个成像过程中记录平均动脉血压（MABP），在功能像数据采集结束时测量动脉血气（paCO2和paO2）排除非生理噪声。所有rsfMRI数据均使用7.0 Tesla MRI扫描仪（Bruker Biospin, Ettlingen）采集，配备有BGA-9梯度线圈 （最高380 mT/m，线性转换速率3,420 T/m/s），72 mm鸟笼发射线圈和四通道接收线圈）。
结构像参数：
快速自旋回波序列，TR/TE = 1200/15 ms，矩阵192*192，FOV = 2*2 cm2, 冠状切片18片，切片厚度 0.60 mm。
功能像参数：
静息态（rsfMRI），平面回波序列（EPI），TR/TE = 1200/15 ms，翻转角30°，矩阵100*100，FOV = 2*2 cm2, 冠状切片18片，层厚0.5 mm，重复次数500帧 （n = 21）和1500 帧（n = 19），对应采集时间10和30分钟。组分析使用500次重复（10 min） 的数据进行分析。单受试者CAP分析使用1500次重复（30 min）的数据进行分析。
数据集2：
为了证实实验的可重复性，实验选取了2014年于意大利理工学院实验室采集的数据作为对照。该实验使用和数据集1同样的麻醉方案和图像参数，每次扫描300帧。相关实验结果发表于2015年NeuroImage。
数据集3：
为了评估该方法在疾病诊断中的作用，实验选取了2017年于意大利理工学院实验室采集的Chd8敲除的自闭症小鼠模型数据进行分析。该实验使用和数据集1同样的麻醉方案和图像参数，每次扫描500帧，相关实验结果发表于2018年Cerebral Cortex。数据集包括15-18周龄小鼠42只（其中Chd8+/+23只，9雄14雌，Chd8+/-19只，6雄13雌）。由于之前的实验中未显示性别对该模型动物行为和功能连接强度有显著影响，本次实验分析未考虑动物性别差异。
静息态功能连接强度计算：
功能像数据依次进行时间片层校正、头动校正、空间标准化配准至小鼠大脑模板，配准后数据空间分辨率0.1*0.1*0.5 mm3（192*192*34）。然后使用头动参数和脑脊液信号进行信号回归，并使用0.01~0.1Hz进行带通滤波，最后使用0.5 mm FWHM高斯核进行空间平滑。静息态功能连接强度的计算先提取种子点区域的时间序列，然后计算全脑逐像素时间序列相关性，针对Fisher-z变换之后的数值进行比较。统计分析阈值T = 7，对应的p < 0.01, Bonferroni 多重比较校正。
基于种子点的共激活特征（CAPs）计算：
根据共激活特征（CAPs）计算方法，先提取感兴趣区域的时间序列，从时间序列中选取了具有高信号强度的帧，然后将这些帧取均值得到基于种子点的共激活特征（seed-based mean CAPs）。根据2012 PNAS文章时间序列中信号强度在前15%的时间帧就可以达到和功能连接强度较高的相关性。本文作者选取了时间序列中信号强度最高的15%取均值（500帧取68±5帧，均值±方差）。统计分析阈值T = 7，对应的p < 0.01, Bonferroni 多重比较校正。
全脑共激活特征（CAPs）计算：
作者首先对个体fMRI数据特征进行了聚类分析。根据文献报道的结果，应先对fMRI图像进行阈值筛选，保留最大10%和最小5%的像素点再进行聚类分析。但本文发现阈值筛选后数据与未进行阈值筛选的数据相比，其聚类结果仅边缘区域有微小差异（< 1.64%），两者相似性可达0.98。故本文后续聚类分析和动态分析均在未经阈值筛选的数据基础上进行。
经预处理的fMRI时间序列先读取成4D数据（包括三个空间维度和一个时间维度ti）。然后通过k-means++算法进行时间维度的聚类，将时间点划分成k个分类，使得分类内部差异D达到最小：
其中Ci代表分类编号，μi代表每一个分类的图像平均值，d（tj，ui）代表时间点j图像和分类图像均值的差异（通过1-相关系数得到）。k-means++ 算法提供了一个优化起始点的方法，让靠后的片层也有机会被选为集群的起始点。该方法在聚类性能和速度两个方面都优于随机种子点的聚类方法。
对整个 rsfMRI 数据集 1（40只 × 500 帧）进行 k 均值聚类，使用不同时间帧中 fMRI 数据之间的 Pearson 相关性作为聚类效果评判标准，重复15 次，每次迭代 500 次。作者尝试将数据分成 2 ～40个类别，并计算增加聚类数聚类算法对整体解释力的影响。
为了评估聚类结果的可重复性，研究中使用的三个独立收集的数据集（n = 40、n = 41 和 n = 23）进行聚类分析，文章分别计算了分类数为k =6~8时三组数据的相关性。相关性计算先使用Hungarian算法对不同数据集的CAP进行匹配，然后进行空间相关性的计算 （Pearson's r），假设 CAP 相关系数 > 0.45，对应p < 10e-5，置换检验。
为了评估这些分类特征在个体水平上的灵敏性，文章使用19个被试的长时程记录数据进行分析（30 分钟，1500 帧）。先分别计算个体水平和组水平CAP模式，然后使用Hungarian算法将个体水平-组水平CAP进行匹配，计算两者之间的相似性。
研究将每个受试者的时间序列从 500 帧依次减少到 72 帧后发现，要保证个体CAP 与组 CAP 模型较高的空间相关性，最小帧数约为 125帧（数据未显示）。
CAP 动力学计算:
为了研究 CAP 变化的形成与消散，文章基于CAP-全脑信号相关性构建了CAP时间序列。先计算CAP与fMRI全时程数据的空间相关性，将相关系数作为CAP时间序列。选择全脑信号与CAP区域最大值时间序列相关性最好的时间点作为CAP代表帧，限制峰值帧数为总帧数的2%，并根据CAP在全脑信号中的占比调整。每个选定的帧都设置为时间 t = 0 的参考事件，然后在 -30 < 30 的时间内观察CAP 组装和拆卸过程（视频 S1、S2、S3、S4、S5 和 S6）。
为了研究CAP 变化的时间特征。文章首先探究了CAP的频谱特征，发现CAPs 1~6 的频谱均不受簇数 k 的影响，并且与fMRI全脑信号特征类似，均在0.01 - 0.03 Hz处富集。为了评估慢速波动（0.01~0.03 Hz）的重要性，文章在时间序列中随机插入CAP1~6，然后计算rsfMRI 时间序列的功率谱。该策略保留了信号幅度分布以及原始信号的相关结构。接下来在 1000 次迭代后进行统计评估，使用单尾双样本 t 检验（FDR 校正）进行统计。进一步的分析还表明，在 k = 20 时获得的 CAP 的频谱特征与 k = 6 观察到的慢速活动峰值相似（数据未显示）。
为了计算CAP在全脑信号波动中的瞬时相位，文章首先对全脑信进行带通滤波（0.01-0.03 Hz），然后使用Hilbert变换将全脑信号分解为具有瞬时相位和幅度特征分析信号。将每个受试者的瞬时全脑信号划分为 [0, 2π] 范围内的周期，探究每个周期内每个 CAP 发生时对应的全脑信号相位。
为了探究 CAP 发生与全脑次慢动力学之间是否存在相位耦合，文章计算了全脑信号周期内给定 CAP 发生与前一个 CAP 发生之间的全脑信号角度差（相位差）。每个 CAP 的全脑相位样本仅在其在该时刻的滤波值在相应实例中高于 1倍SD 时才被考虑。为了排除共线性对与投影方法对相位分析的贡献，文章计算了每对 CAP 的空间相似性 （r） 与它们出现之间的 GS 相位差之间的圆形线性相关性。该分析未显示任何检查的 CAP 对有任何显著关联（p > 0.2，所有 CAP 对，FDR 校正），表明研究中，投影方法不会人为地将具有相似空间分布的集群在相似阶段。
自闭症遗传模型中的 CAP 动力学：
在标准预处理步骤后，将上述相同的分析步骤应用于记录在Chd8+/+ 和Chd8+/- 动物的 rsfMRI 时间序列上，使用 k = 6 在每组中独立计算 CAP。通过三个收敛标准的交集获得组间 CAP 匹配模式。首先分析了每组 CAP 之间的空间相关性，以独立识别三个 CAP 和抗 CAP 对。接下来计算了两个组中已识别的 CAP 对的空间相似性， CAP 对 1-2、3-4 和 5-6 之间的空间相关性在 0.45 到 0.72 之间。对获得的相关系数进行验证，任何其他 Chd8+/- CAP 重新排序都会导致基因型之间的相关匹配低得多。最后，通过观察到匹配的 CAP 的相位出现在不同基因型之间是一致的（图 6C），证实了匹配结果的正确性。
最后，为了评估单个 CAP 在之前描述的海马和运动感觉区域之间基于种子的 rsfMRI 相关差异中的参与，文章以双侧皮层区域为种子点观察了功能连接强度的差异。结果使用双样本 t 检验和FWE校正评估组间差异 （p < 0.05， T（40） = 2.8）。然后重新计算了该种子与海马中过度连接的选定病灶之间的 Pearson 相关性。rsfMRI 相关图为在每个CAP独立回归每个 CAP 的时间过程之前和之后计算。通过重复测量双向方差分析评估 rsfMRI 相关性的差异。
实验结果：
- 基于种子点的共激活模式（CAPs）与静息态功能连接模式（FC）高度相关
文章首先尝试重复文献结果，比较基于种子点的共激活模式与功能连接图之间的关系。文章分别以Acg、Right SS、Left AUD、Right dHCP为种子点，得到的共激活模式与功能连接图均有着较好的一致性。右图以不同比例高信号强度图像计算共激活模式与功能连接图的相关性，结果显示信号强度前15%的图像和功能连接图已经有着较高的相关性（>0.8）。
图S1 基于种子点的共激活模式与功能连接图之间的关系。
Acg 前扣带回，SS 感觉皮层，AUD 听觉皮层，HCP 背侧海马
2.静息态功能成像数据可以使用六种共激活模式来描述
之前的研究表明大脑的共激活模式会在有限个模式中循环出现。基于这种解释框架，作者尝试利用k-means聚类对全时程的静息态数据进行分类，将预设区域、结构特征和强度特征等因素排除在外。这种方式与基于时间点的共激活模式不同，可以同时兼顾高信号强度区域和低信号强度区域，将图像整体较为相似的时间点划分成同一类。
图1
A：基于种子点的静息态功能网络；
B：基于种子点的时间序列（ACg，SSs）和全脑时间序列（Global）以及对应的共激活模式（F1~F4）；
C：不同动物个体可以被划分成相同的几种共激活模式
文章尝试将静息态数据划分成2-40个类别（2<k<40）并探究增加分类数对模式解释的边际效应（图 S2 A）。结果显示分类数在4～8时，边际效应有着明显的变化。当分类数为6时，模型对变异的解释力可以达到59.6%。而当分类数为40时，模型对变异的解释力为63.6%，仅增加3.7%。同时当分类数>6时，增加分类数对变异的解释力提升较小（图S2 B）。文章进一步比较不同数据集（DS1，DS2，Chd8+/+）聚类结果的差异。结果显示，分类数为6时，三个模型整体相似性最好（图 S2 C）。同时作者还比较了分类数为20时各分类之间的相似性，结果显示多个分类具有较高的相关性（>0.5），说明部分分类是冗余的。故文章认为6分类是本数据集最佳的分类模式。
图2
A：分类数k（2-40）增加对变异系数的影响；
B：分类数k逐级增加对变异系数的影响；
C：三个数据集在分类数6～8时共激活模式的相似性
3.六种共激活模式的特征及其代表性
作者将静息态功能连接结果划分成6种共激活模式（CAP）（图）。其中CAP1表现出初级、次级感觉运动皮层的激活，与侧方皮层网络（LCN）类似。同时皮层-边缘区域和海马区域的去激活与默认网络（DMN）类似。而CAP5包含了DMN网络和海马网络。这6种CAP有明显的相关性，其中CAP1-2（r = -0.8）、CAP3-4（r = -0.96）、CAP5-6（r = -0.66）呈现明显的负相关。6种共激活模式占整体的比例差异较小，均在10%～30%之间波动。另外，6种共激活模式持续时间相近，约为5s。
图3
A：6种共激活模式；B：6种共激活模式之间的相关性；C：6种共激活模式的发生率；D：6种共激活模式的持续时间
以这6种共激活模式为标准，文章进一步测试了该分类对个体静息态数据的解析能力。文章对19只动物的静息态数据进行测试，每个数据1500帧数据（30min），根据组分析获得的6种共激活模式计算个体中每一个时间点CAP发生的概率。结果显示这6种CAP在个体水平上同样有着很高的发生概率，并且其模式与CAP模式十分吻合（图 3）。说明这6种共激活模式可以灵敏的反应静息态大脑真实的状态。
图4
A：6种共激活模式在个体水平的出现概率；
B：单个被试的CAP空间分布情况
为了探究信噪比波动对CAP计算的影响，文章提取了脑实质区域信号和空气信号进行比较。结果显示在6种CAP模式中信噪比一致（one-way ANOVA；p = 0.98），说明信噪比的波动对CAP聚类结果的影响较小（图 S5 A）。为了探究运动对CAP计算的影响，文章先对静息态数据进行了筛选去除高运动数据再进行聚类（FD >0.075mm或FD > 0.1mm），结果显示CAP结果与未筛选数据结果相近（相关系数 r > 0.98）（图 S5 B）。另外，高运动数据在6种分类中也呈现均匀分布（图 S5 C），说明运动对CAP聚类结果的影响较小。
图5
A：信噪比在6种CAP中的分布情况；
B：去除高移动数据后的CAP聚类情况；
C、D：高运动数据在6种CAP中的分布情况
4.六种共激活状态的次慢波动力学特征
在确定这6种共激活模式为基本特征之后，作者尝试进一步探究这些模式在整个扫描过程中的动力学特征。文章以这6种特征为基础，依次计算每一种特征和每一个时间点fMRI图像的相关性，构建了6个基于相关性的时间序列，命名为“CAP time courses”。对CAP时间序列进行频谱分析，结果显示6种模式的主要频率均在集中0.01~0.03Hz之间（图 4B），说明这6种CAP在大脑中有规律的循环出现。
图6
A：CAP时间序列和fMRI全脑信号时间序列；
B：CAP时间序列和fMRI全脑信号时间序列的频率分布图
有趣的是，fMRI全脑信号和这6种共激活模式有着类似的频率分布图（图 4B），那么CAP模式和全脑信号有着怎样的关系呢？作者猜测全脑信号的的波动并不仅仅代表大脑整体活动的上升或下降，而是不同特征权重变化的反应。为了验证这个假设，作者先将全脑信号进行滤波，只保留0.01～0.03Hz频率的信号，然后计算每一种CAP在全脑信号波动中发生概率。为了探究全脑信号波动周期的不同部分是否对应不同状态，作者计算了全脑信号对应的相位CAP发生的概率，其中0和π分别代表全脑信号的波峰和波谷。结果显示6种CAP在整个全脑信号波动中并不是均匀分布的，而是集中在全脑信号波动的特定位置（图 5A）。其中CAP 2和CAP 4在波峰区域较为集中，CAP 1和CAP 3在波谷区域较为集中，而CAP 5 和CAP 6在π/2~-π/2区域较为集中（图 5B）。
图7
A：6种CAP模式单个被试fMRI全脑信号中发生的概率；
B：6种CAP模式在一个波动周期中发生的概率；
C：CAP发生的概率及对应的全脑信号相位，其中红色向量代表加权发生率
5.自闭症小鼠模型的共激活模式特征和动力学特征
为了探究这6种共激活模式对脑功能变异的代表性，文章以Chd8基因缺失小鼠为研究对象，之前研究表明这种动物存在自闭症样行为，在海马和运动皮层区域存在明显的功能异常。首先，文章对Chd8+/+正常小鼠的静息态数据进行了分类（图 S2 C），结果表现出与正常动物类似的分类特征。而Chd8+/-单基因缺失小鼠的CAP模式与同批次对照组动物相比有着显著的差异（图6 A、B）。其中CAP1和CAP2的丘脑和背侧海马区域在Chd8+/-单基因缺失小鼠中表现出明显的功能异常。同时，CAP3和CAP4的扣带回和中部丘脑区在Chd8+/-单基因缺失小鼠也表现出明显异常。而CAP5 和CAP6 的感觉运动皮层以及前额叶也出现了共激活强度的变化。说明Chd8+/-单基因缺失小鼠存在多个功能模式的改变，这可能是fMRI自发信号变化的原因。
文章进一步探究了Chd8+/-单基因缺失小鼠CAP动力学的变化（图 6 C）。分别研究对照组和疾病组的6种CAP模式在全脑信号波动中的分布情况。结果显示6种CAP模式全脑信号波动中的分布情况类似，但是除CAP3外，其余所有Chd8+/-自闭症小鼠的CAP模式对应相位均慢于对照组，相位偏差约-30.9°~-46.5°（图6 C）。结果显示Chd8+/-自闭症小鼠的局部特征和动力学特征均有明显变化。
为了探究疾病共激活模式特征和功能连接变化之间的关系，文章以感觉皮层（ss）为种子点分析了Chd8+/-自闭症小鼠模型功能连接强度的变化模式，和文献报道结果一致，多个感觉皮层-海马区域（ss-HCP）的功能连接强度表现出明显上升。为了评估6中CAP模式中哪一种对于功能连接强度的影响较大，文章分别将CAP1~6以及全脑信号（GS）作为协变量进行回归。结果显示当使用CAP1、2作为协变量回归后，实验组和对照组ss-HCP功能连接强度之间的差异消失，而CAP3~6信号回归ss-HCP功能连接强度差异影响有限。说明功能连接强度变化可能和特定的共激活模式有关。同时全脑信号（GS）回归也导致ss-HCP差异减少（这也提示我们全脑信号回归可能导致有效信息被去除）。
图8
A：Chd8+/+正常小鼠和Chd8+/-自闭症模型小鼠的CAP聚类结果；
B：Chd8+/+正常小鼠和Chd8+/-自闭症模型小鼠CAP特征之间的差异；
C：6种特征在全脑信号中出现的概率以及对应的相位；
D：两组基于种子点的静息态功能连接差异；
E：6种CAP信号回归后ss-HCP之间的差异变化
讨论：
rsfMRI数据动态结果显示大脑活动并不是一种固定的状态，而是不断发生的重配置、循环和传播。本文首先通过聚类的方式将大脑划分成了六种状态（CAPs），然后探究了这六种状态的时空特征以及与全脑信号之间的关系，最后利用这用分析方法研究了自闭症小鼠模型的功能网络特征和动力学特征。
文章为fMRI的自发活动提供了解释新的解释框架：
其一，研究证明了不依赖静息态功能连接计算结果，仅通过定位时间序列中的特定帧便可以构建出有效的功能网路特征。
其二，研究证明了大脑活动由多个动态循环的功能网络构成，静息状态下多个功能网络以次慢波形式不断进行着组装和拆卸。
其三，研究基于rsfMRI数据，为自闭症小鼠模型提供了新的解释框架，研究通过6种共激活模式特征更加细化的描述了自闭症小鼠模型的脑网络变化，并且多个共激活模式的动力学特征也发生了明显变化。该解释框架显示的功能异常区海马-感觉皮层与自闭症患者感觉、认知功能障碍的特征一致。
之前研究已经基于脑区之间信号的同步性构建了多个功能网络模型，如默认网络（DMN）、任务响应网络（TPN）、边缘-皮质网络（LCN）等。而本研究通过逐帧聚类得到的6种共激活模式与这些功能网络也存在多个重叠区域。此外，文章还界定了对这些共激活模式在全脑信号的分布规律进行了探究，发现fMRI全脑信号是由多个功能网络有规律的组合叠加形成。尽管这些特征网络变换的驱动力尚不清楚，但本研究中在动物模型上发现的CAP 4基底前脑-皮层特征与人脑基底前脑-皮层特征具有很高的一致性。
神经元层面的研究结果也表明，自发的神经活动具有次慢波特性以及跨越整个背侧皮层的瞬时共激活模式。其中，已有多项研究报道皮质钙活动与空间结构血液动力学波动之间的联系，并且这些波动的空间结构和DMN和LCN网络特征具有较高的一致性。这些结果支撑了全脑 fMRI 波动可能由多个慢波功能网络叠加形成。
实验进一步可以考虑如何从方法学层面更好的比较不同数据集或实验条件下CAP的特征。本文在研究病理模型小鼠时使用的CAP分别从各自的数据集中聚类得到，而同样可以考虑使用组平均的CAP模式探究模型小鼠中不同CAP出现的概率。两者之间的差异则需要进一步研究比较得到。
总结：
本文根据共激活方法的理念，构建了一套新的fMRI信号解释体系，该体系将静息态数据划分成六个特征网络，这些特征网络具有特定的次慢波动力学特征，可以更精细的描述疾病模型的脑网络变化模式。
扫一扫
1338
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
