一類PDGF誘導(dǎo)的腫瘤模型的動力學(xué)性質(zhì)分析

摘要： 考慮一個(gè)以血小板源性生長因子（PDGF）驅(qū)動的反應(yīng)擴(kuò)散神經(jīng)膠質(zhì)瘤數(shù)學(xué)模型. 首先， 對常微分系統(tǒng)給出其平衡點(diǎn)的穩(wěn)定性分析， 以趨化劑產(chǎn)生的速率m作為分支參數(shù)給出正平衡點(diǎn)附近Hopf分支的存在性， 并通過規(guī)范型理論和中心流形定理給出判斷由Hopf 分支產(chǎn)生的周期解穩(wěn)定性公式; 其次， 對反應(yīng)擴(kuò)散系統(tǒng)， 得到當(dāng)擴(kuò)散介入后平衡點(diǎn)不會發(fā)生Turing不穩(wěn)定性； 最后， 通過數(shù)值模擬驗(yàn)證理論分析結(jié)果. 結(jié)果表明， 趨化劑產(chǎn)生的速率m可區(qū)分神經(jīng)膠質(zhì)瘤的類型.

關(guān)鍵詞： 腫瘤模型； 反應(yīng)擴(kuò)散； Hopf分支； 穩(wěn)定性

中圖分類號： O175.21" 文獻(xiàn)標(biāo)志碼： A" 文章編號： 1671-5489（2024）04-0809-12

Dynamical Properties Analysis of a Class ofPDGF-Induced Tumor Models

E Xiqi， WEI Xin， ZHAO Jiantao

（School of Mathematical Science， Heilongjiang University， Harbin 150080， China）

Abstract： We considered a platelet derived growth factor （PDGF） driven reaction-diffusion glioma mathematical model.

Firstly， we gave the stability analysis of the equilibrium point for the ordinary differential system. We took the" rate m generated by chemoattractant as" the bifurcation parameter，

gave the existence of the Hopf bifurcation near the positive equilibrium point， and then gave a formula to judge the stability of the periodic solution produced by the Hopf bifurcation through the gauge

type theory and the central manifold theorem. Secondly， for reaction＼|diffusion systems， we obtained that the equilibrium point" did not occur Turing instability" when diffusion was involved.

Finally， the" theoretical analysis results were verified through numerical simulation. The results show that the rate m generated by chemoattractant can be used to distinguish the types of glioma.

Keywords： tumor model; reaction diffusion; Hopf bifurcation; stability

0 引 言

惡性腫瘤會對人類造成很大危害， 是一類致死率較高的疾病［1］. 通過建立腫瘤數(shù)學(xué)模型可以更好地理解腫瘤增長規(guī)律， 并為腫瘤研究和臨床治療提供科學(xué)依據(jù). 目前， 關(guān)于腫瘤數(shù)學(xué)模型的研究已有很多結(jié)果. 例如， 文獻(xiàn)［2］利用穩(wěn)定性和分支理論分析了相應(yīng)腫瘤模型的動力學(xué)行為. 神經(jīng)膠質(zhì)瘤細(xì)胞是一種惡性腫瘤細(xì)胞， 研究表明， 神經(jīng)膠質(zhì)瘤的增殖和侵襲與一些細(xì)胞生長因子的異常表達(dá)有關(guān)， 其中血小板源性生長因子（platelet derived growth factor， PDGF）是神經(jīng)膠質(zhì)瘤中最常見的驅(qū)動因子之一， 它可以促進(jìn)神經(jīng)膠質(zhì)瘤細(xì)胞的增殖和浸潤.

近年來， 關(guān)于PDGF驅(qū)動的神經(jīng)膠質(zhì)瘤的數(shù)學(xué)模型研究得到廣泛關(guān)注. 例如： 文獻(xiàn)［3］建立了一個(gè)腫瘤細(xì)胞和免疫細(xì)胞相互作用的模型， 模擬了免疫細(xì)胞在不同濃度趨化劑下對腫瘤細(xì)胞的殺傷效果; 文獻(xiàn)［4］建立了一個(gè)考慮PDGF濃度和細(xì)胞自身運(yùn)動的模型， 模擬了PDGF對神經(jīng)膠質(zhì)瘤細(xì)胞增殖和浸潤的影響; 文獻(xiàn)［5］考慮一個(gè)徑向?qū)ΨQ的腫瘤， 建立了由PDGF驅(qū)動的三維膠質(zhì)瘤的數(shù)學(xué)模型， 并通過數(shù)值模擬發(fā)現(xiàn)了趨化劑產(chǎn)生率和趨化系數(shù)這兩個(gè)參數(shù)的重要作用； 文獻(xiàn)［6］研究了一個(gè)關(guān)于PDGF驅(qū)動的神經(jīng)膠質(zhì)瘤的隨機(jī)數(shù)學(xué)模型， 可根據(jù)趨化劑產(chǎn)生率和趨化系數(shù)的參數(shù)范圍區(qū)分兩種膠質(zhì)瘤， 且當(dāng)噪聲強(qiáng)度比趨化系數(shù)大時(shí)， 只會有一種膠質(zhì)瘤， 最后通過數(shù)值模擬進(jìn)行了驗(yàn)證， 并給出了相應(yīng)的醫(yī)學(xué)解釋.

