摘" 要： 主要研究一類病原體與中性粒細胞（PMNs）相互作用的免疫應答模型，詳細分析了其正平衡點的存在條件和類型，并通過MATLAB軟件包對所得到的結論進行了數值模擬。

關鍵詞： 免疫應答模型；正平衡點；尖點

中圖分類號： O175

文獻標識碼： A" 文章編號： 2096-3998（2024）02-0079-07

收稿日期：2023-04-25" 修回日期：2023-06-27

基金項目：國家自然科學基金項目（61763024）

作者簡介：孫得寧（1998—），男，甘肅張掖人，碩士研究生，主要研究方向為微分方程與動力系統；張存華（1972—），女，甘肅武山人，博士，教授，主要研究方向為微分方程與動力系統。

引用格式：孫得寧，張存華.一類免疫應答模型的分支分析.陜西理工大學學報（自然科學版），2024，40（2）：79-85.

人體內主要通過先天免疫系統對病原體作出反應。病原體入侵肺部后，首先巨噬細胞殺死病原體，當巨噬細胞檢測到自身無法遏制的威脅時，中性粒細胞（PMNs）會對病原體作出反應。中性粒細胞與淋巴細胞的比率（NLR）已被視為預測肺炎感染的可靠因子，Wang Dawei等發(fā)現新冠患者在危重期PMNs計數增加且淋巴細胞計數減少（即NLR增加），這表明患者病情可能處于危重狀態(tài)。

Pugliese等提出了一個關于病原體與免疫細胞相互作用的系統，但相互作用項非常具體；D’onofrio對該系統進行了優(yōu)化。從文獻來看，PMNs的活性在正常水平下較低，隨著病原體入侵肺部后迅速增加，因此，Young等構建了如下系統：

=αx-mx1+bx-βxz，

=δx-γxz+η-μz，（1）

其中，x（t）表示t時刻病原體數量的密度，z（t）表示t時刻PMNs肺部中載量的密度，α表示病原體的復制速率，mx1+bx表示通過先天免疫系統清除病原體的數量，-βxz表示PMNs吞噬病原體的數量，-γxz表示單個獨立的免疫細胞抑制病原體的數目有限，η表示正常水平下PMNs的流入量，-μz表示PMNs的自然死亡率，參數α、m、b、β、δ、γ、η、 μ均為正。

Young等發(fā)現當少量病原體入侵人體后，能夠迅速被PMNs清除，即

α＜m，（2）

而他們認為巨噬細胞清除病原體的能力有限，則分析了系統（1）在條件αβ＞δγ下平衡點的類型，且得到邊界平衡點是穩(wěn)定結點，正平衡點是鞍點。然而，參數的實際值因人而異，較大病原體可能不會克服先天免疫系統，所以又考慮了相反的情況

αβ＜δγ，（3）

并分析了系統（1）在條件（3）下平衡點的個數，但沒有提供具體的理論分析。Shi Shujing等分析了系統（1）在條件（3）下平衡點的類型及其穩(wěn)定性，并得到隨著參數值的變化，系統（1）展現了一系列分支現象：包括余維一鞍結點分支和Hopf分支。然而，隨著年齡的增長人體的免疫系統開始衰老和減弱，PMNs的自然死亡率隨之升高，故提出如下系統：

=αx-mx1+bx-βxz，

=δx-γxz+η-μz2。（4）

接下來主要討論系統（4）的正平衡點及其類型。

1" 正平衡點及其類型

首先做如下尺度變換：

x=α2βδx1，" z=αβz1，" t=1αt1，

用x、z、t表示x1、z1、t1，則系統（4）化為

=x-p1x1+p2x-xz，

=x-p3xz+p4-p5z2，（5）

其中，

p1=mα，

p2=α2bβδ，

p3=αγβδ，

p4=βηα2，

p5=μβ。

由式（2）和式（3）知p1gt;1，p3lt;1，而p1、p2、p3、p4、p5均為正參數，因此，僅在下列條件下研究系統（5）：

0lt;p3lt;1

p1gt;1

p2，p4，p5gt;0。（6）

系統的正不變區(qū)域為

Ω={（x，z）|x≥0，z≥0}，

為了討論系統（5）的正平衡點，令

x-p1x1+p2x-xz=0，

x-p3xz+p4-p5z2=0，（7）

由此得

