基于多元線性回歸的真菌分解過程分析

馬秋敏，楊 洪，2，于亦飛

（1.江蘇海洋大學(xué) 理學(xué)院，江蘇 連云港 222005；2.江蘇省海洋資源開發(fā)研究院（連云港），江蘇 連云港 222005）

碳循環(huán)描述了碳在整個地球化學(xué)循環(huán)中的交換過程，是自然界中物質(zhì)和能量循環(huán)的重要組成部分。碳循環(huán)的一部分包括化合物的分解，允許碳被更新并以其他形式使用。這部分過程的一個關(guān)鍵組成部分是植物材料和木質(zhì)纖維的分解，而分解木質(zhì)纖維的關(guān)鍵因素是真菌。2021年美國大學(xué)生數(shù)學(xué)建模競賽A題要求用數(shù)學(xué)建模的方法研究真菌分解木質(zhì)纖維素纖維過程及其相互作用問題。本文擬采用修正Frisch法、最小二乘法、微分方程等研究真菌分解木質(zhì)纖維素纖維的過程，以及不同種類真菌間的相互作用、在環(huán)境影響下的真菌生存能力。

1 “真菌分解速率”模型的建立與分析

1.1 模型建立

基于題中給出的假設(shè)條件，以及影響真菌的外部影響因素（濕度、菌種間的影響、pH值等）和內(nèi)部影響因素（真菌的內(nèi)在性狀，如生長速度、耐濕性等），并且假設(shè)存在下列相關(guān)系數(shù)a，b，c，d，構(gòu)建真菌分解速率y與影響因素間的多元線性回歸模型[1]：

其中，x1為真菌的耐濕性，x2為真菌生長率，x3表示菌種間的影響，x4為pH值，y為真菌分解速率。

1.2 模型求解

由賽題資料可知，真菌分解速率y和生長率x2呈線性關(guān)系，真菌分解速率y與耐濕性x1呈正比關(guān)系。進一步研究真菌耐濕性x1、真菌生長率x2對真菌分解速率y的影響，根據(jù)生長率x2與分解速率y2的數(shù)據(jù)，利用Eviews軟件作x2與y2的散點圖如圖1所示。由圖1可見，x2與y2近似呈對數(shù)關(guān)系。

圖1 真菌生長率x2與分解速率y2關(guān)系

將真菌生長率x2與真菌分解速率y2進行線性回歸，運用最小二乘法得到估計結(jié)果如圖2。

圖2 x2與ln y2最小二乘估計結(jié)果

由圖2可知，

其中，括號內(nèi)的數(shù)字是t值。R2=0.94說明總離差平方和中有94%被樣本回歸直線解釋，只有6%未被解釋。可見回歸方程擬合結(jié)果較好，且t檢驗和F檢驗均符合要求，回歸方程顯著。

將（2）式化簡可得：

同理，得到耐濕性x1與分解速率y1的最小二乘估計結(jié)果如圖3。

圖3 x1與ln y1最小二乘估計結(jié)果

由圖3可得

其中括號內(nèi)的數(shù)字是t值。R2=0.97說明總離差平方和中有97%被樣本回歸直線解釋，只有3%未被解釋。表明回歸方程擬合結(jié)果較好，且t檢驗和F檢驗均符合要求，回歸方程顯著。

將（4）式化簡可得：

同時，在處理真菌耐濕性與生長率關(guān)系時，可發(fā)現(xiàn)不同真菌耐濕性和生長率大致呈線性關(guān)系（這里的線性關(guān)系可看作廣義的二元關(guān)系，只要兩個變量直接有關(guān)聯(lián)即可）。因此，在保持真菌分解速率不變的情況下將兩種因素進行合并，其擬合方程為：

將不同因素影響下的y1、y2進行加總，得出合并處理后的真菌分解速率模型為：

2 真菌間不同關(guān)系模型的建立與分析

2.1 競爭關(guān)系模型的建立與求解

（1）競爭關(guān)系標準模型如下：

其中：x5、x6分別為A、B菌種在t時刻的數(shù)量；rA、rB為A、B菌種的固定增長率；NA、NB是環(huán)境資源容許的菌種最大數(shù)量；σ1表示對于供養(yǎng)A菌種的食物量而言，單位數(shù)量B菌種（相對NB）的消耗為單位數(shù)量A菌種（相對NA）消耗的σ1倍；σ2表示對于供養(yǎng)B菌種的食物量而言，單位數(shù)量A菌種（相對NA）的消耗為單位數(shù)量B菌種（相對NB）消耗的σ2倍[2]。

