牛原玲 陳琳 陳洛南

(1.中南大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,長(zhǎng)沙,410083?2.江西財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與數(shù)據(jù)科學(xué)學(xué)院,南昌,330013?3.中國(guó)科學(xué)院系統(tǒng)生物學(xué)重點(diǎn)實(shí)驗(yàn)室,分子細(xì)胞科學(xué)卓越創(chuàng)新中心(上海生物化學(xué)與細(xì)胞生物學(xué)研究所),上海,200031)

1 引言

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

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

如果隨機(jī)擾動(dòng)出現(xiàn)在微分方程的系數(shù)中,如

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

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

本文主要回顧系統(tǒng)生物學(xué)中的隨機(jī)微分方程模型及相關(guān)的數(shù)值仿真算法.這些隨機(jī)微分方程一般具有高維、高度非線性、真解位于某個(gè)特定的區(qū)域等特點(diǎn),因此對(duì)它們的數(shù)值模擬需要做專門(mén)的研究.鑒于篇幅有限,我們僅從生物化學(xué)(生化)反應(yīng)系統(tǒng)、生態(tài)系統(tǒng)、傳染病模型、群體遺傳學(xué)模型以及細(xì)胞分化模型的數(shù)值仿真來(lái)加以闡述.

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

2 生化反應(yīng)系統(tǒng)的數(shù)值仿真

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

2.1 化學(xué)主方程模型

化學(xué)主方程(Chemical master equation,簡(jiǎn)記為CME)描述了生化反應(yīng)系統(tǒng)在不同時(shí)間點(diǎn)處于每種狀態(tài)的概率:

其中X(t)是狀態(tài)變量,aj(X(t))是傾向函數(shù),aj(X(t))dt表示在下一個(gè)小時(shí)間段[t,t+dt)第j個(gè)反應(yīng)Rj發(fā)生的概率.

化學(xué)主方程似乎是發(fā)現(xiàn)隨機(jī)系統(tǒng)時(shí)間演化規(guī)律的理想方法,它可以返回系統(tǒng)的概率密度函數(shù).然而,化學(xué)主方程只有在非常簡(jiǎn)單的情況下可以求出真解.也就是說(shuō),在大多數(shù)情況下,需要借助數(shù)值方法來(lái)對(duì)其進(jìn)行數(shù)值模擬.此外,由于必須在每個(gè)時(shí)間點(diǎn)對(duì)每個(gè)可能的狀態(tài)進(jìn)行評(píng)估,并且狀態(tài)的數(shù)量可以呈指數(shù)增長(zhǎng),因此計(jì)算量非常大,這使得化學(xué)主方程模型最初只能有效地用于短期模擬小規(guī)?；蛘吆?jiǎn)單系統(tǒng).近年來(lái),雖然化學(xué)主方程的數(shù)值求解方法取得了長(zhǎng)足的進(jìn)步[21,22,23],但當(dāng)分子數(shù)和傾向函數(shù)太大以致于無(wú)法使用化學(xué)主方程計(jì)算時(shí),必須使用新的數(shù)值模擬方法,如Gillespie型算法.

2.2 離散狀態(tài)馬爾可夫鏈模型

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

2.3 郎之萬(wàn)方程模型

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

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

帶邊界條件的生化反應(yīng)模型如下:

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

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

Goldbeter-Koshland 開(kāi)關(guān)是由兩種具有相反作用的酶E1 和E2 促進(jìn)的共價(jià)修飾生化反應(yīng)系統(tǒng),它由以下6 個(gè)反應(yīng)組成:

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

圖1 SSA 與帶反射邊界的隨機(jī)微分方程對(duì)Goldbeter-Koshland 開(kāi)關(guān)數(shù)值模擬的均值和方差(實(shí)線:SSA,虛線:帶反射邊界的隨機(jī)微分方程)[37]

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

3 生態(tài)系統(tǒng)的數(shù)值仿真

