抗病毒藥物治療人類免疫缺陷病毒1型(HIV-1)的分?jǐn)?shù)階模型研究

崔周進(jìn),陳桂東

(1.江蘇海事職業(yè)技術(shù)學(xué)院基礎(chǔ)部,江蘇 南京 211170; 2.解放軍理工大學(xué)理學(xué)院數(shù)學(xué)教研中心,江蘇 南京 211101)

0 引 言

在對(duì)傳染病的諸多研究中，數(shù)學(xué)模型是檢驗(yàn)理論和定量評(píng)估猜想與結(jié)論的實(shí)驗(yàn)工具.自20世紀(jì)80年代首例艾滋病病例被發(fā)現(xiàn)以來，很多學(xué)者在HIV動(dòng)力學(xué)的建模上進(jìn)行研究，絕大多數(shù)的數(shù)學(xué)模型基于整數(shù)階微分方程.例如，文獻(xiàn)[1]中建立了如下HIV動(dòng)力學(xué)模型:

(1)

其中x和y分別表示t時(shí)刻未感染和已感染的CD4+T的數(shù)量，z表示t時(shí)刻血漿中病毒粒子的數(shù)目.系統(tǒng)中所有的系數(shù)均為非負(fù)數(shù)，s是CD4+T細(xì)胞的生成率，μ是人均死亡率，β是CD4+T細(xì)胞被病毒感染率，ε是被感染細(xì)胞的人均消失率，c為感染細(xì)胞與病毒的轉(zhuǎn)化率，而γ是病毒粒子的死亡率.

近幾十年來，由于在科學(xué)和工程眾多領(lǐng)域中的廣泛應(yīng)用，分?jǐn)?shù)階微積分(也即任意實(shí)數(shù)或復(fù)數(shù)階的微積分)受到人們?cè)絹碓蕉嗟年P(guān)注和重視.目前分?jǐn)?shù)階微積分已成為一個(gè)研究熱點(diǎn)，并在多個(gè)領(lǐng)域得到了應(yīng)用，如松弛、振蕩、擴(kuò)散和輸運(yùn)理論、生物組織、高分子材料的解鏈、混沌與湍流、隨機(jī)游走、統(tǒng)計(jì)與隨機(jī)過程、粘彈性力學(xué)及非牛頓流體力學(xué)等．這些領(lǐng)域的應(yīng)用研究又反過來促進(jìn)了分?jǐn)?shù)階微積分理論研究的進(jìn)一步發(fā)展.隨著分?jǐn)?shù)階微分方程的快速發(fā)展，很多研究者嘗試用分?jǐn)?shù)階微分方程建模[2-4].研究表明，在傳統(tǒng)的整數(shù)階方程不能建模的現(xiàn)象中，分?jǐn)?shù)階方程為此提供了可能性.分?jǐn)?shù)階與整數(shù)階模型最大的區(qū)別在于分?jǐn)?shù)階模型擁有記憶，而免疫反應(yīng)的主要特征恰恰包含記憶.為此Arafa等[2]建立如下HIV感染與免疫系統(tǒng)相互作用的分?jǐn)?shù)階動(dòng)力學(xué)模型：

(2)

其中Dα是α階Caputo微分算子，0<α1,α2,α3<1.作者利用廣義Taylor公式與同倫分析法，研究了系統(tǒng)(2)的簡(jiǎn)化模型的解.本文的主要目的是在前期研究的基礎(chǔ)上，探討分?jǐn)?shù)階系統(tǒng)(2)的平衡點(diǎn)的穩(wěn)定性，并通過數(shù)值模擬驗(yàn)證理論分析的結(jié)果.

1 動(dòng)力學(xué)分析

分?jǐn)?shù)階的導(dǎo)數(shù)有許多定義[5-6]，例如Riemann-Liouville定義、Gruunwald-Letnikow定義、Caputo定義、Miller-Ross定義等等.為了方便起見,首先給出本文中采用的Caputo微分算子.

定義1(Caputo定義的分?jǐn)?shù)階導(dǎo)數(shù))

(3)

其中n=[α]+1，Γ為Gamma函數(shù).

Caputo定義由于其Laplace變換公式和整數(shù)階導(dǎo)數(shù)的Laplace變換公式具有相同的形式, 因而在控制理論中廣泛采用.本文研究分?jǐn)?shù)階線性系統(tǒng)的穩(wěn)定性, 采用Caputo分?jǐn)?shù)階導(dǎo)數(shù)定義是合適的.

在平衡點(diǎn)E1=(s/μ,0,0)處，系統(tǒng)(2)的Jacobian矩陣為

那么E1點(diǎn)處對(duì)應(yīng)的特征方程為

在平衡點(diǎn)E2點(diǎn)處對(duì)應(yīng)的特征方程為

λ3+a1λ2+a2λ+a3=0，

定理11)若D(P)>0，則E2局部漸近穩(wěn)定當(dāng)且僅當(dāng)a1>0,a1a2>a3>0.

2)若D(P)<0，a1≥0,a2≥0,a3>0，0.5<αi<2/3，則E2局部漸近穩(wěn)定.

3)若D(P)<0，a1>0,a2>0,a1a2=a3，0.5<αi<1，則E2局部漸近穩(wěn)定.

4)若D(P)<0，a1<0,a2<0，αi>2/3，則E2不穩(wěn)定.

2 數(shù)值模擬

圖1 α1=0.55，α2=0.6，α3=0.65時(shí)的相位圖

