Discovering stochastic dynamical equations from ecologicaltime series data

從生態(tài)時(shí)間序列數(shù)據(jù)中發(fā)現(xiàn)隨機(jī)動(dòng)力學(xué)方程

https://arxiv.org/pdf/2205.02645v6

關(guān)鍵詞：數(shù)據(jù)驅(qū)動(dòng)模型發(fā)現(xiàn)，朗之萬(wàn)動(dòng)力學(xué)，自組織，集體運(yùn)動(dòng)，介觀尺度動(dòng)力學(xué)，數(shù)據(jù)驅(qū)動(dòng)動(dòng)力系統(tǒng)，科學(xué)機(jī)器學(xué)習(xí)，噪聲誘導(dǎo)有序。

摘要：

理論研究表明，隨機(jī)性可以以反直覺(jué)的方式影響生態(tài)系統(tǒng)的動(dòng)力學(xué)行為。然而，若缺乏描述種群或生態(tài)系統(tǒng)演化的控制方程，則難以在真實(shí)數(shù)據(jù)集中準(zhǔn)確判定隨機(jī)性所起的作用。因此，從數(shù)據(jù)中反推控制性隨機(jī)微分方程的逆問(wèn)題具有重要意義。本文提出一種方程發(fā)現(xiàn)方法：以狀態(tài)變量的時(shí)間序列數(shù)據(jù)為輸入，輸出對(duì)應(yīng)的隨機(jī)微分方程。該方法通過(guò)將隨機(jī)微積分中的傳統(tǒng)方法與方程發(fā)現(xiàn)技術(shù)相結(jié)合而實(shí)現(xiàn)。我們通過(guò)若干應(yīng)用案例驗(yàn)證了該方法的普適性：首先，我們刻意選取了若干具有本質(zhì)不同控制方程的隨機(jī)模型，但它們卻產(chǎn)生近乎相同的穩(wěn)態(tài)分布；結(jié)果表明，僅通過(guò)對(duì)時(shí)間序列數(shù)據(jù)的分析，我們即可準(zhǔn)確恢復(fù)其真實(shí)的底層方程，并據(jù)此正確推斷其穩(wěn)定性結(jié)構(gòu)。我們將該方法應(yīng)用于兩個(gè)真實(shí)世界數(shù)據(jù)集——魚(yú)群聚集行為與單細(xì)胞遷移行為，二者在時(shí)空尺度與動(dòng)力學(xué)特性上差異顯著。此外，我們還闡明了該方法的若干局限性與潛在陷阱，并展示了如何通過(guò)診斷性指標(biāo)加以識(shí)別與克服。最后，我們以開(kāi)源軟件包 PyDaDDy（Python Data Driven Dynamics 庫(kù)）的形式公開(kāi)了相關(guān)代碼。

引言
對(duì)復(fù)雜生態(tài)系統(tǒng)的建模，跨學(xué)科普遍采用的核心方法是微分方程（Gotelli 等，2008；Strogatz，2018）。依據(jù)系統(tǒng)的維度及隨機(jī)性的重要性，這些方程可為常微分方程（ODE）、隨機(jī)微分方程（SDE）或偏微分方程（PDE）。即使我們最初僅從細(xì)粒度、局部尺度上的簡(jiǎn)單行為規(guī)則或生態(tài)相互作用出發(fā)，仍可推導(dǎo)出群組、種群乃至生態(tài)系統(tǒng)等宏觀層次的粗粒化動(dòng)力學(xué)描述，其形式即為微分方程（Biancalani 等，2014；Cheng 等，2014；Durrett 與 Levin，1994；Jhawar 等，2019；Loreau，2010；Majumder 等，2021；McKane 與 Newman，2004；Yates 等，2009）。盡管此類基于動(dòng)力系統(tǒng)的方法功能強(qiáng)大，并已帶來(lái)若干關(guān)鍵生物學(xué)洞見(jiàn)，但如何將實(shí)證數(shù)據(jù)與微分方程模型進(jìn)行有意義的整合，仍是一項(xiàng)持續(xù)存在的挑戰(zhàn)。

我們正處大數(shù)據(jù)時(shí)代，高分辨率數(shù)據(jù)對(duì)生態(tài)系統(tǒng)動(dòng)態(tài)的刻畫(huà)日益豐富（Leyk 等，2019；Nathan 等，2022）。此類數(shù)據(jù)覆蓋多個(gè)生物學(xué)組織層次：從個(gè)體（細(xì)胞或動(dòng)物）運(yùn)動(dòng)軌跡（Brückner 等，2019；Nathan 等，2022），到群體屬性（Jhawar 等，2020；Tunstr?m 等，2013；Yates 等，2009），再到種群規(guī)模（Bj?rnstad 與 Grenfell，2001；Stenseth 等，1997）、種群適應(yīng)度（Lenski，2017）及生態(tài)系統(tǒng)狀態(tài)（Carpenter 等，2020；Majumder 等，2019；Xie 等，2008）。尤為重要的是，部分?jǐn)?shù)據(jù)具備高時(shí)間分辨率，有時(shí)甚至兼具空間信息，從而為模型與數(shù)據(jù)的深度融合開(kāi)辟了新路徑。

為準(zhǔn)確刻畫(huà)生物系統(tǒng)的動(dòng)態(tài)，必須將狀態(tài)變量視為非線性且隨機(jī)的，而非僅考慮其線性或平均性質(zhì)。適用于分析此類復(fù)雜隨機(jī)效應(yīng)的合適框架是隨機(jī)微分方程（SDE），它能夠同時(shí)捕捉驅(qū)動(dòng)系統(tǒng)演化的確定性與隨機(jī)性因素。在SDE模型中，噪聲最廣為人知的作用是導(dǎo)致系統(tǒng)圍繞確定性穩(wěn)定態(tài)發(fā)生波動(dòng)——例如，種群受環(huán)境噪聲影響而在承載容量平衡點(diǎn)附近起伏。然而，SDE模型亦預(yù)測(cè)：當(dāng)噪聲強(qiáng)度依賴于系統(tǒng)狀態(tài)本身（即所謂狀態(tài)依賴噪聲）時(shí)，其可生成遠(yuǎn)離確定性穩(wěn)定平衡點(diǎn)的新穩(wěn)態(tài)（Horsthemke 與 Lefever，1984）。例如，在魚(yú)群中，小群體規(guī)模所對(duì)應(yīng)的波動(dòng)反而會(huì)反直覺(jué)地使系統(tǒng)遠(yuǎn)離確定性穩(wěn)定的無(wú)序態(tài)，從而提升群體層面的協(xié)調(diào)性（Biancalani 等，2014；Jhawar 等，2019, 2020）。在干旱生態(tài)系統(tǒng)研究中，數(shù)學(xué)模型預(yù)測(cè)降雨波動(dòng)可在兩個(gè)確定性穩(wěn)定態(tài)之間誘導(dǎo)出一個(gè)新?tīng)顟B(tài)（D’Odorico 等，2005）；這與傳統(tǒng)認(rèn)知——即加性噪聲通常僅引發(fā)系統(tǒng)在多個(gè)穩(wěn)定態(tài)之間躍遷（Guttal 與 Jayaprakash，2007）——形成鮮明對(duì)比。然而，上述關(guān)于隨機(jī)性作用的結(jié)論，僅在已知控制方程的前提下成立；而對(duì)真實(shí)生態(tài)數(shù)據(jù)而言，控制方程往往未知。

