一類兩階段腫瘤-免疫模型的確定性與隨機(jī)分析*

李孟陽 李偉 黃冬梅 楊貴東

（西安電子科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，西安 710071）

引言

癌癥作為一種惡性腫瘤，正嚴(yán)重威脅人類的健康.根據(jù)國(guó)際癌癥研究機(jī)構(gòu)（IARC）發(fā)布的《全球癌癥研究報(bào)告》，2018年全球新增癌癥病例近1800萬，且致死率正在逐年增加[1].眾所周知，盡管臨床上對(duì)癌癥可以進(jìn)行手術(shù)、放療、化療等治療，但不可預(yù)期的復(fù)發(fā)和再次擴(kuò)散，使得攻克癌癥依然是目前醫(yī)學(xué)上最大的難題之一.而生物實(shí)驗(yàn)需要大量的人力、物力、實(shí)驗(yàn)設(shè)備和實(shí)驗(yàn)技術(shù)，各種局限性使得生物數(shù)學(xué)模型以及機(jī)理分析極其重要，成為不可或缺的用來理解腫瘤發(fā)生、演化和擴(kuò)散的重要工具[2].

免疫系統(tǒng)是人體自我保護(hù)的重要組織之一，具有識(shí)別和消滅惡性腫瘤的功能[3].一般來說，如果腫瘤細(xì)胞在某一組織中產(chǎn)生，人體對(duì)腫瘤的免疫反應(yīng)就會(huì)被激活，同時(shí)，效應(yīng)細(xì)胞等多種免疫細(xì)胞會(huì)刺激淋巴細(xì)胞在腫瘤表面進(jìn)行分裂、協(xié)調(diào)，并放置抗體來防御腫瘤.

目前關(guān)于腫瘤-免疫系統(tǒng)相互作用的模型[4-10]主要有兩類.一類是根據(jù)腫瘤細(xì)胞的生化反應(yīng)特征，建立確定性動(dòng)力學(xué)模型，集中討論腫瘤免疫系統(tǒng)解的存在性、穩(wěn)定性、參數(shù)估計(jì)和分岔等[11-16].另一類是考慮細(xì)胞生長(zhǎng)的內(nèi)外微環(huán)境，將隨機(jī)因素引起的隨機(jī)漲落模型化為噪聲，建立隨機(jī)動(dòng)力學(xué)模型，討論隨機(jī)因素對(duì)腫瘤發(fā)生、發(fā)展以及消亡的影響[17-23].

不同類型的模型作用不同，對(duì)理解腫瘤與免疫系統(tǒng)之間的競(jìng)爭(zhēng)都各有幫助.相對(duì)確定性模型而言，隨機(jī)腫瘤-免疫模型的動(dòng)力學(xué)討論文獻(xiàn)積累比較少，變量主要局限于一維或二維情形.這是因?yàn)殡S機(jī)項(xiàng)和非線性項(xiàng)的引入，增加了隨機(jī)腫瘤-免疫模型的理論與數(shù)值求解的難度，特別是噪聲的存在，其分布類型更是直接關(guān)系到隨機(jī)腫瘤-免疫模型是否能夠求解.然而，在腫瘤與免疫系統(tǒng)的斗爭(zhēng)中，對(duì)腫瘤具有抑制作用的效應(yīng)細(xì)胞通常有多種，如巨噬細(xì)胞、細(xì)胞毒性T細(xì)胞、NK細(xì)胞等，它們與腫瘤細(xì)胞之間的相互作用也極其復(fù)雜，三維及更高維的隨機(jī)腫瘤-免疫模型才更符合實(shí)際的腫瘤-免疫關(guān)系.