本文在文獻(xiàn)［5-6］的基礎(chǔ)上， 考慮如下一類由PDGF誘導(dǎo)的反應(yīng)擴(kuò)散腫瘤模型：

Gt=dGΔG+λG1-G+Nk，x∈Ω， tgt;0

At=dAΔA+mGβ+G-γA，x∈Ω， tgt;0，

Nt=dNΔN+αAN-ρN，x∈Ω， tgt;0，

Gν=Aν=Nν=0，x∈Ω， tgt;0，

G（x，0）=G0（x）， A（x，0）=A0（x）， N（x，0）=N0（x），x∈Ω，（1）

其中： Ω瘙綆N且具有光滑邊界Ω; Δ為Laplace算子; G=G（x，t）表示t時(shí)刻x位置神經(jīng)膠質(zhì)瘤細(xì)胞的密度; A=A（x，t）表示t時(shí)刻x位置趨化劑的濃度; N=N（x，t）表示t時(shí)刻x位置浸潤免疫細(xì)胞的密度; dGgt;0， dAgt;0和dNgt;0分別表示膠質(zhì)瘤細(xì)胞、 趨化劑濃度和浸潤免疫細(xì)胞的擴(kuò)散率; λ表示神經(jīng)膠質(zhì)瘤細(xì)胞的自然增長率； k為最大環(huán)境容納量； m和α分別表示神經(jīng)膠質(zhì)瘤細(xì)胞產(chǎn)生趨化劑的速率和趨化系數(shù)； β表示的Michaelis常數(shù)； γ表示趨化劑的降解速率； ρ是浸潤免疫細(xì)胞的清除率. 假設(shè)神經(jīng)膠質(zhì)瘤細(xì)胞的增殖遵循Logistic增長模式， 趨化劑是由神經(jīng)膠質(zhì)瘤細(xì)胞產(chǎn)生的， 用Michaelis-Menten形式mGβ+G模擬趨化劑的產(chǎn)生速率. 模型中參數(shù)均為正數(shù).

目前， Hopf分支是生物模型研究中的一個(gè)熱點(diǎn)問題＼. 本文以趨化劑產(chǎn)生的速率m作為唯一的分支參數(shù)討論常微分系統(tǒng)中正平衡解的穩(wěn)定性和Hopf分支， 并給出其周期解的穩(wěn)定性及分支方向的判別式， 該方法不同于文獻(xiàn)［5-6］以趨化劑產(chǎn)生率和趨化系數(shù)兩個(gè)參數(shù)研究神經(jīng)膠質(zhì)瘤模型的動力學(xué)行為. 本文的主要理論結(jié)果如下. 首先， 對任意正常數(shù)m， 平衡解E0總是不穩(wěn)定的. 其次， 當(dāng)mlt;m0時(shí)， 平衡解E1是局部漸近穩(wěn)定的; 當(dāng)m≥m0時(shí)， E1是不穩(wěn)定的. 此外， 當(dāng)mgt;m0時(shí)， 系統(tǒng)會存在一個(gè)唯一的正平衡點(diǎn)E*， 當(dāng)m∈（m0，m1）時(shí)， E*是局部漸近穩(wěn)定的; 當(dāng)m∈（m1，+∞）時(shí)， E*是不穩(wěn)定的； 當(dāng)m=m1時(shí)， 在E*處會發(fā)生Hopf分支， 并通過規(guī)范型理論和中心流形定理給出判斷由Hopf分支產(chǎn)生的周期解穩(wěn)定性和分支方向的判別公式. 最后， 證明反應(yīng)擴(kuò)散系統(tǒng)（1）的平衡解不會發(fā)生Turing不穩(wěn)定性. 一般地， 神經(jīng)膠質(zhì)瘤的類型可分為wtIDH1型和muIDH1型， 其中wtIDH1型神經(jīng)膠質(zhì)瘤會有較多的免疫細(xì)胞浸潤［5］. 本文研究的理論結(jié)果和數(shù)值模擬在醫(yī)學(xué)上可解釋為： 當(dāng)其他參數(shù)值確定時(shí)， 可通過趨化劑產(chǎn)生的速率m區(qū)分兩種神經(jīng)膠質(zhì)瘤的類型. 即當(dāng)m小于臨界值m0時(shí)， 浸潤免疫細(xì)胞密度最后會消失為零， 從而可判斷神經(jīng)膠質(zhì)瘤的類型為muIDH1型; 當(dāng)m大于臨界值m0時(shí)， 浸潤免疫細(xì)胞會一直存在， 此時(shí)神經(jīng)膠質(zhì)瘤的類型為wtIDH1型. 特別地， 當(dāng)m非常大時(shí)， 神經(jīng)膠質(zhì)瘤中會存在更多的浸潤免疫細(xì)胞. 因此， 本文的研究結(jié)果可在醫(yī)學(xué)上為區(qū)分神經(jīng)膠質(zhì)瘤的類型提供理論依據(jù).

1 常微分系統(tǒng)的動力學(xué)行為

考慮系統(tǒng)（1）對應(yīng)的常微分系統(tǒng)：

dGdt=λG1-G+Nk，

dAdt=mGβ+G-γA，dNdt=αAN-ρN.（2）

通過選擇趨化劑產(chǎn)生的速率m作為系統(tǒng)（2）的主要分支參數(shù)， 討論系統(tǒng)（2）平衡解的穩(wěn)定性和Hopf分支， 以及由Hopf分支產(chǎn)生的周期解的穩(wěn)定性.

