系統生物學中的隨機微分方程數值仿真算法

牛原玲 陳琳 陳洛南

(1.中南大學數學與統計學院,長沙,410083?2.江西財經大學統計與數據科學學院,南昌,330013?3.中國科學院系統生物學重點實驗室,分子細胞科學卓越創新中心(上海生物化學與細胞生物學研究所),上海,200031)

1 引言

隨機微分方程(stochastic differential equation)常常用來描述一些自然社會現象中事物的發展規律,由于考慮了隨機因素的干擾,它往往能比確定性的微分方程更為準確地刻畫變量隨時間的演化規律.隨機微分方程的本質是隨機積分方程,根據所用隨機積分的不同,主要分為It? 型和Stratnovich 型,兩種類型的隨機微分方程是可以相互轉換的.本文我們只關注如下的It? 型隨機微分方程:

其中t∈[0,∞),X(0)=X0,W(t)=(W1(t),W2(t),···,Wd(t))T是d維維納過程,又稱為d維布朗運動.f:R+×Rn→Rn和g:R+×Rn→Rn×d分別稱為漂移項系數和擴散項系數,它們均為[0,∞)上的Borel 可測函數.

如果隨機擾動出現在微分方程的系數中,如

其中ηt為一個隨機過程或者隨機變量,這樣的方程通常叫做帶有隨機輸入的微分方程(random ordinary differential equation).本文不涉及此類微分方程的數值仿真.

作為一個重要的建模工具,隨機微分方程已經在生物、經濟、控制等領域得到廣泛的應用.但是對大部分隨機微分方程而言,一般不能得到其真解.有些隨機微分方程的真解雖然可以得到,但由于形式過于復雜,實際應用中不太方便處理.因此,研究隨機微分方程的數值算法就顯得尤為重要.確定性微分方程的數值算法是一個已經得到長足發展的領域[1–4].近幾十年來,國內外對隨機微分方程數值計算的研究也取得了快速的發展,參見[5–20].常用的求解隨機微分方程的數值方法有Euler-Maruyama 方法、Milstein 方法和Runge-Kutta 方法等.收斂性是數值方法的一個重要評價指標.要求數值方法在一定意義下是收斂的,并且希望其具有較高的收斂階以提高仿真效率.與確定性的情形不同,隨機微分方程數值方法的收斂性有多種定義,其中人們最為關注的分別是強收斂性和弱收斂性.強收斂性要求數值近似的軌道充分接近真實軌道,而弱收斂性關注的則是解過程在矩意義下的收斂性.本文主要涉及前者.此外,我們還希望數值方法能夠保持原系統的一些重要特征,比如保區域,保結構等.

本文主要回顧系統生物學中的隨機微分方程模型及相關的數值仿真算法.這些隨機微分方程一般具有高維、高度非線性、真解位于某個特定的區域等特點,因此對它們的數值模擬需要做專門的研究.鑒于篇幅有限,我們僅從生物化學(生化)反應系統、生態系統、傳染病模型、群體遺傳學模型以及細胞分化模型的數值仿真來加以闡述.

需要說明的是,目前成熟的商業軟件和開源軟件中也存在一些可以用于數值求解隨機微分方程的軟件包:如MATLAB 中的“SDETools”和“SDELab”及R 中的“sde”和“yuima”.其中“SDETools”主要采用的數值格式為Euler–Maruyama 格式和Milstein 格式(要求隨機微分方程具有對角噪聲)?“SDELab”則主要采用θ-Maruyama 和Milstein 格式(不要求隨機微分方程具有對角噪聲)?“sde”僅能數值求解一維的隨機微分方程,主要采用的數值格式為Euler–Maruyama 格式,Milstein 格式和Predictor-Corrector 格式?“yuima”可以數值求解多維的隨機微分方程(可以是分數階噪聲驅動,也可以存在Levy 跳),采用的數值格式為Euler–Maruyama 格式.這些軟件包中采用的數值格式均較為簡單,應用于本文提到的這些復雜模型時效果往往不太理想.

2 生化反應系統的數值仿真

