摘要：

大脑动态功能连接描述了脑区之间短暂的连接。大脑动态特征已与成人的行为情绪和认知功能联系起来，而异常的动态模式与神经发育障碍（如自闭症）相关。尽管在新生儿中已经反复确认了可靠的功能脑网络，但关于动态功能连接的早期发展却鲜有研究。本研究通过功能磁共振成像（fMRI）对足月新生儿（n = 324）和早产儿（n = 66）在出生后几周内的动态功能连接进行了特征描述。研究发现，在人类大脑的出生时已经存在一个动态连接的特征网络，表现为六种短暂的新生儿功能连接状态，并且其动态模式在新生儿期会不断变化。早产儿的动态连接模式较为异常，并与通过18个月时的《学步期儿童自闭症定量筛查表》（Q-CHAT）得分所测量的非典型社交、感知和重复性行为相关。

静息态功能磁共振成像（fMRI）常用于推断大脑不同区域如何协同工作，并构建全脑功能连接网络，这通常基于一个静态的假设。虽然这种方法有助于确定大脑活动的“平均”特征，但大脑的连接本质上是动态的，即非静态的。动态功能连接（dFC）旨在捕捉维持不同脑状态、适应短暂情境以及整合信息所需的持续神经调整。这些动态特性，例如在不同连接模式或“状态”之间的转换，已被与语言、认知和运动功能等过程联系起来。此外，异常的大脑动态特性还与神经发育障碍（如精神分裂症、注意缺陷多动障碍和自闭症谱系障碍）的临床特征及认知功能障碍相关。例如，患有自闭症谱系障碍（ASD）的个体据报道在不同连接模式之间的转换更为直接，而典型发育的个体则通过中间连接模式在这些相同的脑状态之间转换。尽管越来越多的研究表明，神经发育障碍可能起源于围产期，但目前对这一关键阶段的大脑动态特性知之甚少。此外，尚不清楚新生儿的动态功能连接特性是否与儿童期的神经发育相关，尤其是哪些特性提示了后期神经发育障碍的可能性更高。

在本研究中，我们应用了先进的技术，对发育人类连接组计划（dHCP）中的一组足月新生儿（n=324）和早产儿（n=66）在矫正孕龄时的fMRI数据进行了分析，这是目前最大的健康新生儿大脑的公开人群数据集。我们采用了两种方法来研究动态脑功能：第一，我们使用基于Kuramoto Order Parameter（KOP）的测度，即平均同步性和元稳定性，作为评估大脑同步性和灵活性的全局指标；第二，我们通过Leading Eigenvector Analysis（LEiDA）方法识别出与时间相关的“状态”，即短暂的脑同步模式。

我们首先评估了全局动态特征与足月新生儿的扫描时矫正孕龄（PMA）和出生后天数（PND）之间的关联。在足月新生儿中，未观察到扫描时矫正孕龄（PMA）与全局动态特征之间的关联。然而，扫描时的出生后天数（PND）与元稳定性显著相关。与足月新生儿相比，早产儿的平均同步性和元稳定性较低。在整个队列中，元稳定性与任何神经发育指标均无显著关联，但平均同步性与Q-CHAT评分显著相关，即新生儿期的较低平均同步性与18个月时更为异常的社交、感知和重复性行为相关。

我们定义了六种不同的脑状态，通过LEiDA方法从K-Means聚类中启发式获得。简而言之，我们将fMRI时间序列平均到90个皮层和皮下区域，这些区域由适应新生儿大脑的Anatomical Automated Labels（AAL）图谱定义。我们估计了每对区域之间的相位同步性，为每个fMRI时间点获得一个动态功能连接矩阵。我们计算了每个矩阵在每个时间点的第一特征向量，并使用K-Means聚类，K=6。我们识别的六种状态中有三种显示出大脑区域之间的广泛相位一致性，我们称之为全脑全局同步状态，即全局状态A、全局状态B和全局状态C。三种状态更为区域性：一种状态显示枕叶皮层的相位同步；一种状态代表主要在感觉运动皮层的高同步性；一种状态包括额叶皮层、角回和后扣带回皮层的高同步性（为简便起见，我们称之为前顶状态或FP状态）。

我们对足月新生儿（n=324）的六种脑状态的主要动态特征进行了比较，分析指标包括均值同步性、元稳定性、占有率和停留时间。采用三类 ANOVA 模型（Satterthwaite 方法）进行检验，结果显示不同状态在均值同步性（F(5, 1615)=13163；p<0.001）、均值元稳定性（F(5, 1615)=291；p<0.001）、均值占有率（F(5, 1938)=514；p<0.001）和均值停留时间（F(5, 1938)=734；p<0.001）上均存在显著差异。状态间效应表明，与另外五种状态相比，全局状态 A 的均值同步性、均值占有率和均值停留时间更高。此外，全局状态 B 和 C 的均值同步性介于全局状态 A 和其他状态之间。

图 3a 展示了规范的脑状态转换概率分布情况。大多数情况表现为连续多次出现同一种状态，即“停留序列”，六种状态的停留序列概率均超过 89%。若排除这些停留序列，12 种最频繁的转换展现了复杂的动态特征。例如，大脑在全局状态 A（具有高同步性）与其他状态之间转换时，需经过全局状态 B 和 C（具有中等同步性）；而大脑进入前顶状态（FP 状态）时，可通过全局状态 B、枕叶状态和感觉运动状态（区域性状态）实现。