基于淋巴細(xì)胞在免疫應(yīng)答中的生化特點(diǎn)以及腫瘤免疫實(shí)驗(yàn)結(jié)果，DeLisi等[26]提出：只有發(fā)育成熟的淋巴細(xì)胞才能有效殺死腫瘤細(xì)胞，免疫系統(tǒng)對(duì)腫瘤的消除應(yīng)該對(duì)淋巴細(xì)胞按成熟與未成熟兩階段分別進(jìn)行考慮.之后，Liu等[27]對(duì)這一見解給予了肯定，并在文獻(xiàn)[26]的基礎(chǔ)上對(duì)兩階段確定性腫瘤-免疫模型的退化平衡點(diǎn)進(jìn)行了定性分岔分析.再后，龐留勇等[28]對(duì)該模型進(jìn)行了簡(jiǎn)化，重點(diǎn)從理論上研究了Hopf分岔的條件公式.蒙玉倩等[29]進(jìn)一步將該模型進(jìn)行了改進(jìn)，考慮了在分?jǐn)?shù)階情形下該確定性腫瘤-免疫模型的穩(wěn)定性，分析了兩階段淋巴細(xì)胞的動(dòng)態(tài)演變情況.

本文將基于文獻(xiàn)[26-29]的工作，將淋巴細(xì)胞抵御腫瘤的過程分為兩個(gè)階段，即被激活的未成熟淋巴細(xì)胞的生長(zhǎng)階段和能有效消除腫瘤細(xì)胞的成熟淋巴細(xì)胞的防御階段.同時(shí)，本文將在這些確定性模型基礎(chǔ)上繼續(xù)深入探討該模型的相關(guān)性質(zhì)，并將細(xì)胞微環(huán)境產(chǎn)生的隨機(jī)漲落也考慮在內(nèi)，在理論上研究腫瘤細(xì)胞、未成熟淋巴細(xì)胞以及成熟淋巴細(xì)胞之間的三維隨機(jī)模型，探討隨機(jī)漲落對(duì)腫瘤生成以及演變過程的影響.

1 確定性腫瘤-免疫模型

基于上述分析，下面利用Runge-Kutta方法對(duì)確定性模型（2）進(jìn)行數(shù)值模擬來驗(yàn)證結(jié)論的正確性.假設(shè)腫瘤已經(jīng)生成，淋巴細(xì)胞已經(jīng)被激活并進(jìn)行防御，可設(shè)初始點(diǎn)為 P（1，1，1），若參數(shù)值分別為 a=0.9，b=0.6，c=0.6，k=0.8 ，則系統(tǒng)響應(yīng)將漸進(jìn)趨近于一個(gè)穩(wěn)定的無腫瘤平衡點(diǎn)E1（0，1，0），如圖1所示.說明即使腫瘤細(xì)胞的增殖率達(dá)到0.8，免疫系統(tǒng)的60%殺死率也足以清除所有的腫瘤細(xì)胞，最終病人將被治愈，成為無腫瘤的狀態(tài).未成熟的淋巴細(xì)胞在發(fā)育成熟后，不再被激活.

圖 1 E1的相圖，取 a=0.9，b=0.6，c=0.6，k=0.8Fig.1 Phase diagram of E1 for a=0.9，b=0.6，c=0.6，k=0.8

對(duì)于第二個(gè)平衡點(diǎn)E2，與E1的分析過程類似，當(dāng)參數(shù)滿足條件：

圖2給出了條件（8）下的E2穩(wěn)定性區(qū)域劃分.注意到現(xiàn)實(shí)中任何時(shí)候并不期望腫瘤出現(xiàn)或長(zhǎng)期存在，也就意味著E2的穩(wěn)定性對(duì)人類健康沒有好處，它的不穩(wěn)定才會(huì)帶來治療與控制的機(jī)會(huì).因此，圖3考查了初始點(diǎn)為 P（1，1，1）時(shí)不同殺死率下的三類細(xì)胞演化情況.

圖2 E2的穩(wěn)定區(qū)域，取c=0.2，k=0.6Fig.2 The stable region of E2 for c=0.2，k=0.6

