Lévy噪聲和高斯白噪聲共同激勵(lì)的腫瘤細(xì)胞生長模型的平均首次穿越時(shí)間

樊順厚，姚 婷，郭永峰，王琳杰

（天津工業(yè)大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，天津 300387）

癌癥一直是現(xiàn)代醫(yī)學(xué)所關(guān)注的重點(diǎn)問題，科學(xué)家也尋找到了一些治療方法，比如手術(shù)、化療和放射性療法等，但是這些療法也伴有副作用高且治愈率低等缺點(diǎn)。免疫療法是近幾年新提出的一種治療癌癥的方法，利用人體內(nèi)的免疫機(jī)制來對(duì)抗腫瘤細(xì)胞。因此，有必要研究免疫監(jiān)視下的腫瘤生長，以開發(fā)針對(duì)癌癥治療的有效策略。Logistic 系統(tǒng)經(jīng)常被用來描述單一種群的繁殖問題，比如人（蟲）口的增長[1-3]，它也可以來描述生物個(gè)體產(chǎn)生的腫瘤細(xì)胞數(shù)[4-5]。目前人們對(duì)腫瘤細(xì)胞生長模型進(jìn)行了連續(xù)的研究[6-19]，并取得了一系列優(yōu)秀的成果。其中，Zhong 等[6]研究了由純乘性噪聲驅(qū)動(dòng)的抗腫瘤細(xì)胞模型，并發(fā)現(xiàn)了噪聲誘導(dǎo)的隨機(jī)共振現(xiàn)象。Zeng 等[11]研究了免疫反應(yīng)下腫瘤細(xì)胞生長的數(shù)學(xué)模型。文獻(xiàn)[12]研究了時(shí)間延遲對(duì)由乘性噪聲和周期信號(hào)驅(qū)動(dòng)的抗腫瘤模型的影響。文獻(xiàn)[16]研究了由乘性非高斯噪聲與加性高斯色噪聲驅(qū)動(dòng)的抗腫瘤細(xì)胞模型，發(fā)現(xiàn)非高斯噪聲強(qiáng)度和偏離高斯噪聲參數(shù)的增加可以減少腫瘤細(xì)胞的數(shù)量。另外，該系統(tǒng)在噪聲增強(qiáng)穩(wěn)定性（NES）方面得到了研究者的廣泛重視[5，10-11，17]。

平均首次穿越時(shí)間（MFPT）可以用來量化噪聲對(duì)穩(wěn)定點(diǎn)狀態(tài)之間轉(zhuǎn)化行為的影響，表示系統(tǒng)粒子從一個(gè)勢(shì)阱中越過勢(shì)壘進(jìn)入另一個(gè)勢(shì)阱中所需時(shí)間的平均值，可以描述系統(tǒng)的瞬態(tài)性質(zhì)。文獻(xiàn)[17-18]分別研究了色噪聲和周期力對(duì)Logistic 模型中亞穩(wěn)態(tài)穩(wěn)定性的影響，并導(dǎo)出了從亞穩(wěn)態(tài)到穩(wěn)定狀態(tài)的平均首次穿越時(shí)間的表達(dá)式。文獻(xiàn)[20-26]對(duì)噪聲驅(qū)動(dòng)的雙穩(wěn)系統(tǒng)的平均首次穿越時(shí)間進(jìn)行了研究。文獻(xiàn)[22]運(yùn)用最速下降法[27-29]研究了由關(guān)聯(lián)噪聲驅(qū)動(dòng)的雙穩(wěn)系統(tǒng)的平均首次穿越時(shí)間。Li 等[23]研究了由白噪聲激勵(lì)驅(qū)動(dòng)的弱阻尼的一維和多維雙穩(wěn)系統(tǒng)的平均首次穿越時(shí)間，并比較了理論結(jié)果和數(shù)值結(jié)果。文獻(xiàn)[26]研究了乘性非高斯噪聲和加性高斯白噪聲驅(qū)動(dòng)的簡化FHN 神經(jīng)模型的平均首次通過時(shí)間。然而，上述研究中均沒有對(duì)Lévy 噪聲和高斯白噪聲共同激勵(lì)下腫瘤細(xì)胞生長模型的平均首次穿越時(shí)間進(jìn)行研究。