表1 全局动态特征（同步性与元稳定性）与扫描时矫正孕龄（PMA）、出生后天数（PND）的关联，以及与早产关系的统计分析。

足月新生儿脑模块状态的扫描年龄相关特征

扫描时的矫正孕龄（PMA）升高与全局状态 C 的以下指标呈正相关：停留时间（t = 4.4；p < 0.001）、占有率（t = 5.3；p < 0.001）及平均同步性（t = 3.6；p < 0.001）；也与感觉运动状态的占有率（t = 2.8；p < 0.001）和平均同步性（t = 3.1；p = 0.002）呈正相关。相反，PMA 与全局状态 B 的占有率呈负相关（t = −2.8；p = 0.004）。出生后天数（PND）增加则与全局状态 B（t = −2.8；p = 0.004）、枕叶状态（t = −2.5；p = 0.013）及感觉运动状态（t = −2.8；p = 0.006）的停留时间缩短相关。经 FDR 多重比较校正后仍显著的脑状态特征与 PMA/PND 的关联汇总于图 3c、d，并在图 4g 中进行了对比。

足月新生儿脑状态转换与扫描年龄的关联

较高的 PMA 与下列转换概率升高相关：从全局状态 A 到 C（t = 3.2；p = 0.002）、从全局状态 B 到前顶状态（t = 3.1；p = 0.002），以及在全局状态 C 内停留（t = 3.2；p = 0.001）。同时，较高的 PMA 也与以下概率降低相关：从全局状态 C 到枕叶状态（t = −2.8；p = 0.006），以及在全局状态 B 内停留（t = −2.9；p = 0.005）。PMA 与状态转换概率之间的显著关联见图 3e。

早产儿的异常模块脑状态特征

与足月儿相比，早产儿在全局状态 A 的停留时间缩短（t = −4.6；p < 0.001，图 4a）；其占有率亦降低（t = −4.1；p < 0.001）；而全局状态 B（t = 3.4；p = 0.001）、枕叶状态（t = 2.2；p = 0.030）及前顶状态（t = 2.7；p = 0.009）的占有率则升高（图 4b）。早产还伴随以下状态的平均同步性下降：全局状态 A（t = −5.1；p < 0.001）、全局状态 B（t = −3.4；p = 0.001）、全局状态 C（t = −2.8；p = 0.005）、枕叶状态（t = −2.2；p = 0.031）及前顶状态（t = −2.6；p = 0.010）（图 4c）；此外，全局状态 A（t = −2.5；p = 0.014）和前顶状态（t = −4.3；p < 0.001）的元稳定性亦降低（图 4d）。

早产对大脑动态及脑状态转换概率的影响

早产与向枕叶连接模式转换的概率增加相关，具体表现为：从全局状态 A 到全局状态 C 的转换概率升高（t = 4.7；p < 0.001）；从全局状态 C 到枕叶状态的转换概率升高（t = 3.1；p = 0.002）；在全局状态 A 内停留（驻留）的概率降低（t = −4.4；p < 0.001）；详见图 4f。

当以出生时胎龄（连续变量）作为指标评估大脑动态特征变化时，获得了与早产-足月二分类似的结果模式，详见补充图 S2。图 4h 汇总了与早产相关的转换概率显著变化，以及这些变化与扫描时矫正孕龄（PMA）和出生后天数（PND）之间的关系。

图 2 | 新生儿短暂性脑状态

按全局同步性水平从左到右排列。a.使用 AAL 脑区划分得到的新生儿大脑六种脑状态的 LEiDA 向量。b.将 LEiDA 结果映射于脑表面（右侧视图）。c.将 LEiDA 结果映射于脑表面（左侧视图）。Glb. 全局；Occ. 枕叶；SM 感觉运动；FP 前顶叶。原始数据见 Source Data 文件。

图 3 | 足月新生儿（n = 324）的大脑动态

a.包含停留序列在内的全部状态转换。b.排除停留序列后，出现频率最高的 12 种主要状态转换。c.与扫描时矫正孕龄（PMA）显著相关的脑状态特征汇总。d.与扫描时出生后天数（PND）显著相关的脑状态特征汇总。e.与扫描时 PMA 显著相关的状态转换概率汇总。f.与扫描时 PND 显著相关的状态转换概率汇总。统计模型 GLM1（纳入 324 名足月儿）：y -β₀ +β₁PMA +β₂PND +β₃Sex +β₄运动异常帧（FD）。c、d 中数值为 t 统计量。图中所有显著关联均经双侧置换检验，并通过 FDR 多重比较校正（α = 5%）。Glb.全局；Occ. 枕叶；SM 感觉运动；FP 前顶叶；PMA 矫正孕龄；PND 出生后天数。

原始数据见 Source Data 文件。

图 4 | 早产儿（含足月与早产总样本 n = 390）的大脑动态