顯然，若保持成熟淋巴細(xì)胞和腫瘤的增殖率不變，即a=0.9，b=0.6，k=0.8，隨著殺死率的逐漸減小，系統(tǒng)（2）從（a）和（b）中的穩(wěn)定狀態(tài)轉(zhuǎn)變?yōu)椋╟）中的不穩(wěn)定狀態(tài).此外，（a）中的穩(wěn)定點(diǎn)逐漸變成一個(gè)不穩(wěn)定的極限環(huán).這說明免疫系統(tǒng)對(duì)腫瘤的抵抗能力較弱，不能殺死全部的腫瘤細(xì)胞，免疫細(xì)胞與腫瘤細(xì)胞一直處于此起彼伏的較量之中.圖3中z軸為腫瘤細(xì)胞密度.可以看到，隨著參數(shù)c的減小，腫瘤細(xì)胞密度不斷增加，說明腫瘤細(xì)胞一直存在于人體，并正在發(fā)展和惡化.

圖 3 E2的相圖，取 a=0.9，b=0.6，k=0.8Fig.3 Phase diagram of E2 for a=0.9，b=0.6，k=0.8

圖 4 E1的一階矩，取 a=0.9，b=0.6，c=0.6，k=0.8Fig.4 The first-order moments of E1 for a=0.9，b=0.6，c=0.6，k=0.8

2 隨機(jī)腫瘤-免疫模型

細(xì)胞生成的內(nèi)外環(huán)境極其復(fù)雜，腫瘤細(xì)胞和免疫細(xì)胞難免也會(huì)受到環(huán)境隨機(jī)因素的影響，如外部的化學(xué)藥劑、溫度、供氧量和輻射線等，內(nèi)部的溫度、壓力、細(xì)胞大小以及細(xì)胞濃度等.因此，構(gòu)建腫瘤-免疫系統(tǒng)模型時(shí)不可避免地需要考慮隨機(jī)噪聲[31].按照噪聲來源，細(xì)胞內(nèi)微環(huán)境形成的隨機(jī)漲落為加性噪聲，而外部因素形成的噪聲為乘性噪聲.

本節(jié)在確定性腫瘤-免疫系統(tǒng)的基礎(chǔ)上，將微環(huán)境產(chǎn)生的隨機(jī)擾動(dòng)模型化為相互獨(dú)立的高斯白噪聲，考慮如下隨機(jī)情形時(shí)的三維腫瘤-免疫系統(tǒng)：

其中，高斯白噪聲 ξi（t）（i=0，1，2） 滿足均值函數(shù)為 〈ξi（t）〉 =0 ， 相關(guān)函數(shù)為 〈ξi（t）ξi（t′）〉 =2Diiδ（t-t′） ，這里 Dii為 ξi（t） 的噪聲強(qiáng)度.

顯然，隨機(jī)噪聲的存在，使得三類細(xì)胞密度在時(shí)間的演化過程中也具有隨機(jī)性，只能通過概率統(tǒng)計(jì)的方法確定它們的性質(zhì)和規(guī)律.另外，由于隨機(jī)性的存在，本文將主要關(guān)注確定性系統(tǒng)中的兩個(gè)平衡態(tài)，觀察平衡態(tài)附近細(xì)胞密度的隨機(jī)漲落.

2.1 統(tǒng)計(jì)分析

這是一個(gè)三維的確定性偏微分方程，描述了三類細(xì)胞在空間中的聯(lián)合概率分布情況.無論是瞬態(tài)解還是穩(wěn)態(tài)解，方程（11）的求解一般來說都是相當(dāng)困難的.下面將主要討論細(xì)胞變量的統(tǒng)計(jì)矩.

文獻(xiàn)[32]已經(jīng)證明，系統(tǒng)變量xi及其導(dǎo)數(shù)的統(tǒng)計(jì)矩可以分別利用FP方程計(jì)算得出，即

均值函數(shù)以 P（1，1，1）為初始值的解曲線如圖4所示.圖中可以看出，未成熟免疫細(xì)胞與腫瘤細(xì)胞在有腫瘤的初始條件下隨時(shí)間單調(diào)遞減到0，但成熟免疫細(xì)胞先增后減并隨時(shí)間逐漸減小到0.顯然在拐點(diǎn)前，未成熟的淋巴細(xì)胞逐漸成熟，加入抵御腫瘤細(xì)胞的大軍，成熟細(xì)胞數(shù)量增加是必然的.然而在拐點(diǎn)后，單調(diào)的減少與確定性情形截然不同，沒有漸近趨于數(shù)值1.說明噪聲即微環(huán)境改變了成熟淋巴細(xì)胞原有的演變規(guī)律，說明此處微環(huán)境有利于腫瘤的生存和生長(zhǎng).成熟的免疫細(xì)胞在微環(huán)境作用下逐漸死亡并退出防御系統(tǒng).