1.1 平衡解分析

顯然， E0∶=（0，0，0）和E1∶=k，mkγ（β+k），0分別為系統(tǒng)（2）的平凡解和半平凡解， 且二者總存在. 記

m0=ργ（β+k）kα.

直接計(jì)算可知， 當(dāng)0lt;m≤m0時(shí)， 系統(tǒng)（2）無正平衡解; 當(dāng)mgt;m0時(shí)， 系統(tǒng)（2）存在唯一的正平衡解E*∶=（G*，A*，N*）， 其中

G*∶=γρβαm-γρ， A*∶=ρα， N*∶=k-γρβαm-γρ.

1.2 平衡解E0和E1的穩(wěn)定性分析

定理1 對于系統(tǒng)（2）的平衡解E0和E1有如下結(jié)論：

1） 系統(tǒng)（2）的平凡解E0是不穩(wěn)定的;

2） 當(dāng)mlt;m0時(shí)， 系統(tǒng)（2）的平衡解E1是局部漸近穩(wěn)定的； 當(dāng)m≥m0時(shí)， 系統(tǒng)（2）的平衡解E1是不穩(wěn)定的.

證明： 1） 系統(tǒng)（2）在E0處的Jacobi矩陣為

J（E0）∶=λ00mβ-γ00k-ρ，

其特征值為λgt;0，-γlt;0，-ρlt;0. 因此平凡解E0是不穩(wěn)定的.

2） 系統(tǒng)（2）在E1處的Jacobi矩陣為

J（E1）∶=-λ0-λmβ（β+k）2-γ0

00mαkγ（β+k）-ρ，

其特征值為-λ（lt;0），-γ（lt;0），mαkγ（β+k）-ρ. 顯然， 當(dāng)mlt;m0時(shí)， mαkγ（β+k）-ρlt;0， 系統(tǒng)（2）的平衡解E1是局部漸近穩(wěn)定的； 當(dāng)mgt;m0時(shí)， mα

kγ（β+k）-ρgt;0， 系統(tǒng)（2）的平衡解E1是不穩(wěn)定的. 當(dāng)m=m0時(shí)， 有一特征值為0， 此時(shí)不能通過特征值直接判斷該平衡解的穩(wěn)定性. 下證當(dāng)m=m0時(shí)， 系統(tǒng)（2）的平衡解E1是不穩(wěn)定的. 令

=G-k， =A-mαkγ（β+k）， =N，

則系統(tǒng)（2）轉(zhuǎn)化為

ddt=λ（+k）1-++kk，

ddt=m（+k）β+k+-γ-mkβ+k，

ddt=α.（3）

系統(tǒng)（2）的平衡解E1對應(yīng)系統(tǒng)（3）的平衡解1∶=（0，0，0）. 當(dāng)m=m0時(shí)， 系統(tǒng)（3）在1處的Jacobi矩陣為

J（1）=-λ0-λm0β（β+k）2-γ0

000，

其特征值為-λ，-γ，0. 考慮λ≠γ， 通過計(jì)算可得： 特征值-λ對應(yīng)的特征向量為V1=（β+k）2（γ

-λ）m0β，1，0T， 特征值-γ對應(yīng)的特征向量為V2=（0，1，0）T， 特征值0對應(yīng)的的特征向量為V3=

1，m0βγ（β+k）2），-1T.

令T=（V1，V2，V3）， X=（，，）T， 則系統(tǒng)（3）可改寫為

dXdt=J（1）X+F，

其中

F=-λk2-λk，-

mkβ+k+m（+k）β+k+-mβ（β+k）2，αT.

令X=TY， Y=（y1，y2，y3）T， 則有

dYdt=T-1J（1）TY+T-1F，

其中T-1J（1）T=diag（-λ，-γ，0）. 又由X=TY，" Y=（y1，y2，y3）T， 可得

=（β+k）2（γ-λ）mβy1+y3， =y1+y2+mβγ（β+k）2y3， =-y3.

下面計(jì)算T-1F. 直接計(jì)算可得

T=-（β+k）2（γ-λ）mβ，

T*=-10-11-（β+k）2（γ-λ）mβλγ

00（β+k）2（γ-λ）mβ，

T-1=mβ（β+k）2（γ-λ）0mβ（β+k）2（γ-λ）

-mβ（β+k）2（γ-λ）1-λmβγ（β+k）2（γ-λ）00-1，

從而可得

T-1F=mβ（β+k）2（γ-λ）0mβ（β+k）2（γ-λ）

-mβ（β+k）2（γ-λ）1-λmβγ（β+k）2（γ-λ）00-1

-λk2-λk-mkβ+k+m（+k）β+k+-

mβ（β+k）2）α.

記T-1F=（f1，f2，f3）T， 則有

dZdt=BZ+f1f2，

dy3dt=0y3+f3，（4）

其中B=diag（-λ，-γ）， Z=（y1，y2）T. 顯然fi∈C2， fi（0，0，0）=0， Dfi（0，0，0）=0， i=1，2，3. 由中心流形定理可知， 存

在一個(gè)中心流形Z=h（y3）=（h1（y3），h2（y3））T， 滿足

Bh（y3）+f1（h（y3），y3）

f2（h（y3），y3）=Dh（y3）f3（h（y3），y3）.

