復雜生物膜動力系統的顯式有限差分方法研究

孫高峰,牛瑞萍,李 明

(太原理工大學 數學學院,太原 030024)

復雜生物膜動力系統的顯式有限差分方法研究

孫高峰,牛瑞萍,李 明

(太原理工大學 數學學院,太原 030024)

針對刻畫生物膜成長的動力系統設計了一個非負、有界的顯式差分格式,系統是由四個變系數的非線性拋物型偏微分方程組成,描述了營養基質上多種微生物成分間的交互過程；所設計的有限差分格式巧妙地將非線性項轉化成線性形式,確保得到的數值解是非負、有界的。最后通過兩個數值算例檢驗了方法的性能,并給出了宏觀意義下的數值解的誤差。

復雜生物膜系統;非線性反應擴散方程組;顯式有限差分法;非負、有界

在科學計算中經常要用數值求解各類微分方程[1],有限元[2-3]、有限差分法[6-8]和有限體積法[4-5]是最重要的常用方法。本文使用有限差分法求解生物膜系統的數值解.生物膜系統本身物理屬性要求數值解是非負、有界的,因此,設計非負、有界的有限差分技術已成為計算數學的重要研究領域[6-8]。

Eberl et al建立的拋物型偏微分方程組來描述復雜的生物系統與周圍營養基的相互作并證明了該問題在合理的解析約束下,非負、有界解的存在唯一性[9-11]。受Eberl所做工作的啟發,目前的工作主要是制定一種新的顯式有限差分法來有效地求解非線性拋物型偏微分方程組。由于生物膜系統是具有實際意義的體系,其未知變量表示各種微生物及營養成分的密度,是非負、有界的變量。如何設計一個有效的差分格式確保數值解的非負、有界性尤為重要。通常的生物系統的擴散系數是給定非負常數,這樣的系統設計數值計算相對容易,文獻中[6-11]已有許多成功的方法，但Eberl et al的生物系統的擴散系數是未知變量的函數,使得控制方程組是一個非線性的變系數偏微分方程組,這增加了設計數值計算方法的難度;本文的主要內容是在初邊值條件下對復雜生物膜系統顯式有限差分方法進行研究。

1 生物背景及數學模型

復雜生物膜系統由下面四個非線性偏微分方程確定:

(1)

(2)

(3)

(4)

M(X,t)=B(X,t)+I(X,t)+E(X,t) .

(5)

本文,假定μ,kL,kE,kI,kS,YE為非負常數,YH正常數,這些所有的常數均是通過實驗獲得的。和·分別表示空間梯度算子和散度算子。S表示關于最大值標準化后營養濃度;函數M(X,t)表示問題域上的系統生物質量密度。復雜生物膜系統中包含種類型的微生物:B表示活的微生物密度;I表示惰性或死的微生物密度;E表示胞外聚合物密度,這三個變量是關于最大量標準化的得到的量。顯然,對于和M∈[0,1)。從生物學角度看,偏微分方程組(1)-(5)刻畫的是不同生物之間以及生物和營養基之間的相互交互轉化過程。函數D(u)描述生物演化動力學(生長、衰減、繁衍)規律的非線性擴散系數,通常是微生物的函數,因此它是非線性的“根源”,本章使用下列形式的擴散系數:

式中,α≥1,β≥1.

方程(1)刻畫的是生物膜上的營養的消耗和轉化,值得注意的是活的微生物和胞外聚合物是營養的潛在“源”(可以轉化成營養),右端的第三和第四項刻畫的是轉化機制。這里,正常數kL表示活生物轉化為營養物的效率;正常數kE表示胞外聚合物轉化成營養物的效率;右端第一項營養物的擴散項;右端第二項是標準的Monod-type函數,描述的是活生物如何消耗營養物.參數μ表示活生物的最大增長率;參數1/YH表示活生物消耗營養率。

方程(2)刻畫了生物膜上活生物的增長和轉。右端第一項控制活生物的擴散;第二項是標準的Monod-type函數,描述的是營養物如何促進活生物的增長;第三和第四項描述了活生物如何轉化成營養物和惰性生物;正常數kI表示活生物轉化成惰性生物的效率。

方程(3)右端的第一項控制惰性生物的擴散,第二項表示活生物如何促進惰性生物的增長。

方程(4)描述的是胞外聚合物的消耗和轉化。方程右端的第一項控制胞外聚合物的擴散;第二項描述的活生物如何促進胞外聚合物的增長;常數YE表示胞外聚合物消耗活生物的效率。

2 有限差分方法

下面我們給出由方程式(1)-(5)所控制復雜生物膜系統的顯式非負、有界差分逼近形式。設問題域為正規的時-空區域:

(x,y,t)∈[0,X]×[0,Y]×[0,T] .

控制方程(1)-(4)的擴散項可統一表示成下列方程:

·(D(M)f)=D(M)xfx+

