具時滯的腫瘤免疫擴散模型的動力學性質分析

劉高楊, 丁宇婷

(東北林業(yè)大學 理學院, 黑龍江 哈爾濱 150040)

一直以來,腫瘤都是嚴重危害人類健康的重大疾病之一.由于環(huán)境污染的日益嚴重和人口老齡化問題,腫瘤的發(fā)病率越來越高.盡管人類已經(jīng)嘗試多種辦法治療腫瘤,但是如今腫瘤依然是人類難以攻克的難題.相關研究表明,免疫力低下是形成腫瘤的原因之一.隨著分子細胞生物學理論的日益成熟,利用人體的免疫力治療腫瘤成為醫(yī)療工作者的首選方法.

1999年,Sheu等提出了利用免疫逃脫機制治療腫瘤的觀點[1].2000年,Weiner等闡述了放射免疫療法治療腫瘤的原理[2].為尋找提高機體免疫力的方法,Coulon等在實驗中證明樹突狀細胞(dendritic cells,DC)能刺激T細胞的增殖,進而增加機體免疫力.通過刺激輔助性T細胞增殖,達到治療腫瘤的目的,提出利用繼承性細胞免疫療法(adoptive cellular immunotherapy,ACI)刺激免疫系統(tǒng)清除腫瘤細胞[3].2021年,Zhou等指出腫瘤相關巨噬細胞(tumour-associated macrophages,TAM)在治療腫瘤中起著重要的作用[4].如今,“免疫治療法”的相關理論正在日益完善.

近些年來,人們利用數(shù)學模型這一工具研究提高機體免疫力的方法[5-10],并將Hopf分支理論巧妙地應用其中,取得了豐碩的成果[11-17].免疫細胞能夠識別和清除腫瘤細胞,但是由于免疫細胞識別腫瘤細胞具有一定的時間延遲,因此,研究免疫細胞識別腫瘤細胞的時滯對于分析腫瘤細胞與免疫細胞之間的動力學行為是非常重要的.時滯變化能夠改變系統(tǒng)平衡點的穩(wěn)定性,隨著時滯的變化,系統(tǒng)可以展現(xiàn)不同的動力學行為.因此,探索時滯引起機體免疫力變化的規(guī)律成為治療腫瘤的關鍵之一[18-20].

近年來,很多學者構建了各類的時滯微分方程,Kaddar等提出了如下模型用以解釋腫瘤細胞的周期解問題[21]:

式中:E是效應T細胞的數(shù)量;T是腫瘤細胞的數(shù)量;σ是效應T細胞分泌的抗體進入患病部位的速率;ω是效應T細胞的應答速率;δ是效應T細胞的死亡速率;α是腫瘤細胞的最大成長速率;β是環(huán)境最大容納量;τ是免疫系統(tǒng)產(chǎn)生治療腫瘤的化學物質所需的時間;除了ω之外,其他系數(shù)都是正的.

在實際問題中,由于細胞表面糖蛋白等物質的減少,導致腫瘤細胞會發(fā)生擴散,因此,在系統(tǒng)中引入擴散項是必要的.另外,考慮到免疫細胞識別腫瘤細胞存在時間延遲,在模型中引入免疫細胞識別腫瘤細胞的時間延遲τ能夠更準確地刻畫系統(tǒng)的動力學行為.由于腫瘤細胞和免疫細胞只會分布在研究機體范圍內,為此,研究如下具有紐曼邊界條件的腫瘤免疫擴散系統(tǒng):

(1)

式中:u和v分別代表腫瘤細胞和免疫細胞的數(shù)量;d1是腫瘤細胞的擴散速率;d2是免疫細胞的擴散速率;r是腫瘤細胞自身的成長速率;α是腫瘤細胞的最大環(huán)境容納量;m是被免疫細胞殺死的腫瘤細胞的死亡速率;β是免疫細胞的增長速率;d表示免疫細胞的死亡速率;p表示被腫瘤細胞殺死的免疫細胞的死亡速率;τ是免疫細胞識別腫瘤細胞所需的時間;所有的系數(shù)都是正的.

1 系統(tǒng)平衡點的穩(wěn)定性和Hopf分支的存在性

系統(tǒng)(1)的平衡點滿足如下方程:

可得系統(tǒng)(1)的邊界平衡點為E0=(0,0),E1=(α,0).假設

(H1):d<α(β-p).

當系統(tǒng)參數(shù)滿足假設(H1)時,系統(tǒng)存在非平凡平衡點E2=(u*,v*),其中,

(2)

接下來分析系統(tǒng)(1)的平衡點的穩(wěn)定性.

在假設(H1)下,令U(x,t)=(u(x,t),v(x,t))T,系統(tǒng)(1)在平衡點(u0,v0)(這里(u0,v0)=(0,0),(α,0)或(u*,v*))處可寫成如下形式:

式中:

因此,系統(tǒng)(1)在平衡點(u0,v0)處的特征方程為

λ2+Anλ+Bn+C1e-λτ=0,n∈,

(3)