假設(shè)y1=h1（y3）=e2y23+e3y33+o（y33）， y2=h2（y3）=r2y23+r3y33+o（y33），

則有f3（h（y3），y3）=C33y23+o（y23）， 其中C33=αmβ/（γ（β+k）2）gt;0. 系統(tǒng)（4）零解的動力學(xué)行為可由f3（h（y3），y3）

決定. 由于C33gt;0， 所以y3=0是不穩(wěn)定的. 因此， 當(dāng)m=m0時(shí)， 系統(tǒng)（2）的平衡解E1是不穩(wěn)定的. 證畢.

1.3 正平衡點(diǎn)的穩(wěn)定性及Hopf分支分析

系統(tǒng)（2）在E*=（G*，A*，N*）處的Jacobi矩陣為

J（m）∶=-λG*k0-λG*kmβ（G*+β）2-γ0

0αN*0，（5）

J（m）的特征方程為

x3+a2（m）x2+a1（m）x+a0（m）=0，（6）

其中

a2（m）∶=λG*k+γ=λγρβk（αm-γρ）+γ，

a1（m）∶=λγG*k=λγ2ρβk（αm-γρ），

a0（m）∶=αβmλG*N*k（β+G*）2=λγρ［k（αm-γρ）-βγρ］kαm.

當(dāng)mgt;m0時(shí)， E*是系統(tǒng)（2）的唯一正平衡解. 定義

Γ1={m： a1（m）a2（m）-a0（m）gt;0}，Γ2={m： a1（m）a2（m）-a0（m）lt;0}，

Γ3={m： a1（m）a2（m）-a0（m）=0， a′1（m）a2（m）+a1（m）a′2（m）-a′0（m）≠0}.（7）

引理1 假設(shè)mgt;m0， 則有下列結(jié)論：

1） 如果m∈Γ1， 則E*是局部漸近穩(wěn)定的;

2） 如果m∈Γ2， 則E*是不穩(wěn)定的;

3） 如果m∈Γ3， 則系統(tǒng)（2）在E*處產(chǎn)生Hopf分支.

證明： 1） 根據(jù)Routh-Hurwitz準(zhǔn)則知， E*是局部漸近穩(wěn)定的當(dāng)且僅當(dāng)

a0（m）gt;0， a2（m）gt;0， a1（m）a2（m）-a0（m）gt;0.

由mgt;m0， 有a0（m）gt;0和a2（m）gt;0， 因此如果m∈Γ1， 則E* 是局部漸近穩(wěn)定的.

2） 如果m∈Γ2， 則特征方程（6）會存在具有正實(shí)部的根， 因此E*在Γ2中是不穩(wěn)定的.

3） 假設(shè)存在m*∈Γ3， 則有a1（m*）a2（m*）-a0（m*）=0， 特征方程（6）可化為

（x2（m*）+a1（m*））（x（m*）+a2（m*））=0，

它有3個(gè)根： x1=ia1（m*）， x2=-ia1（m*）， x3=-a2（m*）lt;0

. 因此特征方程（6）有一對純虛根和一個(gè)負(fù)根. 對任意的m（mgt;m0）， 假設(shè)特征方程（6）的根為

x1（m）=α（m）+iβ（m）， x2（m）=α（m）-iβ（m）， x3（m）=-a2（m），

則顯然有α（m*）=0和β（m*）=a1（m*）.

下面證明橫截條件α′（m*）≠0成立. 對特征方程（6）關(guān)于m求導(dǎo)， 并分離實(shí)部和虛部， 可得

A1（m）α′（m）-A2（m）β′（m）+A3（m）=0，A2（m）α′（m）+A1（m）β′（m）+A4（m）=0，

其中

A1（m）=3α2（m）+2a2（m）α（m）+a1（m）-3β2（m），A2（m）=6α（m）β（m）+2a2（m）β（m），

A3（m）=α2（m）a′2（m）+a′1（m）α（m）+a′0（m）-a′2（m）β2（m），

A4（m）=2α（m）β（m）a′2（m）+a′1（m）β（m）.

注意到α（m*）=0和β（m*）=a1（m*）， 有

A1（m*）=-2a1（m*），" A2（m*）=2a2（m*）a1（m*），A3（m*）=a′0（m*）-a′2（m）a1（

m*），" A4（m*）=a′1（m）a1（m*）.

經(jīng)計(jì)算可得

α′（m*）=-a′1（m*）a2（m*）+a1（m*）a′

2（m*）-a′0（m*）2（a22（m*）+a1（m*））.

由于a1（m）gt;0， a2（m）gt;0， 因此α′（m*）≠0等價(jià)于

a′1（m*）a2（m*）+a1（m*）a′2（m*）-a′0（m*）≠0.

從而根據(jù)Hopf分支理論［10］可知， 當(dāng)m=m*∈Γ3時(shí)， 系統(tǒng)（2）在E*處產(chǎn)生Hopf分支. 證畢.

根據(jù)引理1， 可得如下關(guān)于E*的穩(wěn)定性和Hopf分支的定理.

定理2 假設(shè)mgt;m0成立， 則：

1） 當(dāng)m∈（m0，m1）時(shí)， 系統(tǒng)（2）的平衡解E*是局部漸近穩(wěn)定的;

2） 當(dāng)m∈（m1，+∞）時(shí)， 系統(tǒng)（2）的平衡解E*是不穩(wěn)定的;

