周期時滯的種群生物學模型研究綜述

樓一均，趙曉強

(1.香港理工大學應用數學系，香港 999077;2.紐芬蘭紀念大學數學與統計系，加拿大 A1C 5S7)

1 引言

在理論生態學研究中，兩類常用的數學模型可以較好地刻畫具有階段結構的種群密度的動力學演化:基于出生，發育，死亡的積分形式(見方程(1)和方程(2))，以及基于狀態(如年齡，尺度等)結構的偏微分方程形式(見方程(3)和方程(4)).本節將對這兩類模型以及它們在不同生物背景假設下的約化作簡要介紹，更多詳細內容可以參閱文獻[1-4].為方便闡述，考慮兩階段的最簡單形式，如在考慮物種生長時，可以將第一階段看作不具有生殖能力的幼年階段，第二階段為具有生殖能力的成年階段.其中物種在第一個階段和第二個階段位于t時刻的種群數量或者密度分別計為N1(t)和N2(t).接下來介紹基于兩種不同思路導出關于N1(t)和N2(t)的方程組.

第一種建模思路是基于對種群中個體的出生，發育及死亡等生態過程進行描述的積分形式.該建模思路可以追溯到Lotka在1907年的工作[5].為了刻畫個體在第一階段的發育率及考慮個體從幼年階段發育到成年階段所需時間的差異性[6]，引入逗留時間概率(sojourn function)函數P(a).該函數描述了在t時刻，具有第一階段停留時間a的個體(也就意味著該個體在t-a時刻進入第一階段)留在第一階段的概率.因此該函數P(a)是單調遞減的且P(0)=1.假設函數P(a)可微，那么關于以逗留時間為隨機變量的概率分布所對應的概率密度函數可以表述為-P′(a)(當P(a)不可微時，相關討論需用黎曼-斯蒂爾杰斯積分(Riemann-Stieltjes integral))[3].根據函數P(a)，可以計算其它相關的生物學指標，如平均逗留時間以及年齡為a的個體剩余逗留時間的期望(expected remaining sojourn times)為[3-4].考慮在t時刻第一階段的種群，其數量包含兩部分:由0時刻引入并在t時刻仍然逗留在該階段的個體以及[0，t]時刻新出生并且在t時刻仍然停留在第一階段的個體.基于此，N1(t)可以表述為

在該表達式中，B(N2(s))代表在s∈[0，t]時刻進入第一階段的輸入率(在此，假設由第二階段的出生率函數B(·)決定).此外，假設在第一階段的存活概率函數為Π1(a).在上式中，由初始0時刻引入，逗留并存活在第一階段的種群數量為

類似地，在第二階段的種群數量可以表述為

方程(2)和方程 (1)相似，但對于第二階段，假設逗留概率為 100%，即個體一旦進入第二階段，將永遠停留在該階段直至死亡(即在第二階段的逗留時間的概率函數為1).在s時刻輸入至第二階段的發育率[4]為

在上述表達式用到了逗留時間的概率密度函數-P′(a)，可以解釋如下:個體在(a，a+h)年齡段離開第一階段的概率為P(a)-P(a+h).取極限h趨向于零時，可得幼年個體在年齡為a時離開第一階段的速率為-P′(a).

另一種常用的模型構建思路是基于個體生理特征，如體型大小(body size)和年齡 (age)，建立相應的數學模型.此時，假設函數u(a，t)和v(η，t)為t時刻第一階段和第二階段的階段年齡(stage-specific age)分別為a和η的種群密度.令第一和第二階段的逗留時間概率函數為P(a)和100%，則t時刻種群數量N1(t)和N2(t)可以分別表述為