D(M)yfy+D(M)(fxx+fyy) ,

式中：函數f依次可代表S,B,I,E。

下面我們使用以下有限差分形式對系統進行近似:

δxxf=

δyyf=

為簡化記號,我們令

下面給出控制方程(1)-(4)式中擴散項的逼近形式為:

·(D(M)f)≈

其中,函數f依次可代表S,B,I,E。為進一步簡化記號,令

因此,復雜的生物膜系統(1)-(4)可由下面有限差分近似:

3 離散有限差分格式

為了表述的方便,下面給出一維模型的有限差分格式。

設a,b為實數,且a

xi=(i-1)Δx.

i=1,2,3,…,M+1,且Δx=(b-a)/M.tk=(k-1)Δt.k=1,2,3,…,P+1,且Δt=T/P.

(6)

(7)

(8)

(9)

(10)

為了說明所提出的有限差分格式的優越性,下面給出另外一種標準中心有限差分格式:

(11)

(12)

(13)

(14)

這里

4 數值算例

下面根據系統定義邊界為Neumann邊界條件的系統宏觀誤差。

聯立系統方程(1)-(5),有:

(15)

將上式方程在域(x,y,t)∈[0,X]×[0,Y]×[0,t]上積分,得到下列積分方程形式:

式中：b1=μc2,b2=(c1-c5)c2,b3=(c1-c5)(c2+c4)-μc2,b4=μc4.

進一步簡化記號,令

這樣方程(15)可以改寫成下列形式:

b1(St-S0)+b2(Bt-B0)+

b3(It-I0)+b4(Et-E0)=0 .

因此，宏觀誤差Em可定義為下列形式:

Em=|b1(SNt-S0)+b2(BNt-B0)+

b3(INt-I0)+b4(ENt-E0)| .

(16)

這里,

SN(x,y,t),BN(x,y,t),IN(x,y,t)及EN(x,y,t)依次表示營養物濃度、活生物密度、死生物密度及胞外聚合物密度在點(x,y,t)處的數值解。

算例1 計算使用的參數μ=0.02,YH=0.003,YE=KE=0.03,KS=0.2,KI=0.04,KL=0.04;α=β=2,問題域為Ω=(0,1),控制方程為(1)-(5),初始條件為:

S0(x)=1 ,

B0(x)=0.15e-200(x-0.3)2+0.2e-150(x-0.65)2,

I0(x)=0.04e-80(x-0.2)2+0.03e-60(x-0.7)2,

E0(x)=0.025 .

在邊界網格點上施加均勻的Neumann邊界條件，使用的網格為Δt=0.001 and Δx=0.01，數值結果如圖 1所示。

圖1 使用有限差分格式(6)-(9)得到的數值解

由式(16)計算出的5個時刻t=0.25,0.5,0.75,1,2的宏觀誤差依次為:

1.589 9×10-4,1.368 3×10-4,

1.131 9×10-4,8.583 0×10-5,3.162 4×10-4.

圖2給出的是4個不同時刻t=5,10,20,30的數值解。相應的宏觀誤差依次為:2.685 6×10-5,2.180 2×10-6,9.097 0×10-5,1.178 6×10-4.在這種長期情況下,營養基質函數是趨于零的,而活微生物擴散到整個空間區域。對于長時間的行為,進一步的實驗證實活的微生物密度最終趨于零,這也是我們所期望的。

綜合圖 1和圖2,我們發現,起初隨著營養基質密度的減小,總的微生物量有所增加,隨著時間的增加,活性微生物轉化成惰性微生物。

從圖3可以觀察到,使用有限差分格式(11)-(14)得到的數值解具有奇異性和震蕩現象,說明該格式不穩定,不能確保數值解是非負、有界的。而本文提出的有限差分格式是穩定的,得到的數值解是非負、有界的。

算例2 計算使用的參數μ=2,YH=0.35,YE=KE=0.03,KS=0.2,KI=0.4,KL=0.001;α=β=4,問題域為Ω=(0,1),控制方程為(1)-(5),初始條件為:

S0(x)=0.8，

B0(x)=0.15e-200(x-0.3)2+0.2e-150(x-0.65)2，

I0(x)=0.02e-80x2+0.01e-60(x-1)2，

E0(x)=0.005 .

在邊界網格點上施加均勻的Neumann邊界條件,使用的網格為Δt=0.001 and Δx=0.01,數值結果見圖4。

圖4給出了4個不同時刻t=1,2,3,10的數值解.相應的宏觀誤差依次為:1.765 1×10-6,6.717 4×10-6,3.037 0×10-6,可以觀察到,在起初階段,隨著營養基質函數的減小,活微生物密度總量增加;當營養基質趨于零時,活微生物密度總減小,逐漸趨于零;而對惰性生物密度總量始終是增加的。這些結果符合系統的特征。

圖2 使用有限差分格式(6)-(9)得到的數值解