有鑒于此，本文探討了Lévy 噪聲和高斯白噪聲共同激勵(lì)下腫瘤細(xì)胞生長模型的平均首次穿越時(shí)間問題。首先給出了Lévy 噪聲和高斯白噪聲共同激勵(lì)下的腫瘤細(xì)胞生長模型，并對(duì)Lévy 噪聲做了簡單介紹；然后給出了Lévy 噪聲序列和高斯白噪聲序列的數(shù)值模擬方法，并利用四階Runge-Kutta 法和Janick-Weron算法計(jì)算出平均首次穿越時(shí)間；最后研究了穩(wěn)定性指標(biāo)α、偏斜參數(shù)β、Lévy 噪聲強(qiáng)度D 和高斯白噪聲強(qiáng)度Q 對(duì)系統(tǒng)平均首次穿越時(shí)間的影響，并進(jìn)一步分析了上述參數(shù)對(duì)腫瘤細(xì)胞生長的生物意義。

1 腫瘤細(xì)胞生長模型及Lévy噪聲

1.1 理論分析

本文考慮的是Logistic 系統(tǒng)下的腫瘤細(xì)胞生長模型[30-37]：

式中：xt為時(shí)間t 關(guān)于腫瘤細(xì)胞最大組織容量的歸一化密度；γ 為免疫速率，它與腫瘤細(xì)胞和免疫細(xì)胞合成復(fù)合物的速率及腫瘤細(xì)胞繁殖速率的比值成正比例；θ 表示腫瘤細(xì)胞和免疫細(xì)胞合成復(fù)合物的速率及復(fù)合物分離成免疫細(xì)胞和死亡腫瘤細(xì)胞的速率的比值的倒數(shù)，所以 θ 為常量，γ 與 θ 的關(guān)系為 θ ＜ 1，0 ＜γ ＜（1+ θ）2/4θ；x 為腫瘤細(xì)胞的初始密度，其勢(shì)函數(shù)為：

勢(shì)函數(shù)U（x）有一個(gè)最大值x2和兩個(gè)最小值x1和x3，如式（3）所示：

式中：x1代表健康狀態(tài)，基本上沒有腫瘤細(xì)胞；x3表示患病狀態(tài)，腫瘤細(xì)胞群處于穩(wěn)定水平；x2表示亞健康狀態(tài)，并且在有擾動(dòng)的情況下會(huì)移動(dòng)到x1或x3。當(dāng)θ=0.1、γ =3.0 時(shí)，x1=0，x2=4，x3=5，如圖1 所示。

圖1 勢(shì)函數(shù)圖像 θ =0.1，γ =3.0，x1=0，x2=4，x3=5Fig.1 Potential function U（x）with parameters，θ =0.1，γ =3.0，x1=0，x2=4，x3=5

模型（1）被用作描述腫瘤細(xì)胞增長的基礎(chǔ)模型，但腫瘤細(xì)胞的增長還會(huì)受到內(nèi)部噪聲和外部噪聲合力的影響。內(nèi)部噪聲主要通過基因突變產(chǎn)生，屬于腫瘤細(xì)胞的自發(fā)組織，可以引入一個(gè)乘性高斯白噪聲表示。外部噪聲可以看成是外部因素的影響，如溫度、藥品、輻射、化學(xué)試劑、寄主的免疫狀態(tài)和營養(yǎng)物質(zhì)的提供等也可以影響腫瘤細(xì)胞的增長，所以可以引入一個(gè)加性Lévy 噪聲來描述。引入噪聲后系統(tǒng)（1）可以寫為[30-37]：

式中：η（t）為高斯白噪聲，具有統(tǒng)計(jì)特性〈η（t）〉=0 和〈η（t）η（t′）〉=2Qδ（t-t′）；Q 為高斯白噪聲的噪聲強(qiáng)度；ξ（t）為穩(wěn)定性指標(biāo)為 α 的 Lévy 噪聲，它是 Lévy 過程 L（t）的時(shí)間導(dǎo)數(shù)。隨機(jī)變量ξ（t）的樣本數(shù)據(jù)服從Lévy 分布Lα，β（ξ；σ，μ），由于Lévy 噪聲的分布函數(shù)和概率密度函數(shù)都不具有顯性表達(dá)式，通?？捎檬剑?）所示的特征函數(shù)Φ（k）來表示[36-43]，即

故當(dāng) α∈（0，1）∪（1，2]時(shí)：

當(dāng) α =1 時(shí)，

