Discovering stochastic dynamical equations from ecologicaltime series data
从生态时间序列数据中发现随机动力学方程
https://arxiv.org/pdf/2205.02645v6
关键词：数据驱动模型发现，朗之万动力学，自组织，集体运动，介观尺度动力学，数据驱动动力系统，科学机器学习，噪声诱导有序。
摘要：
理论研究表明，随机性可以以反直觉的方式影响生态系统的动力学行为。然而，若缺乏描述种群或生态系统演化的控制方程，则难以在真实数据集中准确判定随机性所起的作用。因此，从数据中反推控制性随机微分方程的逆问题具有重要意义。本文提出一种方程发现方法：以状态变量的时间序列数据为输入，输出对应的随机微分方程。该方法通过将随机微积分中的传统方法与方程发现技术相结合而实现。我们通过若干应用案例验证了该方法的普适性：首先，我们刻意选取了若干具有本质不同控制方程的随机模型，但它们却产生近乎相同的稳态分布；结果表明，仅通过对时间序列数据的分析，我们即可准确恢复其真实的底层方程，并据此正确推断其稳定性结构。我们将该方法应用于两个真实世界数据集——鱼群聚集行为与单细胞迁移行为，二者在时空尺度与动力学特性上差异显著。此外，我们还阐明了该方法的若干局限性与潜在陷阱，并展示了如何通过诊断性指标加以识别与克服。最后，我们以开源软件包 PyDaDDy（Python Data Driven Dynamics 库）的形式公开了相关代码。
引言
对复杂生态系统的建模，跨学科普遍采用的核心方法是微分方程（Gotelli 等，2008；Strogatz，2018）。依据系统的维度及随机性的重要性，这些方程可为常微分方程（ODE）、随机微分方程（SDE）或偏微分方程（PDE）。即使我们最初仅从细粒度、局部尺度上的简单行为规则或生态相互作用出发，仍可推导出群组、种群乃至生态系统等宏观层次的粗粒化动力学描述，其形式即为微分方程（Biancalani 等，2014；Cheng 等，2014；Durrett 与 Levin，1994；Jhawar 等，2019；Loreau，2010；Majumder 等，2021；McKane 与 Newman，2004；Yates 等，2009）。尽管此类基于动力系统的方法功能强大，并已带来若干关键生物学洞见，但如何将实证数据与微分方程模型进行有意义的整合，仍是一项持续存在的挑战。
我们正处大数据时代，高分辨率数据对生态系统动态的刻画日益丰富（Leyk 等，2019；Nathan 等，2022）。此类数据覆盖多个生物学组织层次：从个体（细胞或动物）运动轨迹（Brückner 等，2019；Nathan 等，2022），到群体属性（Jhawar 等，2020；Tunstrøm 等，2013；Yates 等，2009），再到种群规模（Bjørnstad 与 Grenfell，2001；Stenseth 等，1997）、种群适应度（Lenski，2017）及生态系统状态（Carpenter 等，2020；Majumder 等，2019；Xie 等，2008）。尤为重要的是，部分数据具备高时间分辨率，有时甚至兼具空间信息，从而为模型与数据的深度融合开辟了新路径。
为准确刻画生物系统的动态，必须将状态变量视为非线性且随机的，而非仅考虑其线性或平均性质。适用于分析此类复杂随机效应的合适框架是随机微分方程（SDE），它能够同时捕捉驱动系统演化的确定性与随机性因素。在SDE模型中，噪声最广为人知的作用是导致系统围绕确定性稳定态发生波动——例如，种群受环境噪声影响而在承载容量平衡点附近起伏。然而，SDE模型亦预测：当噪声强度依赖于系统状态本身（即所谓状态依赖噪声）时，其可生成远离确定性稳定平衡点的新稳态（Horsthemke 与 Lefever，1984）。例如，在鱼群中，小群体规模所对应的波动反而会反直觉地使系统远离确定性稳定的无序态，从而提升群体层面的协调性（Biancalani 等，2014；Jhawar 等，2019, 2020）。在干旱生态系统研究中，数学模型预测降雨波动可在两个确定性稳定态之间诱导出一个新状态（D’Odorico 等，2005）；这与传统认知——即加性噪声通常仅引发系统在多个稳定态之间跃迁（Guttal 与 Jayaprakash，2007）——形成鲜明对比。然而，上述关于随机性作用的结论，仅在已知控制方程的前提下成立；而对真实生态数据而言，控制方程往往未知。
因此，我们自然转向其逆问题：能否仅基于观测到的时间序列数据，构建出相应的SDE模型？答案是肯定的。基于估计所谓“跳跃矩”（jump moments）的方法，原则上可从时间序列中反推随机微分方程（Friedrich 等，2011；Gradišek 等，2000；Rinn 等，2016；Tabar，2019）。此外，在确定性模型领域，近年发展的方程发现（equation discovery）技术，已能从时间序列中推断出简洁、可解释的微分方程模型（Brunton 等，2016；de Silva 等，2020；Rudy 等，2017）。近期，这些技术已被拓展应用于随机动力系统（Boninsegna 等，2018；Brückner 等，2020；Callaham 等，2021；Frishman 与 Ronceray，2020；Huang 等，2022），对生态时间序列数据分析展现出良好前景。
然而，仍存在显著空白：
其一，这些技术散见于物理与工程文献，生物学界对其了解有限；
其二，尚无根本理由确信那些假设噪声结构相对简单的SDE模型，适用于真实生物数据；
其三，验证噪声假设、检验模型正确性的诊断工具至关重要，但在物理与工程文献中新方法研发过程中，此类诊断常被忽视。
因此，目前随机方程发现方法尚难以直接适用于生物学研究。
本文旨在弥合上述鸿沟，提出一个统一框架，用于生物数据集的SDE模型发现与诊断。简言之，本方法允许用户输入一段时间序列，自动推断其底层随机微分方程，并执行诊断以检验所发现模型的有效性。我们认为该方法具备以下新颖性与优势：
第一，可直接从数据出发，获得简洁且可解释的随机动力系统（SDE）模型；
第二，模型推断过程基于输入时间序列在最精细时间尺度上的随机涨落分析，而所发现的动力学模型却能准确复现数据的长时尺度特征——需强调，这些长时特征并未作为建模输入；
第三，方程学习过程几乎无需用户预设函数形式，即可自动发现与给定动力学数据匹配的恰当函数结构，而非仅对用户预定义模型的参数进行拟合。
为验证方法的普适性，我们考察了若干应用场景：
（A）刻意选取确定项与随机项本质不同的SDE，其虽产生高度相似的稳态分布，仍能被准确区分；
（B）经典生态模型（含随机性）；
（C）两类已发表的真实世界数据集，尺度与动力学特性迥异：
(i) 由随机性驱动的鱼群运动（Jhawar 等，2020；见schooling_fish仓库）；
(ii) 由确定性极限环主导、随机性作用微弱的单细胞运动。
我们深入讨论了各类局限性，并详述如何利用诊断工具规避潜在陷阱、检验所发现模型是否满足基本假设并忠实描述数据。最后，为提升方法的可及性，我们以开源软件包PyDaDDy（Python Data Driven Dynamics 库）（Nabeel 等，2024）的形式开放了针对一维/二维数据集的分析代码，见：
https://github.com/tee-lab/PyDaddy
（存档于 https://doi.org/10.5281/zenodo.13777396）。
方法
随机微分方程的数学预备知识
我们的目标是：从对系统状态的测量（表现为时间序列）出发，建模该系统的动力学行为。典型示例包括：单个或多个相互作用种群的种群规模、觅食动物的运动轨迹、景观中的植被覆盖度等。为此，我们采用随机微分方程（SDE）框架。设 x ( t )
为一个 d 维向量，其时间演化可建模为：
逆问题背后的基本原理
我们感兴趣的是从有限采样频率的观测时间序列数据中推断SDE的逆问题。具体来说，我们的目标是找到简单、可解释的解析表达式来描述 f 和 G，而不仅仅是它们的定性形状。
我们解决这个问题的方法包括两个步骤。首先，我们从观测到的时间序列中计算漂移和扩散函数的瞬时估计。接下来，我们使用基于稀疏回归的技术来估计漂移和扩散函数的函数形式。
瞬时漂移和扩散系数
从具有有限采样时间 的观测采样时间序列中，我们可以估计瞬时漂移系数：
漂移系数与扩散系数的方程学习
我们采用一种基于稀疏回归的技术（有时称为“方程学习”），以获得所提取的漂移函数 f f与扩散函数 G G 的可解释解析表达式（Boninsegna 等，2018；Brunton 等，2016）。在无额外假设的前提下，该过程可对 f f 与 G G 的每一分量分别独立进行（即按分量处理）。
模型选择
除基函数外，方程发现算法尚需用户提供额外输入以约束模型选择过程：即稀疏化阈值λ 。该参数值决定了模型复杂性与解释能力之间的权衡——较高的 λ 会剔除更多项，从而生成更简洁的模型，但以拟合数据精度下降为代价。用户可手动设定 λ λ 以针对具体问题优化此权衡；亦可采用自动调参方法。
传统上可使用信息准则（如赤池信息量准则 AIC 或贝叶斯信息量准则 BIC）进行阈值选择；但本文转而采用k 折交叉验证（k-fold cross validation）——这是机器学习文献中常用于依据交叉验证精度选择理想模型的方法（Shalev-Shwartz 与 Ben-David，2014）。
其核心思想是：仅使用部分数据（称为训练集）训练模型，并在与训练集互斥的验证集上评估模型性能。若模型在验证集上表现良好，则说明其未对训练数据中的噪声过拟合，因而具备更强的泛化能力。具体而言，将数据集均分为 k 份；每次选取其中 k − 1
份作为训练集拟合模型，剩余 1 份作为验证集，并计算模型在该验证集上的误差（即验证误差）；重复此过程 k 次，使每一份数据均被用作一次验证集，最终取平均验证误差作为评估指标。
我们选择使交叉验证（CV）误差下降幅度最大的 λ 值作为最优阈值。对许多系统（包括真实世界系统），还可结合对该系统物理/生物学背景的先验知识，人工选定基函数并进行手动模型筛选，从而获得更优、更简洁的模型（参见 Nabeel 等，2023 及补充材料 SI 第 S4 节“真实数据集的模型选择”）。
诊断
数据驱动的随机微分方程（SDE）发现流程，若缺乏充分的诊断检验以合理验证SDE模型的基本假设，则是不完整的。基于 Brückner 等（2019）与 Jhawar 和 Guttal（2020）的思想，我们提出了三类诊断检验，用于评估所发现的SDE模型的可靠性：
- 噪声诊断：方程（1）中的噪声项 η ( t ) 被假设为无关联的高斯白噪声
与传统方法的对比
文献中常用于恢复随机微分方程（SDE）的一种传统方法是：将变量 x x 的取值范围划分为有限个区间（bins），并在每个区间内对 f f 和 G G 进行分箱平均估计。也就是说，漂移函数与扩散函数可由其瞬时形式（式 4、式 6）近似为：
其中， w w 为分箱宽度（bin-width），各区间为 d d 维矩形区间。该方法的缺点在于：每个分箱内的估计相互独立，无法利用任何全局结构信息（例如函数的光滑性），因而估计结果往往较为嘈杂。一种降低估计方差的策略是对时间序列以更大的时间步长 Δ t
进行子采样（Jhawar 与 Guttal，2020），但这可能引入估计偏差（参见补充材料 SI 第 S5.A 节“有限数据下的估计”）。
相比之下，方程学习（equation learning）方法用稀疏回归步骤替代了分箱平均步骤。除能直接输出可解释的解析表达式这一显著优势外，方程学习还具备另外两个优点：
第一，通过合理选择基函数库，可显式纳入全局结构或约束条件（如光滑性、对称性等）；
第二，无需子采样——从而消除了两个任意性较强的参数选择：子采样时间步长 Δ t
与分箱宽度 ε 。
结果
从合成数据集中发现具有差异性的随机微分方程（SDE）
由不同SDE生成的单峰分布
单峰分布广泛见于诸多生物数据集中。然而，此类分布可能源于截然不同的动力学过程，因而可由本质迥异的底层SDE建模。我们旨在表明：即便在此类情形下，应用本文提出的方法仍可从时间序列数据中准确复原原始SDE。
我们考虑以下两个一维随机微分方程：
尽管方程11和12形式迥异，但由它们生成的时间序列及其直方图却极为相似（图2 (A-i, A-ii, B-i, B-ii)）：事实上，x(t) 的稳态分布均为单峰，且对于方程11和12几乎完全相同。从动力系统角度看，方程11在 x* = 0 处存在一个确定性稳定态。在生物学上，x 可被视为种群规模偏离其稳定承载容量的偏差，因此 x 可正可负。加性噪声——代表环境波动——仅使种群动态围绕确定性平衡点扩散，这在时间序列及直方图中均有体现（图2 (A-ii)）。因此，我们可将直方图的峰值视为反映了一个确定性状态。这确实是噪声最广为人知的作用。
另一方面，对于方程12，其确定性稳定平衡点位于 ±1，但直方图的峰值却在 0（图2 (B-ii)）。在此情况下，状态依赖噪声或乘性噪声项彻底改变了系统的稳定性景观，在两个确定性稳定态之间催生出一个新状态，导致峰值出现在 x = 0。这是噪声产生异常效应的一个例子，即“噪声诱导的稳定态”，该概念最初在干旱生态系统模型中被提出（D’Odorico 等，2005）。
针对这些单峰数据集，我们现在提出逆问题：能否仅基于时间序列数据的特征（图2 A-i 与 B-i），推断出正确的底层SDE模型？确实，我们所提出的、结合跳跃矩计算与稀疏回归的方法，能够准确识别出这两个模型的漂移函数与扩散函数的具体形式（图2中，比较A列和B列内的红色线与黑色线；第iii行和第iv行）。至关重要的是，我们能够推断出：图2 A-i所示的时间序列由一个线性漂移项（图2 A-iii）与加性噪声驱动（图2 A-iv）；而图2 B-i所示的时间序列则由一个具有三个根的非线性漂移函数（图2 B-iii）与乘性噪声驱动（图2 B-iv）。此外，我们还能恢复出漂移函数与扩散函数的符号化解析表达式，其结果与我们用于生成合成数据集的原始SDE高度吻合（图2，“Estimated”面板）。
由不同SDE生成的双峰分布
许多生物系统表现出替代稳定态，最简单的情形是双稳态系统，例如干旱生态系统中的草原与林地状态，或湖泊的富营养化与贫营养化状态。具有双稳态的系统，其状态变量通常呈现双峰分布。同样，这种双峰性也可由本质不同的底层过程/SDE生成。为说明这一点，我们考虑以下两个玩具模型：
方程13在 ±√(2/3) 处有两个确定性稳定平衡点，其动力学行为由加性噪声（图2 (C-ii)）扩散至这些平衡点周围；因此，这些是确定性状态。另一方面，方程14仅在 x* = 0 处有一个稳定平衡点；然而，由于乘性噪声项的作用，系统在直方图中呈现出两个远离确定性平衡点的峰值（图2 (D-ii)）；因此，这些是噪声诱导的状态。
我们再次针对这些双峰时间序列数据集（图2 C-i 与 D-i 对比，以及图2 C-ii 与 D-ii 对比）提出逆问题。在此情况下，我们的方法同样能够准确恢复底层的SDE模型，包括其符号函数（参见图2 C和D中的第iii、iv行及“Estimated”面板）。
在理论生态学经典模型中的验证
我们现在将该方法应用于若干理论生态学中的经典模型进行验证。我们考虑了多种单变量和双变量模型的随机化版本。所考虑的单变量模型包括：密度依赖种群增长的逻辑斯蒂模型（不稳定系统）、具有Holling III型功能反应的种群捕捞模型（双稳态系统）（Strogatz, 2018），以及湖泊富营养化模型（Carpenter 等, 1999）（双稳态系统）。所考虑的双变量模型包括：用于种间竞争的Lotka-Volterra竞争模型、具有Holling II型功能反应的非线性捕食者-猎物模型（Alonso 等, 2002），以及Van der Pol振荡器——一个用于描述非线性振荡的最小模型（Strogatz, 2018）。
表1总结了本分析的结果。总而言之，对于所有考虑的模型，SDE发现程序均能从模拟的时间序列中准确地发现正确的模型（详见补充信息SI第S3节——“生态学经典模型验证”）。
有限数据下的不准确模型发现
现实世界的数据存在诸多局限性，例如时间序列较短或采样分辨率不足。理解此类情境如何影响我们所描述的数据发现协议至关重要。此处，我们考察“有限数据”的两种特定子情况：短时间序列和长采样间隔。
为此，我们尝试估计前文所述的方程4中的模型：ẋ = 2x - 3x³ + (1/2)η(t)，该模型表现出双稳态动力学。随后，我们使用一个短时间序列（1000个高分辨率观测值，Δt = 0.01s），其中系统仅探索了状态空间的一小部分。所得模型及其模型诊断结果如补充信息图S11 A、B所示。不出所料，估计的SDE与真实的基础模型并不匹配。当从具有大采样间隔的数据中估计模型时，也会发生类似的错配（补充信息图S11 C、D）。在此情况下，大的采样间隔意味着关于细尺度动力学的信息丢失，因此发现的模型再次偏离实际模型。然而，在这两种情况下，我们都观察到原始时间序列与模型模拟时间序列的自相关函数之间存在巨大差异，这将提醒用户模型估计不准确。另一个因数据有限或嘈杂而可能导致模型误识别的明显迹象是，估计系数中出现较大的误差条，我们也观察到了这一点。
我们在讨论部分详细探讨了这些局限性问题、使用方法时如何避免陷阱、以及如何通过诊断评估所发现的模型，并在补充信息SI第S5节“陷阱与局限性”中对这些局限性进行了定量刻画。
在真实数据集上的应用
起初并不明确这些方法是否适用于复杂的生物数据集。尤其是，随机微分方程（SDE）假设噪声为高斯且无关联的，而实际数据中我们往往预期存在更复杂的噪声结构（例如，非高斯性或具有记忆性）。因此，为验证其适用性与普适性，我们将其应用于两个性质迥异的真实数据集：鱼群聚集运动的群体轨迹数据（Jhawar 等，2020）与受限单细胞迁移数据（Brückner 等，2019）。
这两个数据集在多个关键方面存在显著差异：
第一，鱼群数据刻画了群体层面涌现的秩序性，而细胞数据反映的是个体层面行为；
第二，鱼群的空间尺度相对较大（约1米），但时间动态较快（约0.1秒量级）；而细胞数据空间尺度极小（约10⁻⁴米），运动却极为缓慢，时间尺度达10分钟；
第三，其底层动力学本质不同——鱼群处于噪声诱导态（Jhawar 等，2020），而细胞运动主要由确定性极限环主导，噪声仅起次要作用（Brückner 等，2019）。
鱼群噪声诱导聚集行为的SDE发现
近十年来，动物集体运动在多个生物学领域均取得显著进展。借助前沿成像技术对运动群体进行记录，研究者已获得高时空分辨率的数据，有助于深入探索同步集体运动的机制。集体运动的数学理论常以随机微分方程形式表达其动力学，为这种同步群体行为提供可能解释（Biancalani 等，2014；Jhawar 等，2019；Ramaswamy，2017；Toner 等，2005；Vicsek 与 Zafeiris，2012）。随着高质量动物运动数据的可及性提升，我们现在可提出其逆问题：能否从给定的动物轨迹时间序列中发现随机动力学模型？进一步地，还可探讨随机性在塑造（或破坏）秩序中的根本作用：观测到的集体动力学，究竟与确定性状态一致（即随机性仅在确定性稳定平衡点附近模糊秩序，如图2 A-ii所示），还是与噪声诱导状态一致（即随机性在远离确定性稳定平衡点处生成非平凡状态，如图2 C-ii所示）？
近期一项研究推断，鱼群高度同步的运动行为是一种噪声诱导态（Jhawar 等，2020）。该研究在群体动力学时间序列上采用传统跳跃矩估计方法实现推断。本文中，我们检验能否通过结合跳跃矩计算与稀疏回归的推断流程，复现相同结果；此外，我们还执行了该研究中被忽略的噪声诊断与模型诊断。
我们使用（Jhawar 等，2020）公开发布的Etroplus suratensis（卡拉米恩鱼）群体数据集。实验中，15尾鱼组成的鱼群运动通过高分辨率摄像机记录，并经计算机视觉方法追踪，生成一个30维时间序列（15尾鱼 × 每尾2维位置坐标）。直接在30维空间建模动力学，信息增益有限；相比之下，采用能捕捉群体本质动力学的低维表征更利于理解集体行为。
受物理学文献（Toner 等，2005；Vicsek 等，1995）及其在生物集体运动中的广泛应用（Couzin 等，2002；Jhawar 等，2020；Tunstrøm 等，2013）启发，我们聚焦于群体极化秩序的涌现与演化过程。
该数据集提供了极化向量 m m 的时间序列，采样间隔为0.12秒，持续约1小时，但含若干缺失点（对应追踪失效时刻）。尽管群体极化时间序列展现出显著随机涨落（图3B），其直方图却表明鱼群总体处于有序状态（图3C、D）。此即为我们分析的输入时间序列。
此处，我们借助所开发的PyDaDDy软件包（整合跳跃矩估计与稀疏回归）复现了上述结果（另见补充材料 SI 第 S4.A 节“鱼群数据集的模型选择”）。所发现的方程包含线性漂移项与二次型扩散项，对应的向量随机微分方程形式为：
图4展示了对该发现的SDE模型进行诊断测试的结果。噪声残差 r ( t )
呈高斯分布，符合预期（图4A）。残差噪声 r ( t )中的相关性迅速衰减（图4B）。这些测试在合理程度上支持了关于 η 为高斯白噪声的建模假设。
我们还发现，所发现的方程通过了模型诊断测试。由方程15模拟生成的SDE的直方图与原始时间序列中 m m 的直方图高度吻合（图4C）。模拟时间序列的自相关函数也与原始时间序列表现出合理的吻合（图4D），但有一个显著偏差：数据自相关函数在约10秒的相对较长的时间尺度上呈现负值，之后才收敛至零。鱼群研究论文的作者们已证明，这一特征源于边界效应，对于 Etroplus suratensis 的鱼群动力学而言并不重要。尽管如此，我们强调，本研究所用的SDE发现协议仅利用了极小时间尺度（约0.1秒）下群体极化波动的高度局部信息；它并未使用任何关于群体极化数据频率分布或自相关函数的信息。尽管如此，模拟时间序列在这两个指标上仍与数据表现出良好的一致性。最后但同样重要的是，我们发现该模型具有自洽性——当我们将模拟数据再次输入我们的方程发现协议时，能够恢复出相同的SDE。
受限细胞迁移动力学的SDE发现
从形态发生、伤口愈合到癌症转移，细胞迁移在众多生物学情境中扮演着关键角色，理解细胞如何在复杂环境中移动至关重要。近期一项研究（Brückner 等，2019）探讨了结构化环境中细胞迁移的随机动力学。作者设计了一项实验：一个癌细胞在两个“岛屿”（即“状态”）之间来回迁移（图5A），并使用SDE框架对该二态迁移的动力学进行建模，但未采用显式表达式。他们得出结论：该动力学主要由确定性（漂移）分量中的极限环主导，而随机分量仅起次要作用，影响的是状态间转换的时间尺度。本文中，我们利用所提出的SDE发现协议复现了该研究的关键结果。虽然原研究通过分箱平均、分段表示的方法构建漂移与扩散函数，我们的方法则能够提取出可解释的函数形式。此外，我们还展示了噪声诊断与模型诊断如何引导更精确的模型发现。
该数据集包含149条独立的细胞轨迹重复实验。每条轨迹基于高分辨率图像，每10分钟采集一次，最长持续50小时。图5B展示了细胞在两个岛屿间迁移的一个时间序列示例，图5C-E则描绘了状态变量（位置x和速度v）的直方图。
根据实验设计（图5A）以及细胞位置的双峰直方图（图5C），一个自然的初步尝试是将细胞在两个稳定状态（即着陆垫）之间的跳跃建模为纯粹随机的切换过程。这类似于图2C中描述的模型，其特征是在两个稳定态之间发生随机跃迁。这需要将细胞轨迹 x(t) 建模为如下形式的过阻尼SDE：
然而，我们的模型诊断结果表明，过阻尼模型不足以完整刻画该系统的动力学行为（参见补充材料 SI 第 S4.B 节“细胞迁移数据集的模型选择”）。
因此，我们转而采用欠阻尼模型，引入两个动力学变量——位置（ x ）——对该细胞轨迹进行建模。此类二维描述能够捕捉系统中可能存在的线性或非线性振荡行为；如下文所示，该模型确实最符合实际数据。
该欠阻尼随机动力学遵循如下方程：
其中，我们在应用SDE发现流程前，对位置 x x 和速度 v v 进行了无量纲化缩放。
此缩放步骤对稀疏回归的正确运行是必要的（参见补充材料 SI 第 S4.B 节）。
通过我们的SDE发现流程，我们发现漂移函数 f f（图5E）可很好地用如下形式的三次多项式近似：
我们通过数据发现的细胞迁移SDE模型呈现出可激发流（excitable flow），其定性行为与 Brückner 等（2019）的研究结果一致，表现为弛豫振荡（relaxation oscillations）。此外，我们发现其漂移函数偏离了经典的范德波尔（Van der Pol）振荡器模型——后者是解释弛豫振荡的最简模型。有趣的是，如补充材料第 S4.B 节所示，我们发现扩散项（diffusion term）对准确刻画细胞动力学中的有限边界效应、并满足模型诊断要求至关重要（见下文）。
图6展示了所发现模型的诊断结果：噪声残差在时间上无自相关性（图6B），但其分布呈现比高斯分布更强的中心聚集性（即更尖锐的峰度，图6A）。尽管如此，该模型生成的模拟时间序列在统计特性上（如状态变量直方图与自相关函数）与原始数据总体吻合良好（图6C–F）。此外，该模型具有自洽性（self-consistency）。
此处我们强调模型诊断在准确发现扩散项过程中所起的关键作用：若采用加性噪声（即常数扩散项）配合相同漂移函数，模拟轨迹将频繁突破实验设定的物理边界，导致模拟状态变量直方图与真实数据存在巨大偏差。相比之下，引入四阶扩散项可显著减少此类边界违反行为，从而提升模型自洽性。因此，我们认为该扩散项实质上有效编码了空间约束。
所发现模型的泛化能力验证
为验证SDE发现方法能否获得可泛化模型，我们针对鱼群与细胞迁移两个数据集，均采用留出数据验证（见“诊断”一节中“使用留出数据进行验证”子节）：
- 对鱼群数据集，将4次实验中的2次作为训练集；
- 对细胞迁移数据集，将149条轨迹中的75条用于训练。
在训练集上完成模型发现后，我们分别计算：（i）由所发现模型生成的时间序列的直方图与自相关函数；（ii）留出验证集的对应统计量。
结果表明，对鱼群与细胞两个数据集，模型生成数据的统计特性均与验证集高度匹配，说明所发现的SDE成功捕获了底层动力学的可泛化特征（详见补充材料第 S4.C 节“所发现SDE模型的泛化能力”）。
本文提出了一种方法：以经验观测的时间序列作为输入，发现具备解析可解释性的数据驱动型随机动力学方程。为实现这一目标，我们将传统的跳跃矩法（用于漂移与扩散估计；Tabar, 2019；Van Kampen, 1992）与基于稀疏回归的方程发现技术（Boninsegna 等，2018；Brunton 等，2016）相结合。尤为重要的是，我们强调了对数据驱动模型所依赖假设及其发现结果有效性进行诊断检验的必要性。除在已知真实方程的合成数据集中验证该方法外，我们还通过两个性质迥异的生物时间序列数据集（鱼群聚集与细胞迁移）展示了其普适性。为便于更广泛的生物学研究者使用，我们开发并开源了一个易于使用的 Python 软件包——PyDaDDy（https://pydaddy.readthedocs.io/）。
然而，通过数据驱动的方程发现，我们揭示出：
- 鱼群同步行为实为一种噪声诱导态，其稳定状态远离确定性稳定平衡点；
- 细胞迁移的振荡则源于确定性极限环，随机性仅起次要作用；
- 且该极限环的数学结构不同于经典的范德波尔振子。
换言之，在用户输入极少干预的前提下，我们通过符号化表达的数据驱动动力学模型，成功排除了若干竞争性假说，逼近更真实的机制。
我们强调，该方法具有极高的计算效率。例如，针对软件包中提供的鱼群数据集（二维向量时间序列，约25,000个时间点）：
- 在一台2023款Mac Mini（Apple M2处理器，16GB内存）上，漂移与扩散系数的初始估计仅需约800毫秒，漂移函数拟合约2毫秒，扩散函数拟合约6毫秒；
- 在Google Colab实例上，对应耗时分别约为2.5秒、10毫秒与60毫秒。
陷阱及其规避方法
真实生物系统高度复杂，其“真实”底层动力学可能涉及数十乃至上百个相互作用变量。此外，对此类系统的观测数据常具有以下特征：不完整（仅观测到部分动力学变量）、不精确（含测量噪声，或采样稀疏、时间间隔大）、缺失（若干时间点数据缺失），或代表性不足（数据仅覆盖状态空间的一小部分）。在此情形下，是否仍可用仅含少数变量的随机微分方程（SDE）建模，尚不明确。这就引出一个关键问题：如何确保从经验数据中获得的数据驱动SDE模型确实是系统动力学的有效表征？若盲目将SDE发现流程应用于任意数据集，流程虽总能输出某个SDE，但该模型是否有效、是否忠实反映系统动力学，却无法保证。因此，必须审慎评估SDE发现技术对特定问题或数据集的适用性。
本文方法论（及配套软件包PyDaDDy）的核心理念之一，正是为用户提供模型评估与诊断工具，辅助其作出合理判断。本节及补充材料SI第S5节将讨论在用低维SDE建模生物系统时常见的陷阱、识别方法及（可能的）修正策略。
系统维度问题
原则上，数据驱动方程发现方法适用于任意维度。本文聚焦于一维与二维数据集，旨在阐明随机动力学方程推断的核心原理与应用。为展示方法普适性，我们提供了一个恢复含随机性的三维经典混沌系统——洛伦兹方程的示例代码（见笔记本：Higher dimensions）。
然而高维应用存在若干实际挑战：第一，维度增加导致基函数库规模急剧膨胀，估计易受误差干扰；稀疏回归作为一种正则化手段，可在一定程度上缓解该问题；第二，结合对系统物理、约束与对称性的先验知识，可精心构建小规模基函数库（Nabeel 等，2023）；第三，高维模型的验证本身即具挑战性——需发展不依赖可视化（一/二维可行）的定量统计诊断方法，此为未来研究方向。
需强调的是：当目标是构建基于动力学变量的可解释方程时，天然具备低维表征能力的系统尤为适合数据发现方法。例如，即便原始数据高维（如鱼群数据为 2 N 2N 维， N N 为鱼数），仍常可识别出少数隐变量（latent dimensions）以刻画系统本质动力学。此类降维不仅方便分析，更具理论价值——诸多理论本就建立于对高维生物系统的粗粒化描述之上（Durrett & Levin, 1994），此思路源于物理与应用数学。隐变量可基于粗粒化理论或生物学意义量选定（如鱼群研究中受介观集体运动理论启发，选用二维极化向量）；亦可通过数据驱动技术识别（Greenacre 等，2022；Schmid，2022；Van Der Maaten 等，2009）。当前已有工作探索在同时发现合适潜空间与该空间内动力学方程（针对确定性系统，Champion 等，2019）；将此类方法拓展至随机系统，是值得推进的方向。
然而，某些情形下系统动力学本质高维，或观测变量不足以完备刻画动力学。此时，仅基于观测变量构建的SDE（或任何低维模型）必无法完全捕获系统行为。对此，PyDaDDy 的诊断工具至关重要——可辅助判断模型是否误用。具体而言：
- 若存在未观测变量，残差自相关衰减时间尺度可能显著延长，提示SDE模型未能解释的慢动力学残余；
- 模拟数据与真实数据的直方图或自相关函数若无法匹配，亦表明模型不充分（如以过阻尼SDE建模细胞迁移时失效，见SI S4.B）；
- 另一案例：仅用单物种时间序列拟合两物种相互作用模型，诊断测试即揭示降维模型不完备（见SI S5.B）。
与所有估计方法一样，数据驱动SDE发现亦需质量足够、数量充分的数据。分析表明：时间序列过短或采样频率过低均会导致SDE估计失准。值得庆幸的是，此类情形下模型诊断通常表现不佳，可警示用户谨慎解读结果。总体而言，重建所需数据量取决于底层模型复杂度与测量噪声水平（Fajardo-Fontiveros 等，2023）。尽管难以给出普适数据量阈值，我们在SI S5.A节对PyDaDDy在有限数据下的性能进行了数值探索，发现：方程发现技术在数据受限时，显著优于传统的分箱Kramers–Moyal平均法。
关于测量噪声：
- 当测量噪声较小时，方法仍可有效工作；漂移函数估计在平均意义上基本不受影响；
- 但扩散函数估计会系统性偏高，偏移量约等于测量噪声方差（Böttcher 等，2006）——若目标为构建预测模型，此偏差或可接受；
- 若目标为构建机理模型，则需显式分离测量噪声与生物本征随机性。更广泛地，不同噪声结构（变量间相关、时间相关、测量噪声）如何影响SDE推断，是值得深入研究的开放问题。
最后，时间序列中偶发缺失数据点（如传感器故障所致）影响较小——SDE估计可直接基于可用数据点进行，如鱼群数据集所示。
相关工作与结论性评述
方程学习（equation learning）最初在工程学文献中提出，用于常微分方程与偏微分方程的发现，并已被集成于如PySINDy（de Silva 等，2020）和DataDrivenDiffEq.jl（JuliusMartensen 等，2021）等流行工具包中。然而，这些方法局限于确定性微分方程的发现。虽有少数例外，例如 R 包Langevin（Rinn 等，2016），其可接受一维或二维时间序列输入，并以分箱平均方式估计克雷默–莫亚尔（Kramers–Moyal）系数来获得漂移与扩散项，且提供部分可视化与诊断功能，但其未输出可解释的SDE解析表达式；而且，分箱平均估计本身易受误差影响（Callaham 等，2021；Jhawar 与 Guttal，2020），尤其在数据有限情形下（见补充材料 SI 第 S5.A 节）——而本文所采用的方程学习方法极大缓解了这一问题。
我们强调：绝大多数数据驱动建模工具（包括在动力系统建模领域广泛应用的DiffEqParamEstim.jl（Rackauckas 与 Nie，2017）与Sim.DiffProc（Guidoum 与 Boukhetala，2020）），均采用传统范式——即要求用户预先指定漂移与扩散函数的参数化形式。在此类方法中，用户需先提出一组候选模型，再通过模型选择程序确定最适动力系统。典型地，随机动力系统模型的选择准则依赖于稳态分布与实测数据的匹配程度。然而，此类方法无法识别图2 A-i 与 B-i 所示情形：两个时间序列具有几乎完全相同的稳态分布，却源于本质不同的动力学机制。相比之下，我们的方法不仅可直接从时间序列中发现SDE，更能准确区分本质不同的控制方程。
当前，生物学领域正快速积累跨尺度的大型数据集；随着大数据时代的来临，我们期望：本文提出的这一方法，配合相对易用的开源工具包PyDaDDy，将激励研究者更广泛地将其应用于从数据中刻画控制动力学方程的任务中。借此，我们可在充分考虑系统本征随机性的同时，挖掘具机理基础的预测模型的潜力。
当然，将此类方法拓展至更复杂的现实世界数据集（如高维系统、含时变参数、受观测噪声干扰、具有复杂噪声结构——包括有色噪声与非高斯噪声等），仍将面临诸多挑战；而这些挑战，也正是未来研究的重要机遇所在。
原文链接：https://arxiv.org/pdf/2205.02645v6
热门跟贴