圖3 使用有限差分格式(11)-(14)得到的數值解。使用的控制方程、初始條件、邊界條件、計算參數以及計算網格和圖 1和圖2相同

圖4 使用有限差分格式(6)-(9)得到的數值解

從圖5可以觀察到,對數值算例2使用有限差分格式(11)-(14)得到的數值解也具有奇異性和震蕩現象,說明該格式不穩定,不能確保數值解是非負、有界的.而本章提出的有限差分格式是穩定的,得到的數值解是非負、有界的。

5 結論

本文用新的顯式有限差分法來逼近復雜生物膜系統的正的、有界解,其控制方程是由四個拋物線型偏微分方程構成,偏微分方程中同樣也包含一個與密度相關的非線性擴散反應項.理論上,模型的非負、有界解的存在性和唯一性文獻中已有證明,而獲得該模型的精確解確實是很困難的,因此,我們不得不采用數值模擬的方法來得到近似解。本文的理論和數值研究總結如下:

1) 建立了復雜生物膜系統的顯式有限差分格式;

2) 給出了2個數值算例，結果表明數值結果和實際是相吻合的，能夠反應微生物成長規律;

3) 和標準中心差分方法相比較，本文所提出的方法具有較好的穩定性，不會出現奇異的數值解和震蕩現象;

4) 不足之處，相關收斂性、穩定性等理論體系還沒建立起來。

圖5 使用有限差分格式(11)-(14)得到的數值解。使用的控制方程、初始條件、邊界條件、計算參數以及計算網格和圖4相同

[1] 李治平.偏微分方程數值解講義[M].北京:北京大學出版社,2010.

[2] Liu WK,Liu YL,Farrell D.Immersed finite element method and its applications to biological systems[J].Computer Methods In Applied Mechanics And Engineering,2006,195(13-16):1722-1749.

[3] Liu G R,Quek S S.Finite Element Method:a practical course[M].Burlington MA:BH,2003.

[4] Wang Shuai,Yuan Guangwei,Li Yonghai.A monotone finite volume scheme for advection-diffusion equations on distorted meshes[J].International Journal for Numerical Methods in Fluids,2012,69(7):1283-1298.

[5] Hermeline F.A finite volume method for the approximation of convection-diffusion equations on general meshes[J].International Journal For Numerical Methods In Engineering,2012,91(12):1331-1357.

[6] Macías-Díaz J E.A bounded finite-difference discretization of a two-dimensional diffusion equation with logistic nonlinear reaction[J].Int J Mod Phys C,2011,22:953-966.

[7] Macías-Díaz J E,Jerez-Galiano S,Puri A.Positivity-preserving methods for a linearized Fisher-KPP equation with consistency properties in the energy domain[J].J Differ Equ Appl,2010,16:389-405.

[8] Macías-Díaz J E,Puri A.A boundedness-preserving finite-difference scheme for a dampednonlinear wave equation[J].Appl Numer Math,2010,60:934-948.

[9] Eberl H J,Parker D,van Loosdrecht M C M.A new deterministic spatio-temporal continuum model for biofilm development[J].J Theor Med,2001,3:161-176.

[10] Efendiev M A,Eberl H J,Zelik S V.Existence and longtime behavior ofsolutions of a nonlinear reaction-diffusion system arising in the modeling of biofilms[J].RIMS Kokyuroko,2002,1258:49-71.

[11] Eberl H J.A deterministic continuum model for the formation of EPS in heterogeneous boilarchitectures[J].Proc Biofilms,2004,237-242.

(編輯：朱 倩)

Explicit Finite-difference Method for a ComplexBio-film Dynamics System

SUN Gaofeng,NIU Ruiping,LI Ming

(CollegeofMathematics,TaiyuanUniversityofTechnology,Taiyuan030024,China)

Numerical methods with non-negative and bounded solutions have great significance because of the requirement from many mathematical models with practical background.In this paper,a non-negative and bounded explicit finite-difference scheme is presented for bio-film growth dynamics system.The scheme describes mutual process among different microorganism compositions on nutritional substance and is expressed by four nonlinear parabolic partial differential equations with variable coefficients.The explicit finite-difference scheme neatly transforms the nonlinear form into linear form,which ensures the non-negative and bounded property of solutions.Numerical examples are presented to illustrate the performance of the methods.Error in macro significance is also given.

complex biofilm dynamics system；explicit finite-difference method；nonlinear reaction-diffusion equations；Nonnegative bounded

1007-9432(2015)06-0790-07

2015-04-09

國家自然科學基金資助項目:可商業化光滑有限元分析軟件的研發(11472184)

孫高峰(1976-)， 男，山西臨猗人，博士，主要從事計算數學的研究，(Tel)13453157694

李明(1982-)，男，教授，主要從事無網絡計算方法，偏微方程中的反問題等的研究，(E-mail)liming01@tyut.edu.cn

O241

A

10.16355/j.cnki.issn1007-9432tyut.2015.06.028