因此，我們自然轉(zhuǎn)向其逆問(wèn)題：能否僅基于觀測(cè)到的時(shí)間序列數(shù)據(jù)，構(gòu)建出相應(yīng)的SDE模型？答案是肯定的?；诠烙?jì)所謂“跳躍矩”（jump moments）的方法，原則上可從時(shí)間序列中反推隨機(jī)微分方程（Friedrich 等，2011；Gradi?ek 等，2000；Rinn 等，2016；Tabar，2019）。此外，在確定性模型領(lǐng)域，近年發(fā)展的方程發(fā)現(xiàn)（equation discovery）技術(shù)，已能從時(shí)間序列中推斷出簡(jiǎn)潔、可解釋的微分方程模型（Brunton 等，2016；de Silva 等，2020；Rudy 等，2017）。近期，這些技術(shù)已被拓展應(yīng)用于隨機(jī)動(dòng)力系統(tǒng)（Boninsegna 等，2018；Brückner 等，2020；Callaham 等，2021；Frishman 與 Ronceray，2020；Huang 等，2022），對(duì)生態(tài)時(shí)間序列數(shù)據(jù)分析展現(xiàn)出良好前景。

然而，仍存在顯著空白：
其一，這些技術(shù)散見(jiàn)于物理與工程文獻(xiàn)，生物學(xué)界對(duì)其了解有限；
其二，尚無(wú)根本理由確信那些假設(shè)噪聲結(jié)構(gòu)相對(duì)簡(jiǎn)單的SDE模型，適用于真實(shí)生物數(shù)據(jù)；
其三，驗(yàn)證噪聲假設(shè)、檢驗(yàn)?zāi)Ｐ驼_性的診斷工具至關(guān)重要，但在物理與工程文獻(xiàn)中新方法研發(fā)過(guò)程中，此類診斷常被忽視。
因此，目前隨機(jī)方程發(fā)現(xiàn)方法尚難以直接適用于生物學(xué)研究。

本文旨在彌合上述鴻溝，提出一個(gè)統(tǒng)一框架，用于生物數(shù)據(jù)集的SDE模型發(fā)現(xiàn)與診斷。簡(jiǎn)言之，本方法允許用戶輸入一段時(shí)間序列，自動(dòng)推斷其底層隨機(jī)微分方程，并執(zhí)行診斷以檢驗(yàn)所發(fā)現(xiàn)模型的有效性。我們認(rèn)為該方法具備以下新穎性與優(yōu)勢(shì)：
第一，可直接從數(shù)據(jù)出發(fā)，獲得簡(jiǎn)潔且可解釋的隨機(jī)動(dòng)力系統(tǒng)（SDE）模型；
第二，模型推斷過(guò)程基于輸入時(shí)間序列在最精細(xì)時(shí)間尺度上的隨機(jī)漲落分析，而所發(fā)現(xiàn)的動(dòng)力學(xué)模型卻能準(zhǔn)確復(fù)現(xiàn)數(shù)據(jù)的長(zhǎng)時(shí)尺度特征——需強(qiáng)調(diào)，這些長(zhǎng)時(shí)特征并未作為建模輸入；
第三，方程學(xué)習(xí)過(guò)程幾乎無(wú)需用戶預(yù)設(shè)函數(shù)形式，即可自動(dòng)發(fā)現(xiàn)與給定動(dòng)力學(xué)數(shù)據(jù)匹配的恰當(dāng)函數(shù)結(jié)構(gòu)，而非僅對(duì)用戶預(yù)定義模型的參數(shù)進(jìn)行擬合。

為驗(yàn)證方法的普適性，我們考察了若干應(yīng)用場(chǎng)景：
（A）刻意選取確定項(xiàng)與隨機(jī)項(xiàng)本質(zhì)不同的SDE，其雖產(chǎn)生高度相似的穩(wěn)態(tài)分布，仍能被準(zhǔn)確區(qū)分；
（B）經(jīng)典生態(tài)模型（含隨機(jī)性）；
（C）兩類已發(fā)表的真實(shí)世界數(shù)據(jù)集，尺度與動(dòng)力學(xué)特性迥異：
(i) 由隨機(jī)性驅(qū)動(dòng)的魚(yú)群運(yùn)動(dòng)（Jhawar 等，2020；見(jiàn)schooling_fish倉(cāng)庫(kù)）；
(ii) 由確定性極限環(huán)主導(dǎo)、隨機(jī)性作用微弱的單細(xì)胞運(yùn)動(dòng)。

我們深入討論了各類局限性，并詳述如何利用診斷工具規(guī)避潛在陷阱、檢驗(yàn)所發(fā)現(xiàn)模型是否滿足基本假設(shè)并忠實(shí)描述數(shù)據(jù)。最后，為提升方法的可及性，我們以開(kāi)源軟件包PyDaDDy（Python Data Driven Dynamics 庫(kù)）（Nabeel 等，2024）的形式開(kāi)放了針對(duì)一維/二維數(shù)據(jù)集的分析代碼，見(jiàn)：
https://github.com/tee-lab/PyDaddy
（存檔于 https://doi.org/10.5281/zenodo.13777396）。

方法
隨機(jī)微分方程的數(shù)學(xué)預(yù)備知識(shí)

我們的目標(biāo)是：從對(duì)系統(tǒng)狀態(tài)的測(cè)量（表現(xiàn)為時(shí)間序列）出發(fā)，建模該系統(tǒng)的動(dòng)力學(xué)行為。典型示例包括：?jiǎn)蝹€(gè)或多個(gè)相互作用種群的種群規(guī)模、覓食動(dòng)物的運(yùn)動(dòng)軌跡、景觀中的植被覆蓋度等。為此，我們采用隨機(jī)微分方程（SDE）框架。設(shè) x ( t )
為一個(gè) d 維向量，其時(shí)間演化可建模為：

逆問(wèn)題背后的基本原理

我們感興趣的是從有限采樣頻率的觀測(cè)時(shí)間序列數(shù)據(jù)中推斷SDE的逆問(wèn)題。具體來(lái)說(shuō)，我們的目標(biāo)是找到簡(jiǎn)單、可解釋的解析表達(dá)式來(lái)描述 f 和 G，而不僅僅是它們的定性形狀。

我們解決這個(gè)問(wèn)題的方法包括兩個(gè)步驟。首先，我們從觀測(cè)到的時(shí)間序列中計(jì)算漂移和擴(kuò)散函數(shù)的瞬時(shí)估計(jì)。接下來(lái)，我們使用基于稀疏回歸的技術(shù)來(lái)估計(jì)漂移和擴(kuò)散函數(shù)的函數(shù)形式。

瞬時(shí)漂移和擴(kuò)散系數(shù)

從具有有限采樣時(shí)間 的觀測(cè)采樣時(shí)間序列中，我們可以估計(jì)瞬時(shí)漂移系數(shù)：

漂移系數(shù)與擴(kuò)散系數(shù)的方程學(xué)習(xí)
我們采用一種基于稀疏回歸的技術(shù)（有時(shí)稱為“方程學(xué)習(xí)”），以獲得所提取的漂移函數(shù) f f與擴(kuò)散函數(shù) G G 的可解釋解析表達(dá)式（Boninsegna 等，2018；Brunton 等，2016）。在無(wú)額外假設(shè)的前提下，該過(guò)程可對(duì) f f 與 G G 的每一分量分別獨(dú)立進(jìn)行（即按分量處理）。

模型選擇
除基函數(shù)外，方程發(fā)現(xiàn)算法尚需用戶提供額外輸入以約束模型選擇過(guò)程：即稀疏化閾值λ 。該參數(shù)值決定了模型復(fù)雜性與解釋能力之間的權(quán)衡——較高的 λ 會(huì)剔除更多項(xiàng)，從而生成更簡(jiǎn)潔的模型，但以擬合數(shù)據(jù)精度下降為代價(jià)。用戶可手動(dòng)設(shè)定 λ λ 以針對(duì)具體問(wèn)題優(yōu)化此權(quán)衡；亦可采用自動(dòng)調(diào)參方法。