（2）模型求解

對于兩種有競爭關(guān)系的真菌而言，當(dāng)為了生存爭奪有限的同一食物來源和生活空間時，競爭力較強的種群往往達到環(huán)境允許的最大數(shù)量。任取有競爭關(guān)系的A菌種和B菌種，其中，A菌種是生長緩慢的真菌菌株，更適應(yīng)在變化的環(huán)境下生存和生長，而B菌種生長速度較快，但在相同的環(huán)境變化中不那么茁壯。因此，A菌種固定增長率rA小于B菌種固定增長率rB，而A菌種的競爭能力高于B菌種，即σ1>σ2；在環(huán)境適宜時，取rA=0.2，rB=0.5；σ1=2，σ2=0.5，獲得A、B菌種的數(shù)量關(guān)系與相軌線如圖4。

由圖4可見，最后數(shù)值穩(wěn)定于x6=0，x5=10，即A菌種達到最大數(shù)量，B菌種滅絕。

圖4 A菌種與B菌種數(shù)量關(guān)系及相軌線

環(huán)境變化導(dǎo)致了A、B菌種的自然增長率的快速變化。取rA=0.4，rB=0.3，σ1=2，σ2=0.5，獲取A、B菌種的數(shù)量關(guān)系與相軌線如圖5。

圖5 自然增長率變化后的A、B菌種數(shù)量關(guān)系與相軌線

2.2 捕食關(guān)系模型的建立與求解

任取C菌種和D菌種，假設(shè)其是捕食關(guān)系，其中C菌種是捕食者，D菌種是被捕食者。C菌種主要以木質(zhì)纖維中的營養(yǎng)物質(zhì)和D菌種為食，且當(dāng)捕食D菌種到一定程度則不再捕食；D菌種主要以木質(zhì)纖維中營養(yǎng)物質(zhì)為食。

利用Volterra食餌-捕食者模型來解決以下問題：

其中：t為時間，x（t）、y（t）為D、C菌種在t時刻的數(shù)量；r為D菌種獨立生存時的增長率；d為C菌種獨自存在時的死亡率；a表示C菌種掠取D菌種的能力；b表示D菌種對C菌種的供養(yǎng)能力。

C菌種主要以木質(zhì)纖維中營養(yǎng)物質(zhì)和D菌種為食，當(dāng)環(huán)境快速波動時，仍可以D菌種為食生存；而D菌種主要以木質(zhì)纖維中營養(yǎng)物質(zhì)為食，當(dāng)環(huán)境快速波動時，則極易死亡。

取C菌種初始值為150，D菌種初始值為300，當(dāng)環(huán)境適宜時，取r=2，d=1，a=b=0.01代入上述方程作菌種數(shù)量關(guān)系和相軌線如圖6。

圖6 捕食關(guān)系的C、D菌種數(shù)量關(guān)系及相軌線

由圖6可見：C、D兩菌種相互制約，當(dāng)C菌種逐漸增多時，D菌種逐漸減少；當(dāng)C菌種增加到一定數(shù)量時，由于種群內(nèi)部競爭激烈，又導(dǎo)致其數(shù)量減少，從而使得D菌種天敵減少，進而導(dǎo)致D菌種數(shù)量增加；當(dāng)D菌種增加到一定數(shù)量時，種群間競爭加劇，從而又導(dǎo)致D菌種數(shù)量減少。如此循環(huán)，相互制約發(fā)展。

當(dāng)環(huán)境快速波動時，取r=2，d=1，a=b=0.05，代入方程得圖7。從圖7可以看到，C、D兩菌種仍滿足上述關(guān)系，且在環(huán)境惡劣時，其數(shù)量變化較為顯著。

圖7 D菌對C菌供養(yǎng)能力變化后C、D菌種的數(shù)量關(guān)系及相軌線

2.3 互利共生關(guān)系模型的建立與求解

對于相互依存的兩種真菌而言，相互依存而共生是其普遍現(xiàn)象。任取互利共生關(guān)系的E菌種和F菌種進行研究，其中，E菌種可單獨存在，F(xiàn)菌種為E菌種提供食物，同時不能離開E菌種單獨存活，且有助于E菌種的生長。

建立模型如下：

