紀(jì)元昕，朱家明，王茜瑤，胡學(xué)峰

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

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

1 埃博拉病毒傳染機(jī)制模型

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

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

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

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

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

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

1.2 模型建立

結(jié)合阻滯增長模型和SIS模型進(jìn)行建模。

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

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

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

針對(duì)利比里亞，我們進(jìn)行作圖分析可得，見圖2，3。

圖1 預(yù)測曲線與真實(shí)曲線比較圖

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

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

1.3 結(jié)果綜述

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

表1 各國患病情況信息表

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

2.1 研究思路

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

2.2 模型建立

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

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

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

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

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

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

圖4 登革熱患病率變化圖

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

3 制藥公司的選取與合作

3.1 原理分析

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

我們選取科研攻關(guān)能力a，原材料生產(chǎn)能力b，藥物檢驗(yàn)?zāi)芰，藥物后期跟蹤與反饋能力d。

查閱文獻(xiàn)資料，給出如下四個(gè)指標(biāo)的計(jì)算公式：科研攻關(guān)能力a，原材料生產(chǎn)能力b，藥物檢驗(yàn)試驗(yàn)?zāi)芰，該公司資金流動(dòng)能力d，并研究這四個(gè)指標(biāo)對(duì)降低制藥成本和風(fēng)險(xiǎn)的影響。

3.2 模型建立

對(duì)任一制藥公司i，令x1為該公司科研部門人數(shù)，X是該公司總?cè)藬?shù)，平均研制一款新藥的速度v1，平均每年專利申請(qǐng)數(shù)t，可建立數(shù)學(xué)表達(dá)式：

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

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

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

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

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

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

圖6 ρ變化圖

4 結(jié)論

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

［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.