傳統(tǒng)上可使用信息準(zhǔn)則（如赤池信息量準(zhǔn)則 AIC 或貝葉斯信息量準(zhǔn)則 BIC）進(jìn)行閾值選擇；但本文轉(zhuǎn)而采用k 折交叉驗(yàn)證（k-fold cross validation）——這是機(jī)器學(xué)習(xí)文獻(xiàn)中常用于依據(jù)交叉驗(yàn)證精度選擇理想模型的方法（Shalev-Shwartz 與 Ben-David，2014）。

其核心思想是：僅使用部分?jǐn)?shù)據(jù)（稱為訓(xùn)練集）訓(xùn)練模型，并在與訓(xùn)練集互斥的驗(yàn)證集上評(píng)估模型性能。若模型在驗(yàn)證集上表現(xiàn)良好，則說(shuō)明其未對(duì)訓(xùn)練數(shù)據(jù)中的噪聲過(guò)擬合，因而具備更強(qiáng)的泛化能力。具體而言，將數(shù)據(jù)集均分為 k 份；每次選取其中 k ? 1
份作為訓(xùn)練集擬合模型，剩余 1 份作為驗(yàn)證集，并計(jì)算模型在該驗(yàn)證集上的誤差（即驗(yàn)證誤差）；重復(fù)此過(guò)程 k 次，使每一份數(shù)據(jù)均被用作一次驗(yàn)證集，最終取平均驗(yàn)證誤差作為評(píng)估指標(biāo)。

我們選擇使交叉驗(yàn)證（CV）誤差下降幅度最大的 λ 值作為最優(yōu)閾值。對(duì)許多系統(tǒng)（包括真實(shí)世界系統(tǒng)），還可結(jié)合對(duì)該系統(tǒng)物理/生物學(xué)背景的先驗(yàn)知識(shí)，人工選定基函數(shù)并進(jìn)行手動(dòng)模型篩選，從而獲得更優(yōu)、更簡(jiǎn)潔的模型（參見(jiàn) Nabeel 等，2023 及補(bǔ)充材料 SI 第 S4 節(jié)“真實(shí)數(shù)據(jù)集的模型選擇”）。

診斷
數(shù)據(jù)驅(qū)動(dòng)的隨機(jī)微分方程（SDE）發(fā)現(xiàn)流程，若缺乏充分的診斷檢驗(yàn)以合理驗(yàn)證SDE模型的基本假設(shè)，則是不完整的。基于 Brückner 等（2019）與 Jhawar 和 Guttal（2020）的思想，我們提出了三類診斷檢驗(yàn)，用于評(píng)估所發(fā)現(xiàn)的SDE模型的可靠性：

• 噪聲診斷：方程（1）中的噪聲項(xiàng) η ( t ) 被假設(shè)為無(wú)關(guān)聯(lián)的高斯白噪聲

與傳統(tǒng)方法的對(duì)比
文獻(xiàn)中常用于恢復(fù)隨機(jī)微分方程（SDE）的一種傳統(tǒng)方法是：將變量 x x 的取值范圍劃分為有限個(gè)區(qū)間（bins），并在每個(gè)區(qū)間內(nèi)對(duì) f f 和 G G 進(jìn)行分箱平均估計(jì)。也就是說(shuō)，漂移函數(shù)與擴(kuò)散函數(shù)可由其瞬時(shí)形式（式 4、式 6）近似為：

其中， w w 為分箱寬度（bin-width），各區(qū)間為 d d 維矩形區(qū)間。該方法的缺點(diǎn)在于：每個(gè)分箱內(nèi)的估計(jì)相互獨(dú)立，無(wú)法利用任何全局結(jié)構(gòu)信息（例如函數(shù)的光滑性），因而估計(jì)結(jié)果往往較為嘈雜。一種降低估計(jì)方差的策略是對(duì)時(shí)間序列以更大的時(shí)間步長(zhǎng) Δ t
進(jìn)行子采樣（Jhawar 與 Guttal，2020），但這可能引入估計(jì)偏差（參見(jiàn)補(bǔ)充材料 SI 第 S5.A 節(jié)“有限數(shù)據(jù)下的估計(jì)”）。

相比之下，方程學(xué)習(xí)（equation learning）方法用稀疏回歸步驟替代了分箱平均步驟。除能直接輸出可解釋的解析表達(dá)式這一顯著優(yōu)勢(shì)外，方程學(xué)習(xí)還具備另外兩個(gè)優(yōu)點(diǎn)：
第一，通過(guò)合理選擇基函數(shù)庫(kù)，可顯式納入全局結(jié)構(gòu)或約束條件（如光滑性、對(duì)稱性等）；
第二，無(wú)需子采樣——從而消除了兩個(gè)任意性較強(qiáng)的參數(shù)選擇：子采樣時(shí)間步長(zhǎng) Δ t
與分箱寬度 ε 。

結(jié)果
從合成數(shù)據(jù)集中發(fā)現(xiàn)具有差異性的隨機(jī)微分方程（SDE）

由不同SDE生成的單峰分布
單峰分布廣泛見(jiàn)于諸多生物數(shù)據(jù)集中。然而，此類分布可能源于截然不同的動(dòng)力學(xué)過(guò)程，因而可由本質(zhì)迥異的底層SDE建模。我們旨在表明：即便在此類情形下，應(yīng)用本文提出的方法仍可從時(shí)間序列數(shù)據(jù)中準(zhǔn)確復(fù)原原始SDE。

我們考慮以下兩個(gè)一維隨機(jī)微分方程：

盡管方程11和12形式迥異，但由它們生成的時(shí)間序列及其直方圖卻極為相似（圖2 (A-i, A-ii, B-i, B-ii)）：事實(shí)上，x(t) 的穩(wěn)態(tài)分布均為單峰，且對(duì)于方程11和12幾乎完全相同。從動(dòng)力系統(tǒng)角度看，方程11在 x* = 0 處存在一個(gè)確定性穩(wěn)定態(tài)。在生物學(xué)上，x 可被視為種群規(guī)模偏離其穩(wěn)定承載容量的偏差，因此 x 可正可負(fù)。加性噪聲——代表環(huán)境波動(dòng)——僅使種群動(dòng)態(tài)圍繞確定性平衡點(diǎn)擴(kuò)散，這在時(shí)間序列及直方圖中均有體現(xiàn)（圖2 (A-ii)）。因此，我們可將直方圖的峰值視為反映了一個(gè)確定性狀態(tài)。這確實(shí)是噪聲最廣為人知的作用。

另一方面，對(duì)于方程12，其確定性穩(wěn)定平衡點(diǎn)位于 ±1，但直方圖的峰值卻在 0（圖2 (B-ii)）。在此情況下，狀態(tài)依賴噪聲或乘性噪聲項(xiàng)徹底改變了系統(tǒng)的穩(wěn)定性景觀，在兩個(gè)確定性穩(wěn)定態(tài)之間催生出一個(gè)新?tīng)顟B(tài)，導(dǎo)致峰值出現(xiàn)在 x = 0。這是噪聲產(chǎn)生異常效應(yīng)的一個(gè)例子，即“噪聲誘導(dǎo)的穩(wěn)定態(tài)”，該概念最初在干旱生態(tài)系統(tǒng)模型中被提出（D’Odorico 等，2005）。

