薄興野，謝瑞潔，陳 敏，王艷東，張如意，劉振宇

(1.山西農(nóng)業(yè)大學(xué)信息科學(xué)與工程學(xué)院，山西 晉中 030801；2.山西農(nóng)業(yè)大學(xué)農(nóng)業(yè)工程學(xué)院，山西 晉中 030801)

0 引言

微生物個(gè)體相對(duì)簡(jiǎn)單，但其整個(gè)生長(zhǎng)過程表現(xiàn)出復(fù)雜的行為且數(shù)量龐大，同時(shí)微生物生長(zhǎng)和繁殖表現(xiàn)出種群效應(yīng)[1-2]。當(dāng)前，描述微生物生長(zhǎng)通常使用生長(zhǎng)動(dòng)力學(xué)方程，通常為非線性數(shù)學(xué)模型，是一類確定性的機(jī)理模型[3-4]。這類模型不能表現(xiàn)出微生物生長(zhǎng)過程中存在的多樣性、隨機(jī)性與不確定性，以及對(duì)初始狀態(tài)的敏感性。元胞自動(dòng)機(jī)是解決這些問題的有效工具。

元胞自動(dòng)機(jī)是由無限個(gè)元胞組成的世界，所有元胞都處于某種狀態(tài)并具有相似的形狀，在每次演化過程中，每個(gè)元胞按照相同規(guī)則，根據(jù)本身狀態(tài)和鄰域狀態(tài)改變狀態(tài)。趙峰等[5]使用二維元胞自動(dòng)機(jī)模擬了菌落生長(zhǎng)的分形結(jié)構(gòu)。李向陽(yáng)[6]使用莫諾德方程推導(dǎo)的差分方程模擬了微生物細(xì)胞的生長(zhǎng)。PICIOREANU C等[7]通過使用二維元胞自動(dòng)機(jī)與三維元胞自動(dòng)機(jī)模擬生物膜在凝膠珠生長(zhǎng)，說明三維元胞自動(dòng)機(jī)能更好模擬微生物生長(zhǎng)。目前，由于三維元胞自動(dòng)機(jī)設(shè)計(jì)復(fù)雜、運(yùn)行優(yōu)化難度大，研究者多用二維元胞自動(dòng)機(jī)代替，其不能真實(shí)模擬微生物在生長(zhǎng)中與環(huán)境的相互作用。三維元胞自動(dòng)機(jī)的使用會(huì)將微生物生長(zhǎng)過程中復(fù)雜機(jī)理、動(dòng)態(tài)行為，以及演化的多樣性、隨機(jī)性、不確定性和非線性等特點(diǎn)表現(xiàn)出來。同時(shí)，建立三維元胞自動(dòng)機(jī)模型為研究微生物生長(zhǎng)中與環(huán)境、與個(gè)體之間存在的依賴、競(jìng)爭(zhēng)、關(guān)聯(lián)等復(fù)雜相互作用提供了新思路[8]。

本研究基于Gompertz方程建立的差分方程，使用三維元胞自動(dòng)機(jī)模擬微生物生長(zhǎng)過程，根據(jù)設(shè)計(jì)規(guī)則編寫Python程序，并進(jìn)行仿真試驗(yàn)。仿真試驗(yàn)結(jié)果表明，三維元胞自動(dòng)機(jī)能較好地描述微生物生長(zhǎng)過程，并將生長(zhǎng)過程可視化表現(xiàn)，對(duì)使用三維元胞自動(dòng)機(jī)模擬微生物生長(zhǎng)進(jìn)行了新探索。

1 元胞自動(dòng)機(jī)模型構(gòu)建

采用接近微生物生長(zhǎng)的三維概率性元胞自動(dòng)機(jī)，每個(gè)元胞代表微生物存在所需三維空間，定義了微生物生長(zhǎng)過程中相關(guān)的狀態(tài)。元胞自動(dòng)機(jī)通過制定合理規(guī)則，使各個(gè)元胞在不同狀態(tài)值中不斷轉(zhuǎn)化，模擬微生物在一段時(shí)間內(nèi)的生長(zhǎng)。元胞自動(dòng)機(jī)微生物生長(zhǎng)模型的5個(gè)元祖，有

CA={Cell，State，Space，Neighborhood，Rules}

(1)

式中Cell——元胞，是元胞自動(dòng)機(jī)的基本元素

State——元胞狀態(tài)集合

Cellspace——元胞空間，代表有限空間

Neighborhood——鄰域，代表元胞規(guī)則作用的域

Rules——元胞規(guī)則，是元胞自動(dòng)機(jī)的核心，規(guī)定元胞狀態(tài)的演化

1.1 離散時(shí)間

元胞自動(dòng)機(jī)要對(duì)微生物生長(zhǎng)時(shí)間進(jìn)行離散化，對(duì)于每一個(gè)離散時(shí)間