細胞中發生的大多數生化反應本質上是隨機的.由于傳統實驗成本昂貴,生化反應系統的隨機模擬近年來受到越來越多的關注.假設有一個攪拌良好的生化反應系統,包含N種分子S1,S2,···,SN,這些分子通過多個反應通道在定容恒溫的環境下發生反應.用xi(t)(i=1,···,N)表示t時刻Si的數量,則該系統的狀態可用向量X=(x1(t),x2(t),···,xN(t))T來表示.用矩陣v表示化學計量矩陣,其元素vij表示第j個反應Rj導致分子Si數量的變化.

2.1 化學主方程模型

化學主方程(Chemical master equation,簡記為CME)描述了生化反應系統在不同時間點處于每種狀態的概率:

其中X(t)是狀態變量,aj(X(t))是傾向函數,aj(X(t))dt表示在下一個小時間段[t,t+dt)第j個反應Rj發生的概率.

化學主方程似乎是發現隨機系統時間演化規律的理想方法,它可以返回系統的概率密度函數.然而,化學主方程只有在非常簡單的情況下可以求出真解.也就是說,在大多數情況下,需要借助數值方法來對其進行數值模擬.此外,由于必須在每個時間點對每個可能的狀態進行評估,并且狀態的數量可以呈指數增長,因此計算量非常大,這使得化學主方程模型最初只能有效地用于短期模擬小規?；蛘吆唵蜗到y.近年來,雖然化學主方程的數值求解方法取得了長足的進步[21,22,23],但當分子數和傾向函數太大以致于無法使用化學主方程計算時,必須使用新的數值模擬方法,如Gillespie型算法.

2.2 離散狀態馬爾可夫鏈模型

隨機模擬算法(stochastic simulation algorithm,簡記為SSA)由Gillespie 于1977 年提出[24],該方法可以精確模擬化學主方程給出的完整分布中的各條路徑,是一種離散狀態馬爾可夫鏈方法.由無限多個SSA 模擬建立的概率密度函數將與化學主方程給出的系統真實分布相同.無限次模擬顯然是無法達到的,不過適度多的重復次數可以獲得令人滿意的近似概率密度函數.盡管近年來Gillespie 型方法取得了不小的進展[25–30],但這些算法針對復雜的生化反應系統的計算成本仍然會非常高,尤其是在反應速率和分子數范圍很廣的情況下.

2.3 郎之萬方程模型

在描述生化反應系統時,化學主方程在零分子數時有一個自然邊界條件.Fokker-Planck 方程是對化學主方程的近似[31],近似的過程中可能丟失了一些信息,此時,人們需要給Fokker-Planck方程加上合適的邊界條件,如Neumann 邊界條件[32],反射邊界條件[33],以更好地描述生化反應過程.利用Feynman-Kac 公式,某些偏微分方程的解可以用對應的隨機微分方程的解來表示.與Fokker-Planck 方程數學上等價的隨機微分方程叫做郎之萬方程.該模型是離散狀態馬爾可夫鏈方法的一種有吸引力的替代方法,因為它可以顯著提高計算速度.然而,隨機微分方程模型的解可能會變成負數,甚至是虛數.顯然,負數甚至是虛數作為分子的數目是不合理的.此時,人們在做數值模擬時通常的做法是,當近似值變為負值時,將其設置為0,或者對維納增量重新采樣,直到數值解為非負值.但這些措施會給數值解帶來偏差.文獻[34]探索了一種修改化學郎之萬方程噪聲項的方法,以使變量保持非負.他們只允許物種參與的反應步驟具有隨機性.但是,文獻[34]在去掉乘積噪聲項后得到的改進的保持非負性的郎之萬方程不能保證協方差矩陣的匹配條件.文獻[35]研究了漂移項和擴散項均為隱式的可保非負性的Milstein 方法,用于數值求解化學郎之萬方程.該格式在牛頓迭代中施加了取值非負的約束,計算過程中需要刪除因反應物耗盡而停止的反應通道Rj,即刪除化學計量矩陣v的第j列vj和相應的傾向函數aj,直到反應物濃度再次升高,再將刪除的反應通道Rj恢復,同時修改化學計量矩陣和傾向函數.這些措施使得算法的“保非負性”得以實現.該方案始終滿足協方差矩陣匹配條件.但是因為全隱式方法需要額外的計算,并且計算過程中需要實時修改化學計量矩陣和傾向函數,所以計算成本可能非常高.文獻[36]表明,通過將郎之萬方程的定義域擴展到復空間,得到復數域上的郎之萬方程,該方程得到的反應物的平均濃度落在實數域.但是此時分子數是復數,在生物學上不合理.