式中：穩(wěn)定性指標(biāo) α∈（0，2]，α 決定了分布曲線尾部的平坦程度，當(dāng)α=2 時(shí)，Lévy 分布退化為高斯分布；偏斜參數(shù) β∈[-1，1]，β 決定了分布曲線的對(duì)稱性，β=0說明分布是對(duì)稱的，β ＞ 0 表示分布右偏，β ＜ 0 表示分布左偏；尺度參數(shù)σ∈[0，+∞]，用來表征分布的中心，Lévy 噪聲強(qiáng)度定義為 D= σα；μ∈R。

1.2 數(shù)值模擬

如果ξ 服從α 穩(wěn)定分布，則可以利用Janick-Weron（JW）算法[36-41]進(jìn)行模擬，即當(dāng) α≠1 時(shí)，

當(dāng) α =1 時(shí)，ξ 通過式（9）得到：

式中：V 為（-π/2，π/2）上的均勻分布；W 為均值為 1的指數(shù)分布；V 和W 相互獨(dú)立。

對(duì)高斯白噪聲η（t）的模擬方法[44-45]如下。首先生成一組獨(dú)立的標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)ζi，然后通過式（10），進(jìn)行計(jì)算即可得到高斯白噪聲序列{ηi}，其中Δt 是離散的時(shí)間步長。

由于Lévy 噪聲具有長尾和大面積跳躍的特點(diǎn)，使得系統(tǒng)方程（4）的Fokker-Planck 方程推導(dǎo)和求解十分困難，所以采用數(shù)值方法求解系統(tǒng)方程（4），方法為四階 Runge-Kutta 法和 Janick-Weron（JW）算法相結(jié)合[36-45]，具體如下：

式中：

2 平均首次穿越時(shí)間的計(jì)算

人們所關(guān)注的是腫瘤細(xì)胞是否會(huì)因內(nèi)部或外部的干擾而消亡，這就需要研究系統(tǒng)在引入噪聲后腫瘤細(xì)胞在減少（消亡）和復(fù)發(fā)這2 種狀態(tài)間的轉(zhuǎn)換時(shí)間，而平均首次穿越時(shí)間就可以用來描述這樣的瞬態(tài)性質(zhì)。根據(jù)平均首次穿越時(shí)間的定義[27-29]，運(yùn)用上節(jié)中描述的數(shù)值模擬方法可得到系統(tǒng)響應(yīng)的時(shí)間序列，然后記錄粒子從健康狀態(tài)x1所在勢(shì)阱穿越勢(shì)壘到達(dá)患病狀態(tài)x3所在勢(shì)阱需要的時(shí)間。因?yàn)榇讼到y(tǒng)不是對(duì)稱雙穩(wěn)系統(tǒng)，相反方向的平均首次穿越時(shí)間不相等，所以兩個(gè)相反方向的平均首次穿越時(shí)間都要記錄，最后對(duì)50 000 條響應(yīng)軌道進(jìn)行統(tǒng)計(jì)平均，時(shí)間步長Δt=0.01。

2.1 平均首次穿越時(shí)間作為Lévy噪聲強(qiáng)度D的函數(shù)

圖2 和圖3 分別給出了系統(tǒng)2 個(gè)方向的平均首次穿越時(shí)間T+（x1→x3）和T-（x3→x1）作為Lévy 噪聲強(qiáng)度D 的函數(shù)隨穩(wěn)定性指標(biāo)α 變化的情況。

由圖2 和圖3 可知，在固定參數(shù)Q=0.7 的情況下，無論是 β=0.5，還是 β =-0.5，隨著 Lévy 噪聲強(qiáng)度D 的增加，T+（x1→x3）均逐漸減小且3 條曲線的間距越來越小，而T-（x3→x1）均逐漸增大且3 條曲線的間距越來越大，這說明Lévy 噪聲強(qiáng)度對(duì)系統(tǒng)實(shí)現(xiàn)2 個(gè)狀態(tài)間躍遷的影響不同。此外還可以看出，T+（x1→x3）隨著α 的增大而減小，說明α 的增大使得粒子從x1→x3所需的躍遷時(shí)間減小，而T-（x3→x1）卻隨著α 的增大而增大。由圖2 和圖3 可知，當(dāng)β=0.5 和β=-0.5 時(shí)，兩組圖的差別很小，側(cè)面說明β 的變化對(duì)系統(tǒng)實(shí)現(xiàn)兩個(gè)狀態(tài)間躍遷的影響很小。

