一類免疫應(yīng)答模型的分支分析

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

關(guān)鍵詞： 免疫應(yīng)答模型；正平衡點(diǎn)；尖點(diǎn)

中圖分類號(hào)： O175

文獻(xiàn)標(biāo)識(shí)碼： A" 文章編號(hào)： 2096-3998（2024）02-0079-07

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

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目（61763024）

作者簡(jiǎn)介：孫得寧（1998—），男，甘肅張掖人，碩士研究生，主要研究方向?yàn)槲⒎址匠膛c動(dòng)力系統(tǒng)；張存華（1972—），女，甘肅武山人，博士，教授，主要研究方向?yàn)槲⒎址匠膛c動(dòng)力系統(tǒng)。

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

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

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

=αx-mx1+bx-βxz，

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

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

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

α＜m，（2）

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

αβ＜δγ，（3）

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

=αx-mx1+bx-βxz，

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

接下來主要討論系統(tǒng)（4）的正平衡點(diǎn)及其類型。

1" 正平衡點(diǎn)及其類型

首先做如下尺度變換：

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

用x、z、t表示x1、z1、t1，則系統(tǒng)（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均為正參數(shù)，因此，僅在下列條件下研究系統(tǒng)（5）：

0lt;p3lt;1

p1gt;1

p2，p4，p5gt;0。（6）

系統(tǒng)的正不變區(qū)域?yàn)?/p>

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

為了討論系統(tǒng)（5）的正平衡點(diǎn)，令

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，

則當(dāng)Δgt;0且blt;0時(shí)， f′（x）有兩個(gè)正根，即f（x）存在極值點(diǎn)μ1、 μ2；當(dāng)Δlt;0時(shí)， f′（x）不存在正根，即f（x）是單調(diào)遞增的。

由式（6）得

p22-p22p3gt;0。

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

（a） 三個(gè)不同的正根：E1，E2，E3""" （b） 兩個(gè)正根：重根E*，單根E3""" （c） 兩個(gè)正根：?jiǎn)胃鵈1，重根E*

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

由圖1可知，首先，當(dāng)f（0）lt;0，blt;0，Δgt;0， f（μ1）gt;0且f（μ2）lt;0（μ2gt;μ1gt;0）時(shí)， f（x）有三個(gè)正根（見圖1（a））。

其次，當(dāng)f（0）lt;0，blt;0，Δgt;0， f（μ1）=0或f（μ2）=0時(shí)， f（x）有一個(gè)正的單根和一個(gè)正的重根（見圖1（b）（c））。

再次，當(dāng)f（0）gt;0，blt;0，Δgt;0， f（μ1）=0時(shí)， f（x）有一個(gè)正的重根（見圖1（d））。

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

接下來研究系統(tǒng)（5）正平衡點(diǎn)的類型，系統(tǒng)（5）在E00，p4p5處的Jacobian矩陣為

J（E0）=

1-p1-p4p50

1-p3p4p5-2p4p5

〗，

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

系統(tǒng)（5）在正平衡點(diǎn)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。

根據(jù)Det（J（E））與f′（x）之間的關(guān)系得到

Det（J（E））=

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

令

p*1=（1+p2x）p2x

則當(dāng)p1lt;p*1（p1gt;p*1）時(shí)，Tr（J（E））gt;0（Tr（J（E））lt;0），由以上結(jié)論得到以下定理。

定理1" 當(dāng)（6）中的條件滿足時(shí)，系統(tǒng)（5）總是存在一個(gè)無病平衡點(diǎn)E00，p4p5，它是穩(wěn)定雙曲結(jié)點(diǎn)。此外，

（Ⅰ）若f（0）gt;0、blt;0、Δgt;0、 f（μ2）gt;0，或f（0）gt;0、Δlt;0，則系統(tǒng)（5）不存在正平衡點(diǎn)；

（Ⅱ）若f（0）gt;0，blt;0，Δgt;0， f（μ2）=0，則系統(tǒng)（5）存在唯一正平衡點(diǎn)E*（x*，z*），它是一個(gè)退化平衡點(diǎn)；

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

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

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

定理2" 當(dāng)f（0）gt;0，blt;0，Δgt;0， f（μ2）=0與（6）中的條件滿足時(shí)，系統(tǒng)（5）存在唯一正平衡點(diǎn)E*（x*，z*）。此外，

（Ⅰ）當(dāng)p1lt;p*1（或p1gt;p*1）時(shí)，E*（x*，z*），是包含一個(gè)穩(wěn)定拋物扇形（或不穩(wěn)定拋物扇形）的鞍結(jié)點(diǎn)；

（Ⅱ）當(dāng)p1=p*1時(shí)，E*（x*，z*）是一個(gè)尖點(diǎn)，此外，

（ⅰ）若F≠0，則E*（x*，z*）是一個(gè)余維二的尖點(diǎn)；

（ⅱ）若F=0，則E*（x*，z*）是一個(gè)至少余維三的尖點(diǎn)。

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

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

那么當(dāng)p1lt;p*1（p1gt;p*1）時(shí)，Tr（J（E*））lt;0（Tr（J（E*））gt;0）。

下面來討論E*（x*，z*）的具體類型，令X=x-x*，Z=z-z*，將平衡點(diǎn)移到原點(diǎn)并將該式進(jìn)行泰勒展開，則系統(tǒng)（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表示τ，則系統(tǒng)（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，

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

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

依據(jù)文獻(xiàn)中的定理7.1，結(jié)論（Ⅰ）成立。

令X=x-x*，Z=z-z*，p1=p*1，Det（J（E*））=0，則系統(tǒng)（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表示τ，則系統(tǒng)（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。

根據(jù)文獻(xiàn)中的引理2.4，用x、z表示μ、v，得到系統(tǒng)（13）在原點(diǎn)的小鄰域內(nèi)等價(jià)于系統(tǒng)（14）：

=z，

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

其中，D=d21，F(xiàn)=d22+2d11。經(jīng)計(jì)算得D=d21≠0，則當(dāng)F≠0時(shí)，E*（x*，z*）是一個(gè)余維二的尖點(diǎn)；而當(dāng)F=0時(shí)，E*（x*，z*）是一個(gè)至少余維三的尖點(diǎn)。

2" 數(shù)值模擬

通過MATLAB軟件包對(duì)所獲得的結(jié)論進(jìn)行數(shù)值驗(yàn)證。取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*）是包含一個(gè)穩(wěn)定拋物扇形的鞍結(jié)點(diǎ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*）是包含一個(gè)不穩(wěn)定拋物扇形的鞍結(jié)點(diǎn)（見圖2（b））。

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

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

3" 結(jié)語(yǔ)

本文主要研究了一類描述病原體與宿主免疫反應(yīng)之間相互作用的免疫應(yīng)答模型，并討論了它的正平衡點(diǎn)及其類型，得到該系統(tǒng)至多存在三個(gè)正平衡點(diǎn)，當(dāng)其中兩個(gè)平衡點(diǎn)合為二重平衡點(diǎn)時(shí)，該二重平衡點(diǎn)可能為余維二尖點(diǎn)或余維三尖點(diǎn)。最后通過MATLAB軟件包對(duì)所獲得的結(jié)論進(jìn)行了數(shù)值驗(yàn)證。

［" 參" 考" 文" 獻(xiàn)" ］

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.

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

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

閆夢(mèng)月.腫瘤細(xì)胞與效應(yīng)細(xì)胞相互作用模型的動(dòng)力學(xué)分析.武漢：華中師范大學(xué)，2022.

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

［責(zé)任編輯：魏 強(qiáng)］

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