針對(duì)這些單峰數(shù)據(jù)集，我們現(xiàn)在提出逆問(wèn)題：能否僅基于時(shí)間序列數(shù)據(jù)的特征（圖2 A-i 與 B-i），推斷出正確的底層SDE模型？確實(shí)，我們所提出的、結(jié)合跳躍矩計(jì)算與稀疏回歸的方法，能夠準(zhǔn)確識(shí)別出這兩個(gè)模型的漂移函數(shù)與擴(kuò)散函數(shù)的具體形式（圖2中，比較A列和B列內(nèi)的紅色線與黑色線；第iii行和第iv行）。至關(guān)重要的是，我們能夠推斷出：圖2 A-i所示的時(shí)間序列由一個(gè)線性漂移項(xiàng)（圖2 A-iii）與加性噪聲驅(qū)動(dòng)（圖2 A-iv）；而圖2 B-i所示的時(shí)間序列則由一個(gè)具有三個(gè)根的非線性漂移函數(shù)（圖2 B-iii）與乘性噪聲驅(qū)動(dòng)（圖2 B-iv）。此外，我們還能恢復(fù)出漂移函數(shù)與擴(kuò)散函數(shù)的符號(hào)化解析表達(dá)式，其結(jié)果與我們用于生成合成數(shù)據(jù)集的原始SDE高度吻合（圖2，“Estimated”面板）。

由不同SDE生成的雙峰分布

許多生物系統(tǒng)表現(xiàn)出替代穩(wěn)定態(tài)，最簡(jiǎn)單的情形是雙穩(wěn)態(tài)系統(tǒng)，例如干旱生態(tài)系統(tǒng)中的草原與林地狀態(tài)，或湖泊的富營(yíng)養(yǎng)化與貧營(yíng)養(yǎng)化狀態(tài)。具有雙穩(wěn)態(tài)的系統(tǒng)，其狀態(tài)變量通常呈現(xiàn)雙峰分布。同樣，這種雙峰性也可由本質(zhì)不同的底層過(guò)程/SDE生成。為說(shuō)明這一點(diǎn)，我們考慮以下兩個(gè)玩具模型：

方程13在 ±√(2/3) 處有兩個(gè)確定性穩(wěn)定平衡點(diǎn)，其動(dòng)力學(xué)行為由加性噪聲（圖2 (C-ii)）擴(kuò)散至這些平衡點(diǎn)周?chē)?；因此，這些是確定性狀態(tài)。另一方面，方程14僅在 x* = 0 處有一個(gè)穩(wěn)定平衡點(diǎn)；然而，由于乘性噪聲項(xiàng)的作用，系統(tǒng)在直方圖中呈現(xiàn)出兩個(gè)遠(yuǎn)離確定性平衡點(diǎn)的峰值（圖2 (D-ii)）；因此，這些是噪聲誘導(dǎo)的狀態(tài)。

我們?cè)俅吾槍?duì)這些雙峰時(shí)間序列數(shù)據(jù)集（圖2 C-i 與 D-i 對(duì)比，以及圖2 C-ii 與 D-ii 對(duì)比）提出逆問(wèn)題。在此情況下，我們的方法同樣能夠準(zhǔn)確恢復(fù)底層的SDE模型，包括其符號(hào)函數(shù)（參見(jiàn)圖2 C和D中的第iii、iv行及“Estimated”面板）。

在理論生態(tài)學(xué)經(jīng)典模型中的驗(yàn)證

我們現(xiàn)在將該方法應(yīng)用于若干理論生態(tài)學(xué)中的經(jīng)典模型進(jìn)行驗(yàn)證。我們考慮了多種單變量和雙變量模型的隨機(jī)化版本。所考慮的單變量模型包括：密度依賴種群增長(zhǎng)的邏輯斯蒂模型（不穩(wěn)定系統(tǒng)）、具有Holling III型功能反應(yīng)的種群捕撈模型（雙穩(wěn)態(tài)系統(tǒng)）（Strogatz, 2018），以及湖泊富營(yíng)養(yǎng)化模型（Carpenter 等, 1999）（雙穩(wěn)態(tài)系統(tǒng)）。所考慮的雙變量模型包括：用于種間競(jìng)爭(zhēng)的Lotka-Volterra競(jìng)爭(zhēng)模型、具有Holling II型功能反應(yīng)的非線性捕食者-獵物模型（Alonso 等, 2002），以及Van der Pol振蕩器——一個(gè)用于描述非線性振蕩的最小模型（Strogatz, 2018）。

表1總結(jié)了本分析的結(jié)果。總而言之，對(duì)于所有考慮的模型，SDE發(fā)現(xiàn)程序均能從模擬的時(shí)間序列中準(zhǔn)確地發(fā)現(xiàn)正確的模型（詳見(jiàn)補(bǔ)充信息SI第S3節(jié)——“生態(tài)學(xué)經(jīng)典模型驗(yàn)證”）。

有限數(shù)據(jù)下的不準(zhǔn)確模型發(fā)現(xiàn)

現(xiàn)實(shí)世界的數(shù)據(jù)存在諸多局限性，例如時(shí)間序列較短或采樣分辨率不足。理解此類情境如何影響我們所描述的數(shù)據(jù)發(fā)現(xiàn)協(xié)議至關(guān)重要。此處，我們考察“有限數(shù)據(jù)”的兩種特定子情況：短時(shí)間序列和長(zhǎng)采樣間隔。

為此，我們嘗試估計(jì)前文所述的方程4中的模型：? = 2x - 3x3 + (1/2)η(t)，該模型表現(xiàn)出雙穩(wěn)態(tài)動(dòng)力學(xué)。隨后，我們使用一個(gè)短時(shí)間序列（1000個(gè)高分辨率觀測(cè)值，Δt = 0.01s），其中系統(tǒng)僅探索了狀態(tài)空間的一小部分。所得模型及其模型診斷結(jié)果如補(bǔ)充信息圖S11 A、B所示。不出所料，估計(jì)的SDE與真實(shí)的基礎(chǔ)模型并不匹配。當(dāng)從具有大采樣間隔的數(shù)據(jù)中估計(jì)模型時(shí)，也會(huì)發(fā)生類似的錯(cuò)配（補(bǔ)充信息圖S11 C、D）。在此情況下，大的采樣間隔意味著關(guān)于細(xì)尺度動(dòng)力學(xué)的信息丟失，因此發(fā)現(xiàn)的模型再次偏離實(shí)際模型。然而，在這兩種情況下，我們都觀察到原始時(shí)間序列與模型模擬時(shí)間序列的自相關(guān)函數(shù)之間存在巨大差異，這將提醒用戶模型估計(jì)不準(zhǔn)確。另一個(gè)因數(shù)據(jù)有限或嘈雜而可能導(dǎo)致模型誤識(shí)別的明顯跡象是，估計(jì)系數(shù)中出現(xiàn)較大的誤差條，我們也觀察到了這一點(diǎn)。

我們?cè)谟懻摬糠衷敿?xì)探討了這些局限性問(wèn)題、使用方法時(shí)如何避免陷阱、以及如何通過(guò)診斷評(píng)估所發(fā)現(xiàn)的模型，并在補(bǔ)充信息SI第S5節(jié)“陷阱與局限性”中對(duì)這些局限性進(jìn)行了定量刻畫(huà)。

在真實(shí)數(shù)據(jù)集上的應(yīng)用
起初并不明確這些方法是否適用于復(fù)雜的生物數(shù)據(jù)集。尤其是，隨機(jī)微分方程（SDE）假設(shè)噪聲為高斯且無(wú)關(guān)聯(lián)的，而實(shí)際數(shù)據(jù)中我們往往預(yù)期存在更復(fù)雜的噪聲結(jié)構(gòu)（例如，非高斯性或具有記憶性）。因此，為驗(yàn)證其適用性與普適性，我們將其應(yīng)用于兩個(gè)性質(zhì)迥異的真實(shí)數(shù)據(jù)集：魚(yú)群聚集運(yùn)動(dòng)的群體軌跡數(shù)據(jù)（Jhawar 等，2020）與受限單細(xì)胞遷移數(shù)據(jù)（Brückner 等，2019）。