圖2 T+（x1→x3）作為 Lévy 噪聲強(qiáng)度 D 的函數(shù)隨 α 變化曲線Fig.2 Plot of T+（x1→x3）as function of D with different stability index α

圖3 T-（x3→x1）作為 Lévy 噪聲強(qiáng)度 D 的函數(shù)隨 α 變化曲線Fig.3 Plot of T-（x3→x1）as function of D with different stability index α

圖4 和圖5 分別給出了系統(tǒng)兩個(gè)方向的平均首次穿越時(shí)間T+（x1→x3）和T-（x3→x1）作為Lévy 噪聲強(qiáng)度D 的函數(shù)隨偏斜參數(shù)β 變化的情況。

圖4 T+（x1→x3）作為 Lévy 噪聲強(qiáng)度 D 的函數(shù)隨 β 變化曲線（Q=0.7，α =1.5）Fig.4 Plot of T+（x1→x3）as function of D with different β（Q=0.7，α =1.5）

圖5 T-（x3→x1）作為 Lévy 噪聲強(qiáng)度 D 的函數(shù)隨 β 變化曲線（Q=0.7，α =1.5）Fig.5 Plot of T-（x3→x1）as function of D with different β（Q=0.7，α =1.5）

由圖4 和圖5 可知，在固定參數(shù)Q=0.7，α=1.5的情況下，隨著Lévy 噪聲強(qiáng)度D 的增加，T+（x1→x3）均逐漸減小且曲線的間距越來越小，T-（x3→x1）逐漸增大且曲線的間距越來越大，這也印證了Lévy 噪聲強(qiáng)度D 對(duì)系統(tǒng)實(shí)現(xiàn)兩個(gè)狀態(tài)間的躍遷的影響不同。此外還可以看出，T+（x1→x3）隨著β 的增大而增加，T-（x3→x1）卻隨著β 的增大而減小。 隨著β 的增大，4 條曲線只有微弱的差別，說明β 的變化對(duì)系統(tǒng)實(shí)現(xiàn)兩個(gè)狀態(tài)間躍遷的影響很小。

從生物意義上來看，T+（x1→x3）減小說明腫瘤復(fù)發(fā)的可能性在增加，T-（x3→x1）減小說明腫瘤細(xì)胞消亡的可能性在增加。由圖2 和圖4 可以看出，隨著Lévy 噪聲強(qiáng)度D 增加，T+（x1→x3）逐漸減小，說明腫瘤細(xì)胞越容易復(fù)發(fā)。由圖3 和圖5 可以看出，隨著Lévy 噪聲強(qiáng)度D 增加，T-（x3→x1）在增大，說明腫瘤細(xì)胞消亡的可能性在減小?？傮w來說，Lévy 噪聲強(qiáng)度D的增強(qiáng)增加了粒子處于患病狀態(tài)的時(shí)間，不利于腫瘤細(xì)胞的消亡。

2.2 平均首次穿越時(shí)間作為高斯白噪聲強(qiáng)度Q的函數(shù)

圖6 和圖7 分別給出了系統(tǒng)兩個(gè)方向的平均首次穿越時(shí)間T+（x1→x3）和T-（x3→x1）作為高斯白噪聲強(qiáng)度Q 的函數(shù)隨穩(wěn)定性指標(biāo)α 變化的情況。

圖6 T+（x1→x3）作為高斯白噪聲強(qiáng)度Q 的函數(shù)隨α 變化曲線Fig.6 Plot of T+（x1→x3）as function of Q with different α

圖7 T-（x3→x1）作為高斯白噪聲強(qiáng)度 Q 的函數(shù)隨α 變化曲線Fig.7 Plot of T-（x3→x1）as function of Q with different α

由圖6 和圖7 可知，固定參數(shù)D=0.05，在β=0.5和β=-0.5 的兩種情況下，隨著高斯白噪聲強(qiáng)度Q 的增加，T+（x1→x3）和T-（x3→x1）均逐漸減小且3 條曲線的間距越來越小，說明乘性噪聲強(qiáng)度的增加有利于系統(tǒng)實(shí)現(xiàn)兩個(gè)狀態(tài)間的躍遷。此外還可看出，T+（x1→x3）隨著α 的增大而減小，T-（x3→x1）卻隨著α 的增大而增大。對(duì)于穩(wěn)定性指標(biāo)α 來說，其作用和Lévy 噪聲強(qiáng)度D 類似，也不利于腫瘤細(xì)胞的消亡。