根據McKendrick考慮年齡結構種群密度的建模思想[7](也可參閱文獻[8-12]，得如下偏微分方程

初始時刻的階段密度分布由初值u0(a)和v0(η)給出.在該系統中，μ1和μ2代表第一階段和第二階段的個體死亡率(假設階段年齡為a的個體在每個階段的生存概率Πi(a)符合 Πi(a)=e-μia，則.邊值條件u(0，t) 和v(0，t) 分別代表在t時刻進入第一和第二階段的種群密度，則

當停留函數取特定的分布函數時，基于積分形式(方程(1)和方程(2))和偏微分形式(方程(3)和方程(4))的階段結構模型可以約化為數學上易于處理的常微分方程或者時滯微分方程系統.如假設P(a)服從Erlang分布(即Gamma分布中形狀參數為整數n的情形，兩類模型都可以約化為n+1個方程的階段結構模型[4].其中第一階段種群數量可以看作n個不同子階段的累積.特別地，當n=1以及假設個體在第一階段的平均停留時間為時，逗留函數服從指數分布，即P(a)=e-λa.此時有如下階段結構模型

當在第一階段的停留函數P(a)符合狄拉克(Dirac)分布，即

其中，τ為第一階段的平均停留時間.階段結構模型(方程(1)和方程(2)，或方程(3)和方程(4))可以約化為如下具有時滯τ的時滯微分方程系統

以及

其中u0由積分形式或者偏微分形式的初值決定.當研究種群數量N1(t)和N2(t)的長期動力學行為時，只需要在適當的給定初值意義下研究系統(6).此外，觀察到第二個方程可以從整個系統解耦出來，故可首先著重于了解變量N2(t)的動力學性質，即研究如下方程的動力學行為

根據N2(t)的動力學性質，N1(t)的動力學行為可以由積分形式(1)得到.

當停留函數P(a)符合狄拉克(Dirac)分布時，每個個體在第一階段具有相同的固定停留時間τ(逗留時間分布的方差為零).在很多生物系統中，每個存活個體在特定階段的停留時間τ可能受外界環境影響，而外界環境因素可以由關于時間t的函數給出.由此，發育所需時間用一個關于時間t的函數τ(t)描述更為合適.

在生物學研究中，計算發育所需時間τ(t)也是非常有意義的問題.其中一種方法可以利用發育比率得出[13-14].基于給定時間t所對應的環境條件可以得到所對應的發育比率.如在很多代謝生態學(metabolic theory of ecology)研究中，個體發育比率[15]可以由確定的表達式給出，譬如

其中b0是比率常數，E為活化能常數，k=8.6173×10-5eV K-1為玻爾茲曼常數(Boltzmann constant)，T(t)代表時刻t的熱力學溫度(絕對溫度).假設當個體累積發育比率達至100%時，個體發育至成年階段.假設t時刻成熟的個體在幼年階段的停留時間τ(t)，則這些個體在幼年階段的停留時間段為 [t-τ(t)，t].在此時間段上，個體累積發育比率達至100%.則有

假設任意時刻的發育比率r(ξ)嚴格大于零，則τ(t)是唯一確定的.同時，如果假設t時刻的溫度T(t)是關于t的ω-周期函數，則τ(t)也是周期函數.此外，對該積分形式作微分，得到

在本文中，將展示如何在模型建立時考慮周期性的發育時間τ(t)，以及如何對這些模型進行基本的動力學分析.

2 周期時滯模型的構建及分析

假設t時刻個體發育至第二階段所需的停留時間為τ(t)，則該個體在第一個階段的t-τ(t)時刻出生.由此，可以將(5)式中停留函數P(a)拓展至如下時間依賴形式

將根據年齡結構模型的框架考慮周期性的停留函數P(t，·).為此，引入多元函數u(a，t)表示t時刻具有生理年齡為a的種群密度.考慮周期性的外界環境對個體出生率，死亡率等的影響，則有如下偏微分方程

因此，在t時刻種群數量N1(t)和N2(t)可以分別表述為

進一步，對N1(t)和N2(t)求導可得如下微分方程

由于沒有個體可以永恒存活，故假定具有無窮年齡的種群密度為零，即u(∞，t)=0.假設t時刻的出生率u(0，t)為依賴于時間t和成年物種數量N2(t)的函數，即

此外，需要計算在t時刻的發育率 (1-τ′(t))u(t，τ(t)).此項可以由沿特征線積分法通過u(0，t-τ(t))=B(t-τ(t)，N2(t-τ(t)))表述.為此引入帶參數s的函數Vs(t)=u(t-s，t)，其關于時間變元的導數當t-s≤τ(t)時滿足如下方程

求解該線性常微分方程得到

當t≥:=max{τ(t)}時，設參數值s=t-τ(t)，則

因此，當t≥時，種群動力學可以由下面的時滯微分方程系統描述

觀察到第二個方程可以從整個系統解耦出來

在該模型中，t時刻的發育率為

此外，文獻[16]給出了另一種導出此發育率的思路.為此，作者引入發育水平變量q來刻畫個體在發育過程中的差異性.假設q=0以及q=1分別代表新出生和發育到成年階段時的發育狀態，則在t時刻的種群密度ρ(q，t)的變化率可以表述為[17]

在上述方程右端，J(q，t)描述在t時刻的發育流束(也稱通量).考慮周期性的發育速率，則J(q，t)=r(t)ρ(q，t)，其中r(t)類似于方程 (7)中的發育比率.因此，有

其中該方程的邊界條件由出生率給出ρ(0，t)=B(t，N2(t)).對(11)式沿特征線積分可得

其中τ(q，t)滿足，即在時間段 [t-τ(q，t)，t]上的累積發育水平為q.令τ(t)=τ(1，t)，則有

據此，t時刻的發育率可以表述為

另一方面，由于τ(t)滿足，則根據方程 (7)，發育率 (10)也可寫為

在該模型中，假設出生率函數B(t，N2)，單位死亡率μ1(t)和μ2(t)為正值函數.所有參數為關于時間t的ω-周期函數.此外假設出生率函數B(t，N2)滿足下列條件:

(A1)出生函數B(t，N2)是關于N2變量的利普希茨連續函數，且關于變量N2單調遞增.此外，對于任意t，B(t，0)≡0，以及B(t，N2)有界;

(A2) 函數B(t，N2) 可以表示為B(t，N2)=b(t，N2)N2，這里單位出生率b(t，N2)關于N2是嚴格單調遞減的.

2.1 模型解的基本性態

對于給定初值N2(θ)=ψ(θ)，θ∈[-，0]，方程(9)可以表述為

則關于方程(9)解的存在唯一性可以由時滯微分方程的一般理論得出[18-19].此外，關于 1-τ′(t)的表達式 (7)確保了 1-τ′(t)＞0，則周期時滯模型 (9)的發育率非負，那么根據文獻[47，定理5.2.1和注記5.2.1]可以推導出N2(t)的正性以及有界性(基于出生率的假設(A1)).為了得到系統(9)中N1(t)的正性，將此變量用積分形式表示為

鑒于該積分中被積函數的正性以及有界性，立即可得N1(t)的正性和有界性.具體的證明細節可見文獻[13].

2.2 基本再生數

基本再生數(basic reproduction number)可以計算在無密度制約因素時，單個成年個體在一生中能產生下一代的平均值.一般情況下，可以用基本再生數預測種群是否會滅絕.數學上可以由特定空間上的下一代算子(next generation operator)的主特征值所定義.近二十年來，關于基本再生數的數學理論得到了長足發展[21-28].對于含時滯的周期微分方程系統(9)，可以用文獻[29]的方法定義其基本再生數.為此，先對方程(9)在零解處進行線性化可得

對任意φ∈X，定義

則F(t)和V(t)為關于t的ω-周期函數.此外線性系統(12)可以記為

可以將F(t)看作是出生率相關的泛函，而V(t)可以用來刻畫成年個體的存活概率.定義Z(t，s)，t≥s為由線性方程定義的發展算子，即該算子滿足

尤其在當下社會，隨著社會生產力的提高，人民對物質文化生活的期望也日益提高，會有更多的社會力量去關注“精神食糧”，比如圖書館文獻資源建設的情況。這種具有公益性和慈善性的行為更是對社會的和諧發展，公民的人生觀、價值觀有著積極的指導意義。本文選取在中國歷史上比較具有代表性的階段——民國時期，對圖書館文獻捐贈的歷史進行梳理與探析。

此外，定義Z(t，s)的指數增長界

(C1)對任意t≥0，F(t)X+?R+;

(C2)Ω(Z)＜0.

定義所有從R到R的連續ω周期函數空間為Cω.對此空間賦予最大模范數，則該空間為巴拿赫空間.此外，定義正錐.假設v∈Cω為周期性的初始分布，則F(t-s)vt-s(t≥s≥0)為t-s時刻新發育進入第二階段的種群密度，Z(t，t-s)F(t-s)vt-s為t-s時刻進入第二階段，并存活至t(t≥s)的種群密度.由此

給出了所以初始分布為v(·)在一生中所累積的下一代數目.據此，下一代算子L:Cω→Cω可以定義為

根據文獻 [29]中的理論，可以將基本再生數定義為該算子的譜半徑 (spectral radius)，即R0:=r(L).

引理 2.1R0-1和r(W(ω))-1同號.

2.3 不同的相空間

對于模型(9)，當考慮X為相空間時，系統解的存在唯一性，有界性，以及解對初值的連續依賴性等結論可以由時滯微分方程系統的一般理論直接得出[18-19].但是當考慮相空間為X時，解映射是最終單調(monotone)而非強單調(strongly monotone)的.為了證明全局穩定性等結論，有時需要利用解映射的強單調性(如文獻[30，2.3節]).為此，引入另外一個相空間如下

則 (Y，Y+)為有序巴拿赫空間.對于此新空間中的u:[-τ(0)，+∞)→R和t≥0，記ut∈Y為

則在此新空間下，有如下結論，其具體證明類似于文獻[13]，故此處予以省略.

引理 2.2假設 (A1)成立，對于任意φ∈Y+以及t∈[0，∞)，方程 (9)存在以u0=φ為初值的唯一的非負解u(t，φ).

注 2.1根據解在不同相空間X和Y上解的唯一性可以得知:對于任意ψ∈X+和φ∈Y+，如果ψ1(θ)=φ1(θ)，?θ∈[-τ(0)，0]，則方程 (9)以初值u0=φ的解u(t，φ)和以初值v0=ψ的解符合u(t，φ)=v(t，ψ)，?t≥0.

引理 2.3定義t≥0時刻解映射Qt(φ)=ut(φ)，?φ∈Y+，則Qt是在空間Y+上的ω-周期半流.它符合下面三條性質:(i)Q0=I;(ii)Qt+ω=Qt°Qω，?t≥0;以及 (iii)Qt(φ)關于變量 (t，φ)∈[0，∞)×Y+連續.

進一步，周期半流Qt是最終強單調以及次齊性的，該結論由如下兩個引理給出:

引理 2.4假設 (A1)成立，對于Y+空間上兩個初值φ和ψ，如果滿足φ＞ψ(即φ≥ψ但是)，則方程(9)通過這兩個不同初值u0=φ和v0=ψ的解u(t)和v(t)滿足u(t)＞v(t)(當t＞時).也就意味著當t＞2時，Qt(φ)?Qt(ψ).

引理 2.5假設 (A1)及 (A2)成立，對于空間Y中的任意初值φ?0和任意γ∈(0，1)，當t＞時，解滿足u(t，γφ)＞γu(t，φ).因此對于任意n符合nω＞2，龐加萊映射 (Poincaré map)Qω滿足.

對于給定t≥0，定義G(t)為線性周期系統(12)以Y為相空間的t時刻解映射，即根據系統(12)過初值z0=φ∈Y的解z(t，φ)定義G(t)φ=zt(φ).則下面的結論顯示系統(12)在不同相空間X和Y具有相同的線性穩定性.

引理 2.6在相空間X中的龐加萊映射 (Poincaré map)W(ω)和相空間Y上的龐加萊映射具有相同的譜半徑，即r(W(ω))=r(G(ω)).

根據引理2.1，基本再生數可以決定系統以X為相空間的穩定性.結合上述引理可得基本再生數R0也可以決定線性系統以Y為相空間的穩定性.利用解映射在相空間Y上的最終強單調性(引理2.4)及次齊性(引理2.5)，結合單調動力系統相關理論(見文獻[30，定理2.3.4和引理2.2.1])，可得系統(9)的全局動力學結論:

定理 2.1(i)當R0≤1時，零平衡點在Y+相空間中是全局穩定的.

(ii)當R0＞1時，系統 (9)存在一個唯一的ω-正周期解，該正周期解對于所有初值在Y+{0}中是全局穩定的.

根據N1(t)關于變量N2(t)的表達式

可以進一步得知關于系統(8)中變量N1(t)的動力學行為.從而整個系統的動力學行為可以表述如下:

定理 2.2(i)當R0≤1時，(0，0)是全局穩定的.

(ii)當R0＞1時，系統(8)存在一個唯一的ω-正周期解該正周期解對于所有非平凡解是全局穩定的.

3 不連續周期時滯

具有周期時滯的模型(8)帶有正的修正項1-τ′(t).當時滯關于時間t的導數不存在，即周期時滯函數不可微時，該模型不成立.根據方程(7)，直覺上如果在一個ω周期內某個時刻的發育率為零，則1-τ′(t)不存在.此時基于模型(8)的建模思路不適用.事實上，當生物個體進行滯育(diapause)時，個體的發育比率為零，此時需要用其它方式推導出合理的數學模型.在本節中，主要展示文獻[31]中關于描述滯育現象的數學模型.

滯育是指動物受環境條件的誘導所產生的一種靜止狀態，該生物學策略使動物有能力在不利的季節條件下存活，并與有利的季節條件同步以有效利用資源.文獻[31]以蚊子的滯育現象為研究對象，特別是考慮了兩種不同蚊子物種的滯育策略:庫蚊(Culex)以成蟲形式滯育而伊蚊(Aedes)以卵形式滯育.對每一種滯育策略推導出不同的兩個數學模型.特別地，這兩個模型可以統一為如下一個模型.假設物種在每個時間周期只包含一個壞的環境季節，也即意味著物種在該季節進行滯育.基本的建模思路是將一個周期時間段分為三個時間區間:滯育前時間段 (pre-diapause period)，滯育時間段 (diapause period)，以及滯育后時間段 (post-diapause period).記每個時間段分別為T1:=[n，n+1-τ-τd]，T2:=[n+1-τd-τ，n+1-τ]和T3:=[n+1-τ，t≤n+1].為了闡述方便，假設周期為一年，其中n代表第n年，τ和τd代表正常季節物種所需的發育時滯及壞環境季節的時長(每一年中T2時間段時長).因此可以在每個時間段上分別建立不同的系統來刻畫種群動力學.

在時間段T1:=[n，n+1-τ-τd]上，假設環境因素是不變的，則種群的密度變化仍符合系統(6)的描述，即

在時間段T2:=[n+1-τd-τ，n+1-τ]上，生物個體進行滯育以在極端惡劣環境因素下生存，則幼年個體發育率急速降低，其它的生理活動也大大降低.假設在此時間段上，幼年個體不發育(即從幼年階段到成年階段的發育率為零)以及成年個體不進行生育活動.此外，幼年和成年會受到不同于正常時間段的死亡率影響.基于此，可以用下面系統刻畫種群在該時間段上的動力學:

特別地，當物種以幼年個體的滯育為生存策略時(如伊蚊)，假設在此時間段上成年個體的死亡率極高，即d2?0.而當成年個體滯育時(如庫蚊)，則設此時間段上幼年階段死亡率極高，即d1?0.

在時間段T3:=[n+1-τ，n+1]生物個體能正常生長，發育及生育.但由于T3時間段的長度為τ，則此時間段上所有發育為成年的個體都經歷了T2時間段，也即意味著所有發育個體都經歷了τ+τd發育時滯，同時所對應的生存概率為 e-μ1τ-d1τd.利用系統(6)的建模思路，該時間段上的種群密度變化可以描述為

類似于系統(8)的分析思路，變量N2(t)所對應的方程同樣可以從模型中解耦出來.因此，可先對N2(t)變量進行分析以討論整個物種的種群動力學行為，即先研究如下系統:

在該系統中，n∈N代表第n年.當考慮每個時間段的端點時，考慮的是單邊導數.由于一年內不同時間區間T1及T3上的時滯分別為τ和τ+τd，而在T2時間段上的時滯可以選取任意大于等于零的數值，所以可以將時滯看作關于時間t的周期函數，而該周期函數是不連續的.

對于出生率函數以及死亡率等模型參數，引入下面常用的生物假設[32]:

(H1)出生函數B(N2)是非負的利普希茨連續函數，并且是關于變量N2的單調遞增函數.此外，B(0)=0.存在使得B(N2)e-μ1τ＞μ2N2當，以及B(N2)e-μ1τ＜μ2N2當.

(H2)2τ+τd＜1.

上述條件(H1)為常用的生物假設.此外，蚊子滯育由極端寒冷或者干旱季節所引起，當環境條件變好時，滯育期結束[33].一般地，滯育期對于不同的蚊子物種以及不同地理位置可以維持 3-5個月，即τd在3-5個月取值.同時，幼年蚊子的發育周期較短，基本維持在2-4周[34]，即τ在2-4周取值.從這個意義上，當考慮以年為時間單位時，假設(H2)是合理的.

引入空間Y=C([-τ，0]，R)，對任意φ∈Y，定義范數，則 (Y，‖·‖)為巴拿赫空間.對于給定初值φ∈Y+:=C([-τ，0]，R+)，關于N2(t)變量的系統(13)解的存在唯一性可以通過在不同時間段T1，T2以及T3上分別進行討論.假設系統在t時刻的解為u(t，φ).根據時滯微分方程常用記號ut(θ，φ)=u(t+θ，φ)，?θ∈[-τ，0]，有u0=φ.此外，N2(t)變量的正性可以通過分析積分形式得到

綜合這些分析，并利用文獻[31]的證明細節，可以得到如下關于系統(13)解的存在唯一性以及正性的結論.

定理 3.1假設條件(H1)和條件(H2)存在，則對任意φ∈Y+及初值u0=φ，對所有時刻t∈[0，∞)系統(13)存在唯一的非負及有界的解u(t，φ).

定義 Φt為系統 (13)在集合Y+上的解映射，即如果u(t，φ)是系統 (13)過初值u0=φ∈Y+的解，則 Φt(φ)(θ)=ut(θ，φ)=u(t+θ，φ).

引理 3.1解映射Φt為如下意義下的1-周期半流:(i)Φ0=I，此處I為恒等映射;(ii)Φt+1= Φt°Φ1，?t≥0;(iii)Φt(φ)關于 (t，φ)∈[0，∞)×Y+是連續的.

此外，也可以證明該周期半流是最終強單調及嚴格次齊性的.相關證明可以參閱文獻[31].

引理 3.2對于Y+中給定的兩個初值φ和ψ，假設φ＞ψ(即對于所有s∈[-τ，0]，φ(s)≥ψ(s).此外)，則通過這兩個初值的系統 (13)的解u(t，φ)和v(t，ψ)當t＞τ+τd時滿足u(t，φ)＞v(t，ψ).這意味著當t＞2(τ+τd)時，Φt(φ)?Φt(ψ).

假設出生函數還滿足如下條件:

(H3)出生函數B(N2)可以表示為B(N2)=b(N2)N2以及單位出生率b(N2)關于N2是嚴格單調遞減的.

則在額外假設(H3)下，解映射所定義的周期半流是嚴格次齊性的.

引理 3.3對于Y+中的任意函數φ?0以及λ∈(0，1)，當t＞τ+τd時，有u(t，λφ)＞λu(t，φ).這意味著當整數n滿足n＞2(τ+τd)時，.

基于假設(H1)，系統(13)有一個種群滅亡平衡點0.在此平衡點線性化可得如下線性系統:假設P(t)為線性系統(14)在Y+上的解映射.則可以根據龐加萊映射(Poincaré map)P(1)的譜半徑，記作R，描述線性系統的穩定性.此外，根據單調動力系統理論(如文獻[30，2.3節])以及引理3.2-引理3.3，可以得到如下全局動力學結論.

定理 3.2對于系統(13)，假設條件(H1)-條件(H3)成立，則有如下結論:

(i)當R≤1時，0平衡點在Y+中是全局漸近穩定的.

(ii)當R＞1時，系統存在 1-周期解M*(t)，該周期解相對于Y+{0}是全局漸近穩定的.

通過將N1(t)表示為關于N2(t)的積分形式，則可以得出類似于定理2.2中關于變量N1(t)的全局動力學行為.

在本節最后，介紹一些具有周期時滯的相關研究進展.當系統中的非固定時滯為時間的函數τ(t)，模型往往帶有額外修正項1-τ′(t)[35-37].該項也常見于當時滯由模型變量定義 (state-dependent delay)的系統中[38].關于該項的正性，即 1-τ′(t)＞0，除了可以通過商形式(7)說明外，文獻[35]給出了更多實際問題中的討論.事實上，生物文獻更多地將該項寫成兩個特定時間發育比率的商的形式(7)[39].

這兩種建模思路以及相關模型分析技巧在最近的研究中得到極大的發展.這些想法可以擴展至其它具體物種或者有多個階段的種群動力學研究，如考慮由季節性變化引起蜱蟲各階段發育時滯的周期性變化及其對蜱蟲種群動力學的影響[40]，周期發育時滯的海虱(Salmon lice)生長的階段結構模型[41]，多個周期時滯的階段結構蚊子種群模型[42]等.此外，在疾病傳播研究中，很多疾病的潛伏期(latency period)可能受外界環境影響.譬如溫度變化可以引起蚊子體內瘧原蟲的發育率變化，因此導致的周期性蚊子外潛伏期 (extrinsic incubation period)時滯對瘧疾傳播具有巨大的潛在影響.鑒于此，可以在疾病傳播模型中考慮周期性潛伏期時滯，如周期潛伏期的SEIRS模型[43]，周期性的蚊子外潛伏期時滯的瘧疾傳播模型[16]，具有周期外潛伏時滯的血吸蟲病模型[44]，周期時滯的基孔肯雅熱 (chikungunya fever)模型[45-47]，周期潛伏期的藍舌病 (Bluetongue)傳播模型[48]，周期外潛伏期的西尼羅河病毒傳播模型[49]，周期外潛伏期及垂直傳播的西尼羅河病毒傳播模型[50]，受溫度影響的周期性潛伏期時滯的手足口病 (hand，foot and mouth disease)疾病傳播模型[51]，具有周期性潛伏期的黃龍病(Huanglongbing)模型[52]，以及帶有周期性潛伏期的蝙蝠傳播狂犬病(bat-borne rabies)模型[53]等.

文獻[54]進一步將周期時滯函數推廣至概周期函數，用于研究概周期生態壞境下具有階段結構的種群動力學.此外，也可以在種群生物系統中考慮生物個體的空間移動，從而研究具有周期時滯和空間移動的種群生物系統.

4 考慮個體空間移動的周期時滯模型

將拓展上一節中的具有可微或不連續周期時滯模型，同時考慮生物個體在空間上的移動.主要通過文獻[55-56]的一些結論來展現典型的建模思路以及相關數學結論.

4.1 具有空間移動及可微周期時滯的媒介疾病傳播模型

為了研究空間異質性，媒介(如蚊子)傳播疾病中溫度敏感的內潛伏期(intrinsic incubation period)和外潛伏期(extrinsic incubation period)，以及其它季節性變化等多種因素對疾病傳播的影響，文章[55]提出了一個具有周期時滯的非局部反應擴散系統.該模型主要描述宿主(host)和媒介(vector)的各類倉室種群密度在空間Ω上的動力學變化.將宿主分為易感類 (susceptible)，潛伏類 (exposed)，感染類 (infectious)以及康復類(recovered)，這些類別在時刻t空間位置x上的密度分別為Sh(t，x)，Eh(t，x)，Ih(t，x)以及Rh(t，x).此外，將媒介種群分為易感類，潛伏類，感染類，其密度的時空分布分別為Sv(t，x)，Ev(t，x)和Iv(t，x).假設媒介 (如成年蚊子等)的生命周期較短，故不考慮感染媒介康復的情況.則宿主和媒介總密度分別為

根據病原體在宿主和媒介之間的傳播途徑以及宿主和媒介的種群生長特征，同時考慮宿主和媒介的移動，可以寫出如下系統:

這里的下標h和v分別代表宿主(host)和媒介(vector)相關變量或參數.模型參數的具體意義可見表1.模型(15)中，感染項

表1 模型(15)中參數的生物學意義

中b代表單位時間內每個媒介在宿主上的叮咬次數，h代表每次叮咬從感染媒介到易感宿主的感染概率而v代表從感染宿主到易感媒介的感染概率.空間異質性對物種生長和疾病的傳播影響由依賴于空間變元x的參數體現，而個體的空間移動由宿主擴散系數Dh以及媒介擴散系數Dv描述.此外假設個體不能穿越邊界，故對此偏微分系統，給定諾伊曼邊界條件(Neumann boundary condition).

模型 (15)右端Mh(t，x)和Mv(t，x)分別代表新增感染宿主和媒介的密度變化率，描述所有剛經歷潛伏期從潛伏類轉化為感染類的轉化率.該轉化率類似于模型(8)中的發育率(10)，用于描述相應個體在某一個階段經歷一定時間段至下一個階段的轉化率.由于考慮溫度等季節性變化對潛伏期的影響，假設潛伏期是關于時間t的周期函數.通過考慮宿主在潛伏期[t-τh(t)，t]的空間移動可以導出Mh(t，x)的表達式.為此假設q為刻畫宿主個體感染水平的某種度量，并且q在t時刻的變化率為γh(t)(依賴于t時刻的溫度T(t))，即.則可以用q=qE代表宿主從Sh類剛轉化為Eh類時的感染水平，而q=qI代表將從Eh類轉化為Ih類時個體的感染水平.假設ρ(q，t，x)代表在t時刻x位置具有感染水平q的種群密度，則新增具有感染性的宿主的變化率為Mh(t，x)=γh(t)ρ(qI，t，x).此外，從易感宿主轉化為潛伏類宿主的感染項給出具有qE感染水平的種群密度:βh(t，x)Sh(t，x)Iv(t，x)=γh(t)ρ(qE，t，x).下面將主要展示如何通過計算ρ(qI，t，x)得到Mh(t，x)的表達式.為此，需要刻畫ρ(q，t，x)的變化情況.類似于(11)式中的討論，考慮個體的擴散系數Dh，假設J(q，t，x)為q方向在t時刻x位置的通量，則有如下方程

由于J(q，t，x)=γh(t)ρ(q，t，x)，則有

方程(16)的邊值條件由下式給出

為解上面的方程(16)，引入關于時間t的新變量

由于該函數嚴格遞增，其反函數t=h-1(η)存在.據此定義

方程(16)可以重寫為

設V(s，x)=(s+q-η，s，x)，則有

由于q∈[qE，qI]，η-(q-qE)≤η，從而

在上述表達式中，Γh(t，t0，x，y) 是由方程在諾伊曼邊界條件下定義的格林函數(Green function).因此

假設τ(q，t)代表在t時刻從初始感染水平qE到達感染水平q所需的時間(即個體的感染水平在時間段[t-τ(q，t)，t]上從初始感染水平qE增加到q).考慮到則有

進一步，得到

由于當感染水平為qI時，個體從潛伏類發育到感染類，故在t時刻進入潛伏類的個體所需的潛伏期為τh(t)=τ(qI，t).從而

此外，在方程(17)中設q=qI，則有

對該式求關于t的導數得到

類似地，Mv(t，x)的表達式也可以通過分析媒介在潛伏期[t-τv(t)，t]的空間移動得出:

其中 Γv(t，t0，x，y) 為由方程以及諾伊曼邊界條件所定義的格林函數.將Mh(t，x)以及Mv(t，x)代入方程組 (15)，則得到具有非局部周期時滯的反應擴散系統.在此系統中，周期時滯描述了由季節性溫度條件決定的外潛伏期τv以及內潛伏期τh.

文獻[55]針對模型(15)做了詳細的分析，內容主要包括模型解在給定的生物假設下的全局存在性以及正性，疾病傳播基本再生數的定義，以及由基本再生數刻畫的疾病傳播的閾值結論:當基本再生數小于1時，疾病不能傳播;而當該數大于1時，疾病在空間上持續傳播.此外，當所有參數為常數時，進一步得到了當基本再生數大于1時正平衡點的全局吸引性.該文還利用參數對溫度的依賴關系，擬合出周期參數并進行有意義的數值模擬.更多相關結論見文獻[55].

4.2 具有空間移動及不連續周期時滯的種群模型

對于不連續周期情形(如考慮物種的滯育現象)，也可以將模型(13)進行拓展以考慮個體在空間上的移動.在此，主要總結文獻[56]的模型以展示關于幼年成年兩個階段的種群密度N1(t，x)和N2(t，x)的動力學模型.對于第n年正常生長期，即在時間區間t∈T1:=(n，n+1-τ-τd]時，n=0，1，2，···，種群動力學由下面的模型刻畫

在該系統中，D1和D2代表幼年個體和成年個體的擴散系數，出生率為關于成年種群密度的函數B(N2(t，x))，μ1和μ2代表幼年階段和成年階段的死亡率.其中，發育率M1(t，x;N2)需要通過詳細考慮幼年個體在成熟期的空間移動以及死亡概率等因素.為此，引入t時刻在x∈Ω位置具有年齡a的種群密度p(t，x，a).假設τ為在此時間段上成熟的個體所需要的發育時間，則M1(t，x;N2)=p(t，x，τ).為了得到完整的模型，需要將p(t，x，τ)表述為關于變量I和M的泛函.此外根據對于年齡結構種群密度的McKendrick von-Foerster方程，在T1時間段上該變量可以由下面微分方程描述

由于p(t，x，τ)表示t時刻新發育的個體的密度，這些個體在t-τ時刻出生，故此可以在時間區間 [t-τ，t]上討論變量p(t-τ，y，0)到p(t，x，τ)的演化:時間從t-τ時刻演化到t時刻，個體年齡從a=0發育為a=τ，此外由于個體可能的空間移動，在時間區間 [t-τ，t]上，個體位置可能由y變為x.引入新變量q(s，x):=p(t-τ+s，x，s).根據系統(21)的第一個方程，q(s，x)在年齡區間0＜s≤τ上的演化規則可以由下面方程給出

此演化方程描繪了幼年個體密度的動力學，考慮了個體的死亡和擴散過程.求解初值問題(22)(也可以由系統(21)的第一個方程用沿特征線積分法得出，如文獻[57])可得

其中格林函數Γ(t，x，y)為由偏微分算子[?t-Δ]在給定的邊界條件下所定義的基本解.所以可得發育率

在滯育期T2:=(n+1-τ-τd，n+1-τ]時間段上，由于個體極端的低代謝活動，假設物種沒有出生，發育以及空間移動.在此時間區間上，

此時，幼年和成年階段的自然死亡率d1和d2一般大于T1和T3時間段上的對應的死亡率μ1和μ2.

當進入t∈T3:=(n+1-τ，n+1]時間段時，類似于方程組(20)，有如下方程

在該方程中，主要需要確立發育率M2(t，x;N2).由于T3時間區間長度為個體在正常環境條件下所需發育時間τ，所有在該時間段上新出生的個體不能在該時間區間完成發育.即所有新發育的個體必須在T1時間段出生，并經歷滯育期T2.所以成熟個體經歷的發育時間為τ+τd.此時，M2(t，x;N2)=p(t，x，τ+τd).類似于T1時間段上發育率的推導，可得個體密度從p(t-τ-τd，y，0)演化到p(t，x，τ+τd)的演化規律:時間從t-τ-τd時刻演化到t時刻，個體年齡從a=0發育為a=τ+τd，此外由于個體可能的空間移動，在時間區間[t-τ-τd，t]上，個體位置可能由y變為x.數學上可以描述如下

其中αt=n+1-t是依賴于t的時間刻度.此演化過程考慮了個體的移動和死亡.在壞環境季節(αt，αt+τd]∈T2時間段上，幼年個體不進行空間移動，只有自然死亡率d1.而在 (0，αt]∈T1以及 (αt+τd，τ+τd]∈T3上，幼年個體具有空間移動及自然死亡率μ1.求解該柯西初值問題可得

因此，該時間段上的發育率可以表述為

該發育率考慮了兩個不同時間段的生存概率 e-μ1τ和 e-d1τd，其中在T1和T3時間段上的經歷的時間為τ，在T2時間段上經歷時長為τd.

綜上所述，得出一個具有空間移動及不連續周期時滯的更替模型(20)，模型(23)以及模型(24).對于該模型，文獻[56]考慮了有界區域及無界區域的情形.對于有界區域情形，定義了基本再生數，以及證明了系統的全局穩定性動力學.對于無界區域的情形，討論了系統周期行波解的存在性以及漸近傳播速度等問題.更多的建模過程以及分析細節可以參閱文獻[56].

在本節最后，介紹具有空間移動的周期時滯系統相關研究進展.這些考慮個體空間移動的周期時滯模型研究極大地豐富了種群動力學相關結論，如具有階段結構的蚊子種群生長反應擴散方程模型[58]，具有周期發育時滯的反應擴散階段結構種群模型[59]，帶有周期時滯以及考慮種內競爭的階段結構反應擴散系統[60]，廣義的具有周期時滯的時空媒介傳播傳染病模型[61]，具有周期時滯的反應擴散媒介傳播疾病模型[55]，周期潛伏期的反應擴散多群組傳染病模型[62]，具有非局部擴散及周期性潛伏期的藍舌病傳播模型[63]，周期性潛伏期及蚊子叮咬偏向性的瘧疾傳播反應擴散模型[64].此外，文獻[65]考慮了在時間空間周期的環境中的種群擴散問題，文獻[66]研究了季節更替對具有階段結構種群擴散的影響，文獻[67]討論了帶有周期時滯和非局部擴散周期系統所定義的龐加萊映射的擴散性態.

5 結論

本文主要展示了如何在一個階段結構模型中引入周期發育時滯τ(t)，以及對于此類模型的分析思路.假設幼年個體在任何時刻具有正發育率，則該周期時滯可以通過累積發育率定義，此時，周期時滯為關于時間t的可微函數.但是當幼年個體由于受惡劣生存環境影響導致發育停滯(如物種進行滯育時)，則幼年階段發育時滯為關于時間t的不連續函數.此類種群動力學可以由具有時滯的更替模型(succession model)描述.此外，這些更替模型可以從周期系統的觀點進行研究.

本文中介紹的周期時滯主要用來描述從一個階段到另外一個階段的轉化率(transit rate)，如從幼年到成年階段的發育率，從潛伏階段到感染階段的轉化率等.相較于更為具體的年齡結構或者尺度結構(size-structure)模型，階段結構模型較好地平衡了生物復雜性與數學可分析性等方面的因素.期望此類研究能得到進一步的發展.

致謝 感謝白振國，李志民，羅顏濤和史洋洋對本文初稿提出寶貴意見!