當使用化學郎之萬方程通過連續過程近似化學主方程描述的離散過程時,目前尚不太清楚邊界條件的行為.文獻[37]中,作者通過將邊界條件納入隨機微分方程,得到了確保隨機微分方程模型的生物學上合理的解決方案.由此產生的模型稱為帶反射邊界的隨機微分方程模型.此類方程之前曾用于模擬人類代謝過程[38],離子通道動力學[39]等.

帶邊界條件的生化反應模型如下:

其中W(t)為維納過程.ξ(t)有局部的有界變差,ξ(0)=0.|ξ|(t)非負遞增連續并且滿足(2.3)–(2.4)式.如果X(s)∈?D,那么,n(s)∈N(X(s)),其中N(y)表示當y在邊界上時,指向區域內部的單位向量.(2.3)式表明當且僅當X(t)∈?D時,|ξ|(t)會增長.當X(t)不落在邊界上時,dξ(t)消失,此時帶反射邊界的隨機微分方程(2.2)退化為普通的隨機微分方程

由(2.2)–(2.4)建立的帶反射邊界的隨機微分方程模型似乎是建立生化反應系統模型的理想方法,它的解具有生物意義.但是在大多數情況下,此類隨機微分方程并沒有解析解,因此數值近似是獲得解并研究其性質的重要途徑.這里的定義域D是凸的,因此大多數現有的數值算法可以應用于此處提出的帶反射邊界的隨機微分方程模型的數值模擬.投影法是Euler-Maruyama 方法的簡單擴展,可用于求解帶反射邊界的隨機微分方程:如果Euler-Maruyama 方法獲得的數值解位于D的閉包之外,則將下一時間步的解設置為該點在閉包D上的投影,否則該步的數值逼近等于Euler-Maruyama 方法得到的數值解.數值實驗表明,利用(2.2)–(2.4)確定的帶反射邊界的隨機微分方程模型得到的數值模擬結果與SSA 得到的數值模擬結果非常接近,但是時間成本則大大降低.

Goldbeter-Koshland 開關是由兩種具有相反作用的酶E1 和E2 促進的共價修飾生化反應系統,它由以下6 個反應組成:

取初始狀態(110,100,30,30,100,30)T及參數k1=0.05,k2=0.1,k3=0.1,k4=0.01,k5=0.1,k6=0.1,文獻[37]對上述生化反應系統進行了數值模擬,結果如圖1所示.由圖1可以看到,帶反射邊界的隨機微分方程模型與SSA 得到的數值模擬結果非常接近.而SSA 所花的CPU 時間為51h55min40s,帶反射邊界的隨機微分方程的數值模擬所花的CPU 時間則為3h37min47s,可見后者的數值模擬效率大為提高.

圖1 SSA 與帶反射邊界的隨機微分方程對Goldbeter-Koshland 開關數值模擬的均值和方差(實線:SSA,虛線:帶反射邊界的隨機微分方程)[37]

這一工作的主要貢獻在于提供了一種計算高效的方法來模擬生化反應系統.該方法不僅可以保持離散狀態馬爾可夫模型的隨機行為,并且可以保證真解及數值解都有生物意義,從而使得針對大規模生化反應系統的數值仿真成為可能.值得注意的是,即使帶反射邊界的隨機微分方程模型在許多情況下表現得足夠好,但它并不太適合所有反應通道的傾向函數都非常小的生化反應系統.這種情況下需要用到2.2 節提到的Gillespie 型算法.

3 生態系統的數值仿真

生態數學中,通過“變化率=輸入-輸出”這樣的思路可以建立起常微分方程模型,如在Logistic 模型基礎上得到的Lotka-Volterra 模型.近年來,一些學者對Lotka-Volterra 模型引入了隨機擾動,將其轉化為隨機Lotka-Volterra 模型,其中隨機性通常是由環境的不可預測性和對生態系統的不完全認知引起的.隨機Lotka-Volterra 模型能夠更精確地描述客觀現象.隨機Lotka-Volterra競爭模型,隨機Lotka-Volterra 捕食模型和隨機Lotka-Volterra 合作模型是其中最為典型的代表.