t=kT0

k={0，1，…，N}

式中T0——時(shí)間間隔，即元胞自動(dòng)機(jī)每次演進(jìn)所代表的真實(shí)時(shí)間

N——最大演進(jìn)次數(shù)

1.2 元胞與元胞空間

假設(shè)元胞自動(dòng)機(jī)的元胞是邊長(zhǎng)為L(zhǎng)的微小正方體。隨機(jī)選取一個(gè)的邊長(zhǎng)為1 000L的正方體空間中微生物的生長(zhǎng)。根據(jù)換算可得元胞空間Space={Sijk|i，j，k∈{0，1，2，…，1 000}}，元胞空間中的元胞總數(shù)為1 000×1 000×1 000。

1.3 元胞鄰域

元胞鄰域是元胞規(guī)則的定義域，模型選用三維Moore型鄰域，定義如下

N={Vi=(Vix，Viy，Viz)||Vix-Vox|≤1，|Viy-Voy|≤1，|Viz-Voz|≤1，(Vix，Viy，Viz∈Z3)}

(2)

式中Vox，Voy，Voz——中心元胞坐標(biāo)

Vix，Viy，Viz——鄰居元胞的坐標(biāo)值

每個(gè)元胞有多達(dá)26個(gè)相鄰元胞，遠(yuǎn)超二維元胞自動(dòng)機(jī)。

1.4 元胞狀態(tài)

定義當(dāng)前元胞Cell(x，y，z)在t時(shí)刻的狀態(tài)為Statexyz(t)，每個(gè)元胞只有兩種狀態(tài)，即Statexyz(t)=0表示在t時(shí)刻Cell(x，y，z)為空，Statexyz(t)=i表示在t時(shí)刻Cell(x，y，z)被微生物占據(jù)，i表示元胞中存在的微生物數(shù)量。模型通過各個(gè)元胞狀態(tài)值不斷轉(zhuǎn)化模擬微生物在三維空間中的生長(zhǎng)過程。

2 元胞自動(dòng)機(jī)規(guī)則設(shè)計(jì)

元胞自動(dòng)機(jī)的下一時(shí)刻狀態(tài)與上一時(shí)刻的鄰域元胞和當(dāng)前元胞狀態(tài)有關(guān)，元胞需要根據(jù)設(shè)計(jì)的元胞規(guī)則演化。

2.1 微生物游走規(guī)則

微生物存在一定程度的運(yùn)動(dòng)行為，這是因?yàn)槲⑸镞\(yùn)動(dòng)有助于增強(qiáng)趨化性介導(dǎo)的養(yǎng)分吸收[9]。為更好表現(xiàn)微生物在元胞空間中的運(yùn)動(dòng)，假設(shè)微生物在時(shí)間間隔T0中做b次布朗運(yùn)動(dòng)，模擬微生物在空間中的運(yùn)動(dòng)情況[10]。設(shè)計(jì)隨機(jī)游走規(guī)則為

foriinRange(Sxyz>0)forjinrange(b)Sij(t+1)=Stateij-1，Sij_(t+1)=Stateij_+1

(3)

式中Sxyz——當(dāng)前元胞狀態(tài)，取值≥0，代表當(dāng)前元胞中微生物數(shù)量

i——當(dāng)前元胞正在進(jìn)行布朗運(yùn)動(dòng)的微生物

j——布朗運(yùn)動(dòng)次數(shù)

Sij(t+1)——當(dāng)前元胞第i個(gè)微生物第j次進(jìn)行布朗運(yùn)動(dòng)前微生物所處元胞位置及其狀態(tài)變化規(guī)則

Sij_(t+1)——當(dāng)前元胞第i個(gè)微生物第j次布朗運(yùn)動(dòng)后微生物所處元胞位置及其狀態(tài)變化規(guī)則

2.2 微生物分裂規(guī)則

微生物生長(zhǎng)常用Gompertz方程表示，如下

L(t)=A+Ceap{-exp[-B(t-M)]}

(4)

式中L(t)——菌數(shù)對(duì)數(shù)值，log cfu/mL

A——菌數(shù)初始值，log cfu/mL

C——細(xì)菌最大生長(zhǎng)量，log cfu/mL

M——細(xì)菌在生長(zhǎng)率最大時(shí)所對(duì)應(yīng)的時(shí)間，h

t——生長(zhǎng)時(shí)間，h

B——M時(shí)刻生長(zhǎng)率，log cfu/(mL·h)

推導(dǎo)出差分方程，即

x(t+T0)-x(t)=10L(t+T0)-10L(t)

(5)

其中x(t)即每毫升所含微生物數(shù)量，為元胞自動(dòng)機(jī)的概率演化提供指導(dǎo)。設(shè)計(jì)規(guī)則如下

IfSxyz(t)=1，then

(6)

