基于種群動態膜系統的大熊貓種群模型

田 昊，張葛祥*，榮海娜，Mario J．P REZ-JIM NEZ，Luis VALENCIA-CABRERA，陳 鵬，侯 蓉，齊敦武

(1．西南交通大學電氣工程學院，成都610031; 2．塞維利亞大學計算機科學與人工智能系，塞維利亞41012，西班牙;3．成都大熊貓繁育研究基地，成都610081; 4．四川省瀕危野生動物保護生物學重點實驗室，成都610081)

(*通信作者電子郵箱zhgxdylan@126．com)

0 引言

大熊貓是世界自然基金會的形象大使和當之無愧的世界生物多樣性保護的旗艦物種［1－6］。在研究人員從事與大熊貓活動相關的工作中，需要對大熊貓的生活習性、種群數量、種群年齡結構組成以及引起種群數量變化的各種因素進行分析和預測［7－9］。這些研究對大熊貓的保護工作具有重要意義［10－11］。

在生態系統種群動態性建模中最廣泛使用的方法是常微分方程法和偏微分方程法［12－13］，相關科研人員曾廣泛使用萊斯利矩陣對不同大熊貓種群進行動態預測［14－15］。常微分數學方程通過簡化的數學分析很容易用簡單的公式描述復雜系統，但該方法無法捕捉到空間動態性和隨機效應。偏微分方程建模方法雖克服了常微分方程法建模方法中不能捕獲生態系統空間動態性的缺點，但建模過程中存在非常復雜的數學分析。隨著生態系統復雜性不斷增加，采用常微分方程與偏微分方程方法建立的模型會變得更加難以解決問題，以至于這種建模框架在許多情況下不適用。

目前，關于大熊貓種群動態的研究文獻還比較少。文獻［16－17］利用計算機仿真分析了大熊貓種群生存力;文獻［18－19］分析討論了大熊貓種群動態與它們主要食物來源竹子之間的相互影響作用;文獻［20－21］通過計算機仿真研究了棲息地對大熊貓種群動態的影響;文獻［22－23］研究了汶川地震對大熊貓種群的影響;在文獻［24－25］中，相關研究人員分析了氣候變化對大熊貓種群及其保護政策的威脅和挑戰。這些研究都無法很好地反映大熊貓種群的空間動態性和隨機效應，國外研究人員已經將膜計算數據建模應用于生態系統種群動態性的研究中［26－31］。在國內，本文的實驗室團隊已經使用概率膜系統對成都大熊貓繁育研究基地的圈養大熊貓生態系統進行建模研究［32－33］，該模型粗略地考慮了大熊貓生命周期中的繁殖、死亡、進食、救護行為，繁殖模塊設為每年固定數量的新生大熊貓，并且假定大熊貓為單胞胎繁殖，救護模塊中一年實際救護數量最多為1只個體。根據對譜系數據更深入的研究發現大熊貓雙胞胎現象較為普遍，并且2005—2008年大熊貓的出生率較之后的年份明顯高出很多，查證大熊貓繁殖期后設置合理的大熊貓繁殖年齡，根據處于繁殖期的雌性大熊貓數量作為主要因素之一設計進化規則。救護模塊中加入同一年內出現兩只個體的可能性，最為重要的是隨著野放現象的出現，在計算模型中需要考慮加入野放模塊。而且根據概率膜系統的定義，其無法對兩個大熊貓種群之間的相互作用與交流進行數據建模，如若后期需要對多個種群進行建模研究，此時概率膜系統則無法實現。

在此提出使用種群動態膜系統對中國大熊貓保護研究中心的圈養大熊貓進行數據建模，新模型從個體行為上來模擬中國大熊貓保護研究中心圈養大熊貓的生態系統，更詳細地分析討論了2005—2016年的大熊貓譜系數據，充分考慮了可能會影響該種群動態性的相關因素。該模型更完整地體現出中國大熊貓保護研究中心圈養大熊貓生態系統的特點，更加符合該生態系統的實際情況，同時適用于接下來對中國大熊貓保護研究中心圈養大熊貓生態系統的研究。種群動態膜系統(Population Dynamics P system，PDP system)比微分方程法更復雜，但是使用的語言更接近專家語言(生態學家)。通俗地說，用PDP系統建模方法通過對象集和一些規則來控制它們的進化和相互影響來模擬大熊貓種群現實生活中的進化過程。在復雜系統建模中對象是最重要的組成部分(大熊貓、空間資源、食物等)，它們同時可以根據一系列的重寫規則來描述系統的動態性。每一個相關規則采用概率函數試圖捕捉自然中固有的隨機性。