a 平均停留时间（DT）。b 平均占有率（FO）。c 平均同步性。d 元稳定性。e 与早产显著相关的脑状态特征汇总。f 早产与状态转换概率的关联。g 与早产显著相关的脑状态特征汇总，并与扫描时 PMA、PND 的显著关联进行对比。h与 PMA 升高、PND 增加及/或早产相关的脑状态转换概率汇总。GLM1（324 名足月儿）：y-β₀+β₁PMA+β₂PND+β₃Sex+β₄运动异常帧（FD）。GLM2 及(a)(b)(c)(d)（324 名足月儿 + 66 名早产儿）：y - β₀ + β₁早产 + β₂PMA + β₃Sex + β₄运动异常帧（FD）。p < 0.05，** p < 0.01，*** p < 0.001（双侧置换检验）。e 中数值为 t 统计量。箱线图展示 0、25、50、75、100 百分位数；异常值定义为超出1.5×IQR + 75 百分位数。所有显著关联均经 FDR 多重比较校正（α = 5%）。

Glb. 全局，Occ. 枕叶，SM 感觉运动，FP 前顶叶，PMA 矫正孕龄，PND 出生后天数，DT停留时间，FO 占有率。原始数据见 Source Data 文件。

图 5 | 围产期脑状态特征与18个月校正年龄神经发育结局的关联汇总

a、b分别显示围产期六种脑状态的平均同步性（a）和平均占有率（b）与18个月时 Bayley-III 认知、语言、运动合成得分及 Q-CHAT 得分之间的关联。统计模型 GLM3：Bayley-III 样本：257 名足月儿 + 48 名早产儿Q-CHAT 样本：254名足月儿 + 46 名早产儿

模型公式：y-β₀+β₁GA+β₂PMA+β₃Sex+β₄运动异常帧(FD)+β₅[18个月校正年龄]+β6[评估分项]+β7[多重剥夺指数]

图中数值为 t 统计量；所有显著关联均经双侧置换检验并通过 FDR 多重比较校正（α = 5%）。缩写：Glb. 全局，Occ. 枕叶，SM 感觉运动，FP 前顶叶；Cog. 认知，Lang.语言，Mot. 运动，Comp. 合成得分；Q-CHAT 学步期儿童自闭症定量筛查表。

新生儿脑态与早期儿童神经发育

最后，我们考察了出生后脑态动力学与18个月时通过 Bayley-III 和 Q-CHAT 量表评估的神经发育结局之间的关联。新生儿传感运动（sensorimotor）状态的平均同步化程度越高，与18个月时 Bayley 认知和运动评分表现越差相关（认知：t = −3.2, p = 0.002；运动：t = −3.4, p = 0.002；见图5a）。18个月时更高的 Q-CHAT 得分与前顶叶（frontoparietal）状态的平均同步化程度降低相关（t = −2.7, p = 0.008），并且与传感运动状态的分数占比（fractional occupancy）增加相关（t = 2.5, p = 0.014），以及与全局状态 A 的分数占比降低相关（t = −2.7, p = 0.010；见图5b）。

使用另一种分区方案得到的区域性脑动力学呈现出类似模式

使用替代的 M-CRIB 新生儿脑图谱得到的六种脑状态与 AAL 图谱得到的结果相兼容（补充图 S3），其中包括三个全局同步化状态和三个较为区域限制的状态。用 M-CRIB 图谱得到的状态1、2、3 与 AAL 图谱得到的全局同步化状态一致。M-CRIB 的状态4 包含了一些出现在枕叶（occipital）状态中的结构，状态5 与传感运动状态的结构相似，状态6 与 AAL 得到的前顶叶（FP）状态一致。

此外，对于在足月出生婴儿中使用 M-CRIB 图谱的分析（补充图 S4），校正后发育年龄（PMA）也与状态5（近似于 SM）的分数占比增加（t = 2.9；p = 0.002）和平均同步化程度增加（t = 3.0；p = 0.003）相关；出生后天数（PND）则与状态5 的停留时间（dwell times）减少相关（t = −2.5；p = 0.012）。从状态6（~FP）到状态2（~全局）的转变随 PND 增加也呈上升趋势（t = 3.3；p = 0.002）。早产与 M-CRIB 图谱的状态度量有类似的变化（补充图 S5），表现为状态4（~枕叶）（t = 2.3；p = 0.020）和状态6（~FP）（t = 3.0；p = 0.002）的分数占比增加；以及状态6（~FP）的平均同步化程度降低（t = −2.2；p = 0.027）。关于胎龄（GA）的类似分析见补充图 S6。

在使用 M-CRIB 图谱时，与Bayley-III 的认知（t = −4.0；p < 0.001）和运动（t = −3.4；p < 0.001）分量的发育结局的关联显示，状态5（~SM）的平均同步化程度与这两项呈负相关；更高的 Q-CHAT 得分与状态1 的分数占比降低相关（t = −2.6；p = 0.008）。有关摘要见补充图 S7。使用 M-CRIB 图谱进行分区的结果的详细描述见补充材料。

讨论

我们使用全局及模块化的fMRI信号分析工具研究了大样本足月与早产新生儿的动态功能连接特征。本研究刻画了新生儿大脑状态及其在早期人类发育过程中的动态基本特征。我们发现早产儿的大脑动态受扰，并且早期出生后期的大脑动态特征谱与18个月时的一系列早期儿童神经发育与行为结局相关联。

