基于多元相場理論的細菌生物膜生長動力學建模及其數值模擬

張曙梅,徐向榮,徐浩,周濤,李雙,李妍

(1.安徽工業大學機械工程學院,安徽馬鞍山243032; 2.安徽工業大學數理學院,安徽馬鞍山243032)

基于多元相場理論的細菌生物膜生長動力學建模及其數值模擬

張曙梅1,徐向榮1,徐浩2,周濤1,李雙1,李妍1

(1.安徽工業大學機械工程學院,安徽馬鞍山243032; 2.安徽工業大學數理學院,安徽馬鞍山243032)

細菌生物膜普遍存在于自然界和工業環境中.將細菌生物膜視為由多種組元組成的不可壓縮的黏性流體,結合多元相場理論、物質連續方程、動量守恒方程、溶質場控制方程和生物反應動力方程對細菌生物膜生長動力學進行建模.通過對模型的數值模擬,可以描述細菌和胞外聚合物在靜態環境和動態環境中的生長狀況以及營養物質的消耗規律.

細菌生物膜;相場理論;動力學模型;數值模擬

細菌生物膜是單一或多種細菌為了適應自然環境而形成的微菌落聚集物,其主要成分為多糖蛋白復合物[1-2].細菌生物膜的生長是一個產生、發展、壯大的過程.通過對細菌生物膜生長規律的研究以及建立細菌生物膜生長動力學模型,既可以揭示細菌生物膜的生長機理,又可以有效控制細菌生物膜的形成和發展,抑制其有害作用[3].已有的研究主要是從數學、流體力學、計算機仿真以及實驗的角度對細菌生物膜進行研究,具體可分為低維連續模型、有限擴散凝聚模型、元胞自動機模型和單流體相場模型[4-9].本工作的研究內容是屬于單流體相場模型的范疇.相場理論以Ginzburg-Landau理論為物理基礎,通過微分方程來體現具有特定物理機制的擴散、有序化勢和熱力學驅動的綜合作用,通過計算機編程求解上述方程,從而獲取研究體系在時間和空間上的瞬時狀態[10-11].Cogan等[12-13]提出的模型是將細菌生物膜看作由胞外聚合物和等效溶劑組成的生物凝膠體.由于滲透壓的作用,生物凝膠體吸收或排除溶劑導致其膨脹或收縮.因此,影響細菌生物膜形態的物理化學因素主要是滲透壓.但在細菌生物膜內,細菌和胞外聚合物均屬于高分子,滲透壓不僅作用于胞外聚合物上,而且也作用在細菌上.因此,僅考慮滲透壓對胞外聚合物的作用是不合理的.Zhang等[14-15], Brandon等[16]和Chen等[17]在Cogan研究的基礎上提出了新的模型,將細菌生物膜視為不可壓縮的多組元混合流體,主要考慮細菌生物膜中胞外聚合物、細菌和等效溶劑(營養物質和水)的體積分數的變化.但此模型只考慮了混合自由能密度梯度、流體的靜壓力以及剪切力對細菌生物膜的作用,忽略了細菌生物膜內各組元的摩擦阻力的作用,以及滲透壓對細菌生物膜生長的影響.

鑒于上述模型的不足,本工作提出了一種單流體二組元的細菌生物膜生長動力學模型.本模型既考慮了細菌生物膜內部生物質的分布不均勻性,又考慮了滲透壓、摩擦阻力和流體剪力對細菌生物膜生長的影響.在本工作模型中,仍將整個細菌生物膜系統近似認為是不可壓縮的,并且假設細菌生物膜是由胞外聚合物、細菌和等效溶劑(營養物質和水)組成,側重考慮了滲透壓、流體剪力以及不同組元之間摩擦阻力的作用.通過對建立的多組元細菌生物膜生長動力學模型的數值求解,模擬了靜態和動態環境中細菌生物膜膜內細菌、胞外聚合物的體積分數變化,以及營養物質濃度在空間上的分布規律和流體剪力對細菌生物膜生長的影響.靜態環境,即細菌生物膜系統的體積平均速度為0,vx=0.例如,存有靜態液體的管道和容器、實驗室用于培養細菌生物膜的培養皿和潮濕的載體表面等都屬于細菌生物膜生長的靜態環境.細菌生物膜在此環境中生長,其膜內的胞外聚合物、等效溶劑和營養物質主要通過擴散作用進行物質傳輸以及生物化學反應.動態環境,即細菌生物膜系統的體積平均速度不為0,vx/=0.例如,存有流動液體的管道、污水處理池、受流水沖擊的巖石表面以及任何動物的血管壁等都屬于細菌生物膜生長的動態環境.當細菌生物膜在動態環境中生長時,膜內的聚合物網狀結構、等效溶劑和營養物質主要通過對流-擴散共同作用進行物質傳輸以及生物化學反應.由于在靜態環境和動態環境中細菌生物膜系統的物質傳輸機理是不一樣的,因此其在這兩種條件下的生長是存在差異的.