這兩個(gè)數(shù)據(jù)集在多個(gè)關(guān)鍵方面存在顯著差異：
第一，魚(yú)群數(shù)據(jù)刻畫(huà)了群體層面涌現(xiàn)的秩序性，而細(xì)胞數(shù)據(jù)反映的是個(gè)體層面行為；
第二，魚(yú)群的空間尺度相對(duì)較大（約1米），但時(shí)間動(dòng)態(tài)較快（約0.1秒量級(jí)）；而細(xì)胞數(shù)據(jù)空間尺度極小（約10??米），運(yùn)動(dòng)卻極為緩慢，時(shí)間尺度達(dá)10分鐘；
第三，其底層動(dòng)力學(xué)本質(zhì)不同——魚(yú)群處于噪聲誘導(dǎo)態(tài)（Jhawar 等，2020），而細(xì)胞運(yùn)動(dòng)主要由確定性極限環(huán)主導(dǎo)，噪聲僅起次要作用（Brückner 等，2019）。

魚(yú)群噪聲誘導(dǎo)聚集行為的SDE發(fā)現(xiàn)

近十年來(lái)，動(dòng)物集體運(yùn)動(dòng)在多個(gè)生物學(xué)領(lǐng)域均取得顯著進(jìn)展。借助前沿成像技術(shù)對(duì)運(yùn)動(dòng)群體進(jìn)行記錄，研究者已獲得高時(shí)空分辨率的數(shù)據(jù)，有助于深入探索同步集體運(yùn)動(dòng)的機(jī)制。集體運(yùn)動(dòng)的數(shù)學(xué)理論常以隨機(jī)微分方程形式表達(dá)其動(dòng)力學(xué)，為這種同步群體行為提供可能解釋（Biancalani 等，2014；Jhawar 等，2019；Ramaswamy，2017；Toner 等，2005；Vicsek 與 Zafeiris，2012）。隨著高質(zhì)量動(dòng)物運(yùn)動(dòng)數(shù)據(jù)的可及性提升，我們現(xiàn)在可提出其逆問(wèn)題：能否從給定的動(dòng)物軌跡時(shí)間序列中發(fā)現(xiàn)隨機(jī)動(dòng)力學(xué)模型？進(jìn)一步地，還可探討隨機(jī)性在塑造（或破壞）秩序中的根本作用：觀測(cè)到的集體動(dòng)力學(xué)，究竟與確定性狀態(tài)一致（即隨機(jī)性僅在確定性穩(wěn)定平衡點(diǎn)附近模糊秩序，如圖2 A-ii所示），還是與噪聲誘導(dǎo)狀態(tài)一致（即隨機(jī)性在遠(yuǎn)離確定性穩(wěn)定平衡點(diǎn)處生成非平凡狀態(tài)，如圖2 C-ii所示）？

近期一項(xiàng)研究推斷，魚(yú)群高度同步的運(yùn)動(dòng)行為是一種噪聲誘導(dǎo)態(tài)（Jhawar 等，2020）。該研究在群體動(dòng)力學(xué)時(shí)間序列上采用傳統(tǒng)跳躍矩估計(jì)方法實(shí)現(xiàn)推斷。本文中，我們檢驗(yàn)?zāi)芊裢ㄟ^(guò)結(jié)合跳躍矩計(jì)算與稀疏回歸的推斷流程，復(fù)現(xiàn)相同結(jié)果；此外，我們還執(zhí)行了該研究中被忽略的噪聲診斷與模型診斷。

我們使用（Jhawar 等，2020）公開(kāi)發(fā)布的Etroplus suratensis（卡拉米恩魚(yú)）群體數(shù)據(jù)集。實(shí)驗(yàn)中，15尾魚(yú)組成的魚(yú)群運(yùn)動(dòng)通過(guò)高分辨率攝像機(jī)記錄，并經(jīng)計(jì)算機(jī)視覺(jué)方法追蹤，生成一個(gè)30維時(shí)間序列（15尾魚(yú) × 每尾2維位置坐標(biāo)）。直接在30維空間建模動(dòng)力學(xué)，信息增益有限；相比之下，采用能捕捉群體本質(zhì)動(dòng)力學(xué)的低維表征更利于理解集體行為。

受物理學(xué)文獻(xiàn)（Toner 等，2005；Vicsek 等，1995）及其在生物集體運(yùn)動(dòng)中的廣泛應(yīng)用（Couzin 等，2002；Jhawar 等，2020；Tunstr?m 等，2013）啟發(fā)，我們聚焦于群體極化秩序的涌現(xiàn)與演化過(guò)程。

該數(shù)據(jù)集提供了極化向量 m m 的時(shí)間序列，采樣間隔為0.12秒，持續(xù)約1小時(shí)，但含若干缺失點(diǎn)（對(duì)應(yīng)追蹤失效時(shí)刻）。盡管群體極化時(shí)間序列展現(xiàn)出顯著隨機(jī)漲落（圖3B），其直方圖卻表明魚(yú)群總體處于有序狀態(tài)（圖3C、D）。此即為我們分析的輸入時(shí)間序列。

此處，我們借助所開(kāi)發(fā)的PyDaDDy軟件包（整合跳躍矩估計(jì)與稀疏回歸）復(fù)現(xiàn)了上述結(jié)果（另見(jiàn)補(bǔ)充材料 SI 第 S4.A 節(jié)“魚(yú)群數(shù)據(jù)集的模型選擇”）。所發(fā)現(xiàn)的方程包含線性漂移項(xiàng)與二次型擴(kuò)散項(xiàng)，對(duì)應(yīng)的向量隨機(jī)微分方程形式為：

圖4展示了對(duì)該發(fā)現(xiàn)的SDE模型進(jìn)行診斷測(cè)試的結(jié)果。噪聲殘差 r ( t )
呈高斯分布，符合預(yù)期（圖4A）。殘差噪聲 r ( t )中的相關(guān)性迅速衰減（圖4B）。這些測(cè)試在合理程度上支持了關(guān)于 η 為高斯白噪聲的建模假設(shè)。

我們還發(fā)現(xiàn)，所發(fā)現(xiàn)的方程通過(guò)了模型診斷測(cè)試。由方程15模擬生成的SDE的直方圖與原始時(shí)間序列中 m m 的直方圖高度吻合（圖4C）。模擬時(shí)間序列的自相關(guān)函數(shù)也與原始時(shí)間序列表現(xiàn)出合理的吻合（圖4D），但有一個(gè)顯著偏差：數(shù)據(jù)自相關(guān)函數(shù)在約10秒的相對(duì)較長(zhǎng)的時(shí)間尺度上呈現(xiàn)負(fù)值，之后才收斂至零。魚(yú)群研究論文的作者們已證明，這一特征源于邊界效應(yīng)，對(duì)于 Etroplus suratensis 的魚(yú)群動(dòng)力學(xué)而言并不重要。盡管如此，我們強(qiáng)調(diào)，本研究所用的SDE發(fā)現(xiàn)協(xié)議僅利用了極小時(shí)間尺度（約0.1秒）下群體極化波動(dòng)的高度局部信息；它并未使用任何關(guān)于群體極化數(shù)據(jù)頻率分布或自相關(guān)函數(shù)的信息。盡管如此，模擬時(shí)間序列在這兩個(gè)指標(biāo)上仍與數(shù)據(jù)表現(xiàn)出良好的一致性。最后但同樣重要的是，我們發(fā)現(xiàn)該模型具有自洽性——當(dāng)我們將模擬數(shù)據(jù)再次輸入我們的方程發(fā)現(xiàn)協(xié)議時(shí)，能夠恢復(fù)出相同的SDE。

受限細(xì)胞遷移動(dòng)力學(xué)的SDE發(fā)現(xiàn)