我们发现，在足月出生人群中，早期出生后（扫描时PMA 37–44周）全局动态特征相对稳定，尽管随着出生后在外生活天数的增加，观察到较低的亚稳性（metastability），这表明出生后在外的生命经验会在出生后降低大脑动态的灵活性，并促进更稳定的连接模式。这与我们观察到的另一个现象一致：在足月等效年龄扫描的早产儿（因此在扫描时有更多的出生后外在经历）其亚稳性低于足月出生儿童。然而，早产儿在足月等效年龄还表现出较低的平均同步性（mean synchronisation），这提示与早产相关的大脑动态改变具有独特模式，且该模式独立于出生后外在生活暴露的时长（即扫描时的PND）。先前研究已将较低的亚稳性与成熟个体在外伤性脑损伤后认知灵活性受损联系起来11。然而，本研究中观察到的出生后生活时长与亚稳性之间的负相关不太可能与认知灵活性直接相关，而更可能反映由出生后外在生命经验相关的初级感觉刺激驱动的网络动态精细化。亚稳性在早产儿的足月等效年龄时似乎进一步降低，这可能与一些理论一致，即围产期压力会非典型地加速大脑成熟，但可能带来潜在的长期负面影响34。

我们的研究扩展了此前的工作：先前研究观察到，非常早产出生的儿童在儿童后期表现出亚最佳的神经同步性和改变的全局动态连接模式35。我们的指标表明，早产与大脑动态的关联可能在生命更早阶段就开始存在。早产儿的这些差异可能根源于伴随早产的结构与功能网络框架的改变。例如，以往研究显示早产与脑结构改变36,37、全局功能架构改变19,38以及新生儿期的结构网络改变39相关，这些改变在学龄期仍然存在40并延续到青年期41,42。全局动态特征与后期行为有关：新生儿期较低的平均同步性与18个月时较高的Q-CHAT评分相关。尽管Q-CHAT评分高表示更多自闭症特质，但在本研究中，Q-CHAT衡量的是社会、感觉和重复行为在正常分布上的连续性31,43。因此，尽管全局度量可能有助于提示人类大脑发育基础动态轨迹，但不应过度解读其与Q-CHAT的关联——本研究并非一项自闭症谱系障碍（ASD）研究，且我们也未在神经发育诊断开始正式化的年龄对儿童进行评估。

全局动态的汇总指标本身很可能由更复杂的活动所支撑。因此，我们进一步探讨了模块化脑状态在出生后早期的出现与行为。我们在足月等效年龄的新生儿大脑中表征了六种瞬态状态。在显示广泛一致性的三种状态（全局状态A–C）中，第一种涵盖了几乎所有被分析的皮质区域；第二种表现出对所有感官区（听觉、感觉运动和视觉）的更高贡献；第三种表现出视觉与额顶叶皮质的更高贡献。其余三种状态则跨域范围更受限/区域化，分别由感觉运动、视觉和额顶叶区域有不同贡献。因此，这些状态中的大多数包含了已知比高级网络更早成熟的初级感觉网络16,19,44。这增强了在围产期动态功能连接研究中报告结果的可信度。例如，Ma等人在37名足月新生儿中描述了涵盖默认模式、背侧注意、听觉、感觉运动和视觉网络的四种大脑状态45。在我们使用包含324名足月新生儿的大型队列中，本文建立了六种大脑状态并描述了它们与扫描时年龄和出生后天数的关联。

我们观察到，在新生儿中，全脑同步状态A具有最高的平均同步性、以及最大的占时分数和停留时间；这表明新生儿大量时间处于全局相位同步状态，这与先前使用fMRI45,46和EEG47所描述的主导全脑同步模式相似。结合已有研究，我们的工作支持这样一个观点：大尺度活动在早期大脑发育中起着关键作用。我们的结果亦与这样一种观点一致：这种活动可能支持大尺度皮质网络的形成，并可能促进随后在出生后第一年内成熟的相关长程连接48,49。

我们观察到，随着扫描时校正胎龄（PMA）增加，感知运动皮层内的占用率和平均同步性呈正相关。这支持已有证据，即该系统在出生时在功能和结构上相较于其他系统更为成熟。这可能是因为胎儿末期三个月发生了显著的短程功能重组50,51，以支持对来自足、手和口的感觉信息所需的功能特异性52,53。我们还发现，扫描时PMA增加与进入包含默认模式网络（DMN）前部的前顶叶同步状态的概率增加有关54。这与其他证据一致，尽管该系统在出生时相对不成熟，但在出生后会发生显著变化，前顶叶区域日益被纳入该网络16,19。我们的研究表明，该系统在出生后随着年龄增长被越来越多地招募。

扫描时出生后天数（PND）也与从枕叶到全脑同步状态的转变概率增加显著相关。该发现与神经发育的关键阶段一致。具体而言，来自枕叶皮层状态的转变变化可能有助于视运动能力和早期婴儿在子宫外环境中的感觉整合的成熟55。最后，我们观察到随着出生后生命持续时间的增加，从感知运动向前顶叶结构的转变增加，或许反映了前顶叶系统功能成熟的出现，这与这些区域中大量中间神经元迁移56重合，而感知运动系统相对已较成熟51。