1 模型的建立

本工作將細菌生物膜視為由多種組元組成的黏性流體,即細菌、胞外聚合物和等效溶劑,并運用多元相場理論對細菌生物膜生長動力學進行建模.

1.1 多元細菌生物膜系統的動量方程

設φp為胞外聚合物的體積分數,?b為細菌的體積分數,?s為等效溶劑的體積分數,vp為胞外聚合物的速度,vb為細菌的速度,vs為等效溶劑的速度,ρp為胞外聚合物的密度,ρb為細菌的密度,ρs為等效溶劑的密度.由Navier-Stokes方程得

式中,Πp為作用于胞外聚合物上的滲透壓,Πb為作用于細菌上的滲透壓;ζpb為胞外聚合物和細菌之間的體積摩擦阻力系數,ζp為胞外聚合物和等效溶劑之間的體積摩擦阻力系數,ζb為細菌和等效溶劑之間的體積摩擦阻力系數,ρ為細菌生物膜中胞外聚合物、細菌和等效溶劑的體積平均密度,且ρ=?pρp+?bρb+?sρs,P為作用于整個細菌生物膜上的液體靜壓力,τp為胞外聚合物的應力,τb為細菌的應力,τs為等效溶劑的應力.

將式(1)～(3)相加可得

由?p+?b+?s=1可得令

細菌生物膜的體積平均速度v=?pvp+?bvb+?svs,則式(5)可以改寫為

式中,κp為胞外聚合物的Ginzburg-Landau常數,κb為細菌的Ginzburg-Landau常數.將式(7)代入式(6)可得

由廣義Flory-Huggins混合自由能理論[19]可知:

式中,γ11,γ12分別為胞外聚合物和細菌畸變能系數,γ2為混合自由能系數,T為溫度.

1.2 多組元細菌生物膜系統的連續方程

連續方程是質量守恒定律對于運動流體的表達,由歐拉連續方程的微分形式可得

在黏性流體中,?·v表示流體速度的散度.由式(12)可知,細菌生物膜中胞外聚合物、細菌和等效溶劑的體積平均速度v為無散度平均速度場.

1.3 胞外聚合物、細菌和等效溶劑的輸運方程(相場控制方程)

胞外聚合物、細菌和等效溶劑的傳質過程可以運用多元相場理論來描述.由Hohenberg等[20]提出的相場控制方程可得如下方程.

(1)胞外聚合物的輸運方程為

式中,gp為胞外聚合物的生長率. (2)細菌的輸運方程為

式中,gb為細菌的生長率.

(3)等效溶劑的輸運方程為

式中,ζp為胞外聚合物和等效溶劑之間的體積摩擦阻力系數,ζb為細菌和等效溶劑之間的體積摩擦阻力系數,且二者分別是胞外聚合物、細菌和等效溶劑體積分數的函數,即ζp=ζp(?p,?b,?s),ζb=ζb(?p,?b,?s).分別取ζp=ζ0?p(1-?p),ζb=ζ0?b(1-?b),則胞外聚合物的輸運方程可改寫為

細菌的輸運方程可改寫為

等效溶劑的輸運方程可改寫為

1.4 營養物質的輸運方程(溶質場控制方程)

營養物質的存在是細菌生物膜生長的必要條件,外界營養物質在滲透壓的作用下透過細菌生物膜傳輸到膜內被細菌所消耗,為其提供充足的生長條件.營養物質在細菌生物膜系統中的傳質過程可以運用對流-擴散反應方程[21]來描述,即