生態(tài)數(shù)學(xué)中,通過(guò)“變化率=輸入-輸出”這樣的思路可以建立起常微分方程模型,如在Logistic 模型基礎(chǔ)上得到的Lotka-Volterra 模型.近年來(lái),一些學(xué)者對(duì)Lotka-Volterra 模型引入了隨機(jī)擾動(dòng),將其轉(zhuǎn)化為隨機(jī)Lotka-Volterra 模型,其中隨機(jī)性通常是由環(huán)境的不可預(yù)測(cè)性和對(duì)生態(tài)系統(tǒng)的不完全認(rèn)知引起的.隨機(jī)Lotka-Volterra 模型能夠更精確地描述客觀現(xiàn)象.隨機(jī)Lotka-Volterra競(jìng)爭(zhēng)模型,隨機(jī)Lotka-Volterra 捕食模型和隨機(jī)Lotka-Volterra 合作模型是其中最為典型的代表.

3.1 隨機(jī)Lotka-Volterra 競(jìng)爭(zhēng)模型

Mao 在其專著[40]中研究了如下的隨機(jī)Lotka-Volterra 競(jìng)爭(zhēng)模型:

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

諸如方程(3.1)這樣的隨機(jī)微分方程,具有以下特點(diǎn):高度非線性,正解,多維.這些特點(diǎn)給數(shù)值模擬帶來(lái)了很大的挑戰(zhàn).文獻(xiàn)[41]針對(duì)方程(3.1),在截?cái)嗨惴ǖ幕A(chǔ)上,引入非負(fù)投影來(lái)確保算法的非負(fù)性,給出了一種保正截?cái)郋uler-Maruyama 方法.作者證明了該數(shù)值算法是保正的和強(qiáng)收斂的.數(shù)值算例顯示,該方法具有1/2 階的強(qiáng)收斂階,但作者并沒(méi)有給出相關(guān)的理論證明.

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

文獻(xiàn)[43]借助對(duì)數(shù)變換將方程(3.1)轉(zhuǎn)換成帶有加性噪聲的隨機(jī)微分方程,并對(duì)轉(zhuǎn)換后的方程應(yīng)用顯式Euler-Maruyama 算法進(jìn)行數(shù)值求解.作者通過(guò)對(duì)數(shù)變換使得數(shù)值算法可以保持非負(fù)性,并利用隨機(jī)攝動(dòng)理論證明了算法的1/2 階強(qiáng)收斂性.

3.2 隨機(jī)Lotka-Volterra 捕食模型

除了3.1 節(jié)中介紹的競(jìng)爭(zhēng)模型,捕食模型也是生態(tài)學(xué)中一種重要的模型.文獻(xiàn)[44]研究了這類模型解的存在唯一性.文獻(xiàn)[45,46]分析了隨機(jī)Lotka-Volterra 捕食模型的長(zhǎng)時(shí)間行為.此外,文獻(xiàn)[47]研究了帶時(shí)滯的隨機(jī)Lotka-Volterra 捕食模型的正解的存在性.

文獻(xiàn)[48]考慮了如下的一類隨機(jī)Lotka-Volterra 捕食模型

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

作者證明了模型(3.2)具有唯一的正解,并且具有隨機(jī)辛結(jié)構(gòu),故可以被看作是一個(gè)隨機(jī)哈密爾頓系統(tǒng).同時(shí),作者提出了一種隨機(jī)Runge-Kutta 方法,并且證明了該數(shù)值格式可以保持?jǐn)?shù)值解為正和原系統(tǒng)的隨機(jī)辛結(jié)構(gòu).事實(shí)上,作者是通過(guò)對(duì)數(shù)變換來(lái)實(shí)現(xiàn)保正的.此外,作者還證明了該數(shù)值格式是1 階收斂的.

3.3 隨機(jī)Lotka-Volterra 合作模型