3） 當(dāng)m=m1時(shí)， 系統(tǒng)（2）在平衡解E*處經(jīng)歷Hopf分支.

其中定義m1∶=x1+γραgt;m0.

證明： 經(jīng)計(jì)算可得

a1（m）a2（m）-a0（m）=λγpk2mα（αm-γp）2F（m），

其中

F（m）=αmρβ2γ2λ+kαmβγ2（αm-γρ）-k（αm-γρ）2［（αm-γρ）-γβρ］.

因?yàn)閙gt;m0， 所以有λγρk2mα（αm-γρ）2gt;0. 因此當(dāng)m∈（m0，+∞）時(shí)， a1（m）a2（m）-a0（m）的符號與F（m）的符號相同. 令x=αm-γρ， 則有

F（m）∶=T（x）=-k2x3+kβγ（γ+ρ）x2+（kβγ3ρ+β2 γ2λρ）x+β2γ3λρ2 .

由于mgt;m0等價(jià)于xgt;γρβk， 因此下面討論當(dāng)x∈γρβk，+∞時(shí)T（x）根的情況.

顯然， limx→+∞ T（x）=-∞. 另一方面， 有

Tγρβk=-γ3ρ3β3k+

γ3ρ3β3+γ4ρ2β3k+kγ4ρ2β2+λγ3ρ2β3k+λγ3ρ2β2

=γ3β2ρ2（γ+λ）βk+1gt;0.

由T（x）在γρβk，+∞中的連續(xù)性知， T（x）在γρβk，+∞中至少存在一個(gè)根， 記為x1. 此外， 對于

多項(xiàng)式函數(shù)T（x）， 易見x3的系數(shù)為負(fù)， 而x2，x和x0的系數(shù)均為正. 相鄰的非零系數(shù)符號的變化次數(shù)為1次， 由笛卡爾符號法則可知， T（x）

最多有一個(gè)正根. 根據(jù)上述分析可知， T（x）在γρβk，+∞中存在唯一的正根. 此外， 根據(jù)三次函數(shù)圖像的性質(zhì)， 有

T（x）gt;0，x∈γρβk，x1，lt;0，x∈（x1，+∞），" T′（x1）lt;0.

根據(jù)m1的定義， 有

a1（m1）a2（m1）-a0（m1）=0，a1（m）a2（m）-a0（m）

gt;0，m∈（m0，m1），lt;0，m∈（m1，+∞）.

因?yàn)镕′（m1）=αT′（x1）lt;0， 從而有

a′1（m1）a2（m1）+a1（m1）a′2（m1）-a′0（m1）=λγpk2m1α（αm1-γp）2F′（m1）lt;0.

由式（7）定義可知， Γ1=（m0，m1）， Γ2=（m1，+∞）， Γ3={m1}. 由引理1可得結(jié)論. 證畢.

1.4 Hopf分支周期解的穩(wěn)定性分析

下面討論系統(tǒng)（2）由Hopf分支產(chǎn)生的周期解（Gp（t），Ap（t），Np（t））的穩(wěn)定性.

令y1=G-G*， y2=A-A*， y3=N-N*， 則系統(tǒng)（2）可改寫為

dy1dt=λ（y1+G*）1-y1+y3+G*+N*k，

dy2dt=m（y1+G*）β+y1+G*-γ（y2+A*），

dy3dt=α（y2+A*）（y3+N*）-ρ（y3+N*）.（8）

系統(tǒng)（8）可記為如下形式：

dYdt=J（m）Y+F（Y），（9）

這里J（m）由式（5）定義， Y=（y1，y2，y3）T， F（Y）=（F1（Y），F(xiàn)2（Y），F(xiàn)3（Y））T， 其中

F1（Y）∶=-λky1（y1+y3），F(xiàn)2（Y）∶=mβ+y1+G*y1-mβ

（β+G*）2y1+mG*β+y1+G*-γA*，F(xiàn)3（Y）∶=αy2y3.（10）

如果m=m1， 則J（m）的特征值為

x1（m1）=iω0， x2（m1）=-iω0， x3（m1）=-a2（m1）lt;0.

令ξ1和ξ2分別為J（m1）的特征值iω0和-a2（m1）對應(yīng)的特征向量. 經(jīng)計(jì)算可得

ξ1∶=1γ-ω0iαN*-1-kω0iλG*

，" ξ2∶=λG*kγ-a2（m1）αN*1.

定義

P∶=（Re（ξ1），-Im（ξ1），ξ2）=

10λG*kγγαN*ω0αN*-a2（m1）αN*

-1kω0λG*1，

則

P-1=1detPp11p12p13

p21p22p23p31p32p33，

其中

p11=ω0αN*+kω0a2（m1）λαN*G*， p12=ω0γ，

p13=-λω0G*kγαN*，p21=a2（m1）-γαN*， p22=1+λG*kγ， p23=a2（m1）αN*

+λG*kαN*，p31=kγω0λαN*G*+ω0αN*， p32=-kω0λG*，

p33=ω0αN*，detP=ω0αN*3+

kγλG*+λG*kγ.

令Y=PZ， 其中Z∶=（z1，z2，z3）T， 則系統(tǒng)（9）可寫為

dZdt=P-1J（m1）PZ+G（Z），

其中

P-1J（m1）P=0-ω00ω00000-a2（m1），

