張改梅, 王曉琴, 李建全
(陜西科技大學(xué)文理學(xué)院,西安 710021)
癌癥已成為嚴(yán)重危害人類健康的重大公共衛(wèi)生問題。近期世界衛(wèi)生組織發(fā)表了最新的全球癌癥統(tǒng)計(jì)數(shù)據(jù),2018年全球約有1.81×107例新的癌癥病例和9.6×106例癌癥死亡病例。由于人口的增長和人口老齡化,全球癌癥發(fā)病率和死亡率快速增長[1-2]。研究表明免疫系統(tǒng)能夠識(shí)別和消除惡性腫瘤[3-5],數(shù)學(xué)模型模擬腫瘤免疫動(dòng)力學(xué)不僅有助于理解免疫細(xì)胞和腫瘤細(xì)胞是如何相互作用的,而且也能提供一個(gè)有用的工具來預(yù)測免疫治療的結(jié)果和改善治療策略。
作為腫瘤免疫效應(yīng)階段的執(zhí)行場所,腫瘤微環(huán)境內(nèi)存在許多因素參與腫瘤組織與免疫系統(tǒng)的相互作用。由于腫瘤微環(huán)境的復(fù)雜性,不僅要考慮腫瘤細(xì)胞和效應(yīng)細(xì)胞自身的增殖和死亡、免疫效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的抑制作用,也需要考慮腫瘤對(duì)效應(yīng)細(xì)胞的作用。腫瘤的出現(xiàn)會(huì)刺激效應(yīng)細(xì)胞的增長,同時(shí),腫瘤對(duì)效應(yīng)細(xì)胞也有抑制作用,如腫瘤可通過降低微環(huán)境中葡萄糖的濃度來抑制免疫效應(yīng)細(xì)胞的活性;腫瘤微環(huán)境常呈現(xiàn)酸性,在酸性條件下,免疫效應(yīng)細(xì)胞的活性明顯降低,而且腫瘤細(xì)胞本身也參與腫瘤微環(huán)境中的免疫抑制[6]。
關(guān)于腫瘤免疫相互作用的模型已有許多學(xué)者研究[7-10]。1994年,Kuznetsov等[11]提出了二維腫瘤免疫反應(yīng)的數(shù)學(xué)模型:
(1)
Delisi等[12]、Adam[13]研究了腫瘤細(xì)胞與效應(yīng)細(xì)胞相互作用的常微分方程模型,這兩個(gè)模型均采用Michaelis-Menten型抑制函數(shù)來表示效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的抑制作用。研究結(jié)果均顯示一定范圍內(nèi),效應(yīng)細(xì)胞的增長會(huì)增大腫瘤細(xì)胞的存活率。此外,他們還給出了腫瘤生長不受免疫系統(tǒng)控制變?yōu)閻盒阅[瘤的閾值條件。Kirschner等[14]在文獻(xiàn)[11-12]的基礎(chǔ)上,假設(shè)腫瘤對(duì)效應(yīng)細(xì)胞的刺激項(xiàng)為線性項(xiàng),首次提出了腫瘤細(xì)胞、免疫效應(yīng)細(xì)胞和IL-2之間的相互作用的數(shù)學(xué)模型。研究發(fā)現(xiàn)腫瘤的抗原性對(duì)該模型的動(dòng)力學(xué)性質(zhì)有著非常大的影響,并解釋了腫瘤復(fù)發(fā)的原因。2015年,Yang等[15]在文獻(xiàn)[14]基礎(chǔ)上加入了脈沖免疫治療,建立了一個(gè)與腫瘤-免疫及脈沖免疫治療相關(guān)的數(shù)學(xué)模型,研究發(fā)現(xiàn)效應(yīng)細(xì)胞的初始密度、效應(yīng)細(xì)胞與腫瘤細(xì)胞的比例、免疫治療周期等對(duì)癌癥的治療至關(guān)重要[15]。Parthasakha等[16]提出了具有Monod-Haldane動(dòng)力學(xué)反應(yīng)的腫瘤免疫與IL-2之間相互作用的時(shí)滯數(shù)學(xué)模型,以揭示相關(guān)細(xì)胞間現(xiàn)象的動(dòng)力學(xué)機(jī)制。確定了模型解的正性、有界性和一致持久性,發(fā)現(xiàn)模型在無腫瘤平衡態(tài)下存在跨臨界分支。Duan等[17]對(duì)高斯白噪聲驅(qū)動(dòng)的腫瘤免疫系統(tǒng)的穩(wěn)定性進(jìn)行了研究。研究發(fā)現(xiàn)在這個(gè)系統(tǒng)中,有幾個(gè)穩(wěn)態(tài)。當(dāng)且僅當(dāng)噪聲較弱時(shí),一個(gè)穩(wěn)態(tài)是全局漸近穩(wěn)定的,而當(dāng)噪聲較強(qiáng)或較弱時(shí),另一個(gè)穩(wěn)態(tài)總是不穩(wěn)定的。Li等[18]提出并研究了隨機(jī)切換環(huán)境下腫瘤細(xì)胞與免疫系統(tǒng)之間的競爭模型。得到了腫瘤細(xì)胞滅絕和持續(xù)存在的理論閾值。此外,還進(jìn)行了隨機(jī)模擬來證實(shí)和解釋這些結(jié)論。
由于目前尚未有模型用Michaelis-Menten形式來表示腫瘤對(duì)效應(yīng)細(xì)胞的抑制作用,因此在假設(shè)腫瘤細(xì)胞對(duì)免疫效應(yīng)細(xì)胞具有線性刺激率和Michaelis-Menten型抑制函數(shù)的基礎(chǔ)上,建立了如式(2)所示的腫瘤細(xì)胞與效應(yīng)細(xì)胞相互作用的動(dòng)力學(xué)模型:
(2)
分析了式(2)平衡點(diǎn)的存在性、局部漸近穩(wěn)定性和全局動(dòng)力學(xué)行為,以及模型產(chǎn)生鞍結(jié)點(diǎn)分支的條件,并進(jìn)一步考慮腫瘤細(xì)胞對(duì)效應(yīng)細(xì)胞的抑制率系數(shù)和效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的抑制率系數(shù)對(duì)式(2)動(dòng)力學(xué)行為的影響。最后對(duì)所得結(jié)果進(jìn)行生物學(xué)意義的解釋和數(shù)值模擬來驗(yàn)證所得結(jié)果。
關(guān)于式(2)解的非負(fù)性和有界性有如下結(jié)論。
定理1設(shè)E(t)、T(t)為式(2)滿足非負(fù)初始條件的解, 則E(t)、T(t)是非負(fù)的且最終有界。
證明首先,易知T=0為式(2)第二個(gè)方程的解,根據(jù)解的存在唯一性可知,在非負(fù)初始條件下T(t)≥0成立。當(dāng)E(0)=0時(shí),E′(0)>0,所以存在時(shí)刻t1,當(dāng)t∈(0,t1)時(shí),有E(t)>0;當(dāng)E(0)>0時(shí),假設(shè)存在t2>0使得E(t2)=0且當(dāng)t∈(0,t2)時(shí),有E(t)>0,所以E′(t2)≤0。而由式(2)的第一個(gè)式子可得E′(t2)=s+cT(t2)>0,這與E′(t2)≤0矛盾,所以,式(2)在非負(fù)初始條件下的解保持非負(fù)性。
其次,由式(2)的第二個(gè)方程得:
(3)
s+c(K+ε1)-μE
(4)
(5)
(6)
由式(6)的第二個(gè)方程可以得到:
(7)
當(dāng)T≠K時(shí),將式(7)代入式(6)的第一個(gè)方程得:
(8)
證明由于式(2)有正不變集D,所以可通過判定方程H(T)=0在區(qū)間(0,K)上根的存在性來確定的式(2)平衡點(diǎn)的存在性。直接計(jì)算可得:
(9)
式(9)中:H′(T)為H(T)關(guān)于T求一階導(dǎo)數(shù);H″(T)為H(T)關(guān)于T求二階導(dǎo)數(shù)。
由于式(2)中各參數(shù)都為正數(shù),所以在區(qū)間(0,K)內(nèi)有H″(T)>0,即函數(shù)H(T)在區(qū)間(0,K)的圖形是上凹的。
當(dāng)方程H(T)=0在區(qū)間(0,K)內(nèi)存在一個(gè)單根即T*時(shí),其滿足H′(T*)>0。
(10)
(11)
利用變換:
(12)
將式(11)化為
(13)
式(13)中:l=cαε-pr-pμ。
(14)
將u=av2+o(v2)代入式(14),并比較兩端同次冪系數(shù)可得:
(15)
此時(shí):
(16)
進(jìn)一步,將式(16)代入式(13)的第二個(gè)方程中,得到中心流形上的解滿足:
(17)
綜上所述,對(duì)式(2)邊界平衡點(diǎn)的局部穩(wěn)定性有如下結(jié)論。
(18)
因此
(19)
由式(6)的第一個(gè)方程可得:
(20)
因此:
(21)
進(jìn)一步,由式(7)可得:
(22)
(23)
式(23)在原點(diǎn)附近可表示為
(24)
(25)
于是式(24)可寫為
(26)
(27)
根據(jù)中心流形定理[19],求得式(27)在原點(diǎn)的局部中心流形為
(28)
進(jìn)一步,將式(28)代入式(27)的第一個(gè)方程,可得:
(29)
綜上所述,對(duì)式(2)正平衡點(diǎn)的局部穩(wěn)定性有如下結(jié)論。
對(duì)式(2),記:
(30)
(31)
因此式(2)在區(qū)域D內(nèi)無閉軌線。根據(jù)式(2)平衡點(diǎn)的局部穩(wěn)定性,對(duì)于式(2)平衡點(diǎn)的全局穩(wěn)定性有如下定理。
P0為無腫瘤平衡點(diǎn)圖1 式(2)不存在正平衡點(diǎn)時(shí)的情形Fig.1 The case of formula (2) without the positive equilibrium
P*為模型的唯一有瘤平衡點(diǎn)圖2 式(2)存在唯一正平衡點(diǎn)時(shí)的情形Fig.2 The case of formula (2) with unique positive equilibrium
為模型的有瘤平衡點(diǎn)圖3 式(2)存在正平衡點(diǎn)P**時(shí)的情形Fig.3 The case of formula (2) with positive equilibrium P**
藍(lán)色線表示P*的穩(wěn)定流形圖4 式(2)存在兩個(gè)正平衡點(diǎn)時(shí)的情形Fig.4 The case of formula (2) with two positive equilibrium
為了直觀地展示理論分析的正確性,對(duì)式(2)進(jìn)行數(shù)值模擬。選取如下參數(shù)值:s=1.3×104,c=6.7×10-9,μ=0.0412,ε=1×107,r=1.5×10-1。
(2)取參數(shù)值p=5,α=1.041×10-7,K=570 000。此時(shí)R0=4.57>1,式(2)有唯一正平衡點(diǎn)P*(4.3×104,5.5×105),且該點(diǎn)在區(qū)域D內(nèi)是全局漸近穩(wěn)定的(圖2)。
將分析腫瘤細(xì)胞對(duì)效應(yīng)細(xì)胞的抑制率系數(shù)p和效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的抑制率系數(shù)α對(duì)模型動(dòng)力學(xué)行為的影響。
首先,當(dāng)s=1.3×104,c=6.7×10-9,μ=0.041 2,ε=1×107,r=1.5×10-1,α=1.041×10-7,K=570 000時(shí),取p分別為0、4、8、12、16。由圖5可以看到,當(dāng)p=0時(shí),腫瘤細(xì)胞的數(shù)量處于最低水平,而效應(yīng)細(xì)胞的數(shù)量處于最高水平。當(dāng)p增大時(shí),腫瘤細(xì)胞的數(shù)量會(huì)隨著腫瘤對(duì)效應(yīng)細(xì)胞的抑制率系數(shù)p的增大而增加,而效應(yīng)細(xì)胞的數(shù)量會(huì)隨著p的增大而減少。其次,當(dāng)s=1.3×104,c=6.7×10-9,μ=0.041 2,ε=1×107,r=1.5×10-1,p=5,K=570 000時(shí),由圖6可以看到,隨著效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的抑制率系數(shù)α的增大,腫瘤細(xì)胞的數(shù)量逐漸減少,而效應(yīng)細(xì)胞的數(shù)量逐漸增加,當(dāng)α=15×10-7時(shí),腫瘤細(xì)胞的數(shù)量減少到0,也即此時(shí)腫瘤會(huì)滅絕。
圖5 腫瘤細(xì)胞對(duì)效應(yīng)細(xì)胞的抑制率系數(shù)p對(duì)模型動(dòng)力學(xué)行為的影響Fig.5 Effects of the inhibition coefficient p of tumor cells to effector cells on the dynamics of the model
圖6 效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的抑制率系數(shù)α對(duì)模型動(dòng)力學(xué)行為的影響Fig.6 Effects of the inhibition coefficient α of effector cells to tumor cells on the dynamics of the model
通過觀察腫瘤細(xì)胞對(duì)效應(yīng)細(xì)胞抑制率系數(shù)與效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞抑制率系數(shù)的變化對(duì)模型的影響(圖5,圖6)可得到如下結(jié)論。
(1)可以通過一些方法來降低腫瘤細(xì)胞對(duì)效應(yīng)細(xì)胞的抑制率系數(shù),逆轉(zhuǎn)免疫抑制狀態(tài),提高免疫效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的殺傷力。
(2)增大效應(yīng)細(xì)胞對(duì)腫瘤細(xì)胞的抑制率對(duì)腫瘤的清除有良好的效果。