式中x、y、z——當(dāng)前元胞坐標(biāo)值

Statexyz(t+1)——元胞在t+1時(shí)刻狀態(tài)

Statexyz——當(dāng)前元胞在t時(shí)刻狀態(tài)值

N(x，y，z)——t時(shí)刻元胞鄰域中微生物的數(shù)量

Pt0——當(dāng)前元胞在t時(shí)刻分裂的概率，為0～1隨機(jī)數(shù)

Pt1——當(dāng)前元胞在t時(shí)刻的概率閾值

概率閾值定義如下

Pt1=dx/(x(t)×(1-Pc))

(7)

式中dx——t至t+1時(shí)刻微生物的理論最大增長(zhǎng)量

x(t)——t時(shí)刻理論生物量

假設(shè)當(dāng)前元胞鄰域中元胞狀態(tài)非0數(shù)量為Nnum，則Pc為模型中每個(gè)元胞>Nnum或元胞狀態(tài)>1時(shí)，元胞中的微生物數(shù)量與模型中當(dāng)前時(shí)刻微生物數(shù)量之間的比值。

3 試驗(yàn)與仿真

3.1 試驗(yàn)材料和方法

(1)菌種：試驗(yàn)所用大腸桿菌由山西農(nóng)業(yè)大學(xué)生命科學(xué)學(xué)院提供。

(2)藥品：酵母提取物(山東玉寶生物科技有限公司)；胰蛋白胨(武漢吉業(yè)升化工有限公司)；氯化鈉[西格瑪奧德里奇(上海)貿(mào)易有限公司]；瓊脂(廣州宏程生物科技有限公司)、蛋白胨(PEP，山東玉寶生物科技有限公司)。

(3)蛋白胨溶液：使用BS224S電子分析天平稱取0.2 g蛋白胨放入300 mL三角錐形瓶中，用100 mL量筒經(jīng)兩次量取200 mL蒸餾水加入錐形瓶后，手持玻璃棒進(jìn)行攪拌，直至觀察到瓶中液體清澈，無固體粉末沉淀，再進(jìn)行封口后放入高壓滅菌鍋滅菌，取出后放入雙人單面凈化工作臺(tái)中備用。

(4)微生物檢測(cè)：通過平板計(jì)數(shù)法對(duì)微生物進(jìn)行檢測(cè)。待檢測(cè)的大腸桿菌稀釋到適合濃度時(shí)，將其涂于培養(yǎng)基上，培養(yǎng)出細(xì)菌菌落后，計(jì)算生長(zhǎng)出的細(xì)菌菌落數(shù)，得到單位體積細(xì)菌原始菌落數(shù)為計(jì)算出的菌落數(shù)乘以稀釋倍數(shù)。為了計(jì)數(shù)準(zhǔn)確，需要將待測(cè)細(xì)菌稀釋濃度進(jìn)行嚴(yán)格控制，確保平板上的細(xì)菌菌落數(shù)在30～300范圍[11]。自然條件下在蛋白胨溶液中接種大腸桿菌培養(yǎng)12 h。

3.2 仿真結(jié)果

仿真試驗(yàn)利用Python語(yǔ)言，按照設(shè)計(jì)的結(jié)構(gòu)編制元胞自動(dòng)機(jī)程序，模擬以上試驗(yàn)中大腸桿菌在蛋白胨溶液中生長(zhǎng)。演化程序包括以下功能。利用Gompertz方程計(jì)算出每一時(shí)刻的微生物數(shù)量，以及每次迭代的概率閾值，根據(jù)比例計(jì)算出在三維空間中微生物初始數(shù)量，并隨機(jī)生成初始位置。在三維空間中，自動(dòng)描繪微生物生長(zhǎng)的狀態(tài)圖。計(jì)算每一時(shí)刻模型中微生物數(shù)量，并按比例計(jì)算出1 mL中的數(shù)量，描繪曲線圖。

根據(jù)實(shí)際測(cè)得數(shù)據(jù)計(jì)算得出Gompertz方程的參數(shù)：A=7.32，C=2.585，B=0.417 76，M=2.096 9。設(shè)置時(shí)間間隔t=0.5 h，鄰域中微生物閾值Nnum=15，布朗運(yùn)動(dòng)次數(shù)b=2。將值賦給程序，運(yùn)行結(jié)果如圖1所示。

圖1 元胞自動(dòng)機(jī)運(yùn)行結(jié)果Fig.1 Cellular automata running results

其中圖1a為元胞自動(dòng)機(jī)初始轉(zhuǎn)態(tài)，圖1b為迭代24次的程序運(yùn)行結(jié)果，可以看到微生物在經(jīng)過布朗運(yùn)動(dòng)與分裂后，微生物有明顯數(shù)量變化并向周圍進(jìn)行擴(kuò)散，基本可以模擬微生物生長(zhǎng)過程中的分裂與運(yùn)動(dòng)行為。