圖8 和圖9 分別給出了系統(tǒng)兩個(gè)方向的平均首次穿越時(shí)間T+（x1→x3）和T-（x3→x1）作為高斯白噪聲強(qiáng)度Q 的函數(shù)隨偏斜參數(shù)β 變化的情況。

圖8 T+（x1→x3）作為高斯白噪聲強(qiáng)度 Q 的函數(shù)隨 β 變化曲線（D=0.05，α =1.5）Fig.8 Plot of T+（x1→x3）as function of Q with different β（D=0.05，α =1.5）

由圖8 和圖9 可看出，在固定參數(shù)D=0.05，α=1.5的情況下，隨著高斯白噪聲強(qiáng)度Q 的增加，T+（x1→x3）和T-（x3→x1）均逐漸減小，這也印證了乘性噪聲強(qiáng)度的增加有利于系統(tǒng)實(shí)現(xiàn)兩個(gè)狀態(tài)間的躍遷。此外還可以看出，T+（x1→x3）隨著β 的增大而增加，T-（x3→x1）卻隨著β 的增大而減小。

由圖4、圖5、圖8 和圖9 可以看出，隨著偏斜參數(shù)β 的增加，T+（x1→x3）逐漸增大，T-（x3→x1）逐漸減小，說明偏斜參數(shù)β 的增加有利于粒子向無腫瘤細(xì)胞狀態(tài)轉(zhuǎn)變，腫瘤細(xì)胞消亡的可能性在增加。但是隨著β的增大，曲線之間只有微弱的差別，這說明β 的增加雖有利于腫瘤細(xì)胞的消亡，但從整體來說，加性噪聲的增加還是不利于人體恢復(fù)健康。這可能與不正確的外部環(huán)境助長了腫瘤細(xì)胞的生長有關(guān)。由圖7 和圖9 可看出，隨著高斯白噪聲強(qiáng)度Q 的增加，T-（x3→x1）逐漸減小，說明當(dāng)腫瘤細(xì)胞的數(shù)量處于高水平時(shí)，增加高斯白噪聲的強(qiáng)度有利于腫瘤細(xì)胞的消亡。

3 結(jié) 論

本文研究了在腫瘤細(xì)胞的治療過程中，環(huán)境波動(dòng)對(duì)系統(tǒng)平均首次穿越時(shí)間的影響。為了更好地模擬現(xiàn)實(shí)情況，本文引入Lévy 噪聲和高斯白噪聲來模擬環(huán)境擾動(dòng)。利用四階Runge-Kutta 法和Janick-Weron 算法，分別對(duì)平均首次穿越時(shí)間作為Lévy 噪聲強(qiáng)度和高斯白噪聲強(qiáng)度的函數(shù)進(jìn)行研究。結(jié)果表明：

（1）隨著Lévy 噪聲強(qiáng)度D 的增加，T+（x1→x3）逐漸減小且3 條曲線的間距越來越小，而T-（x3→x1）逐漸增大且3 條曲線的間距越來越大，說明Lévy 噪聲強(qiáng)度對(duì)系統(tǒng)實(shí)現(xiàn)兩個(gè)狀態(tài)間的躍遷的影響不同。但從生物意義上來看，這種影響卻是一致的，Lévy 噪聲強(qiáng)度D 的增加，不利于腫瘤細(xì)胞的消亡，在實(shí)際處理中應(yīng)該避免。

（2）隨著高斯白噪聲強(qiáng)度Q 的增加，T+（x1→x3）和T-（x3→x1）均逐漸減小，說明高斯白噪聲強(qiáng)度的增加有利于系統(tǒng)實(shí)現(xiàn)兩個(gè)狀態(tài)間的躍遷。當(dāng)腫瘤細(xì)胞的數(shù)量處于高水平時(shí)，增加高斯白噪聲的強(qiáng)度有利于腫瘤細(xì)胞的消亡。

（3）對(duì)于穩(wěn)定性指標(biāo)α 來說，T+（x1→x3）隨著α 的增大而減小，T-（x3→x1）卻隨著α 的增大而增大，其作用和Lévy 噪聲強(qiáng)度D 類似；而偏斜參數(shù)β 起到了相反的作用。

（4）綜上所述，選擇適當(dāng)?shù)膮?shù)可以控制腫瘤細(xì)胞的增長，進(jìn)而達(dá)到理想的治療效果。