式中,c為營養物質的濃度,Ds為營養物質的擴散系數,gc為營養物質的消耗率.

1.5 多組元細菌生物膜中各組元的生物反應動力學

Monod方程常被用來描述生物的生長速率與基質濃度的關系[22].因此,運用Monod方程可以得到胞外聚合物、細菌生長率以及營養物質消耗率的表達式.

(1)胞外聚合物的生長率為

式中,r1為胞外聚合物的最大生長率,r2為胞外聚合物的最大消耗率,K3,K4,K5均為半飽和常數,ε1和ε2均為反應動力學方程的尺度參數.

(2)細菌的生長率為

式中,Cb為細菌的最大生長率,ε為反應動力學方程的尺度參數,?b,max為細菌體積分數的最大容量,Kc為半飽和常數.

(3)營養物質的消耗率為

式中,Ac為營養物質的最大消耗率,K1為半飽和常數.

1.6 本構方程

本構方程是反映物質宏觀性質的數學模型,主要描述應力和應變率或應力張量與應變張量之間的函數關系.本工作將細菌生物膜中的胞外聚合物、細菌和等效溶劑均視為不可壓縮的廣義黏性牛頓流體,因此運用廣義牛頓流體的本構方程可以分別描述胞外聚合物、細菌和等效溶劑的應力.由于本工作主要考慮細菌生物膜在平均速度場中的生長,因此胞外聚合物、細菌和等效溶劑的應變率張量將運用平均應變率張量表示.由廣義牛頓流體的本構方程可得

2 數值求解與模擬

2.1 一維動力學模型

由于細菌生物膜的形狀比較復雜,大體呈“蘑菇”狀,這對細菌生物膜生長動力學方程邊界條件的確定以及方程的數值求解帶來困難.經實驗測定,大部分細菌生物膜的厚度在1～2 mm范圍[23]內,因此,細菌生物膜生長動力學方程的生長尺度較小,可以將其截面形狀等效為Hz×L矩形,如圖1所示.

圖1 多組元細菌生物膜的截面形狀等效圖F ig.1 Equivalent figure ofmulticomponent bacterial biofi lm’s cross section shape

本工作主要研究多組元細菌生物膜沿載體垂直向上方向(z方向)的生長,即確定多組元細菌生物膜的生長動力學方程在z∈[0,Hz]區域的數值解.細菌生物膜內生物質在x和y方向的生長和物質傳輸可以認為是均勻進行的,因此多組元細菌生物膜生長動力學方程的所有變量均是(z,t)的函數,即?p=?p(z,t),?b=?b(z,t),vx=vx(z,t),c=c(z,t).

基于此,多組元細菌生物膜生長動力學方程可以簡化為一維動力學方程,即

多組元細菌生物膜一維生長動力學方程的邊界條件為

多組元細菌生物膜一維生長動力學方程的初始條件為

表1列出了多組元細菌生物膜生長動力學數值模擬所需的物性參數.

表1 多組元細菌生物膜系統的物性參數Tab le 1 Physical parameters ofmulticomponent bacterial biofi lmsystem

2.2 多組元細菌生物膜生長動力學方程的求解

在對多組元細菌生物膜生長動力學進行模擬時,關鍵步驟是對動力學方程進行數值求解.數值計算的實質就是將動力學方程對時間和空間進行一定形式的離散化處理,然后結合方程的初始條件及邊界條件對離散后的方程進行迭代求解.本工作基于Crank-Nicolson有限差分格式求解細菌生物膜生長動力學方程.該方法是屬于對偏微分方程求解的加權隱式方法中的一種[24].多組元細菌生物膜生長動力學模型求解系統的前處理和數值計算程序實現的流程如圖2所示.

圖2 多組元細菌生物膜生長動力學模型求解系統的前處理和數值計算程序流程圖Fig.2 Flow chart of themulticomponent bacterial biofi lmgrow ing dynamicalmodel’s solving systemof pretreatment processing and numerical calculation program

3 多組元細菌生物膜在靜態環境和動態環境中的生長

本工作結合多組元細菌生物膜生長動力學模型的數值模擬結果,分析比較了靜態和動態環境中細菌生物膜胞外聚合物的體積分數和營養物質濃度在24～96 h(1～4 d)的分布情況,以及兩種條件之間的差異.