早产与非典型神经发育有关22，包括自闭症诊断率更高57–62。我们组先前的研究表明，早产儿其功能结构存在改变16,19,38。在此，我们扩展该研究，报告早产儿的动态功能连接亦被改变，表现为枕叶和额叶状态的分数占用率增加、从全局到枕叶状态的转变增加；以及高同步性全脑状态的驻留时间减少。仅有一项其他研究探讨了早产与非典型动态功能连接的关联45。他们报告了在一个传感运动与听觉皮层之间连接更强的状态中平均驻留时间显著更短，以及全局状态的平均驻留时间显著更长45。然而，直接将其与我们的发现比较具有挑战性，一方面因为我们的研究在ROI数量上分辨率更高，另一方面我们的样本量也更大，统计功效更强。我们在新生儿中观察到多个全局状态，这些状态由枕叶和前顶叶区域以不同方式贡献，并且发现早产与状态A驻留时间呈负相关。总之，早产儿在动态功能连接上表现为向枕叶和前顶叶同步特征的转变，并抑制全脑同步模式。早产已知影响终生认知63,64。我们的结果提示，出生后不久即存在的脑动态改变可能具有功能后果。

我们的工作与其他研究一致，均确认新生儿期是感知运动皮层发育的关键时期50,52,53，随后向更高阶网络的连接转变16,19。我们在此将该观察扩展，显示改变的脑动态与18个月时的一般发育结局以及更特定的社交、感觉和重复行为相关联。围绕足月期，感知运动状态的同步性较低与18个月时更好的认知和运动结局（Bayley-III）呈正相关。该关联在使用AAL和M-CRIB两个图谱的分析中均可见。然而，当在出生附近高全脑同步状态A的分数占用率较低且感知运动状态的分数占用率增加时，18个月时Q-CHAT量表反映出更多非典型的社交、感觉和重复特征（在AAL图谱下）。因此，我们的研究表明，大脑动态与自闭症相关特质的联系不仅限于成年期的状态转变15,65，也可能包含在早期发育阶段确立的状态占用、整体同步性及转变的改变。一种可能解释感知运动皮层动态与后期社交、感觉和重复特质关联的机制是，新生儿期的脑动态改变可能使个体更容易表现出对感官刺激的异常反应66。感知运动状态较高的分数占用率与较高Q-CHAT评分的正相关，也可能代表在新生儿期对该网络的过度依赖——这可能影响更高阶网络的发展。

我们强调并未对Q-CHAT分数较高的儿童进行18个月之后的随访，因此本研究并未评估已确诊自闭症的预测指标，且研究对象并非已诊断的自闭症人群。然而，最终被诊断为自闭症的儿童以及具有更广泛表型的个体，可能在更早阶段就已出现相关特质。这也是开发Q-CHAT的出发点31,43。此外，越来越多的不同模态证据表明，未来被确诊为自闭症的婴儿从6–9个月大起就在神经生物学和生理学上显示差异67–72。因此我们的工作将与自闭症相关的早期社会认知特征，与出生时的动态功能连接，特别是感觉网络的动态，联系起来。这一在捕捉早期与自闭症相关特征（尽管未必用于诊断）的量表分数与感官系统动态之间的相关性，符合在自闭症个体终生中感觉过程的重要性。纵向研究显然是必要的，但我们的工作为感觉系统沿着神经多样性轨迹与自闭症特征跨越生命期的联系增加了证据。感觉差异常是预示自闭症诊断的早期特征70,73，并且持续存在，成为诊断的核心特征74。例如，我们先前报道具有更高自闭症家族风险的新生儿在感知运动皮层呈现更高的区域同质性75，这与我们记录到该区域更高的分数占用率相吻合。已诊断个体的研究报告了运动皮层异常激活76，这可能影响将视觉输入转化为运动理解的能力，进而对社交互动有潜在影响。这些感觉系统的功能性脑差异被认为源于兴奋-抑制通路特别是GABA神经传递的改变。相关证据包括感觉处理与感觉皮层GABA水平差异之间的紧密关系77，我们也曾在成人中提供直接实验证据，显示自闭症中视觉处理差异依赖于GABA78。