定理1 對于任意τ≥0,系統(tǒng)(1)在平衡點E0處經(jīng)歷B -T分支.

證明 根據(jù)方程(3)可得,在平衡點E0(即u0=v0=0)處的特征方程為

λ2+(d1+d2)n2λ+d1d2n4=0,n∈.

(4)

令λ1和λ2是方程(4)的根,則λ1+λ2=-(d1+d2)n2,λ1λ2=d1d2n4.

當n=0時,λ1=λ2=0,因此,系統(tǒng)(1)在平衡點E0處經(jīng)歷B -T分支.

定理2 對于任意τ≥0,系統(tǒng)(1)在平衡點E1處經(jīng)歷不動點分支.

證明 根據(jù)方程(3)可得,在平衡點E1(即u0=α,v0=0)處的特征方程為

λ2+(d1n2+d2n2+r)λ+d2n2(d1n2+r)=0,n∈.

(5)

因此,

λ1+λ2=-(d1n2+d2n2+r),λ1λ2=d2n2(d1n2+r).

當n=0時,λ1+λ2<0,λ1λ2=0,所以,當n=0時,特征方程(5)有一個單零根,因此,系統(tǒng)(1)在平衡點E1處經(jīng)歷不動點分支.

下面分析系統(tǒng)(1)在平衡點E2處的穩(wěn)定性.

定理3 當τ=0時,在假設(H1)下,平衡點E2=(u*,v*)是局部漸近穩(wěn)定的.

證明 根據(jù)方程(3)可得,當τ≠0時,系統(tǒng)(1)在平衡點E2=(u*,v*)(即u0=u*,v0=v*)處的特征方程為

λ2+Anλ+Bn+C2e-λτ=0,n∈,

(6)

當τ=0時,系統(tǒng)(1)在平衡點E2=(u*,v*)處的特征方程為

(7)

令λ1和λ2是方程(7)的根,于是,

因此,λ1和λ2具有嚴格的負實部.所以,當τ=0時,平衡點E2是局部漸近穩(wěn)定的.

為了方便說明,引入如下記號:

當τ≠0時,設λ=iω(ω>0)是系統(tǒng)(1)在平衡點E2處特征方程(6)的根,分離實部和虛部,得

(8)

兩邊平方相加,得

(9)

由式(8)可得

(10)

令z=ω2,可得

(11)

當0≤n

(12)

由方程(10)可得

(13)

(14)

因此,當τ=τ*時,系統(tǒng)(1)在平衡點E2處經(jīng)歷Hopf分支.結合上面分析,可得定理4.

2 Hopf分支的規(guī)范型及分支周期解的穩(wěn)定性

在本節(jié),利用多尺度方法推導系統(tǒng)(1)在平衡點E2=(u*,v*)處的Hopf分支的規(guī)范型.

(15)

下面用多尺度方法推導系統(tǒng)(15)的Hopf分支規(guī)范型.將式(15)的解設為如下形式:

(17)

式中,

此時,對t的導數(shù)可化為

令

式中,j∈.

由式(17)可以得

(18)

將時滯部分進行泰勒展開,可得:

(19)

式中,uj,1=uj(T0-1,T1,T2,…),vj,1=vj(T0-1,T1,T2,…),j∈.

將式(17)～式(19)代入到式(15),對應ε前的系數(shù),可得:

(20)

因此,式(20)的解為

(21)

式中,c.c.代表前面項的復共軛.

對應ε2前的系數(shù),可得如下等式:

(22)

(23)

(24)

式中,

假設

(25)

(26)

(27)

對應ε3前的系數(shù),可得到如下表達式:

(28)

(29)

式中,

(30)

式中,

因此,系統(tǒng)(15)(或系統(tǒng)(1))截斷到三階的Hopf分支規(guī)范型為

(31)

式中,M和χ分別由式(24)和式(30)給出.

1) 當Re(M)μ<0時,約化在中心流形上的分支周期解不穩(wěn)定;

2) 當Re(M)μ>0時,約化在中心流形上的分支周期解穩(wěn)定.

3 數(shù)值模擬

根據(jù)文獻[22]所述,取如下參數(shù):r=0.18 d,α=5.0×106個活細胞,d2=0.2 d,d=0.015 2 d,p=3.422×10-10個活細胞/d,d1=0.02 d,m=1.101×10-7個活細胞/d,β=6.2×10-9個活細胞/d.由式(2)可得,系統(tǒng)平衡點E2=(2.594 8×106,7.864 3×105),由式(12)得

由定理1可知,任意τ≥0,系統(tǒng)(1)在平衡點E0處經(jīng)歷B -T分支.以(0.000 000 001,0.000 000 1)為初值,取τ=0,平衡點E0是局部漸近穩(wěn)定的,如圖1所示;取τ=0.1,平衡點E0是不穩(wěn)定的,如圖2所示.這說明在B-T分支附近,當參數(shù)在不同區(qū)域時,系統(tǒng)(1)可產(chǎn)生不同的動力學現(xiàn)象.

(a) u(x,t)(b) v(x,t)