3.1 隨機Lotka-Volterra 競爭模型

Mao 在其專著[40]中研究了如下的隨機Lotka-Volterra 競爭模型:

其中x(t)=(x1(t),x2(t),···,xd(t))T是生態系統中d個相互作用的種群在t時刻的狀態,b=(b1,···,bd)T∈Rd,σ=(σ1,···,σd)T∈Rd,A=(aij)d×d∈Rd×d為系統的參數.文獻[40]指出,當矩陣A的元素均為非負,即aij≥0(1≤i,j≤d)時,對任意的初始值x(0)∈,方程(3.1)有唯一的強解x(t),并且滿足x(t)∈,其中={(x1,···,xd)T∈Rd:xi≥0,1≤i≤d}.作者還研究了帶延遲的隨機Lotka-Volterra 競爭模型真解的性質.

諸如方程(3.1)這樣的隨機微分方程,具有以下特點:高度非線性,正解,多維.這些特點給數值模擬帶來了很大的挑戰.文獻[41]針對方程(3.1),在截斷算法的基礎上,引入非負投影來確保算法的非負性,給出了一種保正截斷Euler-Maruyama 方法.作者證明了該數值算法是保正的和強收斂的.數值算例顯示,該方法具有1/2 階的強收斂階,但作者并沒有給出相關的理論證明.

文獻[42]則在矩陣A的所有元素非負且對角線元素為正的條件下,針對方程(3.1)提出了一種線性隱式Milstein 方法.該方法不僅能保正,并且能保證解的長時間有界性.數值實驗表明,該方法具有1 階的強收斂階,但作者理論證明方法的強收斂階僅為1/2.

文獻[43]借助對數變換將方程(3.1)轉換成帶有加性噪聲的隨機微分方程,并對轉換后的方程應用顯式Euler-Maruyama 算法進行數值求解.作者通過對數變換使得數值算法可以保持非負性,并利用隨機攝動理論證明了算法的1/2 階強收斂性.

3.2 隨機Lotka-Volterra 捕食模型

除了3.1 節中介紹的競爭模型,捕食模型也是生態學中一種重要的模型.文獻[44]研究了這類模型解的存在唯一性.文獻[45,46]分析了隨機Lotka-Volterra 捕食模型的長時間行為.此外,文獻[47]研究了帶時滯的隨機Lotka-Volterra 捕食模型的正解的存在性.

文獻[48]考慮了如下的一類隨機Lotka-Volterra 捕食模型

其中X(t)=(x1(t),···,xd(t))T和Y(t)=(y1(t),···,yd(t))T分別是t時刻被捕食者和捕食者的數目,η(1)=表示在沒有食物的情況下捕食者的自然死亡率,而η(2)=表示沒有捕食者時被捕食者的自然增長率.

作者證明了模型(3.2)具有唯一的正解,并且具有隨機辛結構,故可以被看作是一個隨機哈密爾頓系統.同時,作者提出了一種隨機Runge-Kutta 方法,并且證明了該數值格式可以保持數值解為正和原系統的隨機辛結構.事實上,作者是通過對數變換來實現保正的.此外,作者還證明了該數值格式是1 階收斂的.

3.3 隨機Lotka-Volterra 合作模型

Lotka-Volterra 合作模型描述了自然界以及人類社會中物種或群體之間的互惠互利關系,如蜜蜂和花.文獻[49] 將白噪聲引入到一類時滯Lotka-Volterra 模型的內稟增長率中,建立了持續性的充分條件,并獲得了平穩分布的存在唯一性.文獻[50] 研究了食物有限情況下的隨機Lotka-Volterra 合作模型,并給出了持續和滅絕的充分條件.文獻[51]考慮了如下帶有年齡結構的隨機Lotka-Volterra 合作模型:

其中X(t,a) 和Y(t,a) 分別是兩個物種在t時刻、a年齡時的密度,γ1(a) 和γ2(a) 分別是物種的內稟增長率,α11(a)和α22(a)表示種內競爭率,α12(a)和α21(a)表示種間合作率,γ1(a),γ2(a),α11(a),α22(a),α12(a)和α21(a)均為正數.

與前面提到的其他模型不同,帶年齡結構的隨機Lotka-Volterra 合作模型(3.3)是一個隨機偏微分方程.事實上,考慮到年齡結構的種群模型都是偏微分方程.文獻[51]證明了方程(3.3)全局解的存在唯一性,并且分析了Euler-Maruyama 方法用于方程(3.3)的數值求解時的收斂率.但[51]并沒有證明數值解的保正性.

值得注意的是,這里談到的Lotka-Volterra 模型不僅可用于描述種群動力學,還可用于物理學和經濟學等方面的建模.

4 隨機傳染病模型的數值仿真

近年來,世界經濟和人類健康均受到新型冠狀病毒的影響,該病毒已經嚴重影響了人們正常的生產生活.歷史上也出現過許多類似的傳染病,例如黑死病,“非典”等.數理流行病模型已經成為了解疫情,掌握疫情變化的重要工具.數學家們建立了很多著名的傳染病數學模型,例如SIS,SIQS,SIR,SEIR 模型等.

事實上,在現實生活中,任何一個系統都會受到或多或少的隨機干擾.正因如此,對具有隨機擾動的傳染病模型的研究是非常重要且是有意義的.據此我們可以得到隨機版本的SIS,SIQS,SIR,SEIR 模型等.這些模型都是具有非線性系數的隨機微分方程,它們的解析解通常很難求出.因此需要使用計算機來對其進行數值模擬,從而達到對傳染病模型的研究.系數不滿足單邊Lipschitz 條件,數值解需要保持真解位于特定區域的性質,這是對隨機傳染病模型進行數值模擬的兩大難點.

4.1 隨機SIR 模型

1927 年數學家Kermack 與McKendrick 在研究傳染病時提出了SIR 倉室模型[52].此模型是最為經典的傳染病數學模型,當時該模型的提出主要針對歐洲出現的黑死病.模型將人群分為3 類,分別為易感者(S 類)、感染者(I 類)和康復者(R 類).該模型適用于治愈后康復者具有很強的免疫力,不會再次被感染的傳染病,如水痘等.此外,該模型還適用于一些致死性的傳染病,因為死亡的病人也可以歸入“康復者(R)”.這里的“康復”可以理解為是退出了系統.

文獻[53]將所提出的數值算法成功運用于如下的隨機SIR 模型:

其中參數Λ,β,μ和γ為正常數,Λ 表示易感個體自然出生率,μ表示每個個體類別自然死亡率,β為疾病接觸率,γ表示受感染個體的恢復率.在初值假設S(0)+I(0)+R(0)=Λ/μ下,模型(4.1)滿足群體總人口不變,即在任意時刻t都有S(t)+I(t)+R(t)=Λ/μ.文獻[53]提出了一種增量馴服的Euler 算法,并得到了模型(4.1)數值近似的強收斂性.

文獻[54]針對模型(4.1)構造了一種維納截斷的線性隱式算法,在確保數值近似取值非負的基礎上,得到了算法是強1/2 階收斂的,并討論了算法的數值動力學行為.

文獻[55]則在另一類隨機擾動下研究了如下隨機SIR 模型:

其中受感染個體以ε的速率變得易感并轉移到易感類別,其他參數與模型(4.1)一致.但由于模型中康復者也受到隨機因素影響,模型(4.2)的人口總數帶有隨機性.針對模型(4.2),作者提出了一種基于勒讓德多項式的勒讓德譜配置方法來求解.

4.2 隨機SIS 模型

在隨機SIR 模型的數值分析中,如何得到算法的收斂階仍然有一些困難.但對相對簡單的隨機SIS 傳染病模型,在其數值近似方面有著豐富的研究成果.

在文獻[56]中,作者考慮了如下的隨機SIS 模型:

其中S(t),I(t)表示t時刻的易感者和感染者數目,N是疾病在其中傳播的人口總數,μ是人均死亡率,1/γ為平均感染期,β是疾病傳播系數.SIS 模型可以看作SIR 模型的特殊情形,適用于感染后不產生永久免疫力的疾病建模,如一些性傳播疾病和細菌性疾病.對于這些疾病,個體起初易感,在某個階段感染該疾病,并在短的感染期后再次易感.文獻[56]證明了模型(4.3)具有唯一的全局正解I(t),并建立了I(t)滅絕和持續的條件.

文獻[57]提出了一種顯式的數值格式,稱為Lamperti 平滑截斷格式(LST),用于對隨機SIS 模型的數值模擬.該方案將Lamperti 型變換(實質上為對數變換)與顯式截斷方法相結合,所得的數值近似結果保持了原始隨機微分方程的解位于特定區域的性質,并且具有1 階均方收斂階.這里使用的Lamperti 型變換只能用于一維的隨機微分方程,變換后的隨機微分方程的系數滿足單邊Lipschitz 條件.由于(4.3)式中的變量滿足S(t)+I(t)=N,所以Lamperti 型變換在這里是可以使用的.此時(4.3)式變為:

當βN-μ-γ=5,β=0.8,σ=0.1,N=10 時,作者給出了用Lamperti 平滑截斷格式求解(4.4)式的收斂階,如圖2所示.其中LST 即為這里的Lamperti 平滑截斷格式.由圖2 可知,該方法具有1 階收斂階.雖然LBE,即Lamperti 向后Euler 格式[63],也具有1 階收斂階,但由于LBE 是隱式的,需要求解超越方程,所以LST 算法的運行時間明顯低于LBE 算法.顯然,LST 算法比LBE 算法更有效.

圖2 LST 算法求解隨機SIS 模型的均方誤差[57]

在文獻[58]中,作者通過結合對數變換和Euler-Maruyama 方法,為隨機SIS 模型構造了一種保正的數值方法,證明了該算法不僅保留了原始隨機微分方程的值域,而且在所有p>0 的情況下,在有限時間間隔內具有1 階Lp收斂階.

文獻[59]考慮了隨機SIS 模型的帶截斷維納過程的分裂步θ方法,證明了該數值格式具有1/2 階均方收斂階,數值解保正且有界.雖然收斂階比[57]提出的LST 算法要低,但是不需要做Lamperti 型變換,可以直接近似求解原方程.事實上,該方法通過將維納過程截斷來保正,因此無需做變換.

此外,文獻[60]構造了一種帶有兩個隨機擾動項的SIS 模型,證明了該模型具有唯一且有界的全局解,分析了模型中感染人群I(t)持續、滅絕的條件.之后文獻[61]通過在平方根內取絕對值的方式,針對該模型提出了一種保取值范圍的數值算法,并得到了該算法的收斂性,但沒有給出算法的收斂階.

5 群體遺傳學模型的數值仿真

群體遺傳學中一個非常重要的模型是Wright-Fisher 模型,于20 世紀30 年代由Wright 和Fisher 獨立引入,它描述了種群遺傳結構中的隨機進化.

文獻[62]考慮了如下的Wright-Fisher 模型:

其中參數的具體意義見[62].Wright-Fisher 模型最初用于模擬有限群體內的遺傳漂移,在文獻[62]中被用作心臟和神經細胞內離子通道動力學的模型.該方程的解落在區間[0,1]內,但多數數值方法無法在數值模擬時保證此類取值范圍.[62]結合隱式平衡算法和分裂步算法,提出了一種平衡隱式分裂步算法,可使得數值解始終以概率1 保持在[0,1]內,并證明了該方法的強收斂性.數值實驗表明,該格式具有強1/2 階收斂性.將該方法推廣到多維情況,數值試驗表明該算法仍然具有1/2 階的強收斂性.取A=1,B=2,C=0.2462,將平衡分裂步算法應用于Wright-Fisher 模型(5.1)的數值實驗結果如圖3所示.

圖3 平衡隱式分裂步算法求解Wright-Fisher 模型(5.1)的收斂性.紅色虛線是斜率為1/2 的參考線[62]