z=1-p11+p2x，

（p22-p22p3）x3+（2p2-2p2p3+p1p2p3+p22p4-p22p5）x2+

（1-p3+p1p3+2p2p4-2p2p5+2p1p2p5）x+p4-p5+2p1p5-p21p5=0，

令

f（x）=（p22-p22p3）x3+（2p2-2p2p3+p1p2p3+p22p4-p22p5）x2+

（1-p3+p1p3+2p2p4-2p2p5+2p1p2p5）x+p4-p5+2p1p5-p21p5，

則

f′（x）=3（p22-p22p3）x2+（4p2-4p2p3+2p1p2p3+2p22p4-2p22p5）x+

1-p3+p1p3+2p2p4-2p2p5+2p1p2p5，

令

a=3p22-3p22p3，

b=4p2-4p2p3+2p1p2p3+2p22p4-2p22p5，

c=1-p3+p1p3+2p2p4-2p2p5+2p1p2p5，

Δ=b2-4ac，

則當Δgt;0且blt;0時， f′（x）有兩個正根，即f（x）存在極值點μ1、 μ2；當Δlt;0時， f′（x）不存在正根，即f（x）是單調遞增的。

由式（6）得

p22-p22p3gt;0。

那么得到f（x）根的5種情況（見圖1）。

（a） 三個不同的正根：E1，E2，E3""" （b） 兩個正根：重根E*，單根E3""" （c） 兩個正根：單根E1，重根E*

（d） 一個正根：重根E*""""""""""""" （e） 不存在正根

由圖1可知，首先，當f（0）lt;0，blt;0，Δgt;0， f（μ1）gt;0且f（μ2）lt;0（μ2gt;μ1gt;0）時， f（x）有三個正根（見圖1（a））。

其次，當f（0）lt;0，blt;0，Δgt;0， f（μ1）=0或f（μ2）=0時， f（x）有一個正的單根和一個正的重根（見圖1（b）（c））。

再次，當f（0）gt;0，blt;0，Δgt;0， f（μ1）=0時， f（x）有一個正的重根（見圖1（d））。

最后，當f（0）gt;0，blt;0，Δgt;0， f（μ2）gt;0或f（0）gt;0，Δlt;0時， f（x）不存在正根（見圖1（e））。

接下來研究系統（5）正平衡點的類型，系統（5）在E00，p4p5處的Jacobian矩陣為

J（E0）=

1-p1-p4p50

1-p3p4p5-2p4p5

〗，

該矩陣有兩個特征值λ1=1-p1-p4p5和λ2=-2p4p5，再由（6）知p1gt;1，則λ1lt;0，因此E0是穩(wěn)定雙曲結點。

系統（5）在正平衡點E（x，z）處的Jacobian矩陣為

J（E）=

p1p2x（1+p2x）2-x

1-p3+p1p31+p2x-p3x-2p5+2p1p51+p2x

〗，

由Jacobian矩陣得到J（E）的行列式為

Det（J（E））=

2p21p2p5x-（p1p2p3x2+2p1p2p5x）（1+p2x）

（1+p2x）3

，

J（E）的跡為

Tr（J（E））=

p1p2x（1+p2x）2+

2p1p51+p2x-

p3x-2p5。

根據Det（J（E））與f′（x）之間的關系得到

Det（J（E））=

xf′（x）（1+p2x）2。（8）

令

p*1=（1+p2x）p2x

則當p1lt;p*1（p1gt;p*1）時，Tr（J（E））gt;0（Tr（J（E））lt;0），由以上結論得到以下定理。

定理1" 當（6）中的條件滿足時，系統（5）總是存在一個無病平衡點E00，p4p5，它是穩(wěn)定雙曲結點。此外，

（Ⅰ）若f（0）gt;0、blt;0、Δgt;0、 f（μ2）gt;0，或f（0）gt;0、Δlt;0，則系統（5）不存在正平衡點；

（Ⅱ）若f（0）gt;0，blt;0，Δgt;0， f（μ2）=0，則系統（5）存在唯一正平衡點E*（x*，z*），它是一個退化平衡點；