從形態(tài)發(fā)生、傷口愈合到癌癥轉(zhuǎn)移，細(xì)胞遷移在眾多生物學(xué)情境中扮演著關(guān)鍵角色，理解細(xì)胞如何在復(fù)雜環(huán)境中移動(dòng)至關(guān)重要。近期一項(xiàng)研究（Brückner 等，2019）探討了結(jié)構(gòu)化環(huán)境中細(xì)胞遷移的隨機(jī)動(dòng)力學(xué)。作者設(shè)計(jì)了一項(xiàng)實(shí)驗(yàn)：一個(gè)癌細(xì)胞在兩個(gè)“島嶼”（即“狀態(tài)”）之間來(lái)回遷移（圖5A），并使用SDE框架對(duì)該二態(tài)遷移的動(dòng)力學(xué)進(jìn)行建模，但未采用顯式表達(dá)式。他們得出結(jié)論：該動(dòng)力學(xué)主要由確定性（漂移）分量中的極限環(huán)主導(dǎo)，而隨機(jī)分量?jī)H起次要作用，影響的是狀態(tài)間轉(zhuǎn)換的時(shí)間尺度。本文中，我們利用所提出的SDE發(fā)現(xiàn)協(xié)議復(fù)現(xiàn)了該研究的關(guān)鍵結(jié)果。雖然原研究通過(guò)分箱平均、分段表示的方法構(gòu)建漂移與擴(kuò)散函數(shù)，我們的方法則能夠提取出可解釋的函數(shù)形式。此外，我們還展示了噪聲診斷與模型診斷如何引導(dǎo)更精確的模型發(fā)現(xiàn)。

該數(shù)據(jù)集包含149條獨(dú)立的細(xì)胞軌跡重復(fù)實(shí)驗(yàn)。每條軌跡基于高分辨率圖像，每10分鐘采集一次，最長(zhǎng)持續(xù)50小時(shí)。圖5B展示了細(xì)胞在兩個(gè)島嶼間遷移的一個(gè)時(shí)間序列示例，圖5C-E則描繪了狀態(tài)變量（位置x和速度v）的直方圖。

根據(jù)實(shí)驗(yàn)設(shè)計(jì)（圖5A）以及細(xì)胞位置的雙峰直方圖（圖5C），一個(gè)自然的初步嘗試是將細(xì)胞在兩個(gè)穩(wěn)定狀態(tài)（即著陸墊）之間的跳躍建模為純粹隨機(jī)的切換過(guò)程。這類似于圖2C中描述的模型，其特征是在兩個(gè)穩(wěn)定態(tài)之間發(fā)生隨機(jī)躍遷。這需要將細(xì)胞軌跡 x(t) 建模為如下形式的過(guò)阻尼SDE：

然而，我們的模型診斷結(jié)果表明，過(guò)阻尼模型不足以完整刻畫(huà)該系統(tǒng)的動(dòng)力學(xué)行為（參見(jiàn)補(bǔ)充材料 SI 第 S4.B 節(jié)“細(xì)胞遷移數(shù)據(jù)集的模型選擇”）。

因此，我們轉(zhuǎn)而采用欠阻尼模型，引入兩個(gè)動(dòng)力學(xué)變量——位置（ x ）——對(duì)該細(xì)胞軌跡進(jìn)行建模。此類二維描述能夠捕捉系統(tǒng)中可能存在的線性或非線性振蕩行為；如下文所示，該模型確實(shí)最符合實(shí)際數(shù)據(jù)。

該欠阻尼隨機(jī)動(dòng)力學(xué)遵循如下方程：

其中，我們?cè)趹?yīng)用SDE發(fā)現(xiàn)流程前，對(duì)位置 x x 和速度 v v 進(jìn)行了無(wú)量綱化縮放。
此縮放步驟對(duì)稀疏回歸的正確運(yùn)行是必要的（參見(jiàn)補(bǔ)充材料 SI 第 S4.B 節(jié)）。

通過(guò)我們的SDE發(fā)現(xiàn)流程，我們發(fā)現(xiàn)漂移函數(shù) f f（圖5E）可很好地用如下形式的三次多項(xiàng)式近似：

我們通過(guò)數(shù)據(jù)發(fā)現(xiàn)的細(xì)胞遷移SDE模型呈現(xiàn)出可激發(fā)流（excitable flow），其定性行為與 Brückner 等（2019）的研究結(jié)果一致，表現(xiàn)為弛豫振蕩（relaxation oscillations）。此外，我們發(fā)現(xiàn)其漂移函數(shù)偏離了經(jīng)典的范德波爾（Van der Pol）振蕩器模型——后者是解釋弛豫振蕩的最簡(jiǎn)模型。有趣的是，如補(bǔ)充材料第 S4.B 節(jié)所示，我們發(fā)現(xiàn)擴(kuò)散項(xiàng)（diffusion term）對(duì)準(zhǔn)確刻畫(huà)細(xì)胞動(dòng)力學(xué)中的有限邊界效應(yīng)、并滿足模型診斷要求至關(guān)重要（見(jiàn)下文）。

圖6展示了所發(fā)現(xiàn)模型的診斷結(jié)果：噪聲殘差在時(shí)間上無(wú)自相關(guān)性（圖6B），但其分布呈現(xiàn)比高斯分布更強(qiáng)的中心聚集性（即更尖銳的峰度，圖6A）。盡管如此，該模型生成的模擬時(shí)間序列在統(tǒng)計(jì)特性上（如狀態(tài)變量直方圖與自相關(guān)函數(shù)）與原始數(shù)據(jù)總體吻合良好（圖6C–F）。此外，該模型具有自洽性（self-consistency）。

此處我們強(qiáng)調(diào)模型診斷在準(zhǔn)確發(fā)現(xiàn)擴(kuò)散項(xiàng)過(guò)程中所起的關(guān)鍵作用：若采用加性噪聲（即常數(shù)擴(kuò)散項(xiàng)）配合相同漂移函數(shù)，模擬軌跡將頻繁突破實(shí)驗(yàn)設(shè)定的物理邊界，導(dǎo)致模擬狀態(tài)變量直方圖與真實(shí)數(shù)據(jù)存在巨大偏差。相比之下，引入四階擴(kuò)散項(xiàng)可顯著減少此類邊界違反行為，從而提升模型自洽性。因此，我們認(rèn)為該擴(kuò)散項(xiàng)實(shí)質(zhì)上有效編碼了空間約束。

所發(fā)現(xiàn)模型的泛化能力驗(yàn)證

為驗(yàn)證SDE發(fā)現(xiàn)方法能否獲得可泛化模型，我們針對(duì)魚(yú)群與細(xì)胞遷移兩個(gè)數(shù)據(jù)集，均采用留出數(shù)據(jù)驗(yàn)證（見(jiàn)“診斷”一節(jié)中“使用留出數(shù)據(jù)進(jìn)行驗(yàn)證”子節(jié)）：

• 對(duì)魚(yú)群數(shù)據(jù)集，將4次實(shí)驗(yàn)中的2次作為訓(xùn)練集；
• 對(duì)細(xì)胞遷移數(shù)據(jù)集，將149條軌跡中的75條用于訓(xùn)練。

在訓(xùn)練集上完成模型發(fā)現(xiàn)后，我們分別計(jì)算：（i）由所發(fā)現(xiàn)模型生成的時(shí)間序列的直方圖與自相關(guān)函數(shù)；（ii）留出驗(yàn)證集的對(duì)應(yīng)統(tǒng)計(jì)量。