G（Z）∶=（G1（Z），G2（Z），G3（Z））T=P-1F（PZ）.

記F（PZ）=（f1（Z），f2（Z），f3（Z））T， 則由式（10）， 可得

f1（Z）=-λkz1+λG*kγz3kω0λG*z2+1+λG*kγ

z3，f2（Z）=m1G*+z1+λG*kγ

z3β+G*+z1+λG*kγz3-γA*-m1β（β+G*）2z1+λG*kγz3，

f3（Z）=1N*（γz1+ω0z2-a2（m1）z3）-z1+kω0λG*z2+z3.

進(jìn)一步計(jì)算可得

G1（Z）=p11f1（Z）+p12f2（Z）+p13f3（Z）detP，

G2（Z）=p21f1（Z）+p22f2（Z）+p23f3（Z）detP，

G3（Z）=p31f1（Z）+p32f2（Z）+

p33f3（Z）detP.

為研究Hopf分支周期解的穩(wěn)定性， 需要在m=m1和（z1，z2，z3）=（0，0，0）時(shí)計(jì)算以下變量［10-11］：

g11=14［（G1）z1z1+（G1）z2z2+i（（G2）z1z1+（G2）z2z2）］，

g02=14［（G1）z1z1-（G1）z2z2-2（G2）z1z2+i（（G2）z1z1-（G2）z2z2+2（G1）z1z2）］，

g20=14［（G1）z1z1-（G1）z2z2+2（G2）z1z2

+i（（G2）z1z1-（G2）z2z2-2（G1）z1z2）］，

τ11=14a2（m1）（（G3）z1z1+（G3）z2z2），

τ20=14（a2（m1）+2iω0）-1（（G3）z1z1-（G3）z2z2-2i（G3）z1z2），

G110= "12［（G1）z1z3+（G2）z2z3+i（（G2）z1z3-（G1）z2z3）］，

G101= "12［（G1）z1z3-（G2）z2z3+i（（G2）z1z3+（G1）z2z3）］，

G21= "18［（G1）z1z1z1+（G1）z1z2z2+（G2）z1z1z2+（G2）z2z2z2

+ "i（（G2）z1z1z1+（G2）z1z2z2-（G1）z1z1z2-（G1）z2z2z2）］，

g21= "G21+2G110τ11+G101τ20.

由文獻(xiàn)［10］， 定義第一Lyapunov系數(shù)c1（m1）為

c1（m1）=i2ω0g11g20-2g112

-g0223+g212.

因此

Re（c1（m1））=12Re（g21）-12ω0（Re（g11）Im（g20

）+Im（g11）Re（g20））.（11）

由定理2的證明， 可知

a1（m1）a′2（m1）+a′1（m1）a2（m1）-a′0（m1）lt;0，

故α′（m1）gt;0. 綜上可得如下定理.

定理3 對于系統(tǒng)（2）的Hopf分支周期解， 有如下結(jié)論：

1） 如果Re（c1（m1））lt;0， 則Hopf分支周期解是軌道漸近穩(wěn)定的， 且分支方向是向前的;

2） 如果Re（c1（m1））gt;0， 則Hopf分支周期解是不穩(wěn)定的， 且分支方向是向后的.

2 反應(yīng)擴(kuò)散系統(tǒng)的動力學(xué)行為

文獻(xiàn)［12］研究表明， 擴(kuò)散能使一個(gè)系統(tǒng)的穩(wěn)態(tài)解失去穩(wěn)定性， 這種不穩(wěn)定性即為由擴(kuò)散引起的不穩(wěn)定性或Turing不穩(wěn)定性. 下面討論在常

微分系統(tǒng)（2）中穩(wěn)定的常數(shù)平衡解， 當(dāng)擴(kuò)散介入后其穩(wěn)定性是否發(fā)生變化， 即研究系統(tǒng)（1）常數(shù)穩(wěn)態(tài)解是否會發(fā)生Turing不穩(wěn)定性.

定理4 對于系統(tǒng)（1）， 有如下結(jié)論：

1） 如果m∈（0，m0）， 則系統(tǒng)（1）的平衡解E1是局部漸近穩(wěn)定的;

2） 如果m∈（m0，m1）， 則系統(tǒng)（1）的平衡解E*是局部漸近穩(wěn)定的.

證明： 系統(tǒng)（1）在E1處線性化系統(tǒng)系數(shù)矩陣為

（E1）∶=-λ+dGΔ0-λmβ（β+k）2-γ+dAΔ0

00mαkγ（β+k）-ρ+dNΔ，

根據(jù)文獻(xiàn)［13-14］可知， （E1） 的特征值可由n確定：

n∶=-λ-dGμn0-λmβ（β+k）2-γ-dAμn0

00mαkγ（β+k）-ρ-dNμn，

其中0=μ0lt;μ1lt;…μnlt;…為-Δ在區(qū)域Ω上滿足齊次Neumann邊界條件的特征值. 當(dāng)m∈（0，m0）時(shí)， 對任意n∈瘙綃， n的特征值為

-λ-dGμnlt;0， -γ-dAμnlt;0， mαk-kγα-βγργ（β+k）-dNμnlt;0.

所以系統(tǒng)（1）的平衡解E1是局部漸近穩(wěn)定的.

系統(tǒng)（1）在E*處線性化系統(tǒng)系數(shù)矩陣為