（Ⅲ）若f（0）lt;0，blt;0，Δgt;0， f（μ1）=0或f（μ2）=0，則系統（5）存在兩個不等的正平衡點：一個退化平衡點E*（x*，z*）和一個雙曲平衡點Ei（i=1，3）。當p1lt;p*1時，Ei（i=1，3）為穩(wěn)定的雙曲結點或焦點；當p1gt;p*1時，Ei（i=1，3）為不穩(wěn)定的雙曲結點或焦點；當p1=p*1時，Ei（i=1，3）為弱焦點或中心；

（Ⅳ）若f（0）lt;0，blt;0，Δgt;0， f（μ1）gt;0且f（μ2）lt;0，則系統（5）存在三個不等的正平衡點：E2為雙曲鞍點，Ei（i=1，3）為雙曲結點或焦點；當p1lt;p*1時，Ei為穩(wěn)定的雙曲結點或焦點；當p1gt;p*1時，Ei為不穩(wěn)定的雙曲結點或焦點；當p1=p*1時，Ei為弱焦點或中心。

證明" 根據等式（8）與圖1可知Det（J（Ei））gt;0（i=1，3），Det（J（E2））lt;0，Det（J（E*））=0，則E1、E3為雙曲結點或焦點，E2為雙曲鞍點，E*為退化平衡點。

定理2" 當f（0）gt;0，blt;0，Δgt;0， f（μ2）=0與（6）中的條件滿足時，系統（5）存在唯一正平衡點E*（x*，z*）。此外，

（Ⅰ）當p1lt;p*1（或p1gt;p*1）時，E*（x*，z*），是包含一個穩(wěn)定拋物扇形（或不穩(wěn)定拋物扇形）的鞍結點；

（Ⅱ）當p1=p*1時，E*（x*，z*）是一個尖點，此外，

（?。┤鬎≠0，則E*（x*，z*）是一個余維二的尖點；

（ⅱ）若F=0，則E*（x*，z*）是一個至少余維三的尖點。

證明" 令Tr（J（E*））=0得到

p1=p*1=（1+p2x）p2x，

那么當p1lt;p*1（p1gt;p*1）時，Tr（J（E*））lt;0（Tr（J（E*））gt;0）。

下面來討論E*（x*，z*）的具體類型，令X=x-x*，Z=z-z*，將平衡點移到原點并將該式進行泰勒展開，則系統（5）化為

dXdt=a11X+a12Z+a13X2-XZ，

dZdt=a21X+a22Z-p3XZ，

（9）

其中，

a11=p1p2x*（1+p2x*）2，

a12=-x*，

a13=p1p2（1+p2x*）4，

a21=p1p31+p2x*-p3+1，

a22=2p1p51+p2x*-p3x*-2p5。

令

X=-x*μ+p1p2x*（1+p2x*）2v，

Y=-p1p2x*（1+p2x*）2μ+1-p3+p1p31+p2x*v，

τ=p1p2x*（1+p2x*）2+2p1p51+p2x*-p3x*-2p5t。

用t表示τ，則系統（9）化為

=b11μ2+b12μv+b13v2，

v·=v+b21μ2+b22μv+b23v2，（10）

其中，

b11=a12a13a21a22+a11a21a22-p3a211a22（a11a21+a21a22）（a11+a22），

b12=（2a11a12a13-a12a21+a211+p3a11a12-p3a211a12）（a11+a22）a12a22（a11+a22），

b13=a211a13a22-a11a21a22+p3a211a22（a12a22+a11a12）（a11+a22），

b21=p3a11a12a11a21+a21a22+a11a12a13a21a22+a211a21a22-p3a311a22（a11a221+a221a22）（a11+a22），

b22=p3a211-p3a12a21a21（a11+a22）+2a211a12a13-a11a12a21+a311+p3a211a12-p3a311a12a12a21a22，

b23=a311a13a22-a211a21a22+p3a311a22（a12a21a22+a11a12a21）（a11+a22）-p3a11a11+a22，

經計算判斷b11≠0，依據文獻中11.3節(jié)的中心流形理論，得到限制在中心流形上的方程為

dμdt=b11μ2+o（|μ|3），（11）

依據文獻中的定理7.1，結論（Ⅰ）成立。

令X=x-x*，Z=z-z*，p1=p*1，Det（J（E*））=0，則系統（5）化為

dXdt=c11X+c12Z+c13X2-XZ，

dZdt=c21X+c22Z-p3XZ，（12）

其中，

c11=p1p2x*（1+p2x*）2，

c12=-x*，

c13=p1p2（1+p2x*）4，