結(jié)果表明，對(duì)魚(yú)群與細(xì)胞兩個(gè)數(shù)據(jù)集，模型生成數(shù)據(jù)的統(tǒng)計(jì)特性均與驗(yàn)證集高度匹配，說(shuō)明所發(fā)現(xiàn)的SDE成功捕獲了底層動(dòng)力學(xué)的可泛化特征（詳見(jiàn)補(bǔ)充材料第 S4.C 節(jié)“所發(fā)現(xiàn)SDE模型的泛化能力”）。

本文提出了一種方法：以經(jīng)驗(yàn)觀測(cè)的時(shí)間序列作為輸入，發(fā)現(xiàn)具備解析可解釋性的數(shù)據(jù)驅(qū)動(dòng)型隨機(jī)動(dòng)力學(xué)方程。為實(shí)現(xiàn)這一目標(biāo)，我們將傳統(tǒng)的跳躍矩法（用于漂移與擴(kuò)散估計(jì)；Tabar, 2019；Van Kampen, 1992）與基于稀疏回歸的方程發(fā)現(xiàn)技術(shù)（Boninsegna 等，2018；Brunton 等，2016）相結(jié)合。尤為重要的是，我們強(qiáng)調(diào)了對(duì)數(shù)據(jù)驅(qū)動(dòng)模型所依賴假設(shè)及其發(fā)現(xiàn)結(jié)果有效性進(jìn)行診斷檢驗(yàn)的必要性。除在已知真實(shí)方程的合成數(shù)據(jù)集中驗(yàn)證該方法外，我們還通過(guò)兩個(gè)性質(zhì)迥異的生物時(shí)間序列數(shù)據(jù)集（魚(yú)群聚集與細(xì)胞遷移）展示了其普適性。為便于更廣泛的生物學(xué)研究者使用，我們開(kāi)發(fā)并開(kāi)源了一個(gè)易于使用的 Python 軟件包——PyDaDDy（https://pydaddy.readthedocs.io/）。

然而，通過(guò)數(shù)據(jù)驅(qū)動(dòng)的方程發(fā)現(xiàn)，我們揭示出：

• 魚(yú)群同步行為實(shí)為一種噪聲誘導(dǎo)態(tài)，其穩(wěn)定狀態(tài)遠(yuǎn)離確定性穩(wěn)定平衡點(diǎn)；
• 細(xì)胞遷移的振蕩則源于確定性極限環(huán)，隨機(jī)性僅起次要作用；
• 且該極限環(huán)的數(shù)學(xué)結(jié)構(gòu)不同于經(jīng)典的范德波爾振子。

換言之，在用戶輸入極少干預(yù)的前提下，我們通過(guò)符號(hào)化表達(dá)的數(shù)據(jù)驅(qū)動(dòng)動(dòng)力學(xué)模型，成功排除了若干競(jìng)爭(zhēng)性假說(shuō)，逼近更真實(shí)的機(jī)制。

我們強(qiáng)調(diào)，該方法具有極高的計(jì)算效率。例如，針對(duì)軟件包中提供的魚(yú)群數(shù)據(jù)集（二維向量時(shí)間序列，約25,000個(gè)時(shí)間點(diǎn)）：

• 在一臺(tái)2023款Mac Mini（Apple M2處理器，16GB內(nèi)存）上，漂移與擴(kuò)散系數(shù)的初始估計(jì)僅需約800毫秒，漂移函數(shù)擬合約2毫秒，擴(kuò)散函數(shù)擬合約6毫秒；
• 在Google Colab實(shí)例上，對(duì)應(yīng)耗時(shí)分別約為2.5秒、10毫秒與60毫秒。

陷阱及其規(guī)避方法

真實(shí)生物系統(tǒng)高度復(fù)雜，其“真實(shí)”底層動(dòng)力學(xué)可能涉及數(shù)十乃至上百個(gè)相互作用變量。此外，對(duì)此類系統(tǒng)的觀測(cè)數(shù)據(jù)常具有以下特征：不完整（僅觀測(cè)到部分動(dòng)力學(xué)變量）、不精確（含測(cè)量噪聲，或采樣稀疏、時(shí)間間隔大）、缺失（若干時(shí)間點(diǎn)數(shù)據(jù)缺失），或代表性不足（數(shù)據(jù)僅覆蓋狀態(tài)空間的一小部分）。在此情形下，是否仍可用僅含少數(shù)變量的隨機(jī)微分方程（SDE）建模，尚不明確。這就引出一個(gè)關(guān)鍵問(wèn)題：如何確保從經(jīng)驗(yàn)數(shù)據(jù)中獲得的數(shù)據(jù)驅(qū)動(dòng)SDE模型確實(shí)是系統(tǒng)動(dòng)力學(xué)的有效表征？若盲目將SDE發(fā)現(xiàn)流程應(yīng)用于任意數(shù)據(jù)集，流程雖總能輸出某個(gè)SDE，但該模型是否有效、是否忠實(shí)反映系統(tǒng)動(dòng)力學(xué)，卻無(wú)法保證。因此，必須審慎評(píng)估SDE發(fā)現(xiàn)技術(shù)對(duì)特定問(wèn)題或數(shù)據(jù)集的適用性。

本文方法論（及配套軟件包PyDaDDy）的核心理念之一，正是為用戶提供模型評(píng)估與診斷工具，輔助其作出合理判斷。本節(jié)及補(bǔ)充材料SI第S5節(jié)將討論在用低維SDE建模生物系統(tǒng)時(shí)常見(jiàn)的陷阱、識(shí)別方法及（可能的）修正策略。

系統(tǒng)維度問(wèn)題

原則上，數(shù)據(jù)驅(qū)動(dòng)方程發(fā)現(xiàn)方法適用于任意維度。本文聚焦于一維與二維數(shù)據(jù)集，旨在闡明隨機(jī)動(dòng)力學(xué)方程推斷的核心原理與應(yīng)用。為展示方法普適性，我們提供了一個(gè)恢復(fù)含隨機(jī)性的三維經(jīng)典混沌系統(tǒng)——洛倫茲方程的示例代碼（見(jiàn)筆記本：Higher dimensions）。

然而高維應(yīng)用存在若干實(shí)際挑戰(zhàn)：第一，維度增加導(dǎo)致基函數(shù)庫(kù)規(guī)模急劇膨脹，估計(jì)易受誤差干擾；稀疏回歸作為一種正則化手段，可在一定程度上緩解該問(wèn)題；第二，結(jié)合對(duì)系統(tǒng)物理、約束與對(duì)稱性的先驗(yàn)知識(shí)，可精心構(gòu)建小規(guī)?；瘮?shù)庫(kù)（Nabeel 等，2023）；第三，高維模型的驗(yàn)證本身即具挑戰(zhàn)性——需發(fā)展不依賴可視化（一/二維可行）的定量統(tǒng)計(jì)診斷方法，此為未來(lái)研究方向。

需強(qiáng)調(diào)的是：當(dāng)目標(biāo)是構(gòu)建基于動(dòng)力學(xué)變量的可解釋方程時(shí)，天然具備低維表征能力的系統(tǒng)尤為適合數(shù)據(jù)發(fā)現(xiàn)方法。例如，即便原始數(shù)據(jù)高維（如魚(yú)群數(shù)據(jù)為 2 N 2N 維， N N 為魚(yú)數(shù)），仍常可識(shí)別出少數(shù)隱變量（latent dimensions）以刻畫(huà)系統(tǒng)本質(zhì)動(dòng)力學(xué)。此類降維不僅方便分析，更具理論價(jià)值——諸多理論本就建立于對(duì)高維生物系統(tǒng)的粗?；枋鲋希―urrett & Levin, 1994），此思路源于物理與應(yīng)用數(shù)學(xué)。隱變量可基于粗?；碚摶蛏飳W(xué)意義量選定（如魚(yú)群研究中受介觀集體運(yùn)動(dòng)理論啟發(fā)，選用二維極化向量）；亦可通過(guò)數(shù)據(jù)驅(qū)動(dòng)技術(shù)識(shí)別（Greenacre 等，2022；Schmid，2022；Van Der Maaten 等，2009）。當(dāng)前已有工作探索在同時(shí)發(fā)現(xiàn)合適潛空間與該空間內(nèi)動(dòng)力學(xué)方程（針對(duì)確定性系統(tǒng)，Champion 等，2019）；將此類方法拓展至隨機(jī)系統(tǒng)，是值得推進(jìn)的方向。