在文獻[63]中,針對這類具有非Lipschitz 漂移或擴散系數且取值在區域D=(l,r)內的標量隨機微分方程,作者提出了一種Lamperti-Euler 格式(LBE).首先借助Lamperti 變換x(t)=F(y(t)),使用向后Euler 格式近似變換過程x(t),再用Lamperti 逆變換進行反向變換,得到原始隨機微分方程的近似值.該數值格式具有1 階收斂性,并且可以保持原始隨機微分方程的取值在D內.

這里的LBE 格式是較早使用Lamperti 變換實現數值解保正的格式.隨機SIS 模型對應的Lamperti 變換是對數變換,而Wright-Fisher 模型對應的則是三角變換.需要注意的是,變換后的隱式數值格式在每一步都需要求解超越方程,計算成本較高.為了克服隱式算法計算成本較大的難點,最近文獻[64]在Lamperti 變換的基礎上提出一種新的顯式傾斜截斷算法,不僅得到了算法的強1 階收斂性,還得到了數值近似的平穩分布收斂性.

6 細胞分化模型的數值仿真

細胞分化過程形成了形態、結構、功能各異的細胞類型,使得生物世界豐富多彩.對基因調控網絡的建模對理解細胞的分化至關重要.構建生物系統的表觀遺傳景觀是研究細胞分化的重要工具.

目前主要有以下兩種方法可以用來構造勢能景觀:一種是基于數據的方法,即從實驗數據,特別是單細胞RNA 測序(scRNA-seq) 數據中識別細胞類型、分化軌跡和細胞多能性,如文獻[65,66,67]等?另一種方法是基于模型的方法,即構建基因調控網絡的動力學方程,分析系統的動力學行為,并使用數值模擬構建勢能景觀,如文獻[68,69]等.

在文獻[68]中,作者用如下帶生滅項的Fokker-Planck 方程作為種群水平上的細胞分化模型:

其中x表示基因表達水平,b(x)表示基因間的相互作用,ε代表噪聲的強度,R(x)表示x處的細胞的出生與死亡率,即生滅項.同時,作者用如下的隨機微分方程作為單細胞水平的細胞分化模型:

其中Xt(ω)為細胞ω在t時刻的基因表達水平,ρt(ω)代表細胞ω的時變權重.由It? 公式可知,(6.1),(6.2)兩個模型的等價關系可以由下式得到

這里,δ是狄拉克函數,E 表示對所有可能的軌道ω取期望.

在構建勢能景觀的過程中,針對低維模型,作者使用有限元方法或者有限差分方法直接對種群水平上的細胞分化模型(6.1)進行數值求解直至穩態?針對高維模型,則對單細胞水平的細胞分化模型(6.2)進行平均場近似,然后再根據所獲得的穩態概率密度函數及作者提出的勢能景觀分解理論,得到細胞類型勢能景觀U(x)和體現細胞干性的多潛能性勢能景觀V(x):U(x)=-εlogPU(x),V(x)=-εlog (P0(x)/PU(x)),其中PU(x)表示當R(x)給定時,由(6.1)或(6.2)得到的穩態概率密度函數,P0(x)表示當R(x)≡0 時,由(6.1)或(6.2)得到的穩態概率密度函數.

7 結論

隨著近年來交叉學科的蓬勃發展,涌現了大量與系統生物學相關的隨機微分方程模型.這些模型由于來源于現實世界,有一定的生物背景,通常具有以下一個或者若干個特點:解位于特定區域,高維,系數不滿足Lipschitz 條件等.以上這些特點為數值算法的構造及理論分析帶來了挑戰.常用的數值算法用來做這類模型的數值仿真時要么不收斂,要么不能保持真解的某些性質,如位于特定的區域等.Lamperti 型變換算法,隱式高階算法或者平衡隱式算法可能是針對這類問題構造數值格式的備選方案.目前已有的保取值域的方法主要包括:投影(文獻中也稱截斷),相關文獻有[41]等?對數變換,相關文獻有[47,57,58]等?平衡隱式方法,相關文獻有[62]等?維納增量截斷,相關文獻有[59]等.

對于引言中提到的帶有隨機輸入的微分方程(random ordinary differential equation)的數值模擬,可參考文獻[70].系統生物學家通常使用這類方程對生態系統進行建模,如文獻[71,72]等.

鑒于篇幅和時間有限,本文所綜述的工作難免掛一漏萬,敬請讀者海涵.