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

樓一均，趙曉強

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

以及

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

觀察到第二個方程可以從整個系統(tǒng)解耦出來

在該模型中，t時刻的發(fā)育率為

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

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

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

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

據(jù)此，t時刻的發(fā)育率可以表述為

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

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

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

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

2.1 模型解的基本性態(tài)

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

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

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

2.2 基本再生數(shù)

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

對任意φ∈X，定義

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

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

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

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

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

(C2)Ω(Z)＜0.

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

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

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

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

2.3 不同的相空間

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

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

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

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

注 2.1根據(jù)解在不同相空間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+連續(xù).

進一步，周期半流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)為線性周期系統(tǒng)(12)以Y為相空間的t時刻解映射，即根據(jù)系統(tǒng)(12)過初值z0=φ∈Y的解z(t，φ)定義G(t)φ=zt(φ).則下面的結論顯示系統(tǒng)(12)在不同相空間X和Y具有相同的線性穩(wěn)定性.

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

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

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

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

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

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

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

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

3 不連續(xù)周期時滯

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

滯育是指動物受環(huán)境條件的誘導所產生的一種靜止狀態(tài)，該生物學策略使動物有能力在不利的季節(jié)條件下存活，并與有利的季節(jié)條件同步以有效利用資源.文獻[31]以蚊子的滯育現(xiàn)象為研究對象，特別是考慮了兩種不同蚊子物種的滯育策略:庫蚊(Culex)以成蟲形式滯育而伊蚊(Aedes)以卵形式滯育.對每一種滯育策略推導出不同的兩個數(shù)學模型.特別地，這兩個模型可以統(tǒng)一為如下一個模型.假設物種在每個時間周期只包含一個壞的環(huán)境季節(jié)，也即意味著物種在該季節(jié)進行滯育.基本的建模思路是將一個周期時間段分為三個時間區(qū)間:滯育前時間段 (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代表正常季節(jié)物種所需的發(fā)育時滯及壞環(huán)境季節(jié)的時長(每一年中T2時間段時長).因此可以在每個時間段上分別建立不同的系統(tǒng)來刻畫種群動力學.

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

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

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

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

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

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

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

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

(H2)2τ+τd＜1.

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

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

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

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

定義 Φt為系統(tǒng) (13)在集合Y+上的解映射，即如果u(t，φ)是系統(tǒng) (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+是連續(xù)的.

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

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

假設出生函數(shù)還滿足如下條件:

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

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

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

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

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

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

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

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

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

這兩種建模思路以及相關模型分析技巧在最近的研究中得到極大的發(fā)展.這些想法可以擴展至其它具體物種或者有多個階段的種群動力學研究，如考慮由季節(jié)性變化引起蜱蟲各階段發(fā)育時滯的周期性變化及其對蜱蟲種群動力學的影響[40]，周期發(fā)育時滯的海虱(Salmon lice)生長的階段結構模型[41]，多個周期時滯的階段結構蚊子種群模型[42]等.此外，在疾病傳播研究中，很多疾病的潛伏期(latency period)可能受外界環(huán)境影響.譬如溫度變化可以引起蚊子體內瘧原蟲的發(fā)育率變化，因此導致的周期性蚊子外潛伏期 (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]進一步將周期時滯函數(shù)推廣至概周期函數(shù)，用于研究概周期生態(tài)壞境下具有階段結構的種群動力學.此外，也可以在種群生物系統(tǒng)中考慮生物個體的空間移動，從而研究具有周期時滯和空間移動的種群生物系統(tǒng).

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

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

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

為了研究空間異質性，媒介(如蚊子)傳播疾病中溫度敏感的內潛伏期(intrinsic incubation period)和外潛伏期(extrinsic incubation period)，以及其它季節(jié)性變化等多種因素對疾病傳播的影響，文章[55]提出了一個具有周期時滯的非局部反應擴散系統(tǒng).該模型主要描述宿主(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).假設媒介 (如成年蚊子等)的生命周期較短，故不考慮感染媒介康復的情況.則宿主和媒介總密度分別為

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

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

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

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

模型 (15)右端Mh(t，x)和Mv(t，x)分別代表新增感染宿主和媒介的密度變化率，描述所有剛經歷潛伏期從潛伏類轉化為感染類的轉化率.該轉化率類似于模型(8)中的發(fā)育率(10)，用于描述相應個體在某一個階段經歷一定時間段至下一個階段的轉化率.由于考慮溫度等季節(jié)性變化對潛伏期的影響，假設潛伏期是關于時間t的周期函數(shù).通過考慮宿主在潛伏期[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)式中的討論，考慮個體的擴散系數(shù)Dh，假設J(q，t，x)為q方向在t時刻x位置的通量，則有如下方程

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

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

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

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

方程(16)可以重寫為

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

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

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

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

進一步，得到

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

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

對該式求關于t的導數(shù)得到

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結論

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

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

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