本研究的一些优点和局限性应予以说明。我们采用了先进的新生儿fMRI采集和专用的新生儿多频带EPI流程，具有高时间分辨率（TR = 392 ms），但也承认这在深部灰质结构（包括丘脑、基底节和脑干）中信噪比较低79。我们选择了LEiDA处理流程，因为它提供了大脑状态的时变指标，并已成功用于成人脑动态研究28,29。大多数先前关于动态脑功能连接的研究采用滑动窗口方法，这需要在分析中选择任意参数，例如窗口和步长大小。对这些参数尚无共识，因为使用小窗口可能放大假性变化，而大窗口可能淡化脑动态中的尖锐变化80–82。因此，选择像LEiDA这样的时变方法是我们方法学的一大优点，尽管像隐马尔可夫模型（HMM）等其他技术也可用83。我们的分析将脑状态发生视为相互独立，即不带记忆。未来研究可以受益于开发考虑先前出现状态如何影响当前状态的度量。我们的结果也可能受限于使用AAL图谱定义感兴趣区（ROI）。尽管该图谱已为新生儿使用做过调整32，但基于结构的图谱可能与功能边界在脑内（尤其在新生儿期）相关性较差2,84。未来工作的一个可能解决方案是为新生儿fMRI研究采用多模态生成的表面图谱85。然而，我们也使用了另一种分区方案（M-CRIB图谱）重复了我们的流程，获得了与新生儿AAL图谱相似的结果。dHCP队列侧重于表征典型发育，因此本研究中足月与早产婴儿的数量不平衡。纳入的早产婴儿大多为中度或晚期早产，且大多“健康”，未发现临床上重要的影像学异常。尽管这可能更能代表一般早产人群（80%的早产婴儿为中度或晚期早产86），但我们不能将结果外推至极早或极重度早产儿，或有显著白质损伤的婴儿。在我们的分析中，除将多重剥夺指数（Index of Multiple Deprivation）作为协变量外，也未考虑社会人口学因素和早期神经发育结局中的社会剥夺影响。然而，多项研究表明，家庭心理社会和社会人口学因素对儿童期脑发育有显著影响；例如母亲压力、抑郁、低教育水平、母亲移民身份、母亲年龄大于35岁、父亲年龄超过38岁以及低家庭收入等因素均与较差的发展结局相关87–91。因此，未来研究可受益于评估这些更广泛社会人口学指标与早期脑动态之间的联系92。

总之，在本研究中我们评估了新生儿的大规模功能脑动力学和瞬态脑状态。我们的方法使我们能够定义一组由早期脑发育阶段已被证实建立的结构组成的六个基本瞬态脑状态。我们强调了脑成熟对脑动力学的影响，以及与早产相关的非典型模式。出生时的大脑状态动力学似乎具有功能相关性，因为它们与早期儿童期的一系列神经发育结局相关。这鼓励进一步研究以了解这些指标的预测价值及其调节机制，从而在适当情况下指导支持和干预。

方法

参与者

我们分析了来自 dHCP（发布版本 3）23 的共 809 份数据中的 390 份数据，样本包括足月出生（n = 324）和早产（出生时胎龄 (GA) <37 周；n = 66）婴儿。完整的纳入和排除标准以及每一步被排除的样本详见补充图 S8。足月出生婴儿出生时中位胎龄为 40.14 周（四分位距 IQR = 1.86 周），并在出生后不久进行扫描（扫描时的中位围产期年龄（PMA）= 41.57 周，IQR = 2.43 周）。早产儿出生时的中位胎龄为 33.36 周（IQR = 5.86 周），在等同足月的时间点接受扫描（扫描时中位 PMA = 40.5 周，IQR = 2.71）。表 2 显示了样本的人口统计学资料。纳入本研究个体的扫描时 PMA 与出生时胎龄的分布见补充图 S9。所有儿童均被邀请到位于伦敦圣托马斯医院的发育大脑中心（Centre for the Developing Brain）在预期出生日期后 18 个月由有经验的小儿科医生或心理学家进行神经发育评估。Bayley 婴幼儿发展量表第三版（Bayley-III）30 用于评估总样本中 305 名个体（其中 257 名足月出生、48 名早产）在运动、语言和认知领域的一般发育结果（分数越高表示能力越强）。18 个月校正年龄时的量化儿童自闭症清单（Q-CHAT）31 在 300 名个体（254 名足月、46 名早产）中可用，作为测量在人群中呈连续分布的非典型社会、感觉和重复行为的指标31。尽管较高的 Q-CHAT 得分可能指示更多的临界或亚临界自闭症特征，但我们强调使用该量表的目的是捕捉 Bayley-III 无法覆盖的行为特征，而非作为 ASD 的筛查工具。我们还从每位受试者母亲出生时的地址获取了多重剥夺指数（IMD）排名——一种基于地理位置的综合剥夺程度测量93——并在评估新生儿脑动力学与随后神经发育及行为结果关系的模型中将其作为协变量。

MRI 采集

我们评估了 dHCP 在Evelina 新生儿成像中心（Evelina London Children’s Hospital）获取的 fMRI 扫描，使用 3 特斯拉Philips Achieva 系统（Philips Medical Systems）。研究获得英国国家研究伦理委员会（UK National Research Ethics Authority，批准号14/LO/1169）的伦理批准，并在数据采集前获得所有参与家庭的书面知情同意。扫描在无镇静的情况下进行，采用专用的新生儿设定及优化的转运系统、定位装置、听力保护以及定制的 32 通道接收头线圈和声学罩94。扫描期间由新生儿护士或儿科医生监督，持续监测心率、血氧饱和度和体温。血氧依赖（BOLD）功能磁共振成像采用多切片回波平面成像序列并结合多波束激发（multiband 因子 9）采集（重复时间 TR = 392 ms，回波时间 TE = 38 ms，体素大小 = 2.15 × 2.15 × 2.15 mm，翻转角 = 34°，45 层，总时长 = 15分 3 秒，体积数 = 2300）23。用于脑形态测量和临床报告的结构像也一并采集23。T1 加权图像的重建空间分辨率为 0.8 × 0.8 × 0.8 mm，视野 = 145 × 122 × 100 mm，TR = 4795 ms。T2 加权图像的重建空间分辨率为 0.8 × 0.8 × 0.8 mm，视野 = 145 × 145 × 108 mm，TR = 12 s，TE = 156 ms。具有过度运动（超过 10% 的运动异常帧19）或临床显著的偶发 MRI 发现（放射学评分为 4 或 5，表示白质、皮质、基底节或小脑存在重大病变——如 dHCP 数据库所述23）的 fMRI 数据集被排除。在双胞胎/三胞胎扫描的情况下，仅纳入运动异常帧最少的一个婴儿。