1 譜系數據分析處理

本文以中國動物園協會已發布的全球圈養大熊貓譜系數據為依據，從中搜集統計中國大熊貓保護研究中心大熊貓種群數據，并以其作為研究對象。根據對大熊貓譜系數據中年齡結構、繁殖率和死亡率的分析考量，以及與大熊貓專家研討得到的相關意見，在大熊貓從出生到死亡的整個生命過程中，依照其不同時期的行為特點將其分為幼年期、亞成年期、青年期、中年期、中老年期和老年期共6個年齡階段，并確定繁殖期所處的年齡范圍。以6個年齡階段為基礎，追蹤每只大熊貓的行為，統計確定大熊貓在各階段的食物種類、進食數量、死亡率以及在繁殖期內的大熊貓繁殖率。

同時，由于人為因素的干擾，在此需要考慮野外大熊貓救護行為對大熊貓種群動態性產生影響的因素，根據往年救護大熊貓的數量、年齡和性別的統計和分析，確定每年可能救護的大熊貓數量和年齡分布。野外的雌性大熊貓由于生存環境等原因無法同時撫養雙胞胎幼體，之前都是考慮大熊貓是單胞胎的哺乳動物，但在圈養大熊貓生態系統中隨著生存環境的改善和人工繁殖技術的發展，雙胞胎個體都能得到良好的成長環境，故在此需要把雙胞胎的可能性納入模型。通過與中國大熊貓保護研究中心專家討論研究發現，不同時間段的大熊貓出生率有較大差別，2005—2008年出生率較高，之后隨著大熊貓種群數量的增長、年齡結構的穩定，出生率也隨之下降并穩定下來，故本文首次在種群動態膜系統的計算模型中加入分段式的進化規則，對2008年前后兩個時期的大熊貓種群參數進行不同的參數設置，使其更符合中國大熊貓保護研究中心圈養大熊貓生態系統的種群特點。圈養大熊貓除了保護種群這個重要原因，更是希望通過野放讓大熊貓回歸自然，復壯野外小種群。故隨著野化技術的成熟，從2012年開始選擇合適的大熊貓進行野化放歸。

模型中所需的參數如下所示:

i表示圈養大熊貓的性別，其中i=1表示雄性大熊貓，i=2表示雌性大熊貓。

ki，1表示圈養大熊貓年齡達到亞成年階段;ki，2表示圈養大熊貓年齡達到成年階段;ki，3表示圈養大熊貓年齡達到中年階段;ki，4，1表示圈養大熊貓年齡達到中老年階段;ki，4，2表示圈養大熊貓年齡達到老年階段;ki，5表示圈養大熊貓的最大年齡;ki，6表示圈養大熊貓處于幼年階段的死亡率;ki，7表示圈養大熊貓處于亞成年階段的死亡率;ki，8表示圈養大熊貓處于成年階段的死亡率;ki，9表示圈養大熊貓處于中年階段的死亡率;ki，10表示圈養大熊貓處于中老年階段的死亡率;ki，11表示圈養大熊貓處于老年階段的死亡率;ki，12表示圈養大熊貓進入繁殖期的年齡;ki，13表示圈養大熊貓結束繁殖期的年齡。

g1表示繁殖期大熊貓生育單胞胎的概率;g2表示繁殖期大熊貓生育雙胞胎的概率;g3表示一年中供給圈養大熊貓的竹筍量(kg);g4表示一年中供給圈養大熊貓的竹子量(kg);g5表示一年中供給圈養大熊貓的其他食物(例如水果、牛奶等)量(kg)。