然而，某些情形下系統(tǒng)動(dòng)力學(xué)本質(zhì)高維，或觀測(cè)變量不足以完備刻畫(huà)動(dòng)力學(xué)。此時(shí)，僅基于觀測(cè)變量構(gòu)建的SDE（或任何低維模型）必?zé)o法完全捕獲系統(tǒng)行為。對(duì)此，PyDaDDy 的診斷工具至關(guān)重要——可輔助判斷模型是否誤用。具體而言：

• 若存在未觀測(cè)變量，殘差自相關(guān)衰減時(shí)間尺度可能顯著延長(zhǎng)，提示SDE模型未能解釋的慢動(dòng)力學(xué)殘余；
• 模擬數(shù)據(jù)與真實(shí)數(shù)據(jù)的直方圖或自相關(guān)函數(shù)若無(wú)法匹配，亦表明模型不充分（如以過(guò)阻尼SDE建模細(xì)胞遷移時(shí)失效，見(jiàn)SI S4.B）；
• 另一案例：僅用單物種時(shí)間序列擬合兩物種相互作用模型，診斷測(cè)試即揭示降維模型不完備（見(jiàn)SI S5.B）。

與所有估計(jì)方法一樣，數(shù)據(jù)驅(qū)動(dòng)SDE發(fā)現(xiàn)亦需質(zhì)量足夠、數(shù)量充分的數(shù)據(jù)。分析表明：時(shí)間序列過(guò)短或采樣頻率過(guò)低均會(huì)導(dǎo)致SDE估計(jì)失準(zhǔn)。值得慶幸的是，此類情形下模型診斷通常表現(xiàn)不佳，可警示用戶謹(jǐn)慎解讀結(jié)果?？傮w而言，重建所需數(shù)據(jù)量取決于底層模型復(fù)雜度與測(cè)量噪聲水平（Fajardo-Fontiveros 等，2023）。盡管難以給出普適數(shù)據(jù)量閾值，我們?cè)赟I S5.A節(jié)對(duì)PyDaDDy在有限數(shù)據(jù)下的性能進(jìn)行了數(shù)值探索，發(fā)現(xiàn)：方程發(fā)現(xiàn)技術(shù)在數(shù)據(jù)受限時(shí)，顯著優(yōu)于傳統(tǒng)的分箱Kramers–Moyal平均法。

關(guān)于測(cè)量噪聲：

• 當(dāng)測(cè)量噪聲較小時(shí)，方法仍可有效工作；漂移函數(shù)估計(jì)在平均意義上基本不受影響；
• 但擴(kuò)散函數(shù)估計(jì)會(huì)系統(tǒng)性偏高，偏移量約等于測(cè)量噪聲方差（B?ttcher 等，2006）——若目標(biāo)為構(gòu)建預(yù)測(cè)模型，此偏差或可接受；
• 若目標(biāo)為構(gòu)建機(jī)理模型，則需顯式分離測(cè)量噪聲與生物本征隨機(jī)性。更廣泛地，不同噪聲結(jié)構(gòu)（變量間相關(guān)、時(shí)間相關(guān)、測(cè)量噪聲）如何影響SDE推斷，是值得深入研究的開(kāi)放問(wèn)題。

最后，時(shí)間序列中偶發(fā)缺失數(shù)據(jù)點(diǎn)（如傳感器故障所致）影響較小——SDE估計(jì)可直接基于可用數(shù)據(jù)點(diǎn)進(jìn)行，如魚(yú)群數(shù)據(jù)集所示。

相關(guān)工作與結(jié)論性評(píng)述

方程學(xué)習(xí)（equation learning）最初在工程學(xué)文獻(xiàn)中提出，用于常微分方程與偏微分方程的發(fā)現(xiàn)，并已被集成于如PySINDy（de Silva 等，2020）和DataDrivenDiffEq.jl（JuliusMartensen 等，2021）等流行工具包中。然而，這些方法局限于確定性微分方程的發(fā)現(xiàn)。雖有少數(shù)例外，例如 R 包Langevin（Rinn 等，2016），其可接受一維或二維時(shí)間序列輸入，并以分箱平均方式估計(jì)克雷默–莫亞爾（Kramers–Moyal）系數(shù)來(lái)獲得漂移與擴(kuò)散項(xiàng)，且提供部分可視化與診斷功能，但其未輸出可解釋的SDE解析表達(dá)式；而且，分箱平均估計(jì)本身易受誤差影響（Callaham 等，2021；Jhawar 與 Guttal，2020），尤其在數(shù)據(jù)有限情形下（見(jiàn)補(bǔ)充材料 SI 第 S5.A 節(jié)）——而本文所采用的方程學(xué)習(xí)方法極大緩解了這一問(wèn)題。

我們強(qiáng)調(diào)：絕大多數(shù)數(shù)據(jù)驅(qū)動(dòng)建模工具（包括在動(dòng)力系統(tǒng)建模領(lǐng)域廣泛應(yīng)用的DiffEqParamEstim.jl（Rackauckas 與 Nie，2017）與Sim.DiffProc（Guidoum 與 Boukhetala，2020）），均采用傳統(tǒng)范式——即要求用戶預(yù)先指定漂移與擴(kuò)散函數(shù)的參數(shù)化形式。在此類方法中，用戶需先提出一組候選模型，再通過(guò)模型選擇程序確定最適動(dòng)力系統(tǒng)。典型地，隨機(jī)動(dòng)力系統(tǒng)模型的選擇準(zhǔn)則依賴于穩(wěn)態(tài)分布與實(shí)測(cè)數(shù)據(jù)的匹配程度。然而，此類方法無(wú)法識(shí)別圖2 A-i 與 B-i 所示情形：兩個(gè)時(shí)間序列具有幾乎完全相同的穩(wěn)態(tài)分布，卻源于本質(zhì)不同的動(dòng)力學(xué)機(jī)制。相比之下，我們的方法不僅可直接從時(shí)間序列中發(fā)現(xiàn)SDE，更能準(zhǔn)確區(qū)分本質(zhì)不同的控制方程。

當(dāng)前，生物學(xué)領(lǐng)域正快速積累跨尺度的大型數(shù)據(jù)集；隨著大數(shù)據(jù)時(shí)代的來(lái)臨，我們期望：本文提出的這一方法，配合相對(duì)易用的開(kāi)源工具包PyDaDDy，將激勵(lì)研究者更廣泛地將其應(yīng)用于從數(shù)據(jù)中刻畫(huà)控制動(dòng)力學(xué)方程的任務(wù)中。借此，我們可在充分考慮系統(tǒng)本征隨機(jī)性的同時(shí)，挖掘具機(jī)理基礎(chǔ)的預(yù)測(cè)模型的潛力。

當(dāng)然，將此類方法拓展至更復(fù)雜的現(xiàn)實(shí)世界數(shù)據(jù)集（如高維系統(tǒng)、含時(shí)變參數(shù)、受觀測(cè)噪聲干擾、具有復(fù)雜噪聲結(jié)構(gòu)——包括有色噪聲與非高斯噪聲等），仍將面臨諸多挑戰(zhàn)；而這些挑戰(zhàn)，也正是未來(lái)研究的重要機(jī)遇所在。

原文鏈接：https://arxiv.org/pdf/2205.02645v6