埃博拉病毒傳染機制與研藥策略的數(shù)模分析

紀元昕，朱家明，王茜瑤，胡學峰

（安徽財經(jīng)大學統(tǒng)計與應用數(shù)學學院；安徽蚌埠233030）

埃博拉病毒經(jīng)歷了大爆發(fā)穿過非洲大陸在2014年，感染人數(shù)之多，死亡率之高引起了全球廣泛關注。同時，埃博拉病毒的傳染居然周期爆發(fā)性，主要針對非洲一帶地區(qū)。因此研究埃博拉的傳染機制與內(nèi)在規(guī)律是進一步抵制并消除病毒的必要前提［1-3］；之后，是否應研制藥物，研發(fā)周期與研藥量均是需要解決的問題；最后，政府如何選取制藥公司進行低成本高水平的制藥開發(fā)是我們在疾病理論基礎上展開實際解決方案的手段。

1 埃博拉病毒傳染機制模型

1.1 數(shù)據(jù)分析

查找到Guinea，Liberia，Sierra Leone等國家感染人數(shù)各天累計的感染人數(shù)加和作為總體的患病情況，通過相關數(shù)據(jù)我們可以觀察到的是在該病毒從爆發(fā)到傳染的較長一段時間內(nèi)該病毒，該病毒滿足指數(shù)增長模型，該疫情被發(fā)現(xiàn)的前中期，因為缺乏必要的醫(yī)療手段和醫(yī)療常識，該病毒的繁衍速度可近似看作是自然條件下種群的增長［4-6］。

為了有效將時間進行排序，定義2014年3月25日為發(fā)現(xiàn)疫情的第一天，以后分別按照天數(shù)間隔進行排序，根據(jù)散點圖的大概走勢，猜測在前期及0～220 d之間，感染人數(shù)的增長趨勢呈指數(shù)增長，令：x（t）=x0ert對函數(shù)進行擬合。

兩邊同時取 ln，y=rt+a，y=ln x，a=ln x0，利用 Matlab 輸入相關數(shù)據(jù)處理得擬合函數(shù)：y=0.02x+4.7，x0=109.95；所以有：

做出在0～300內(nèi)的預測曲線和在真實變化趨勢。

從圖1中可以看出基本上在0～255 d內(nèi)，感染人數(shù)基本上呈指數(shù)增長型趨勢，但在255天以后真實的感染人數(shù)較指數(shù)增長明顯減緩，這從宏觀上說明埃博拉病毒的傳播速度在一定程度上受到了阻礙。

1.2 模型建立

結合阻滯增長模型和SIS模型進行建模。

在人類還沒有有效措施能夠制止埃博拉疫情的傳播并治愈已經(jīng)感染的病人時，將埃博病毒的傳播看成是一個病毒種群在自然環(huán)境下的種群數(shù)量增長情況。因此可以運用阻滯增長模型—logistic模型對其發(fā)展規(guī)律進行研究。

求出阻滯模型變化率曲線，得出在253 d時非洲三國埃博拉病毒的感染者增長率最大為320，大概從500 d開始，病人總數(shù)趨于穩(wěn)定，最大感染人數(shù)為40000人。

現(xiàn)在針對三個國家的疫情在自然增長情況下分別進行建模。從三個國家的患病人數(shù)與時間變化圖中，我們可以發(fā)現(xiàn)三個國家的感染情況基本符合阻滯增長的形式，這一點與總的趨勢也相同。三個國家的患病人數(shù)初始階段呈現(xiàn)指數(shù)增長模式，后來增加速度在250 d之后普遍出現(xiàn)變慢，這說明疫情在隔離，治療等相關措施下受到了一定控制。

針對利比里亞，我們進行作圖分析可得，見圖2，3。

圖1 預測曲線與真實曲線比較圖

圖2 利比內(nèi)亞感染人數(shù)變化圖

圖3 利比內(nèi)亞患病增加率

1.3 結果綜述

從表中可看出，最早達到增長高峰的是幾內(nèi)亞為191 d，最遲的是塞拉利昂在256 d，按照患病前在人數(shù)的多少獲取權重進行求和，研發(fā)期限最長為為244.35 d。

表1 各國患病情況信息表

2 藥物研發(fā)必要性判別

2.1 研究思路

為了研究是否有其他疾病需要新藥研制，我們選取該地區(qū)較為顯著的疾病種類：肺結核，艾滋病與馬來熱。根據(jù)這三種疾病在該階段內(nèi)的感染人數(shù)，將其與埃博拉感染人數(shù)進行比較判斷，進行方差檢驗，若有顯著差異說明該種疾病需要新藥的救治。在對疾病進行判斷之后，我們對各疾病每年感染人數(shù)提取異常值，即疾病未得到有效控制有進一步爆發(fā)趨勢的年份，將其作為研制新藥的最遲開始時間，制藥期限以病情下次大爆發(fā)為界，進行控制。

2.2 模型建立

由于在上一題中我們估測出了埃博拉病毒的感染人數(shù)，且該估測是在無其他新藥開發(fā)的情況下得到的。因此在本題的求解中，我們可依據(jù)埃博拉病毒感染數(shù)作為基準來研究其他疾病感染的波動情況。由于疾病類型的復雜多樣，為了簡便計算，我們選取了烏干達地區(qū)較為顯著的三大類傳染性疾?。悍谓Y核，艾滋病以及登革熱。肺結核與艾滋病數(shù)據(jù)較為完整齊全，而登革熱數(shù)據(jù)有所缺失，因此我們采用多項式擬合近似求的其他年份的感染數(shù)據(jù)。