治療藥物會(huì)影響產(chǎn)生病毒顆粒的感染細(xì)胞的死亡率. 但是因?yàn)樗幬锟赡懿皇?00%有效，因此只有部分受感染的細(xì)胞將恢復(fù)到未感染類.文獻(xiàn)[2]給出模型(2)中典型的參數(shù)取值，取s=0.272，μ=0.001 36，ε=0.33，γ=2，β=0.000 27，c=50，帶有如下初值條件x(0)=200,y(0)=0,z(0)=1，則平衡點(diǎn)E2=(48.89,0.6228,15.57)，a1=2.34，a2=0.012 96，a3=0.002 77，D(P)<0，a1>0，a2>0，a3>0.根據(jù)定理1可得，當(dāng)0.5<αi<2/3時(shí)，平衡點(diǎn)E2局部漸近穩(wěn)定.相位圖被用來確定系統(tǒng)的動(dòng)力學(xué)行為，借助于Adams-Bashforth-Moulton預(yù)估校正算法，利用Matlab軟件可以得到模型(2)在不同階導(dǎo)數(shù)下的相位圖.圖1是模型(2)中分?jǐn)?shù)階導(dǎo)數(shù)分別取成α1=0.55，α2=0.6，α3=0.65時(shí)的相位圖；圖2是模型(2)中分?jǐn)?shù)階導(dǎo)數(shù)取成α1=α2=α3=0.6時(shí)的相位圖；圖3是模型(2)中導(dǎo)數(shù)取為α1=α2=α3=1.0時(shí)的相位圖，通過圖形可以驗(yàn)證上述結(jié)論的正確性.

圖2 α1=α2=α3=0.6時(shí)的相位圖

圖3 α1=α2=α3=1.0時(shí)的相位圖

3 總 結(jié)

傳染病動(dòng)力學(xué)是對(duì)傳染病進(jìn)行理論性定量研究的一種重要方法，是根據(jù)疾病的發(fā)生及在種群內(nèi)的傳播、發(fā)展規(guī)律，以及與之有關(guān)的社會(huì)等因素，建立能反映傳染病動(dòng)力學(xué)特性的數(shù)學(xué)模型.通過對(duì)模型動(dòng)力學(xué)性態(tài)的定性、定量分析和數(shù)值模擬，來分析疾病的發(fā)展過程，揭示流行規(guī)律，預(yù)測(cè)變化趨勢(shì)，分析疾病流行的原因和關(guān)鍵因素，尋求預(yù)防和控制的最優(yōu)策略，為防治決策提供理論依據(jù).在生物學(xué)中，研究人員發(fā)現(xiàn)，生物細(xì)胞膜具有分?jǐn)?shù)階的電子傳導(dǎo)性，因而歸類為分?jǐn)?shù)階模型；另外，需要特別強(qiáng)調(diào)的是分?jǐn)?shù)階微分算子與記憶有著重大關(guān)聯(lián)，已經(jīng)證實(shí)由分?jǐn)?shù)階微分方程建立的某些生物學(xué)模型比整數(shù)階更有優(yōu)勢(shì)．隨著分?jǐn)?shù)階微分方程的發(fā)展，分?jǐn)?shù)階模型更加向?qū)嶋H靠攏.本文中所利用的分?jǐn)?shù)階的Routh-Hurwitz準(zhǔn)則，Adams-Bashforth-Moulton預(yù)估校正算法可以推廣到類似的模型的研究中去.關(guān)于分?jǐn)?shù)階系統(tǒng)的動(dòng)力學(xué)性質(zhì)，例如分?jǐn)?shù)階傳染病模型平衡點(diǎn)的存在性和穩(wěn)定性，分?jǐn)?shù)階微分系統(tǒng)階數(shù)如何影響模型的動(dòng)力學(xué)性態(tài)，分?jǐn)?shù)階時(shí)滯傳染病模型的最優(yōu)控制問題，并給出分?jǐn)?shù)階動(dòng)力系統(tǒng)控制的條件，這些都是下一步研究所關(guān)注的方向.

[1]Tuckwell H C, Wan F Y M. Nature of equilibria and effects of drug treatments in some simple viral population dynamical models[J]. IMA J Math Appl Med Biol,2000，17:311-327.

[2]Arafa A A M, Rida S Z, Khalil M. The effect of anti-viral drug treatment of human immunodeficiency virus type 1(HIV-1) described by a fractional order model[J]. Applied Mathematical Modelling,2013,37:2189-2196.

[3]Cui Z J, Mao Z S, Yang S J,etal. Approximate analytical solutions of fractional perturbed diffusion equation by reduced differential transform method and the homotopy perturbation method[J]. Mathematical Problems in Engineering,2013,Article ID 186934.

[4]Cui Z J, Yu P N, Zong W. Dynamical behaviors and chaos in a new fractional-order financial system: proceeding of the Fifth International Workshop on Chaos-fractals Theories and Applications[C].Dalian:5th IWCFTA,2012:109-113.

[5]Podlubny I. Fractional differential equations: vol.198 of Mathematics in Science and Engineering[M]. San Diego, Calif：Academic Press,1999.

[6]Hilfer R. Applications of fractional calculus in physics[M]. Singapore：World Scientific,2000.

[7]Ahmed E, E1-Sayed A M A, E1-Saka H A A. On some routh-hurwitz conditions for fractional order diferential equations and their applications in lorenz，rssler，Chua and Chen systems[J]. Physics Letters A,2006,358:1-4．

[8]Faieghi M R, Delavari H. Chaos in fractional-order Genesio-Tesi system and its synchronization[J]. Commun Nonlinear Sci Numer Simulat,2012,17(2):731-741．