c21=p1p31+p2x*-p3+1，

c22=2p1p51+p2x*-p3x*-2p5。

令

X=-x*μ-v，

Z=-p1p2x*（1+p2x*）2μ，

τ=（1+p2x*）2p1p2x*（1+p2x*）t。

用t表示τ，則系統（12）化為

=v+d11μ2+d12μv，

v·=d21μ2+d22μv+d23v2，（13）

其中，

d11=p3c11c12c21，

d12=p3c11c21，

d21=p3c11c212-c11c212c13-c211c12c21，

d22=c211-p3c11+2c11c12c13c21，

d23=-c11c13c21。

根據文獻中的引理2.4，用x、z表示μ、v，得到系統（13）在原點的小鄰域內等價于系統（14）：

=z，

=Dx2+Fxz+o（|x，z|3），（14）

其中，D=d21，F=d22+2d11。經計算得D=d21≠0，則當F≠0時，E*（x*，z*）是一個余維二的尖點；而當F=0時，E*（x*，z*）是一個至少余維三的尖點。

2" 數值模擬

通過MATLAB軟件包對所獲得的結論進行數值驗證。取p1=3 462/2 500+0.1，p2=23，p3=7/10，p4=194/12 000，p5=13/25，滿足條件p1lt;p*1，得到E*（x*，z*）是包含一個穩(wěn)定拋物扇形的鞍結點（見圖2（a））；取p1=57 598/237 862+0.1，p2=58/91，p3=78/62，p4=86/91-0.046，p5=65/93，滿足條件p1gt;p*1，得到E*（x*，z*）是包含一個不穩(wěn)定拋物扇形的鞍結點（見圖2（b））。

（a） p1lt;p*1時，E*（x*，z*）是包含一個穩(wěn)定拋物扇形的鞍結點

（b） p1gt;p*1時，E*（x*，z*）是包含一個不穩(wěn)定拋物扇形的鞍結點

3" 結語

本文主要研究了一類描述病原體與宿主免疫反應之間相互作用的免疫應答模型，并討論了它的正平衡點及其類型，得到該系統至多存在三個正平衡點，當其中兩個平衡點合為二重平衡點時，該二重平衡點可能為余維二尖點或余維三尖點。最后通過MATLAB軟件包對所獲得的結論進行了數值驗證。

［" 參" 考" 文" 獻" ］

WANG Dawei，HU Bo，HU Chang，et al.Clinical characteristics of 138 hospitalized patients with 2019 novel coronavirus-infected pneumonia in Wuhan，China.Jama，2020，323（11）：1061-1069.

PUGLIESE A，GANDOLFI A.A simple model of pathogen-immune dynamics including specific and nonspecific immuni-ty.Mathematical Biosciences，2008，214（1/2）：73-80.

D’ONOFRIO A.On the interaction between the immune system and an exponentially replicating pathogen.Mathematical Biosciences And Engineering，2010（7）：579-602.

YOUNG T R，BUCKALEW R，MAY A K，et al.A low dimensional dynamical model ofthe initial pulmonary innate reponse to infection.Mathematical Biosciences，2012，235（2）：189-200.

SHI Shujing，HUANG Jicai，WEN Jing.Bifurcation Analysis of a Dynamical Model for the Innate Immune Response to Initial Pulmonary Infections.International Journal of Bifurcation and Chaos，2020，30（16）：1-22.

張錦炎，馮貝葉.常微分方程幾何理論與分支問題.北京：北京大學出版社，1997：329-339.

張芷芬，丁同仁，黃文灶，等.微分方程定性理論.北京：科學出版社，1997：125-140.

閆夢月.腫瘤細胞與效應細胞相互作用模型的動力學分析.武漢：華中師范大學，2022.

PERKO L.Differential Equations and Dynamical Systems.New York： Springer，2001：147-163.

［責任編輯：魏 強］

Bifurcation analysis of a class of immune response models

SUN Dening，" ZHANG Cunhua

School of Mathematics and Physics， Lanzhou Jiaotong University， Lanzhou 730070， China

Abstract：" We mainly study the immune response model of the interaction between a class of pathogens and neutrophils（PMNs）， and analyze the conditions and type of their positive equilibrium point in detail， and carry out numerical simulation of the conclusions obtained by MATLAB software package.

Key words：" immune response model; positive equilibrium point; sharp points