fi，1表示每只幼年大熊貓一年中需要消耗的竹筍總量;fi，2表示每只幼年大熊貓一年中需要消耗的竹子總量;fi，3表示每只幼年大熊貓一年中需要消耗的其他食物總量;fi，4表示每只亞成年大熊貓一年中需要消耗的竹筍總量;fi，5表示每只亞成年大熊貓一年中需要消耗的竹子總量;fi，6表示每只亞成年大熊貓一年中需要消耗的其他食物總量;fi，7表示每只成年大熊貓和老年大熊貓一年中需要消耗的竹筍總量;fi，8表示每只成年大熊貓和老年大熊貓一年中需要消耗的竹子總量;fi，9表示每只成年大熊貓和老年大熊貓一年中需要消耗的其他食物總量。

cmin表示每年救護大熊貓的最少數量;cmax表示每年救護大熊貓的最大數量;cmaxage表示救護大熊貓的最大年齡;pcc表示一年中救護c只大熊貓的概率;pgi表示救護大熊貓性別為i的概率;paj表示救護大熊貓年齡為j的概率;pdd表示一年中野放d只大熊貓的概率;psi表示野放大熊貓性別為i的概率;pyj表示野放大熊貓年齡為j的概率。

2 種群動態膜系統

2．1 種群動態膜系統定義

膜計算是自然計算的新分支，由歐洲科學院院士、羅馬尼亞科學院院士Paun提出［34］。膜計算受生物細胞的啟發，旨在從細胞組織的結構和功能中抽象出計算模型［35－36］。膜系統(計算模型)由膜結構、初始對象集和進化規則集三要素構成，膜系統中膜結構分隔了對象和規則所處的區域，對象按照所在區域的規則不斷進化，進化規則以極大并行的方式執行。膜系統運行開始前的狀態稱為初始格局，膜系統自初始格局開始根據不同的進化規則從一個格局進化到下一個格局，直到沒有規則可以執行時系統終止，此時的膜系統狀態稱為終止格局，同時終止格局時可以指定膜系統里不同區域中的對象數目作為膜系統的輸出結果。多環境概率膜系統是在膜系統的基礎上引入概率函數而得到的一種新的膜系統概念，種群動態膜系統是多環境膜系統的一個變種［37］。

一個度為(q，m)，同時滿足q，m≥1且T≥1種群動態膜系統可形式化定義為:

其中:

G=(V，S) 是一個 V={e1，e2，…，em} 的有向圖。

Γ和Σ是字母表。Γ是包含有限對象的非空字母表，由種群動態膜系統所有區域中的所有對象組成;Σ代表字母表Γ中可以出現在所有環境中的對象組成的字母表。

Τ表示生態系統進化過程被仿真的年數。

RE表示有限規則集，形如

其中:x，y1，y2，…，yh∈ Σ，(ej，ejl) ∈ S，1 ≤l≤ h，1 ≤ k≤ n，且p，p'是定義域為{1，2，…，T}的可計算的函數。

u是集合{1，2，…，q}×{0，+，－}中元素單射標號的有q個節點的根樹，如果膜標號為(i，α)那么這個膜將表示為同時該膜的標號為i控制電荷為α，樹的根節點設為1作為關聯標簽。

對于每一個r∈R和1≤j≤m，fr，j是一個定義域為{1，2，…，T}的可計算函數。

對于每一個 i，j(1 ≤ i≤ q，1 ≤ j≤ m)，Mi，j是字母表 Γ中的有限多重集。

對于每一個j，1≤j≤m，Ej是字母表Σ中的有限多重集。

2．2 計算模型設計

整個模型由初始模塊、繁殖模塊、救護模塊、進食模塊、死亡模塊、野放模塊和更新模塊6個部分組成。

計算模型考慮使用度為(2，1)時間單元為T的種群動態膜系統模型

定義如下:

G是一個空表。

工作字母表Γ為:

其中:Xi，j，y表示繁殖模塊之前性別為i年齡為j的大熊貓個體，Yi，j，y表示進入死亡模塊性別為i年齡為j的大熊貓個體，Zi，j，y表示存活性別為 i年齡為 j的大熊貓個體，Wi，j，y表示進食模塊之后性別為i年齡為j的大熊貓個體，Ci，j，y表示救護性別為i年齡為j的大熊貓個體，y表該個體目前所處的仿真周期。對象S，B，O表示不同類型的食物:分別為竹筍、竹子和其他食物。對象F，A1，G1是輔助符號分別表示如下:F是一個符號用來作為在每個周期初始產生新的食物，A1是一個符號用來作為觸發救護個體的行為;G1是一個符號用來作為觸發野放個體的行為，下標1表示目前為第一個仿真周期。