*∶=a11（m）+dGΔ0a13（m）

a21（m）a22（m）+dAΔ00a32（m）dNΔ.

同理， （E*） 的特征值可由*n確定：

*n∶=a11（m）-dGμn0a13（m）a21（m）a22（m）-dAμn0

0a32（m）-dN μn，

*n的特征方程為

x3+a2，n（m）x2+a1，n（m）x+a0，n（m）=0，

其中

a2，n（m）∶= "（dG+dA+dN）μn+a2（m），a1，n（m）∶= "（dGdA+dGdN+dAdN）μ2n

- "［a22（m）dG+a11（m）dA+（a11（m）+a22（m））dN］μn+a1（m），

a0，n（m）∶= "dGdAdNu3n-（dGdNa22（m）+dAdNa11（m））μ2n+（a11（m）a22（m）dN）μn+a0（m），

這里a2（m），a1（m）和a0（m）與式（6）相同. 當(dāng)m∈（m0，m1）時(shí)，

a0（m）gt;0， a1（m）gt;0， a2（m）gt;0， a1（m）a2（m）-a0（m）gt;0.

注意到dGgt;0， dAgt;0， dNgt;0， a11（m）lt;0， a22（m）lt;0， 可得

a2，n（m）gt;0， a1，n（m）gt;0， a0，n（m）gt;0.

直接計(jì)算可得

a2，n（m）a1，n（m）-a0，n（m）∶=h1μ3n+h2μ2n+h3μn+h4，

其中

h1= "（dG+dA）（dG+dN）（dA+dN），h2=-（dA+dN）（2dG+dA+dN）a11（m）-a22（m）（dG+dN）（dG+2dA+dN），

h3= "［2a11（m）a22（m）+a222（m）］dG+［2a11（m）a22（m）+a211（m）］dA+［（a11（m）+a22（m））2］dN，

h4= "a2（m）a1（m）-a0（m）.

當(dāng)m∈（m0，m1）時(shí)， 易得higt;0（i=1，2，3，4）， 從而a2，n（m）a1，n（m）-a0，n（m）gt;0. 所以當(dāng)m∈（m0，m1）時(shí)， 系統(tǒng)（1）的平衡解E*是局部漸近穩(wěn)定的. 證畢.

3 數(shù)值模擬

下面通過數(shù)值模擬驗(yàn)證上述理論分析結(jié)果. 選取如下參數(shù)值：

λ=0.48， β=0.1， γ=2.185， α=0.6， k=1， ρ=0.9.

經(jīng)計(jì)算得m0=3.602 5， m1=5.589 1.

對于常微分系統(tǒng)（2）， 根據(jù)定理1和定理2知， 若m∈（0，3.602 5）， 則系統(tǒng)（2）的平衡解E1是

局部漸近穩(wěn)定的; 若m∈（3.602 5，5.589 1）， 則系統(tǒng)（2）的平衡解E*是局部漸近穩(wěn)定的; 若mgt;5.589 1， 則系統(tǒng)（2）的平衡解E*是不穩(wěn)定的. 當(dāng)m=5.589 1時(shí)， 系統(tǒng)（2）會產(chǎn)生Hopf

分支， 在E*附近分支出一族周期解， 此時(shí)E*=（0.141 8，1.5，0.858 2）. 利用1.4節(jié)中給出的計(jì)算公式， 計(jì)算可得：

g11=-0.102 0+0.036 5i，" g02=0.211 6-0.332 5i，g20=-0.041 0+0.141 1i，" g21=-0.042 5-0.007 4i，

τ11=0.161 0，" τ20=-0.150 8+0.050 0i，G21=-0.052 9+0.091 3i，" G110=-0.331 2+0.095 4i，

G101=-0.233 6+0.000 036i.

從而Re（c1（m1））=-0.041 6lt;0. 由定理3知， 分支周期解是軌道漸近穩(wěn)定的， 且分支方向是向前的. 取m=5.6， 可觀察到系統(tǒng)（2）存在周期解， 如圖1所示.

選取dG=10， dA=1， dN=20. 當(dāng)m=2∈（0，m0）時(shí)， 根據(jù)定理4， 系統(tǒng)（1）的平衡解E1=（1，0.832 1，0）是局部漸近穩(wěn)定的， 如圖2所示.

如果選取m=4∈（m0，m1）， 根據(jù)定理4知， 系統(tǒng)（1）的平衡解E*=（0.453 6，1.5，0.546 4）是局部漸近穩(wěn)定的， 如圖3所示.

圖1 當(dāng)m=5.6時(shí)系統(tǒng)（2）產(chǎn)生的穩(wěn)定周期解

Fig.1 Stable periodic solutions generated by system （2） when m=5.6

圖2 當(dāng)m=2時(shí)系統(tǒng)（1）的平衡解E1

=（1，0.832 1，0）的局部漸近穩(wěn)定性Fig.2 Local asymptotic stability of equilibrium solution E1=（1，0.832 1，0） of system （1） when m=2

圖3 當(dāng)m=4時(shí)系統(tǒng)（1）平衡解E*=（0.453 6，1.5，0.546 4）的局部漸近穩(wěn)定性

Fig.3 Local asymptotic stability of equilibrium solution E*=（0.453 6，1.5，0.546 4） of system （1） when m=4

