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

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

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

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

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

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

1.1 理論分析

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

式中：xt為時間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ù)U（x）有一個最大值x2和兩個最小值x1和x3，如式（3）所示：

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

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

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

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

當(dāng) α =1 時，

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

1.2 數(shù)值模擬

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

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

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

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

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

式中：

2 平均首次穿越時間的計算

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

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

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

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

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

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

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

圖4 T+（x1→x3）作為 Lévy 噪聲強度 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 噪聲強度 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 噪聲強度D 的增加，T+（x1→x3）均逐漸減小且曲線的間距越來越小，T-（x3→x1）逐漸增大且曲線的間距越來越大，這也印證了Lévy 噪聲強度D 對系統(tǒng)實現(xiàn)兩個狀態(tài)間的躍遷的影響不同。此外還可以看出，T+（x1→x3）隨著β 的增大而增加，T-（x3→x1）卻隨著β 的增大而減小。 隨著β 的增大，4 條曲線只有微弱的差別，說明β 的變化對系統(tǒng)實現(xiàn)兩個狀態(tài)間躍遷的影響很小。

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

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

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

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

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

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

圖8 T+（x1→x3）作為高斯白噪聲強度 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的情況下，隨著高斯白噪聲強度Q 的增加，T+（x1→x3）和T-（x3→x1）均逐漸減小，這也印證了乘性噪聲強度的增加有利于系統(tǒng)實現(xiàn)兩個狀態(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 可看出，隨著高斯白噪聲強度Q 的增加，T-（x3→x1）逐漸減小，說明當(dāng)腫瘤細(xì)胞的數(shù)量處于高水平時，增加高斯白噪聲的強度有利于腫瘤細(xì)胞的消亡。

3 結(jié) 論

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

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

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

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

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