考慮細(xì)胞變量的二階矩，取n=2，則得到一個(gè)關(guān)于細(xì)胞密度方差和協(xié)方差的6個(gè)常微分方程：

求解該方程組，則細(xì)胞密度 xi（i=0，1，2） 的方差和協(xié)方差如圖5所示.

圖5 不同噪聲強(qiáng)度D22時(shí)E1的二階矩.a=0.9，b=0.6，c=0.6，k=0.8，D00=0.001，D11=0.001Fig.5 The second-order moments of E1 with different noise intensities D22.a=0.9，b=0.6，c=0.6，k=0.8，D00=0.001，D11=0.001

圖5給出了三種細(xì)胞的六個(gè)二階矩在腫瘤細(xì)胞微環(huán)境噪聲強(qiáng)度下的時(shí)間演變情況.可以看出，當(dāng)噪聲強(qiáng)度較小時(shí)，比較而言，成熟的淋巴細(xì)胞前期波動(dòng)較大，但隨著腫瘤-免疫系統(tǒng)的不斷演化，三類細(xì)胞的數(shù)量波動(dòng)逐漸減緩，趨于平緩到一個(gè)較小的值.然而，當(dāng)腫瘤細(xì)胞微環(huán)境的噪聲強(qiáng)度較大時(shí)，三類細(xì)胞波動(dòng)不斷持續(xù)增大，腫瘤與免疫細(xì)胞之間的斗爭(zhēng)更加激烈，并始終保持較大的數(shù)值.說明微環(huán)境此時(shí)嚴(yán)重影響細(xì)胞間的演變，腫瘤與免疫系統(tǒng)間的競(jìng)爭(zhēng)此起彼伏，一直處于膠著的狀態(tài).

均值函數(shù)解曲線如圖6所示.參數(shù)值相同時(shí)，三種細(xì)胞的均值函數(shù)前期呈現(xiàn)出較大的振蕩，隨著時(shí)間的增加最終趨于0.其中，腫瘤細(xì)胞受微環(huán)境影響產(chǎn)生的振幅較大，它更易受到微環(huán)境的干擾.

圖 6 E2的一階矩，取 a=0.9，b=0.6，c=0.2，k=0.8Fig.6 The first-order moments of E2 for a=0.9，b=0.6，c=0.2，k=0.8

將方程（21）代入式（13）得到變量 yi（i=0，1，2）的6個(gè)二階矩的微分方程組：

圖7給出了E2附近的二階矩隨噪聲強(qiáng)度的變化.這里出現(xiàn)一個(gè)非常有趣的結(jié)果，即腫瘤細(xì)胞無論噪聲強(qiáng)度大小，始終呈現(xiàn)出很大的波動(dòng)，而免疫細(xì)胞幾乎沒有波動(dòng).這種情況說明，在腫瘤一直存在的狀態(tài)下，細(xì)胞與免疫細(xì)胞長(zhǎng)期處于競(jìng)爭(zhēng)狀態(tài)，微環(huán)境對(duì)細(xì)胞間的演變影響較小.

圖 7 E2的二階矩，取 a=0.9，b=0.6，c=0.08，k=0.8，D00=0.001，D11=0.001Fig.7 The second-order moments of E2 for a=0.9，b=0.6，c=0.08，k=0.8，D00=0.001，D11=0.001

注意到當(dāng) bk-c=0時(shí)，平衡點(diǎn) E2（x2*，y2*，z2*）的x2*與z2*均變?yōu)?.此時(shí)，未成熟的免疫細(xì)胞與腫瘤細(xì)胞都不存在，平衡點(diǎn)E2可轉(zhuǎn)化為E1.因此，接下來主要研究有腫瘤平衡點(diǎn)E2的有關(guān)性質(zhì).