Lotka-Volterra 合作模型描述了自然界以及人類社會(huì)中物種或群體之間的互惠互利關(guān)系,如蜜蜂和花.文獻(xiàn)[49] 將白噪聲引入到一類時(shí)滯Lotka-Volterra 模型的內(nèi)稟增長(zhǎng)率中,建立了持續(xù)性的充分條件,并獲得了平穩(wěn)分布的存在唯一性.文獻(xiàn)[50] 研究了食物有限情況下的隨機(jī)Lotka-Volterra 合作模型,并給出了持續(xù)和滅絕的充分條件.文獻(xiàn)[51]考慮了如下帶有年齡結(jié)構(gòu)的隨機(jī)Lotka-Volterra 合作模型:

其中X(t,a) 和Y(t,a) 分別是兩個(gè)物種在t時(shí)刻、a年齡時(shí)的密度,γ1(a) 和γ2(a) 分別是物種的內(nèi)稟增長(zhǎng)率,α11(a)和α22(a)表示種內(nèi)競(jìng)爭(zhēng)率,α12(a)和α21(a)表示種間合作率,γ1(a),γ2(a),α11(a),α22(a),α12(a)和α21(a)均為正數(shù).

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

值得注意的是,這里談到的Lotka-Volterra 模型不僅可用于描述種群動(dòng)力學(xué),還可用于物理學(xué)和經(jīng)濟(jì)學(xué)等方面的建模.

4 隨機(jī)傳染病模型的數(shù)值仿真

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

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

4.1 隨機(jī)SIR 模型

1927 年數(shù)學(xué)家Kermack 與McKendrick 在研究傳染病時(shí)提出了SIR 倉(cāng)室模型[52].此模型是最為經(jīng)典的傳染病數(shù)學(xué)模型,當(dāng)時(shí)該模型的提出主要針對(duì)歐洲出現(xiàn)的黑死病.模型將人群分為3 類,分別為易感者(S 類)、感染者(I 類)和康復(fù)者(R 類).該模型適用于治愈后康復(fù)者具有很強(qiáng)的免疫力,不會(huì)再次被感染的傳染病,如水痘等.此外,該模型還適用于一些致死性的傳染病,因?yàn)樗劳龅牟∪艘部梢詺w入“康復(fù)者(R)”.這里的“康復(fù)”可以理解為是退出了系統(tǒng).

文獻(xiàn)[53]將所提出的數(shù)值算法成功運(yùn)用于如下的隨機(jī)SIR 模型:

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

文獻(xiàn)[54]針對(duì)模型(4.1)構(gòu)造了一種維納截?cái)嗟木€性隱式算法,在確保數(shù)值近似取值非負(fù)的基礎(chǔ)上,得到了算法是強(qiáng)1/2 階收斂的,并討論了算法的數(shù)值動(dòng)力學(xué)行為.

文獻(xiàn)[55]則在另一類隨機(jī)擾動(dòng)下研究了如下隨機(jī)SIR 模型:

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

4.2 隨機(jī)SIS 模型

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

在文獻(xiàn)[56]中,作者考慮了如下的隨機(jī)SIS 模型:

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

文獻(xiàn)[57]提出了一種顯式的數(shù)值格式,稱為L(zhǎng)amperti 平滑截?cái)喔袷?LST),用于對(duì)隨機(jī)SIS 模型的數(shù)值模擬.該方案將Lamperti 型變換(實(shí)質(zhì)上為對(duì)數(shù)變換)與顯式截?cái)喾椒ㄏ嘟Y(jié)合,所得的數(shù)值近似結(jié)果保持了原始隨機(jī)微分方程的解位于特定區(qū)域的性質(zhì),并且具有1 階均方收斂階.這里使用的Lamperti 型變換只能用于一維的隨機(jī)微分方程,變換后的隨機(jī)微分方程的系數(shù)滿足單邊Lipschitz 條件.由于(4.3)式中的變量滿足S(t)+I(t)=N,所以Lamperti 型變換在這里是可以使用的.此時(shí)(4.3)式變?yōu)?

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

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

在文獻(xiàn)[58]中,作者通過(guò)結(jié)合對(duì)數(shù)變換和Euler-Maruyama 方法,為隨機(jī)SIS 模型構(gòu)造了一種保正的數(shù)值方法,證明了該算法不僅保留了原始隨機(jī)微分方程的值域,而且在所有p>0 的情況下,在有限時(shí)間間隔內(nèi)具有1 階Lp收斂階.

