用生物物理启发的神经ODE推断单细胞转录组随机动态.
0. TL; DR
单细胞RNA测序技术（scRNA-seq）极大地推动了我们对细胞异质性和命运决定的理解，但从静态的“快照”数据中推断随机动态过程仍然是一个根本性挑战。现有方法面临一个关键的权衡：机理模型通常假设过于严格，限制了生物学真实性；而纯数据驱动的方法则牺牲了可解释性，难以进行深入的机理探索。
本文介绍了一种名为DynNet的深度学习方法，它将神经ODE（Neural ODEs）与生物物理模型以及基因表达动态的先验知识相结合，旨在学习基因调控系统在细胞命运决定过程中的随机动态。通过在合成数据集上的基准测试，DynNet展示了其推断稳定细胞状态、重构动态轨迹以及表征多稳态细胞命运转变的能力。在一个肝细胞分化数据集上，DynNet成功推断了发育轨迹和潜在的细胞命运景观，揭示了不同细胞状态的稳定性及转变概率。最后，在应用于上皮-间质转化（Epithelial-mesenchymal transition, EMT）数据集时，DynNet进一步捕捉了EMT过程中的关键基因调控和状态转换路径。
1. 背景介绍
随机动态转换是生物、物理和化学系统中的基本现象，对于理解细胞命运决定尤为关键。多稳态（multi-stability）的存在允许细胞在响应内外信号时，在不同的稳定命运之间进行切换，这在干细胞分化、免疫细胞激活和癌细胞表型转换等过程中起着核心作用。
scRNA-seq技术使我们能够在单细胞分辨率上研究细胞功能，但其破坏性的测量方式意味着我们无法直接追踪单个细胞随时间变化的分子轨迹。这一局限性催生了大量的计算方法，旨在从离散的时间点快照中重构连续的动态过程。
现有方法大致可分为两类：
- 马尔可夫模型（Markov-model-based）方法：这类方法，如CellRank，专注于重构细胞状态轨迹和排列伪时间，但它们缺乏一个与底层基因调控机制明确关联的、可定量的、可解释的动态模型。
- 深度学习方法：神经ODE（Neural ODEs）利用神经网络和微分方程，能够从数据中连续地重构基因表达动态。然而，经典的确定性Neural ODEs无法捕捉多稳态动态，因为在任何一个点上，动态方向都是唯一的。此外，这类纯数据驱动的框架难以融入先验知识（如已知的基因调控关系），并且通常在低维隐空间中推断动态，限制了在单个基因层面的直接可解释性。
为了弥合机理模型与数据驱动方法之间的鸿沟，作者提出了DynNet（Stochastic Dynamics Inference Network）。DynNet的核心思想是利用随机微分方程（Stochastic Differential Equations, SDEs）来解耦细胞动态，将其分解为一个描述基因调控的确定性漂移项（drift）和一个代表随机波动的扩散项（diffusion）。通过将生物物理模型（如Hill函数）嵌入Neural ODE框架，并利用先验的基因调控网络（GRN）作为约束，DynNet旨在从稀疏的时间序列scRNA-seq数据中，学习一个既符合数据又具有生物学可解释性的连续随机动态模型。
2. DynNet 方法
DynNet是一个混合建模框架，它基于动力系统理论，旨在从多个离散时间点的scRNA-seq快照数据中重构连续的细胞命运转换动态。
2.1 核心方程：随机微分方程 (SDE)
DynNet使用随机微分方程（SDEs）来描述基因表达的动态过程：
\[d\mathbf{x} = \mathbf{f}(\mathbf{x})dt + \mathbf{\sigma}(\mathbf{x})d\mathbf{W}\]其中$\mathbf{x} \in \mathbb{R}^M$ 是一个细胞在某个时间的$M$维基因表达状态。$\mathbf{f}(\mathbf{x}) \in \mathbb{R}^M$ 是drift项，描述了由非线性基因调控相互作用驱动的确定性动态，即基因表达的变化趋势。$\mathbf{\sigma}(\mathbf{x})d\mathbf{W}$ 是diffusion项，描述了系统中的随机波动，$\mathbf{W}(t)$是标准布朗运动。
2.2 DynNet框架的构成
DynNet由两个核心的神经网络构成：
速度估计网络 (Velocity Estimation Network)
速度估计网络用于估计漂移项$\mathbf{f}(\mathbf{x})$。为了融入生物学机理，该网络的结构被设计为显式地编码基因调控关系。对于第$i$个基因，其动态由以下公式描述：
\[\mathbf{f}_i(\mathbf{x}) = \sum_{N_{ij}=1} \frac{a_{ij}x_j^{n_h}}{s_j^{n_h} + x_j^{n_h}} + \sum_{N_{ij}=-1} \frac{b_{ij}s_j^{n_h}}{s_j^{n_h} + x_j^{n_h}} - k_i x_i\]该公式基于经典的Hill函数，分别对来自其他基因$j$的激活和抑制作用进行建模。$a_{ij}$和$b_{ij}$分别代表基因$j$对基因$i$的激活和抑制强度。$s_j$是基因$j$的半激活/抑制阈值。$n_h$是Hill系数，控制调控的陡峭程度（本文中固定为4）。$k_i$是基因$i$的降解速率。$N$是一个先验的GRN邻接矩阵，用于指导模型的学习。
噪声估计网络 (Noise Estimation Network)
这是一个Multi-Layer Perceptron (MLP)，用于估计扩散矩阵$\mathbf{D}(\mathbf{x})$（其中$2\mathbf{D}(\mathbf{x}) = \mathbf{\sigma}(\mathbf{x})\mathbf{\sigma}(\mathbf{x})^T$）。为了简化计算，该网络估计一个对角扩散矩阵，主要用于建模内在的生化波动。
2.3 约束与优化策略
为了从稀疏的时间点数据中学到可靠的动态，DynNet引入了两种约束机制，并采用交替优化的训练策略。
约束机制
- 硬约束 (Hard constraints)：直接“屏蔽”掉先验GRN中不存在的调控连接，只优化存在的连接参数。
- 软约束 (Soft constraints)：通过正则化项惩罚与先验GRN不符的连接，使其参数趋向于零，但保留了发现新连接的可能性。
交替优化 (Alternating Optimization)
优化速度估计器：
首先固定噪声，优化速度估计器的参数$\theta = (\mathbf{A}, \mathbf{B}, \mathbf{S})$（分别对应激活矩阵、抑制矩阵和阈值向量）。这一步的损失函数$\mathcal{L}_v$是一个复合损失：
\[\mathcal{L}_v = \mu_1 \mathbb{E}_T[|\mathbf{f}(\mathbf{x}_{t_n})|^2] + \mu_2 \frac{1}{T} \sum_{t=1}^T \mathbb{E}_t[|\hat{\mathbf{w}}_t - \mathbf{w}_t|^2] + \lambda_1 \mathcal{C}(\theta)\]第一项动态一致性损失确保在最后一个时间点（假设为稳态），系统的漂移项趋近于零。第二项权重匹配损失通过高斯混合模型（GMM）估计真实数据和模型模拟数据在各个稳态上的细胞分布权重，并最小化二者的差异，以确保模型能正确捕捉多稳态的全局分布。第三项是基于先验GRN的软约束或硬约束。
优化噪声估计器：
固定优化好的速度估计器，转而优化噪声估计器的参数$\vartheta$。这一步的损失函数$\mathcal{L}_{\mathcal{D}}$是最小化模型模拟的边际基因表达分布与真实数据分布之间的Kullback-Leibler (KL)散度：
\[\mathcal{L}_{\mathcal{D}} = \frac{1}{T} \sum_{t=1}^T \mathbb{E}_t[\text{KL}(\hat{p}_{j,t} | p_{j,t})] + \lambda_2 \mathcal{R}_2(\vartheta)\]2.4 下游分析：能量景观与动态路径
训练完成后，DynNet可以提供丰富的下游分析：
- 能量景观构建：通过求解福克-普朗克方程（Fokker-Planck equation）的稳态解，可以得到系统的稳态概率分布$P_{ss}(\mathbf{x})$。能量景观$U(\mathbf{x})$则定义为$U(\mathbf{x}) = -\ln P_{ss}(\mathbf{x})$。景观中的“山谷”（局部最小值）对应稳定的细胞表型（吸引子），而“山脊”和“鞍点”则揭示了状态转变的能垒和关键路径。
- 状态转换分析：利用景观和动力学模型，可以计算最可能的转换路径（Minimum Action Path）和平均首达时间（First Passage Time, FPT），定量地分析不同细胞命运间的转换概率和速度。
- In Silico扰动：通过在模型中修改特定基因的表达或调控参数，可以模拟基因敲除或过表达的效应，预测其对细胞命运的影响。
3. 实验分析
作者通过一系列合成数据集和真实生物学数据集，系统地评估了DynNet的性能。
实验一：MISA模型的双稳态动态恢复
作者首先在一个经典的双基因互抑制自激活（Mutual Inhibition and Self-Activation, MISA）模型上测试DynNet。该模型能产生三个稳态。
- 图A, B：展示了MISA模型的网络结构和其生成的三稳态数据（S1, S2, S3）。
- 图C, D：DynNet准确地识别了三个稳态，其预测的稳态中心位置、协方差区域以及基因X1和X2的边际分布均与基准数据高度吻合。
- 图E：重构的能量景观清晰地展示了三个吸引子（“山谷”），与稳态位置一致，速度场（箭头）也正确地指向了吸引子。
- 图F：DynNet预测的决策边界（Decision Boundaries）与基准数据推断的边界高度一致，准确地划分了三个吸引盆，证明其能可靠地预测未定命运细胞的最终归宿。
- 图G, H：在速度场重构的准确性上，无论训练数据的时间点是2个、3个还是5个，DynNet的均方误差（MSE）始终低于经典的Neural ODE方法，并且能更准确地捕捉到鞍点等关键动态特征，显示了其在稀疏时间数据下的鲁棒性。
实验二：发育分叉过程的动态恢复
作者构建了一个包含7个基因的调控网络，模拟了从一个干细胞样状态（I）分叉到两个不同分化状态（D1, D2）的过程。
- 图A-C：展示了网络结构、PCA可视化的分叉轨迹以及两个终末状态下差异表达的基因。
- 图D, E：DynNet成功地从模拟数据中恢复了正确的稳态边际分布和指向两个终末状态的速度场。
- 图F：模型准确地捕捉到了下游标志基因G6和G7在分叉点后的表达分离，重现了细胞命运决定的动态过程。
- 图G, H：比较了软约束和硬约束的效果。
- 图H显示，硬约束收敛更快。
- 图G显示，软约束虽然收敛较慢，但能够发现一些由间接调控（如G5对G7）引起的“伪”正相关连接，这既是挑战也是其探索未知调控的潜力所在。对于已知的连接，其权重预测是准确的。这表明软约束更适合作为捕捉功能依赖性的工具，而非精确的结构推断工具。
实验三：肝细胞分化的分叉动态解码
作者将DynNet应用于一个真实的肝细胞分化数据集，该过程涉及从iPSCs分化到成熟肝细胞样（MH）或肝芽（LB）两种不同命运。
- 图A-C：展示了分化路径示意图、基于20个核心基因构建的先验GRN，以及PCA展示的分叉轨迹。
- 图D, E：DynNet的预测结果与真实数据高度一致。
- 图D显示，模型预测的稳态基因表达边际分布与真实MH和LB状态的数据吻合良好。
- 图E重构的能量景观清晰地显示了两个吸引子（MH和LB），速度场在肝内胚层（HE）阶段出现分叉，与实验观察一致。
- 图F, G：动态路径分析。
- 图F的状态转换网络显示，通过中间态的路径（如DE→HE→early LB）比直接跳跃（DE→LB）具有更低的“转换代价”（action），表明中间态是促进状态转换的。
- 图G的平均首达时间（FPT）分析显示，到达MH状态的时间比LB状态更短，暗示了向MH分化的偏好。
- 图H, I：In silico基因敲除实验。
- 图H显示，敲除MH特异性基因（如TOB1）会使细胞命运偏向LB状态，而敲除LB特异性基因（如RANBP1）则会偏向MH状态。
- 图I的FPT分析进一步量化了这一效应，例如敲除TOB1显著缩短了到达LB状态的时间。这些模拟结果与已知的生物学功能一致，验证了模型预测的可靠性。
实验四：EMT景观转换的解码
最后，作者将DynNet应用于一个包含EMT（上皮-间质转化）和MET（间质-上皮转化）过程的时间序列数据集。
- 图A, B：展示了从数据中直接估计的能量景观。在TGF-β撤除后，间质态（M）的能量升高，上皮态（E）的能量降低，驱动了MET的发生。
- 图C-E：DynNet在5个核心基因上训练后，准确地重构了MET过程。
- 图C的能量景观与数据分布吻合。
- 图D的稳态边际分布与真实数据吻合。
- 图E在成对基因表达空间中展示了不同状态的分布。
- 图F, G：
- 图F的速度场显示了向E和M两个稳态的流动，并暗示了M态在MET后仍保留一定的稳定性。
- 图G推断出的加权GRN与已知的生物学知识一致：同类型标志物基因间相互激活，不同类型间相互抑制。
- 图A, B, C：转换路径分析。
- 图A, B显示MET（M→E）更倾向于通过中间态（IM），而EMT（E→M）则更直接。
- 图C的转换网络定量地显示，通过中间态的路径具有更高的概率和更低的转换代价。
- 图D, E (EMT过程模拟)：通过在已训练好的MET模型中加入TGF-β对CDH1的抑制作用，DynNet能够准确地模拟EMT过程，其预测的稳态分布和能量景观与真实数据一致。
- 图F, G, H：关键基因识别。
- 图F的参数敏感性分析表明，增强M态标志物（如CD44）的相互激活会稳定M态，反之亦然。
- 图G显示了沿MET路径关键基因的表达动态。
- 图H的模拟敲除实验显示，敲除CD44会使M细胞回到E态，而敲除CDH1会促进E细胞向M态转变，这与它们在EMT中的核心作用一致。