2.2 敏感性分析

顯然，免疫系統(tǒng)對(duì)腫瘤細(xì)胞的殺死率、腫瘤細(xì)胞的自我增殖率來說都是十分重要的參數(shù)，它們的取值將直接影響到腫瘤是否能夠得到控制.因此這部分將引入敏感系數(shù)的概念，來討論腫瘤-免疫系統(tǒng)中各類細(xì)胞的穩(wěn)態(tài)值對(duì)參數(shù)取值的敏感依賴性.敏感系數(shù)定義如下：

其中，Mi（i=0，1，2） 表示腫瘤 -免疫系統(tǒng)中第 i種細(xì)胞的穩(wěn)態(tài)值，p表示系統(tǒng)中某參數(shù).

在腫瘤與免疫系統(tǒng)的競(jìng)爭(zhēng)關(guān)系中，成熟淋巴細(xì)胞的殺死率大小是一個(gè)至關(guān)重要的參數(shù)，決定了免疫系統(tǒng)是否能夠有效清除腫瘤細(xì)胞，并控制腫瘤的發(fā)展和擴(kuò)散.因此，這里主要考查殺死率參數(shù)c對(duì)第二個(gè)穩(wěn)態(tài)值E2的影響，即M0=x*2，M1=y*2，M2=z*2.由公式（29）有

因此，圖8給出了第二個(gè)穩(wěn)態(tài)值E2附近三種細(xì)胞表型對(duì)殺死率的敏感性曲線.可以看出，當(dāng)參數(shù)a，b，k固定時(shí)，隨著c的增加，第一階段免疫細(xì)胞和腫瘤細(xì)胞都敏感依賴于參數(shù)c，并呈非線性增加.第二階段免疫細(xì)胞的敏感系數(shù)為0，即對(duì)成熟淋巴細(xì)胞的密度沒有影響.觀察敏感系數(shù)曲線的發(fā)展趨勢(shì)，不難判斷，當(dāng)殺死率進(jìn)一步增大時(shí)，腫瘤細(xì)胞密度和未成熟淋巴細(xì)胞的密度將強(qiáng)烈地依賴參數(shù)c，因?yàn)樾枰罅康奈闯墒炝馨图?xì)胞轉(zhuǎn)化為具有殺傷力的成熟淋巴細(xì)胞，并保持成熟淋巴細(xì)胞的恒定，去抵抗腫瘤細(xì)胞.這個(gè)過程中，免疫細(xì)胞與腫瘤細(xì)胞在數(shù)量上變化較大，故依賴于參數(shù)c.

圖8 三種細(xì)胞表型對(duì)殺死率c的敏感性依賴.a=0.9，b=0.6，k=0.78Fig.8 The sensitivity of three cell phenotypes on the parameter of killing rate c.a=0.9，b=0.6，k=0.78

3 結(jié)論

本文研究了高斯白噪聲作用下的兩階段腫瘤-免疫細(xì)胞模型.對(duì)確定性模型中穩(wěn)態(tài)解的動(dòng)力學(xué)性質(zhì)進(jìn)行了分析和仿真.考慮到細(xì)胞微環(huán)境對(duì)系統(tǒng)的影響，將高斯白噪聲引入模型，基于朗之萬理論和FP方程給出了隨機(jī)腫瘤免疫模型細(xì)胞的一階矩和二階矩方程，并討論了細(xì)胞的方差和敏感系數(shù)等性質(zhì).研究表明噪聲會(huì)改變免疫系統(tǒng)對(duì)腫瘤細(xì)胞的防御規(guī)律.當(dāng)噪聲強(qiáng)度較大時(shí)，腫瘤細(xì)胞更易受到微環(huán)境影響，其密度在短時(shí)間內(nèi)發(fā)生劇烈變化，腫瘤與免疫細(xì)胞之間的斗爭(zhēng)更加激烈.敏感系數(shù)表明腫瘤細(xì)胞對(duì)免疫系統(tǒng)的殺死率參數(shù)c很敏感.最后通過數(shù)值方法，驗(yàn)證了分析結(jié)果的準(zhǔn)確性.