数据处理

单个 fMRI 数据集按照 dHCP 专用的新生儿处理流程79 进行预处理。简要来说，使用 topup95 校正由场非均匀引起的局部畸变；进行体内和体间运动校正；并使用刚性体重对齐和切片到体积的涡流校正96 校正相关的动态畸变。使用 FSL FIX97 回归剔除残余运动、多波束采集及心肺生理伪影。

头部运动已被证明会在静息态 fMRI 中产生伪影和系统性相关98。除了 dHCP 管道中为最小化运动伪影对 BOLD 信号影响而实施的特定处理步骤79，我们还使用帧间位移（FD）98 评估每个数据集的运动幅度。为最小化运动（及任何可能相关差异）对脑状态判定的影响，仅选择运动异常帧少于 10% 的参与者（运动异常定义为 FD > 第 75 百分位 + 1.5*IQR）进入最终分析子样本。运动异常帧总数也被作为控制协变量纳入后续统计模型。

我们使用专用的新生儿组织分割流程99 将 T2 加权体积分割为包括白质、灰质和脑脊液在内的九类组织。我们使用针对新生儿大脑映射的自动解剖标注（AAL）图谱100 将每个受试者的 T2加权体积分划为 90 个皮质及皮下结构分区，该图谱已在dHCP 高分辨率模板上进行适配并手动校正101。我们基于T2 加权对比和预先获得的分割，采用基于微分同胚的对称图像配准方法（SyN）102 将 AAL 图谱从模板空间非线性配准到每一受试者的原始空间。灰质分割和分区随后使用 dHCP 功能处理流程中提供的基于边界的线性配准方法传播到每个受试者的 fMRI 空间。然后在每个 AAL 分区与相应灰质、深灰质或基底节分割掩模的交集区域内计算平均 BOLD 时间序列。我们还采用了另一种分区方案，遵循相同流程但使用来自 M-CRIB 图谱103 的 80 个皮质区域，M-CRIB 是 Desikan–Killiany 图谱104 的新生儿适配版。

表2 | 足月出生组与早产组的人口学详情

我们对 BOLD 时序进行了带通滤波，使用二阶巴特沃斯滤波器，频段为 0.02–0.10 Hz29，并通过希尔伯特变换对每个分区在每个时间点获得相位φj(t)。对于给定的实值信号 s(t)，我们构造了复值信号z(t)33,105，如下所示：

其中 H[s(t)] 表示作用于实值信号s(t)的希尔伯特变换，其定义如下，p.v. 为柯西主值105：

每个分区的相位 φj(t) 可以直接从z(t) 计算得出：

Kuramoto 有序参数

Kuramoto 有序参数（KOP）用于衡量多个振荡器的全局同步水平，其定义见下式，其中φj(t)表示振荡器j在给定时间的信号相位。在我们的研究中，每个脑区（AAL区域）被视为一个独立的振荡器。

一旦在每个时间点得到了KOP，我们计算了平均KOP（随时间的同步性聚合，称为“平均同步化”）和KOP的标准差（称为“亚稳定性”）25,26。

平均KOP提供了对全脑同步水平的总体度量，而亚稳定性则衡量不同振荡器之间的同步随时间的波动程度，即大脑的灵活性11,25,27。

主特征向量分析（LEiDA）

KOP分析能够揭示所有振荡器（脑区）层面的全局动态特性，但无法指出哪些具体脑结构可能参与了这些变化。为评估此类模块化（局部）特性，我们应用了LEiDA方法——该方法使我们能够研究不同脑区集合间的相位相干性。为此，我们首先在每个时刻（TR）计算脑区i与脑区j之间的相位差，使用余弦距离：

对于每个fMRI体积，我们首先构建相位差矩阵。随后我们使用LEiDA 方法28,29 得到一个低维表示，其中LEiDA向量对应于矩阵Δφij(t)分解的第一特征向量。先前已有研究表明该方法能够揭示网络和图的社区结构信息106。一旦得到所有LEiDA向量，我们使用 K-Means对它们进行聚28,29，最优簇数K（为6）通过Calinski–Harabasz 和 Davies–Bouldin 方法启发式确定（补充图S10）。

每个簇代表一组LEiDA 向量，我们将每一簇称为一种脑态。可以用三类主要指标研究这些脑态的动态特性：分数占有率（fractional occupancy）——表示在某一脑态上花费的总时间比例或该脑态的出现概率；驻留时间（dwell time）——表示在该脑态上平均连续停留的时间；以及脑态之间转移的马尔可夫转移概率28,29。此外，我们还通过对属于某一簇的体积对应的体积计算值取平均，得到每一脑态的平均同步化和亚稳定性值。

统计分析