3.3 試驗(yàn)驗(yàn)證

微生物生長(zhǎng)過程中活菌總數(shù)直接表現(xiàn)了微生物生長(zhǎng)狀態(tài)。圖2中帶點(diǎn)虛線為實(shí)際測(cè)出大腸桿菌菌數(shù)值(cfu/mL)，虛線為根據(jù)Gompertz方程繪制出的曲線圖，帶點(diǎn)實(shí)線為元胞自動(dòng)機(jī)模擬仿真大腸桿菌在蛋白胨溶液中菌數(shù)值(cfu/mL)。由圖2可以初步看出，元胞自動(dòng)機(jī)的設(shè)計(jì)與參數(shù)選擇基本合理，可以與試驗(yàn)中實(shí)際大腸桿菌的生長(zhǎng)菌數(shù)值吻合，并提供了三維可視化的參考。同時(shí)模型在初始時(shí)間與Gompertz方程基本相似，隨著時(shí)間的進(jìn)行模型活菌總數(shù)雖與Gompertz方程出現(xiàn)了較大的分離，但還是較為接近真實(shí)值。

圖2 仿真過程中大腸桿菌總數(shù)隨時(shí)間變化Fig.2 Total number of Escherichia coli changes over time during simulation

由于微生物數(shù)量龐大，通常會(huì)使用活菌對(duì)數(shù)值來表示微生物生長(zhǎng)，大腸桿菌與元胞自動(dòng)機(jī)仿真活菌數(shù)對(duì)數(shù)值如表1所示。以培養(yǎng)時(shí)間t(h)為橫坐標(biāo)，活菌數(shù)值log(cfu/mL)為縱坐標(biāo)，大腸桿菌生長(zhǎng)曲線擬合結(jié)果如圖3所示。其中帶點(diǎn)虛線為數(shù)據(jù)真實(shí)值表示，帶點(diǎn)實(shí)線為仿真生長(zhǎng)擬合曲線。

表1 大腸桿菌與元胞自動(dòng)機(jī)仿真活菌數(shù)對(duì)數(shù)Tab.1 Logarithm of number of live bacteria simulated by Escherichia coli and cellular automata

圖3 大腸桿菌擬合結(jié)果Fig.3 Escherichia coli fitting results

進(jìn)行擬合計(jì)算，其精度為R2=0.985 68?？梢娖鋽M合精度較高，元胞自動(dòng)機(jī)對(duì)大腸桿菌的生長(zhǎng)仿真較好。在培養(yǎng)時(shí)間為0時(shí)，數(shù)據(jù)差值相對(duì)較大，這是由于元胞自動(dòng)機(jī)隨機(jī)生成的初始狀態(tài)導(dǎo)致，隨著時(shí)間推移，仿真得到的大腸桿菌對(duì)數(shù)值接近真實(shí)值。

4 結(jié)論

基于微生物生長(zhǎng)機(jī)理建立了模擬微生物生長(zhǎng)的三維概率元胞自動(dòng)機(jī)，并得到以下結(jié)論。

(1)采用布朗運(yùn)動(dòng)方法可有效地在元胞自動(dòng)機(jī)中模擬微生物的運(yùn)動(dòng)特性，微生物的運(yùn)動(dòng)能力隨離散時(shí)間t內(nèi)的布朗運(yùn)動(dòng)次數(shù)增大而增強(qiáng)，通過合理設(shè)置運(yùn)動(dòng)次數(shù)能較好模擬實(shí)際微生物運(yùn)動(dòng)。

(2)使用Gompertz方程推導(dǎo)的差分方程制定合理元胞規(guī)則來模擬微生物的分裂，在各個(gè)離散時(shí)間內(nèi)元胞自動(dòng)機(jī)仿真得到的微生物數(shù)量與實(shí)驗(yàn)數(shù)據(jù)接近。

(3)本次仿真通過設(shè)置元胞自動(dòng)機(jī)時(shí)間間隔為t=0.5 h，鄰域中微生物閾值Nnum=15，布朗運(yùn)動(dòng)次數(shù)b=2，Gompertz方程參數(shù)A=7.32，C=2.585，B=0.417 76，M=2.096 9，模擬大腸桿菌在蛋白胨溶液生長(zhǎng)，擬合精度較高，并表現(xiàn)了微生物在三維空間中的分布與微生物生長(zhǎng)的隨機(jī)性與不確定性。

本研究未考慮微生物死亡狀態(tài)并假定生長(zhǎng)過程中營(yíng)養(yǎng)物質(zhì)充分，研究接下來會(huì)將微生物死亡狀態(tài)加入探討，同時(shí)引入營(yíng)養(yǎng)物質(zhì)不均衡狀況并建立基于基因表達(dá)的元胞規(guī)則。