由定理2可知,任意τ≥0,系統(tǒng)(1)在平衡點E1處經(jīng)歷不動點分支.以(4.99×106,0.01)為初值,取τ=0,平衡點E1是不穩(wěn)定的,如圖3所示;取τ=0.009,平衡點E1是不穩(wěn)定的,如圖4所示.比較圖3和圖4,雖然平衡點E1都是不穩(wěn)定的,但是第一個分量從穩(wěn)定變成不穩(wěn)定,說明腫瘤細胞發(fā)生了擴散,證明了即使免疫細胞識別腫瘤細胞的時滯發(fā)生微小的變化,卻能影響腫瘤細胞能否發(fā)生擴散.

(a) u(x,t)(b) v(x,t)圖2 當τ=0.1時,系統(tǒng)(1)的平衡點E0是不穩(wěn)定的Fig.2 The equilibrium E0 of system(1) is unstable when τ=0.1

(a) u(x,t)(b) v(x,t)圖3 當τ=0時,系統(tǒng)(1)的平衡點E1是不穩(wěn)定的Fig.3 The equilibrium E1 of system(1) is unstable when τ=0

(a) u(x,t)(b) v(x,t)圖4 當τ=0.009時,系統(tǒng)(1)的平衡點E1是不穩(wěn)定的Fig.4 The equilibrium E1 of system(1) is unstable when τ=0.009

由定理3和定理4可知,當τ∈[0,τ*)時,平衡點E2是局部漸近穩(wěn)定的;當τ>τ*時,平衡點E2是不穩(wěn)定的.

當τ=0時,以(2.5 947×106,7.8 642×105)為初值,平衡點E2是局部漸近穩(wěn)定的,如圖5所示.

(a) u(x,t)(b) v(x,t)圖5 當τ=0時,系統(tǒng)(1)的平衡點E2是局部漸近穩(wěn)定的Fig.5 The equilibrium E2 of system(1) is locally asymptotically stable when τ=0

當τ=3<τ*時,以(2.594 7×106,7.864 2×105)為初值,平衡點E2是局部漸近穩(wěn)定的,如圖6所示.

(a) u(x,t)(b) v(x,t)圖6 當τ=3<τ*時,系統(tǒng)(1)的平衡點E2是局部漸近穩(wěn)定的Fig.6 The equilibrium E2 of system (1) is locally asymptotically stable when τ=3<τ*

取上述參數(shù),當τc=τ*=92.331 2時,由式(24)和式(30),經(jīng)過MATLAB計算可得,Re(M)<0,Re(χ)>0,根據(jù)定理5,可得系統(tǒng)(15)(或系統(tǒng)(1))在μ>0時具有不穩(wěn)定的周期解,并且Hopf分支方向是向前的(μ>0).

從生物學的意義上講,當免疫細胞識別腫瘤細胞需要時間較長時,免疫系統(tǒng)的免疫能力極低,導致腫瘤細胞生長活躍,容易癌變,有向身體其他部位擴散和轉移的危險.

4 結 語

本文基于KADDAR和ALAOUI的模型,在模型中引入免疫細胞識別腫瘤細胞的時間延遲τ和擴散項,改進得到了一個偏泛函微分方程(partial function differential equation,PFDE),討論了該系統(tǒng)的動力學性質.由于腫瘤細胞會發(fā)生擴散和免疫細胞識別腫瘤細胞存在時間延遲這兩個現(xiàn)象都是影響腫瘤免疫系統(tǒng)動力學性質的重要因素,因此,相對于常微分方程(ordinary differential equation,ODE)和泛函微分方程(function differential equation,FDE),研究結果對于腫瘤的控制和治療更具有指導意義.研究了時滯對于腫瘤免疫擴散系統(tǒng)平衡點的穩(wěn)定性以及Hopf分支周期解穩(wěn)定性的影響.選擇時滯作為分支參數(shù),建立了系統(tǒng)非負平衡點穩(wěn)定的條件.時滯可以影響平衡點的穩(wěn)定性,當時滯在一定范圍內,系統(tǒng)非平凡平衡點是局部漸近穩(wěn)定的,一旦時滯穿過臨界值,系統(tǒng)平衡點失去了穩(wěn)定性,并且系統(tǒng)產(chǎn)生了不穩(wěn)定的周期解.當時滯取臨界值時,系統(tǒng)經(jīng)歷了Hopf分支.這意味著,當免疫細胞能較快的識別腫瘤細胞時,人體具有較強的免疫力,能抑制腫瘤細胞的增長和擴散;反之,腫瘤細胞迅速增殖,大量聚積在人體,對人體健康產(chǎn)生很大的危害.系統(tǒng)臨界值τ*的發(fā)現(xiàn),有助于生物學家對腫瘤細胞的成長進行預測、判斷和控制,進而采取有效的治療手段,最大限度的發(fā)揮細胞免疫療法在治療腫瘤中的作用.該研究只考慮了免疫細胞識別腫瘤細胞的時間延遲對治療腫瘤的影響,在以后的工作中,將致力于研究其它因素對控制腫瘤細胞擴散的影響,找到抑制腫瘤細胞惡性生長的條件.