首先，我们将样本限制为足月出生个体（n = 324），评估脑成熟度与出生后经历（分别以扫描时的PMA 与PND 表示）与全局和模块化脑动态的关联。其次，我们在整个 390 人样本中评估早产出生与脑动态的关联。早产以二元变量编码：若出生时胎龄（GA）小于 37 周则为1（早产），若出生时胎龄≥ 37 周则为0（足月）。最后，我们评估全局与模块化动态特征与18个月校正龄的后续神经发育结局的关联（n = 305 — Bayley-III；n = 300 — Q-CHAT）。

全局动态

我们通过拟合线性模型GLM1（324 名足月婴儿）来表征年龄与出生后经历（扫描时的 PMA与PND）之间的关联：y-β0 + β1PMA + β2PND + β3Sex + β4[运动异常值(FD)]。为表征早产对脑动态的影响，我们拟合了线性模型 GLM2（324 名足月与 66 名早产婴儿）：y - β0 + β1Preterm-born + β2PMA + β3Sex + β4[运动异常值 (FD)]。我们还在模型 GLM3 中评估了脑全局动态与 18 个月的认知与行为结局（即 Bayley-III 与 Q-CHAT）之间的关联（用于 Bayley-III 分数的样本为 257 名足月与 48 名早产婴儿；用于 Q-CHAT 的样本为 254 名足月与 46 名早产婴儿）：y - β0 + β1GA + β2PMA + β3Sex + β4[运动异常值 (FD)] + β5[评估时的校正龄] + β6[被评估的量表分项] + β7[多重剥夺指数]（其中“被评估的量表分项”包括 Bayley 的认知、语言、运动或 Q-CHAT 总分）。

模块化动态：脑态

首先，我们通过类型 III 方差分析（使用 Satterthwaite 方法107）和线性混合效应模型 GLME1（包含 324 名足月）检验本研究定义的脑态在平均同步化、亚稳定性、分数占有率与个体驻留时间方面的差异：y- β0 + β1State + (1 | Subject ID)，其中Subject ID 作为随机效应。随后，通过拟合 GLM1（324名足月婴儿）：y - β0 + β1PMA + β2PND + β3Sex + β4[运动异常值 (FD)]，我们表征了扫描时的年龄（PMA）与出生后经历（PND）对六个脑态的分数占有率、驻留时间、平均同步化与亚稳定性的影响。

其次，为表征早产对脑态及脑态间转移概率的影响，我们拟合了 GLM2（324 名足月与 66 名早产婴儿）：y-β0 + β1Preterm-born + β2PMA + β3Sex + β4[运动异常值 (FD)]，以及 GLM3（324 名足月与 66 名早产婴儿）：y- β0 + β1GA + β2PMA + β3Sex + β4[运动异常值(FD)]。

第三，我们在模型 GLM4 中评估了脑动态与神经发育结局（即 18 个月的 Bayley-III 与 Q-CHAT）之间的关联（用于 Bayley-III 的样本为 257 名足月与 48 名早产婴儿；用于 Q-CHAT 的样本为 254 名足月与 46 名早产婴儿）：y-β0 + β1GA + β2PMA + β3Sex + β4[运动异常值 (FD)] + β5[评估时的校正龄] + β6[被评估的量表分项] + β7[多重剥夺指数]（其中“被评估的量表分项”包括 Bayley 的认知、语言、运动或 Q-CHAT 总分）。

统计显著性与重复测量

我们使用双侧置换检验（permutation tests），对所有广义线性模型（GLMs）进行了10,000次重复，以评估每个感兴趣变量的统计显著性。P值未做校正地报告，并强调那些在状态间多重比较校正后仍显著的结果，校正方法为Benjamini–Hochberg FDR，显著性水平（α）取5%108。

报告摘要

有关研究设计的更多信息可在与本文链接的 Nature Portfolio Reporting Summary 中查阅。

数据可用性

已将适配至 dHCP 模板空间的AAL-UNC 大脑图谱以及本研究生成的预处理 BOLD 时间序列数据存放于 Zenodo 数据库，访问号为 10.5281/zenodo.7053984。fMRI 数据集和临床数据因 dHCP 数据发布条款而受限访问，访问可通过 https://data.developingconnectome.org 申请获得。随本文提供了源数据。

代码可用性

BOLD 信号波动的动态特性使用 dynFC：CoDe-Neuro 的动态功能连接工具（Dynamic Functional Connectivity Tools）109 评估，该工具由 Python v3.7 编写的脚本组成（https://code-neuro.github.io/dynfc/），并依赖以下支持库：Numpy、SciPy、Scikit-learn、pickle、h5py、pandas、os、sys 和 feather。

统计分析和图表在 R 语言及辅助包中完成，所用包包括 ggplot2、tidyr、dplyr、cowplot、purrr、RColorBrewer、knitr、janitor、ggExtra、stringr、rjson、broom、Tidymodels、coin、shadowtex、effsize、modelr、ggimage、ggpubr、patchwork、ggbeeswarm、ggrepel、ggtext、MetBrewer、lmerTest、forcats、stateR、p-testR、broom.-mixed、lme4 和 DiagrammeR。

参考文献

解读：万家利

审核：薛晓萌