綜上， 本文在齊次Neumann邊界條件下研究了一類PDGF誘導(dǎo)的腫瘤反應(yīng)擴(kuò)散模型. 首先， 通過以趨化劑產(chǎn)生的速率m為分支參數(shù)， 給出了該模型在正平衡點(diǎn)附近Hopf分支的存在性

， 并通過規(guī)范型理論和中心流形定理給出了判斷由Hopf分支產(chǎn)生的周期解穩(wěn)定性的公式. 其次， 對反應(yīng)擴(kuò)散系統(tǒng)的研究表明， 當(dāng)擴(kuò)散介入后穩(wěn)態(tài)解不會產(chǎn)生Turing不穩(wěn)定性.

因此經(jīng)典的自擴(kuò)散不能導(dǎo)致Turing斑圖的出現(xiàn)， 數(shù)值模擬也證實(shí)了該結(jié)論.

參考文獻(xiàn)

［1］ 吳菲， 林國楨， 張晉昕. 我國惡性腫瘤發(fā)病現(xiàn)狀及趨勢 ［J］. 中國腫瘤， 2012， 21（2）： 81-85. （WU F， LIN G Z， ZHANG J X. An Overview

of Cancer Incidence and Trend in China ［J］. China Cancer， 2012， 21（2）： 81-85.）

［2］ 楊艷紅， 伏升茂. 一類帶時(shí)滯的腫瘤免疫模型的Hopf分支 ［J］. 數(shù)學(xué)進(jìn)展， 2021， 50（3）： 38

3-398. （YANG Y H， FU S M. Hopf Bifurcation of a Tumor Immune Model with Time Delay ［J］. Advances in Mathematics， 2021， 50（3）： 383-398.）

［3］ EFTIMIE R， BRAMSON J L， EARN D J D. Interactions between the Immune System and C

ancer： A Brief Review of Non-spatial Mathematical Models ［J］. Bulletin of Mathematical Biology， 2011， 73： 2-32.

［4］ LEDER K， PITTER K， LAPLANT Q， et al. Mathematical Modeling of PDGF-Driven Gliob

lastoma Reveals Optimized Radiation Dosing Schedules ［J］. Cell， 2014， 156（3）： 603-616.

［5］ NIU B， ZENG X Y， PHAN T A， et al. Mathematical Modeling of PDGF-Driven Glioma Rev

eals the Dynamics of Immune Cells Infiltrating into Tumors ［J］. Neoplasia， 2020， 22（9）： 323-332.

［6］ PHAN T A， NGUYEN H D， TIAN J P. Deterministic and St

ochastic Modeling for PDGF-Driven Gliomas Reveals a Classification of Gliomas ［J］. Journal of Mathematical Biology， 2021， 83（2）： 22-1-22-51.

［7］ 劉琪， 常笑源. 帶有Allee效應(yīng)的擴(kuò)散時(shí)滯單種群模型的分支分析 ＼. 東北師大學(xué)報(bào)（自然科學(xué)版）， 2023， 55（2）： 11＼|15.

（LIU Q， CHANG X Y. Bifurcation Analysis of a Diffusive Single Population Model with Allee Effect" and Time Delay ＼. Journal of Northeast Normal

University （Natural Science Edition）， 2023， 55（2）： 11＼|15.）

［8］ 高鶴， 李秀玲. 二維具時(shí)滯捕食＼|食餌共生模型的Hopf分支 ＼. 東北師大學(xué)報(bào)（自然科學(xué)版）， 2024， 56（1）： 23＼|28.

（GAO H， LI X L. Hopf Bifurcation of a Two＼|Dimensional Predator＼|Prey Symbiotic Model with Time Delay ＼. Journal of Northeast Normal University

（Natural Science Edition）， 2024， 56（1）： 23＼|28.）

［9］ 趙浛弛， 李杰梅." 一類腫瘤＼|免疫模型的穩(wěn)定性與Hopf分支分析 ＼. 吉林大學(xué)學(xué)報(bào)（理學(xué)版）， 2024， 62（2）： 189＼|196.

（ZHAO H C， LI J M. Stability and Hopf Bifurcation Analysis of a Class of Tumor＼|Immune Models ＼. Journal of Jilin University （Science Edition）， 2024， 62（2）： 189＼|196.）

［10］ HASSARD D D， KAZARINOFF N D， WAN Y H. Theory and App

lications of Hopf Bifurcation ［M］. Cambridge： Cambridge University Press， 1981： 86＼|91.

［11］ WANG M， YI F. On the Dynamics of the Diffusive Field

-Noyes Model for the Belousov-Zhabotinskii Reaction ［J］. Journal of Differential Equations， 2022， 318（13）： 443-479.

［12］ TURING A M. The Chemical Basis of Morphogenesis ［J

］. Philosophical Transactions of the Royal Society of London， Series B： Biological Sciences， 1952， 237： 37-72

.［13］ JANG J， NI W M， TANG M. Global Bifurcation and Structure of Turing Patterns in the 1-D Lengyel-Epstein Model ［J］. Journal of Dy

namics and Differential Equations， 2004， 16（2）： 297-320.

［14］ YI F， WEI J， SHI J. Bifurcation and Spatiotemporal Patterns in a Homogenous Diff

usive Predator-Prey System ［J］. Journal of Differential Equations， 2009， 246（5）： 1944-1977.

（責(zé)任編輯： 李 琦）