圖3和4是不同時刻胞外聚合物和細菌的體積分數在靜態和動態環境中的分布曲線,其中初始時刻t=0.由圖可知,胞外聚合物和細菌的體積分數及其分布空間在兩種環境中均隨著時間的增加而不斷增大,但在載體附近,靜態環境中胞外聚合物和細菌的體積分數均大于動態環境,且二者的差值隨著時間的增加而增大.這是由于動態環境中胞外聚合物和細菌在對流-擴散的共同作用下傳質速度較快、傳質阻力小,胞外聚合物和細菌能迅速地向無膜區域傳輸,從而使得動態環境中載體附近的胞外聚合物和細菌的體積分數較小.動態環境中胞外聚合物和細菌的體積分數的分布空間均大于靜態環境,且二者的差值也是隨著時間的增加而增大.during diff erent time

圖3 不同時刻,胞外聚合物體積分數在靜態環境和動態環境中的分布對比曲線Fig.3 Distribution curves of the extracellu lar polymer’s volume fraction in the static and dynamic environments during diff erent time

圖4 不同時刻,細菌體積分數在靜態環境和動態環境中的分布對比曲線Fig.4 Distribution curves of the bacteria’s volume fraction in the static and dynamic environments

這說明動態環境中細菌生物膜厚度增長較快.通過對比可知,動態環境更適合細菌生物膜的生長,并且膜內物質的對流傳質過程是促進細菌生物膜生長的主要因素.

圖5不同時刻營養物質濃度在靜態環境和動態環境中的分布對比曲線Fig.5 Distribution curves of the nutrients’concentration in the static and dynamic environments during diff erent time

圖5 為不同時刻營養物質濃度在靜態環境和動態環境中的分布曲線,其中初始時刻t=0.由圖可知,兩種環境中載體附近營養物質濃度的分布曲線呈非線性,而遠離載體處的分布曲線可以近似認為呈線性,且直線的斜率隨著時間的增加而增大.同一時刻動態環境中營養物質消耗較多,營養物質的濃度下降較快,濃度差值隨著時間的增加而增大,且遠離載體處營養物質濃度梯度較大.通過對比可知,對流-擴散傳質過程促進營養物質的傳輸和消耗,使得細菌生物膜獲取的營養物質更多,生長速率更快.

4 結論

本工作主要運用多元相場理論控制方程,并結合物質連續方程和動量守恒方程建立了耦合流場的相場模型,并且添加溶質場控制方程和生物反應動力方程,對細菌生物膜生長動力學進行了綜合建模.對建立的動力學模型進行了一維數值模擬,得到如下結論.

(1)細菌生物膜在靜態和動態環境中生長時,胞外聚合物和細菌的體積分數,以及營養物質濃度在空間上的分布規律均是相似的.載體附近營養物質濃度呈非線性分布,而遠離載體處營養物質濃度呈近似線性分布.

(2)動態環境中胞外聚合物、細菌和營養物質通過對流-擴散作用進行傳質,加快了細菌生物膜的生長.動態環境更適合細菌生物膜的生長,并且膜內物質的對流傳質過程是促進細菌生物膜生長的主要因素.

(3)本工作模型在已有模型的基礎上進行改善,考慮了滲透壓、摩擦阻力等多種因素對生物膜生長的影響,且與現實生活中生物膜生長情況以及國內外的實驗結果相一致.這也充分反映了本模型的可行性.

[1]楊朵,張正.細菌生物膜及其相關研究進展[J].中國實驗診斷學,2007,11(10):1416-1422.

[2]MARTIN K,MARTINAH,H ANS-J OACHMI B,et al.Automated con focal laser scanningmicroscopy and semia-automatedmiage processing for analysisofbiofi lms[J].Appl Environ Microbiol,1998, 64:4115-4124.

[3]陳黎明,柴立和.生物膜動力學的研究現狀與展望術[J].力學進展,2005,35(3):411-416.

[4]尚倩倩,方紅衛,何國建.水利工程中的生物膜研究進展[J].中國科技論文在線,2010,5(7): 563-568.

[5]F LEMMING H C,W INGENDER J,MAYER C,et al.Cohesiveness in biofi lmsmatrix polymers[J]. J Theor Biol,2009,253:689-708.