根據(jù)現(xiàn)有的2000到2008年的烏干達登革熱患病人數(shù)，用多項式擬合出2009-2012年的患病人數(shù)。具體患病數(shù)見圖4。我們將登革熱患病人數(shù)與肺結核和艾滋病感染人數(shù)進行簡單加權作為總體患病人數(shù)進行綜合分析。此時，為了判斷該階段是否需要研制新藥，可將整體的疫情數(shù)據(jù)與之前估測的埃博拉感染數(shù)據(jù)進行方差分析檢驗，若有顯著差異，說明整體的疾病未如同埃博拉一般得到一定控制，而有持續(xù)爆發(fā)加重的趨勢，則需要研發(fā)新藥。若無顯著差異，則無需研制新藥。

通過方差分析，得到：F=529.97，F(xiàn) 值較大，且 P=0＜0.05，通過了顯著性檢驗，認為兩個數(shù)據(jù)之間存在顯著性差異，即總體疾病疫情上有存在其他疾病肆虐的情況，需要研發(fā)新藥物進行治療。

對于藥物研制的最遲期限與可能的總產(chǎn)量，通過對每年的總體疾病感染人數(shù)進行觀測，找出數(shù)據(jù)的異常波動點，作為研發(fā)開端的依據(jù)，并將藥物研制周期控制在疾病的下一次異常爆發(fā)之前，使之得以投入市場，阻止疫情的下一次爆發(fā)。

通過Matlab做出總體患傳染病人數(shù)的波動圖，見圖5。

由圖可看出，在2004年以前總體疾病人數(shù)趨于平緩，但在2004年開始突然大幅度增長，說明此時有新的疫情的爆發(fā)，并為需研制新藥的最佳開始時間。

圖4 登革熱患病率變化圖

圖5 總患病人數(shù)波動圖

3 制藥公司的選取與合作

3.1 原理分析

考慮到公司聯(lián)合降低制藥要有效降低成本和風險，就需要在一定公司范圍內(nèi)選取關于制藥過程中幾個花費和風險最低的部分進行組合，這些部分來源的公司便是最優(yōu)組合的公司。

我們選取科研攻關能力a，原材料生產(chǎn)能力b，藥物檢驗能力c，藥物后期跟蹤與反饋能力d。

查閱文獻資料，給出如下四個指標的計算公式：科研攻關能力a，原材料生產(chǎn)能力b，藥物檢驗試驗能力c，該公司資金流動能力d，并研究這四個指標對降低制藥成本和風險的影響。

3.2 模型建立

對任一制藥公司i，令x1為該公司科研部門人數(shù)，X是該公司總人數(shù)，平均研制一款新藥的速度v1，平均每年專利申請數(shù)t，可建立數(shù)學表達式：

令y1為該公司近三年平均每年生產(chǎn)藥物量，Y為同年該型藥物研制總量，v2為平均生產(chǎn)一噸新藥速度：

令ys為該公司該年生產(chǎn)實際銷售該種藥品總數(shù)，y為該年該藥物生產(chǎn)量，v2為該公司檢驗藥品科研人數(shù)，

令g，c分別為該年該公司近三年平均每年營業(yè)額和所有成本，r%為該公司營業(yè)額增長率，R為

建立一種限制數(shù)量上不斷增長的方法：在選定的公司中如果對于一定范圍內(nèi)i=1，2…m…k（k≥m），隨著i的不斷增長T也不斷增長。

首先對任意公司i求出其綜合能力：T=λ1ai+λ2bi+λ3ci+λ4di然后按照綜合能力排序 Tm＞Tk＞Tl＞…＞Tmin，按照從大到小的順序依次作出構造一個在數(shù)量上隨著i增長 T 不斷增長的數(shù)列 y1，y2，y3…yn+k。求其增長速度 ρ，以及閾值η，在增長速度小于閾值時我們認為此時增加公司的數(shù)目已經(jīng)沒有必要。

因此，在滿足條件時，選用公司組合進行聯(lián)合制藥的效率最高。

圖6 ρ變化圖

4 結論

以上各模型經(jīng)過各軟件的檢驗，具有一定的合理性。但同時模型的建立必須依靠必要的假設，在本文中我們假設埃博拉病毒與其他疾病無交叉作用，并在第二題運用第一題的埃博拉傳染病模型時忽略該病毒的周期性，因此與實際的疾病數(shù)會有偏差。綜上來說，我們依據(jù)埃博拉現(xiàn)有的感染病例數(shù)建立了較符合它傳播特征的傳染病模型。同時對藥物研發(fā)方面進行了研究與推測，最后對制藥公司的評估與選擇提供了較新穎的理論方案。

［1］Xinxin Wang,Shengqiang Liu.An epidemic model with different distributed latencies ［J］.Applied Mathematics and Computation 2014,241:259-266.

［2］Li D,Gui C,Luo X.Impulsive Vaccination SEIR Model with Nonlinear Incidence Rate and Time Delay ［J］.MathematicalProblemsin Engineering,2013.

［3］Buonomo B,d’Onofrio A,Lacitignola D.Modeling of pseudo-rational exemption to vaccination for SEIRdiseases［J］.Journal of Mathematical Analysisand Applications,2013,404（2）:385-398.

［4］Zhang J,Jia J,Song X.Analysis of an SEIREpidemic Model with Saturated Incidence and Saturated Treatment Function［J］.The Scientific World Journal,2014.

［5］Liu M,Bai C,Wang K.Asymptotic stability of a two-group stochastic SEIRmodel with infinite delays［J］.Communications in Nonlinear Scienceand Numerical Simulation,2014,19（10）:3444-3453.

［6］Wang X,Wei L,Zhang J.Dynamical analysisand perturbation solution of an SEIR epidemic model［J］.Applied Mathematics and Computation,2014,232:479-486.