T表示在進化生態系統中仿真的年數。

u=［［］1］2表示采用兩層膜結構。

M1，1和M2，1是Γ中的有限多重集，描述區域u中的初始對象:

qi，j表示性別為i年齡為j的初始大熊貓數量

M2，1={F A1G1}

1)初始規則。

為大熊貓提供食物

2)繁殖規則。

未達到繁殖期年齡的個體

3 實驗仿真結果

為了驗證上述設計的中國大熊貓保護研究中心圈養大熊貓種群動態膜系統計算模型的有效性和正確性，以2005—2016年全球大熊貓譜系數據為依據，分別統計了各年份中國大熊貓保護研究中心圈養大熊貓的數量，用實際數據和實驗仿真結果進行對比。把2005年各年齡段的大熊貓數量作為初始參數輸入，2006—2016年各年齡段的大熊貓數量作為結果輸出。仿真實驗所用計算機為戴爾 Inspiron N5110，2．10 GHz，3．16 GB，使用 MeCoSim 仿真軟件作為仿真平臺。

實驗中的參數設置為:仿真周期cycles設置為11，實驗中一個仿真周期代表大熊貓生態系統中的一個自然年，周期運行步數 steps為5，為保證實驗結果的可靠性仿真次數simulations設置為1000次。

大熊貓生態系統通常涉及參數由于自然變化以及許多因素影響大熊貓個體的健康、死亡率和生物學現象，故所研究的過程本質上是隨機的。在種群定投膜系統中，進化規則以概率方式執行，可以有效處理含隨機性、不確定性的數據。部分進化規則的關鍵參數設置如表1～3所示。

表2 大熊貓出生率及其繁殖期年齡Tab．2 Birth ratio and reproductive age of giant panda

表3 野放大熊貓數量及年齡Tab．3 Quantity and age of released giant panda

不同年份雄性大熊貓和雌性大熊貓數量變化的實驗仿真結果與實際數據對比分別如圖1所示。

圖1 雌雄大熊貓數量實驗結果與實際數據對比Fig．1 Comparison of experimental and real results with respect to the numbers of female and male giant pandas

不同年份大熊貓總體數量變化實驗仿真結果與實際數據進行對比，結果如圖2所示。

圖2 大熊貓總體數量變化Fig．2 Changes in the overall number of giant pandas

從圖1～2的仿真結果可以看出，本文設計的種群動態膜系統模型仿真結果與真實數據相比較，雄性大熊貓仿真結果與實際數據的相對誤差最大不超過±8．95%，大部分控制在±2．5%以內;雌性大熊貓仿真結果與實際數據的相對誤差最大不超過±7．81%，大部分控制在±2．4%以內;全體大熊貓仿真結果與實際數據的相對誤差最大不超過±4．13%，大部分控制在±2．7%以內。同時由于種群基數較小導致相對誤差較為明顯，但依然遠小于±10%的相對誤差標準。同時與之前采用概率膜系統設計的簡單模型相比，該模型實驗結果誤差更小，擬合度更高，故本文設計的大熊貓種群動態膜系統較好地模擬了中國大熊貓保護研究中心大熊貓生態系統種群動態變化趨勢。

4 結語

本文提出了一種研究大熊貓種群動態的新方法，種群動態膜系統模型能夠有效預測種群未來的變化趨勢，可以及時評估種群發展的合理性，設計出一個基于大熊貓譜系數據的具有雙層膜結構、對象集和一系列進化規則的種群動態膜系統模擬實際生態系統中大熊貓固有的隨機進化過程。這是首次搜集和分析如此完整的譜系數據，并以此為基礎研究大熊貓的種群動態。同時有專門開發的仿真軟件MeCoSim來協助設計該計算模型并仿真驗證所提出的計算模型。它還可以進行基于仿真實驗結果的決策，診斷種群異常的原因，分析外界因素對種群動態性的影響，為大熊貓基地的工作人員提供決策參考，避免盲目決策。