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

關鍵詞： 腫瘤模型； 反應擴散； Hopf分支； 穩(wěn)定性

中圖分類號： O175.21" 文獻標志碼： 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 引 言

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

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

本文在文獻［5-6］的基礎上， 考慮如下一類由PDGF誘導的反應擴散腫瘤模型：

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時刻x位置神經膠質瘤細胞的密度; A=A（x，t）表示t時刻x位置趨化劑的濃度; N=N（x，t）表示t時刻x位置浸潤免疫細胞的密度; dGgt;0， dAgt;0和dNgt;0分別表示膠質瘤細胞、 趨化劑濃度和浸潤免疫細胞的擴散率; λ表示神經膠質瘤細胞的自然增長率； k為最大環(huán)境容納量； m和α分別表示神經膠質瘤細胞產生趨化劑的速率和趨化系數(shù)； β表示的Michaelis常數(shù)； γ表示趨化劑的降解速率； ρ是浸潤免疫細胞的清除率. 假設神經膠質瘤細胞的增殖遵循Logistic增長模式， 趨化劑是由神經膠質瘤細胞產生的， 用Michaelis-Menten形式mGβ+G模擬趨化劑的產生速率. 模型中參數(shù)均為正數(shù).

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

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

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

dGdt=λG1-G+Nk，

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

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

1.1 平衡解分析

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

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

直接計算可知， 當0lt;m≤m0時， 系統(tǒng)（2）無正平衡解; 當mgt;m0時， 系統(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有如下結論：

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

2） 當mlt;m0時， 系統(tǒng)（2）的平衡解E1是局部漸近穩(wěn)定的； 當m≥m0時， 系統(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）-ρ. 顯然， 當mlt;m0時， mαkγ（β+k）-ρlt;0， 系統(tǒng)（2）的平衡解E1是局部漸近穩(wěn)定的； 當mgt;m0時， mα

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

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

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

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

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

ddt=α.（3）

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

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

000，

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

-λ）m0β，1，0T， 特征值-γ對應的特征向量為V2=（0，1，0）T， 特征值0對應的的特征向量為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.

下面計算T-1F. 直接計算可得

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. 由中心流形定理可知， 存

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

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

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

假設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）零解的動力學行為可由f3（h（y3），y3）

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

1.3 正平衡點的穩(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.

當mgt;m0時， 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 假設mgt;m0， 則有下列結論：

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

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

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

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

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）會存在具有正實部的根， 因此E*在Γ2中是不穩(wěn)定的.

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

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

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

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

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

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

下面證明橫截條件α′（m*）≠0成立. 對特征方程（6）關于m求導， 并分離實部和虛部， 可得

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*）.

經計算可得

α′（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等價于

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

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

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

定理2 假設mgt;m0成立， 則：

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

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

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

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

證明： 經計算可得

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

其中

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

因為mgt;m0， 所以有λγρk2mα（αm-γρ）2gt;0. 因此當m∈（m0，+∞）時， 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等價于xgt;γρβk， 因此下面討論當x∈γρβk，+∞時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，+∞中至少存在一個根， 記為x1. 此外， 對于

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

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

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，+∞）.

因為F′（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可得結論. 證畢.

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

下面討論系統(tǒng)（2）由Hopf分支產生的周期解（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）對應的特征向量. 經計算可得

ξ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.

進一步計算可得

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）時計算以下變量［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.

由文獻［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分支周期解， 有如下結論：

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

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

2 反應擴散系統(tǒng)的動力學行為

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

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

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

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ù)文獻［13-14］可知， （E1） 的特征值可由n確定：

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

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

其中0=μ0lt;μ1lt;…μnlt;…為-Δ在區(qū)域Ω上滿足齊次Neumann邊界條件的特征值. 當m∈（0，m0）時， 對任意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）相同. 當m∈（m0，m1）時，

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.

直接計算可得

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）.

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

3 數(shù)值模擬

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

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

經計算得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)定的. 當m=5.589 1時， 系統(tǒng)（2）會產生Hopf

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

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. 當m=2∈（0，m0）時， 根據(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 當m=5.6時系統(tǒng)（2）產生的穩(wěn)定周期解

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

圖2 當m=2時系統(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 當m=4時系統(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誘導的腫瘤反應擴散模型. 首先， 通過以趨化劑產生的速率m為分支參數(shù)， 給出了該模型在正平衡點附近Hopf分支的存在性

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

因此經典的自擴散不能導致Turing斑圖的出現(xiàn)， 數(shù)值模擬也證實了該結論.

參考文獻

［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］ 楊艷紅， 伏升茂. 一類帶時滯的腫瘤免疫模型的Hopf分支 ［J］. 數(shù)學進展， 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效應的擴散時滯單種群模型的分支分析 ＼. 東北師大學報（自然科學版）， 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］ 高鶴， 李秀玲. 二維具時滯捕食＼|食餌共生模型的Hopf分支 ＼. 東北師大學報（自然科學版）， 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分支分析 ＼. 吉林大學學報（理學版）， 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.

（責任編輯： 李 琦）