[6]任偉.不同清洗劑及其不同作用時間對內鏡生物膜清除的實驗研究[D].廣州:南方醫科大學,2011.

[7]李鵬.尾端縫線雙J管細菌學檢測的臨床意義[D].蘇州:蘇州大學,2012.

[8]ANGUIGE K,K ING J R.Mathematical modeling of therapies targeted at bacterial quorumsensing[J].Mathematical Biosciences,2004,19:239-283.

[9]李亮亮,徐向榮,周濤,等.生物膜數學建模的發展[J].廣東化工,2013,40(19):64-65.

[10]MOELANS N,B LANpAIN B,W OLLANTS P．An introduction to phase-field modeling of microstructure evolution[J]．Computer Coupling of Phase Diagrams and Thermochemistry,2008,32: 268-294．

[11]呂冬蘭.流場作用下合金枝晶生長的相場法數值模擬[D].南昌:南昌航空大學,2009.

[12]C OGAN N G,K EENER J P.The role of the biofi lmmatrix in structural development[J].Math Med Biol,2004,21:147-166.

[13]C OGAN N G,K EENER J P.Channel formation in gels[J].SIAMJ Appl Math,2004,65:1839-1854.

[14]Z HANG T,C OGAN N C,W ANG Q.Phase-field models for biofi lmsⅠ.Theory and one dimensional simulations[J].SIAMJ Appl Math,2008,69:641-669.

[15]Z HANG T,C OGAN N C,W ANG Q.Phase-field models for biofi lmsⅡ.2D numerical simu lations of biofi lm-flow interactions[J].CommComput Phys,2008(4):72-101.

[16]B RANDON L,W ANG Q,Z HANG T Y.Amulticomponent model for biofi lm-d rug interaction[J]. DCDS-B,2012,15(2):417-456.

[17]C HEN C,W ANG Q.3-D numericalsimulations ofbiofi lmflow s[J].East Asian Journalon Applied Mathematics,2011(1):197-214.

[18]AKIRAO.Phase transition dynamics[M].London:Cambridge University Press,2008.

[19]C HARLES W,E GBERT H,D ALE K,et al.How myxobacteria glide[J].Current Biology,2002, 12(6):369-377.

[20]H OHENBERG P C,H ALpERIN B I.Theory of dynamic critical phenomena[J].Rev Mod Phys, 1977,49:435-456.

[21]周光炯,嚴宗毅,許世雄,等.流體力學[M].北京:高等教育出版社,2000.

[22]W ANG Q,Z HANG T Y.Review ofmathematicalmodels for biofi lms[J].Solid State Communications,2010,150:1009-1022.

[23]AKIRAO.Theory of phase transition in polymer gels[D].Kyoto:Kyoto University,1992:6-10.

[24]J OHN D A.計算流體力學基礎及其應用[M].吳頌平,劉趙淼,譯.北京:機械工業出版社,1995.

Bacterial biofi lmgrow th dynamicalmodeling and numerical simu lation based on mu ltiple-phase field theory

ZHANG Shumei1,XU Xiangrong1,XU Hao2,ZHOU Tao1,LIShuang1,LIYan1
(1.School of Mechanical Engineering,Anhui University of Technology,Maanshan 243032, Anhui,China; 2.School of Mathematical Science and Engineering,Anhui University of Technology, Maanshan 243032,Anhui,China)

Bacteria biofi lms are w idespread in nature and industrial environments.Abacteria biofi lmis regarded as an incompressible viscous fluid composed ofmultiplecomponents.The theory of multiple-phase field,material continuity equation,momentumconservation equation,field control equations and biological reaction dynamic equation are combined tomodelbacterialbiofi lmgrow th dynamics.Anumericalsolution of themultiplecomponent bacterial biofi lmgrow th dynamical model is used to describe the grow th of bacteria,the law of nutriment’s consumption,and extracellular polymers in static and dynamic environments.

bacteria biofi lm;phase field theory;dynamicalmodel;numerical simu lation

Q 66

A

1007-2861(2017)04-0563-12

DO I:10.12066/j.issn.1007-2861.1690

2015-06-29

國家自然科學基金資助項目(31300125)

徐浩(1982—),男,講師,研究方向為生物膜、生物計算.E-mail:ahaxuo@ahut.edu.cn