文獻(xiàn)[59]考慮了隨機(jī)SIS 模型的帶截?cái)嗑S納過(guò)程的分裂步θ方法,證明了該數(shù)值格式具有1/2 階均方收斂階,數(shù)值解保正且有界.雖然收斂階比[57]提出的LST 算法要低,但是不需要做Lamperti 型變換,可以直接近似求解原方程.事實(shí)上,該方法通過(guò)將維納過(guò)程截?cái)鄟?lái)保正,因此無(wú)需做變換.

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

5 群體遺傳學(xué)模型的數(shù)值仿真

群體遺傳學(xué)中一個(gè)非常重要的模型是Wright-Fisher 模型,于20 世紀(jì)30 年代由Wright 和Fisher 獨(dú)立引入,它描述了種群遺傳結(jié)構(gòu)中的隨機(jī)進(jìn)化.

文獻(xiàn)[62]考慮了如下的Wright-Fisher 模型:

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

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

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

這里的LBE 格式是較早使用Lamperti 變換實(shí)現(xiàn)數(shù)值解保正的格式.隨機(jī)SIS 模型對(duì)應(yīng)的Lamperti 變換是對(duì)數(shù)變換,而Wright-Fisher 模型對(duì)應(yīng)的則是三角變換.需要注意的是,變換后的隱式數(shù)值格式在每一步都需要求解超越方程,計(jì)算成本較高.為了克服隱式算法計(jì)算成本較大的難點(diǎn),最近文獻(xiàn)[64]在Lamperti 變換的基礎(chǔ)上提出一種新的顯式傾斜截?cái)嗨惴?不僅得到了算法的強(qiáng)1 階收斂性,還得到了數(shù)值近似的平穩(wěn)分布收斂性.

6 細(xì)胞分化模型的數(shù)值仿真

細(xì)胞分化過(guò)程形成了形態(tài)、結(jié)構(gòu)、功能各異的細(xì)胞類型,使得生物世界豐富多彩.對(duì)基因調(diào)控網(wǎng)絡(luò)的建模對(duì)理解細(xì)胞的分化至關(guān)重要.構(gòu)建生物系統(tǒng)的表觀遺傳景觀是研究細(xì)胞分化的重要工具.

目前主要有以下兩種方法可以用來(lái)構(gòu)造勢(shì)能景觀:一種是基于數(shù)據(jù)的方法,即從實(shí)驗(yàn)數(shù)據(jù),特別是單細(xì)胞RNA 測(cè)序(scRNA-seq) 數(shù)據(jù)中識(shí)別細(xì)胞類型、分化軌跡和細(xì)胞多能性,如文獻(xiàn)[65,66,67]等?另一種方法是基于模型的方法,即構(gòu)建基因調(diào)控網(wǎng)絡(luò)的動(dòng)力學(xué)方程,分析系統(tǒng)的動(dòng)力學(xué)行為,并使用數(shù)值模擬構(gòu)建勢(shì)能景觀,如文獻(xiàn)[68,69]等.

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

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

其中Xt(ω)為細(xì)胞ω在t時(shí)刻的基因表達(dá)水平,ρt(ω)代表細(xì)胞ω的時(shí)變權(quán)重.由It? 公式可知,(6.1),(6.2)兩個(gè)模型的等價(jià)關(guān)系可以由下式得到

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

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

7 結(jié)論

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

對(duì)于引言中提到的帶有隨機(jī)輸入的微分方程(random ordinary differential equation)的數(shù)值模擬,可參考文獻(xiàn)[70].系統(tǒng)生物學(xué)家通常使用這類方程對(duì)生態(tài)系統(tǒng)進(jìn)行建模,如文獻(xiàn)[71,72]等.

鑒于篇幅和時(shí)間有限,本文所綜述的工作難免掛一漏萬(wàn),敬請(qǐng)讀者海涵.