其中：x7、x8為E、F菌種在t時刻的數(shù)量；rE、rF為E、F菌種固定增長率；NE、NF是環(huán)境資源容許的E、F菌種最大數(shù)量；σ3表示對于供養(yǎng)E菌種的食物量而言，單位數(shù)量F菌種（相對NF）的消耗為單位數(shù)量E菌種（相對NE）消耗的σ3倍；σ4表示對于供養(yǎng)F菌種的食物量而言，單位數(shù)量E菌種（相對NE）的消耗為單位數(shù)量F菌種（相對NF）消耗的σ4倍。

在溫度適宜時，取rE=1.8，rF=1.5；σ3=0.1，σ4=3，作菌種數(shù)量關(guān)系與相軌線如圖8。

圖8 互利共生關(guān)系的E、F菌種數(shù)量關(guān)系與相軌線

從圖8可看出，當(dāng)t趨于無窮大時，兩種真菌的最大數(shù)量均增大且相互依存，在無外界因素劇烈變化的條件下，兩種真菌數(shù)量趨于穩(wěn)定。

當(dāng)環(huán)境快速波動時，取rE=1.8，rF=1.5；σ3=0.3，σ4=3，得到E菌和F菌的數(shù)量關(guān)系及相軌線如圖9。

圖9 σ3變化后E、F菌種數(shù)量關(guān)系及相軌線

從圖9可看出，在環(huán)境快速波動時，隨著σ3的增加，兩種真菌間共生率增大，表明其促進作用也增大，且兩種真菌數(shù)量隨時間增加達到平衡位置的值也增大。

3 不同環(huán)境影響下真菌分解速率模型

3.1 模型的建立

不同種類真菌的耐濕性和生活的最佳溫度都和真菌生長率有著密切關(guān)系，也決定了真菌最適宜的生存環(huán)境條件。針對這兩種因素建模，分析其最適宜的生存環(huán)境條件，最終得到真菌生長率與真菌所處環(huán)境的溫度和耐濕性的關(guān)系式。

3.2 模型求解

由賽題所給條件可得生長率x2、溫度x3與耐濕性x1的數(shù)據(jù)，利用Eviews軟件作x2與x1的散點圖如圖10。由圖10可見，x2與x1近似呈對數(shù)關(guān)系。

圖10 x1與x2散點圖

將真菌生長率x2與真菌分解速率x1進行線性回歸，運用最小二乘法得到估計結(jié)果如圖11。

圖11 ln x21與x1的最小二乘估計結(jié)果

由圖11可得

R2=0.99，Rˉ2=0.99，F(xiàn)=593.54，DW=1.69其中，括號內(nèi)的數(shù)字為t值。R2=0.99說明總離差平方和中有99%被樣本回歸直線解釋，只有1%未被解釋。可見回歸方程擬合結(jié)果較好，且t檢驗和F檢驗均符合要求，回歸方程顯著[3]。

將（11）式化簡可得：

同理可得圖12的散點圖。由圖12可見，x2與x3近似呈線性關(guān)系。

圖12 x2與x3散點圖

擬合得x22=-2.194952+0.397 525x3

由t檢驗可知，常數(shù)項不顯著，故舍去。再次擬合得：

按上述步驟檢驗，得知回歸方程顯著。

綜上可得，真菌生長率與溫度及耐濕性的關(guān)系式為：

假設(shè)溫度為5~35℃，且x3的系數(shù)為正，由式（15）可知，溫度對生長率始終起促進作用。為指數(shù)函數(shù)，且當(dāng)-0.500434+2.755192x1>0時，即x1>0.181633時，為單調(diào)遞增函數(shù)，即真菌的生長率隨著溫度與耐濕性的增大而增加。反之，若-0.500 434+2.755 192x1<0，則為單調(diào)遞減，且當(dāng)x1→∞時，x21→0，即真菌生長率只與溫度有關(guān)。

分析表明，干旱地區(qū)與半干旱地區(qū)真菌的耐濕性最弱，且干旱地區(qū)與半干旱地區(qū)晝夜溫差較大，因此真菌生長率較低。喬木林和熱帶雨林地區(qū)真菌的耐濕性最強，因此可視為只受溫度的影響，其生長率與溫度成正比關(guān)系。溫帶地區(qū)真菌耐濕性較強，且溫度適宜，生長率明顯高于干旱地區(qū)、半干旱地區(qū)及喬木林和熱帶雨林地區(qū)，適宜真菌生長[4]。