隨著計算生物學和系統生物學研究應用范圍的不斷拓展,復雜的生物過程需要不同的分析方法去應對。針對典型的流行病免疫防治過程的動態模型,進行分析評價。系統模型經離散化后,應用非線性規劃法求解,獲得流行病傳播免疫過程中各項指標的動態仿真結果;依據生化系統理論方法將動態模型轉換為S系統,其特征值表明該系統局部穩定,并存在振蕩特性;接著分析了速率常數對模型狀態影響的靈敏度以及三個關鍵參數對系統指標的對數增益影響,通過靈敏度分析,評價了系統的魯棒性,為今后模型的建立與修正工作做出指導;同時,此項工作也拓寬了生化系統理論的應用范圍。
引用本文： 潘多濤, 史洪巖, 黃明忠, 袁德成. 流行病防治過程SIR模型的穩定性分析. 生物醫學工程學雜志, 2015, 32(5): 1013-1018. doi: 10.7507/1001-5515.20150180 復制
引言
在當下的信息時代，生物學家已經意識到傳統通過實驗的研究方式面臨成本高、效率低的問題。在計算機應用領域不斷拓展的背景下，計算生物學與系統生物學得到快速發展。依據長期積累的實驗數據，利用數學建模的方法去研究并預測不同生物過程的動態特性進而有效指導實際實驗研究方向，已成為生物研究領域的關注重點。隨著研究的不斷深入，不同類型生物過程的復雜程度也同時提高，與之對應的模型評價與分析工作也越來越困難，這勢必需要新的建模手段和分析方法去應對[1-4]。
流行病是在外部環境下，于短時間內在人群中廣泛傳播并極可能面臨大面積感染的傳染病。病原體的抗原變異在流行病傳播過程中是常見的現象，對流行病的防治是一個巨大的困難。疫苗可使機體產生特異性抗體，眾所周知，通過疫苗的接種可使人體獲得對某種病原體的免疫力，比如通過種痘的方式已在全世界范圍徹底抑制了天花。但目前，仍然存在很多沒有疫苗可應對的傳染病，這對流行病防治來說是極大的挑戰。通過建立流行病的數學模型，可較準確地評價和預測流行病的傳播、發展趨勢以及評估疫苗的有效性，有利于提前試行有效的防治措施，目前在這方面已有大量的工作報道[5-6]，但多集中在模型的建立和參數辨識上[7-8]，而模型的有效性與評價方法研究相對較少，由于流行病學的實際數據需要長期的積累，因此模型的可靠性的提高能夠保證數據應用的時效性。
針對流行病防治過程的“易感-感染-恢復”(Susceptible-Infective-Recovered，SIR)模型，在前期工作基礎上[9-10]，采用生化系統理論(Biochemical Systems Theory，BST)的方法對其進行評價，對模型特性進行分析并評價其魯棒性。
1 方法
不同類型生物過程的動態模型大多采用微分方程或微分代數方程來描述，針對難以分析的復雜模型，可利用BST理論方法建立與之對應的S系統來幫助分析。現有許多工具可用來求解，除MATLAB等常規軟件外，還可選用開源免費的軟件包，如微分/代數解題包SUNDIALS[11]和非線性規劃解題器IPOPT[9, 12]等。在使用計算機對生物數學模型進行分析時，為確保相關信息的正確描述，一般需要一種通用的數據描述、儲存規范。目前系統生物學標記語言(Systems Biology Markup Language，SBML)使用最為廣泛，當前已有多種公用模型數據庫[13]以SBML為主要的儲存格式。
1.1 動態模擬
不同規模的生物過程通常描述為復雜程度不同的非線性系統，其動態特性采用前期工作中開發的通用可擴展集成規劃工具GIEPT[9]求解。當生物系統的復雜性達到一定程度后，其數學模型不能采用常規的數值積分法求解，須將其轉化為非線性規劃問題(nonlinear programming，NLP)后采用動態優化的方法進行計算。通常NLP的標準形式如下：
| $ \begin{array}{l} min\int\limits_0^{{t_f}} {\Phi \left( {{z_d}\left( t \right),{z_a}\left( t \right),u\left( t \right),\theta ,\varepsilon } \right){\rm{d}}} t\\ {\rm{s}}{\rm{.t}}{\rm{.:}}\frac{{{\rm{d}}{z_d}}}{{{\rm{d}}t}} = {f_d}\left( {{z_d}\left( t \right),{z_a}\left( t \right),u\left( t \right),\theta ,\varepsilon } \right)\\ \;\;\;\;\;\;\;\;0 = {f_a}\left( {{z_d}\left( t \right),{z_a}\left( t \right),u\left( t \right),\theta ,\varepsilon } \right)\\ \;\;\;\;\;\;\;\;\;z_d^L \le {z_d}\left( t \right) \le z_d^U\\ \;\;\;\;\;\;\;\;\;z_a^L \le {z_a}\left( t \right) \le z_a^U\\ \;\;\;\;\;\;\;\;\;{u^L} \le u\left( t \right) \le {u^U}\\ \;\;\;\;\;\;\;\;\;z\left( 0 \right) = {z_0} \end{array} $ |
其中：zd(t)、za(t)和u(t)分別表示系統的微分狀態變量、代數狀態變量以及控制變量；θ表示系統內可變參數(為降低復雜度，只考慮經簡化的初始狀態)，ε描述了系統內具有固定值的參數。選取N個離散點；采用文獻[14-16]報道的方法將模型離散化，并構造非線性規劃的標準數學描述。
1.2 BST模型
BST[17-18]是一種可對系統進行模擬分析的數學計算方法。最初，BST主要是針對生化過程設計提出的，本文將其應用到相近的流行病傳播過程的分析中。
BST是基于動力學的方法，其最為常用的數學描述稱為S系統，它是由一種冪律函數描述的。取X代表目標系統中的某一狀態量，將對X有貢獻(流入)的因素合計為正向流量(記為V+)，使X有損失(流出)的因素合計為負向流量(記為V-)，則S系統可由下面的數學式表示：
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{i}}=V_{i}^{+}-V_{i}^{-}={{\alpha }_{i}}\underset{j=1}{\overset{n+m}{\mathop \prod }}\,X_{j}^{{{g}_{i,j}}}-{{\beta }_{i}}\underset{j=1}{\overset{n+m}{\mathop \prod }}\,X_{j}^{{{h}_{i,j}}},i=1,2,\cdots n $ |
其中Xi是系統的狀態變量(亦稱之為非獨立變量)；Xj是包括狀態變量在內的所有影響系統響應速率的因素；n為狀態變量的個數，m為獨立變量的個數；αi為流入(正向)速率常數，βi為流出(負向)速率常數；gi,j和hi,j分別為各因素的動力階。速率常數和動力階分別由下式給出：
| $ {g_{i,j}} = {\left. {\frac{{\partial V_i^ + }}{{\partial {X_j}}} \cdot \frac{{{X_j}}}{{V_i^ + }}} \right|_{op}},{h_{i,j}} = {\left. {\frac{{\partial V_i^ - }}{{\partial {X_j}}} \cdot \frac{{{X_j}}}{{V_i^j}}} \right|_{op}} $ |
| $ {\alpha _i} = {\left. {\frac{{V_i^ + }}{{_{\mathop \prod \limits_{j = 1}^{n + m} X_j^{{g_{i,j}}}}}}} \right|_{op}},{\beta _i} = {\left. {\frac{{\partial V_i^ - }}{{\mathop \prod \limits_{j = 1}^{n + m} X_j^{{h_{i,j}}}}}} \right|_{op}} $ |
其中OP為工作點，一般取自系統的穩態值。
由于S系統具有其他動態系統所不具備的特性，使得借助BST方法，能夠對生物系統進行數學分析。尤其體現在，當生物系統的規模逐步拓展，采用微分方程描述的機理模型的復雜性會無限提高，但與之對應的S系統的基本結構仍會保持不變。這就使得S系統能夠應對不同規模的系統[18]。
2 應用
當前研究較多的流行病模型大多以Kermack和McKendrick提出的經典倉室模型[5]為基礎，本文研究的模型即來源于典型的SIR流行病學模型[19]。該模型描述了疫苗接種防治過程中的免疫逃避動態過程：在某一群體中，通過接種疫苗使得原有土著致病菌株(百日咳)維持較低的發病率，引入突變菌株后，假定兩種菌株具有相同繁殖率，則這兩種菌株可達到共存的相對穩定狀態。文獻中原有數學模型狀態變量(或參數)與S系統因素的對應關系如表 1所示。
參數值來源于在線數據庫BioModels[13]。其中三個參數：疫苗覆蓋率(即新生兒接種率)p，接種后易感人群數量的下降率τ以及接種后傳染周期的縮減率η是接下來模型分析的考察目標。
表1
系統原參量與S系統因素的對應關系
Table1.
Correspondence between original variables and elements of S-system
2.1 動態過程仿真
使用GIEPT調用IPOPT進行求解。在利用IPOPT求解之前，如前所述，需將動態模型進行離散化，轉化成NLP。由于原有數學模型為單純的微分方程組，并未有顯式的目標函數，因此必須構造一個恰當的目標函數形式。這里，以離散點處斜率估算誤差總和最小為目標函數，如式(5)所示。
| $ \min :\sum\limits_{i}^{n}{\sum\limits_{j}^{N}{\left\| {{\overset{\centerdot }{\mathop{x}}\,}_{i,j}}-\frac{{{x}_{i,j+1}}-{{x}_{i,j}}}{{{t}_{i,j+1}}-{{t}_{i,j}}} \right\|}} $ |
式(5)中，x代表SIR系統數學模型中的狀態變量；n為狀態變量的個數；N為離散分段數。至此，形如式(1)的標準形式已經具備，可以采用非線性規劃的方法進行求解。為了避免數據點過多產生堆積，每隔6個時刻繪制一次數據點，如圖 1所示。
圖1
SIR系統的動態過程
子圖為
sub-figure describes the dynamic process of
從結果可以看出，該過程是一個非常復雜的非線性過程，伴隨有不同程度的振蕩，圖示中各指標經過振蕩后均能達到穩定，這種具有減幅振蕩的暫態過程表明系統的特征值具有非零的虛部。考察SIR系統的主要指標S(易感比率)、I(感染比率)以及R(恢復比率)，三種指標特別是土著菌株與突變菌株的感染率(I1、I2)的變化情況符合文獻[19]中的描述：兩種菌株達到共存的穩定狀態(圖 1中子圖)。
圖 1所示的SIR系統動態過程事實上是一種理想狀態，其描述了雙菌株在長期的致病及接種疫苗防治過程中達到穩定的情形，而在實際情況中，若再出現新的突變致病菌株，則有流行病大爆發的危險；另一方面，多菌株病原體的出現，也給疫苗研究工作帶來很大困難。因此，通過模型預測流行病的防治效果及潛在威脅就顯得尤為重要，而這其中的關鍵問題是模型是否可靠。
2.2 模型穩定性與靈敏度分析
利用前述介紹的BST理論方法，將原有的SIR模型微分方程中的因式分為正向和負向部分，依照表 1的參量對應關系，將其轉化為BST模型(即S系統)，如式(6)所示。
| $ \begin{align} & {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{1}}=0.0775X_{1}^{-0.0091}X_{4}^{0.0305}X_{5}^{0.0525}X_{6}^{0.0459}X_{10}^{0.8803}X_{12}^{-0.2459}- \\ & 1298.2{{X}_{1}}X_{2}^{0.2424}X_{3}^{0.3898}X_{7}^{0.0985}X_{8}^{0.1309}X_{9}^{0.1385} \\ \end{align} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{2}}=539.6934{{X}_{1}}X_{2}^{0.711}X_{7}^{0.289}-17.4{{X}_{2}} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{3}}=772.1643{{X}_{1}}X_{3}^{0.5914}X_{8}^{0.1986}X_{9}^{0.2101}-17.4{{X}_{3}} $ |
| ${{\overset{\centerdot }{\mathop{\text{X}}}\,}_{4}}=17.38{{X}_{2}}-243.0366X_{3}^{0.5391}X_{4}^{0.181}X_{9}^{0.1915} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{5}}=17.38{{X}_{3}}-113.3209X_{2}^{0.5987}{{X}_{5}}X_{7}^{0.2434} $ |
| $ \begin{align} & {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{6}}=0.0111X_{6}^{-0.2296}X_{12}^{-0.2459}- \\ & 17.7973X_{3}^{0.5194}{{X}_{6}}X_{8}^{0.1744}X_{9}^{0.1845}X_{11}^{-7.9048} \\ \end{align} $ |
| ${{\overset{\centerdot }{\mathop{\text{X}}}\,}_{7}}=269.8467X_{2}^{0.711}{{X}_{5}}X_{7}^{0.289}-34.78{{X}_{7}} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{8}}=386.0822X_{3}^{0.5914}{{X}_{4}}X_{8}^{0.1986}X_{9}^{0.2101}-34.78{{X}_{8}} $ |
| $ \begin{align} & {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{9}}=29.9152X_{3}^{0.5914}{{X}_{6}}X_{8}^{0.1986}X_{9}^{0.2101}X_{11}^{-9}- \\ & 69.5323{{X}_{9}}X_{13}^{0.9994} \\ \end{align} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{10}}=133.9756X_{7}^{0.2678}X_{8}^{0.3558}X_{9}^{0.3764}X_{13}^{0.3764}-0.07{{X}_{10}} $ |
首先求解模型的特征值：－29.75，－34.78，－27.70，－0.44+3.51i，－0.44－3.51i，－0.14+1.69i，－0.14－1.69i，－0.39，－0.02，－0.09。由于所有特征值的實部均為負值，這表明該系統具有局部穩定的特點；而有4個特征值的虛部非零，這說明系統可能存在振蕩的特性，這點在圖 1所示的動態過程中已有體現。
模型的穩定性分析揭示了系統在非獨立變量(狀態變量)受到擾動時能否回到先前穩態的特性；而靈敏度分析則能夠預測獨立變量(參數)在發生微小改變后，系統所達到的新穩態與原穩態的差異，而這種相對擾動量的差異程度即可表征模型的魯棒性，通常情況下，魯棒的模型具有較低的參數靈敏度。借助S系統所具有的性質，可以很方便地分析系統參量的穩態靈敏度[20-21]。
2.2.1 速率常數引起的狀態變量的相對變化
靈敏度定義為：
| $ S\left( {{X}_{i,{{\alpha }_{j}}}} \right)=\partial \left( \ln {{X}_{i}} \right)/\partial \left( {{\alpha }_{j}} \right) $ |
負向速率常數的靈敏度S(Xi,βj)具有相同的定義，可以證明S(Xi,αj)與S(Xi,βj)絕對值相同，但符號相反[20]。因此只討論正向速率常數的靈敏度即可，計算結果如圖 2所示。
圖2
速率常數的靈敏度對比
Figure2.
Sensitivities of rate constant
從圖 2可以看出，在總共100個靈敏度中，絕對值均未超過5。其中對X1、X5和X9(對應的物理量分別為S、R2和IV2)的影響最小，靈敏度小于1；影響最大的為X2和X8(I1和J2)，靈敏度在4~5之間。參考常規的生化反應過程分析經驗[18]，靈敏度絕對值在0~5范圍內，表明其影響并不顯著，系統比較穩定，但對擾動仍會響應。
2.2.2 關鍵參數引起的狀態變量的相對變化
除了速率常數，模型的可靠性與參數的選取有密不可分的關系。為了考察參數的靈敏度，取原SIR模型中的三個關鍵參數τ、p和η作為S系統中的獨立變量X11、X12和X13來分析它們的可靠性。獨立變量的靈敏度一般稱之為對數增益(logarithmic gain)，其定義如下：
| $ L\left( {{X}_{i}},{{X}_{j}} \right)=\partial \left( \ln {{X}_{i}} \right)/\partial \left( {{X}_{j}} \right) $ |
可以看出，式(8)與式(7)有類似的數學描述。計算結果以柱狀圖給出，如圖 3所示。
圖3
關鍵參數的靈敏度對比
Figure3.
Sensitivities of key parameters
從結果中可以看出p和η對狀態變量的影響有限(靈敏度均不超過2)，相對地，參數τ對系統的影響作用較為顯著，尤其是對X6(V，已接種并對病原體1免疫)的靈敏度接近8，未超過10。文獻[18]指出，若對數增益在10以上，則模型對參數的擾動會有顯著的響應，文獻[19]中對τ取兩個不同值(0.7和0.9)后的結果顯示，兩者收斂狀態未見顯著差異。雖然參數τ并未導致系統穩態的顯著變化，但相對于另外兩個參數p和η來說，其對數增益明顯過高，這在實際模型評價及預測分析中應引起足夠的重視。通過對數增益分析可知，參數的設定對系統的可靠性至關重要，這對模型的建立與應用可起到有價值的指導作用。
通過分析系統特征值、速率常數的靈敏度以及關鍵參數對系統狀態影響的對數增益，該系統是局部穩定的，參量靈敏度也在可接受范圍，系統是魯棒的；而通過分析系統存在的振蕩特性可為疫苗的供給與接種以及流行病的防治工作給予指導。需要強調的是，前文的靈敏度分析參考了常規生化反應過程的經驗，而在流行病SIR模型中的分析評價指標需要在實踐中進一步驗證。
3 結語
通過對流行病的傳播過程進行數學建模分析，可對流行病的防治工作提供科學的指導，而流行病模型的常規分析方法需要應用復雜的數學理論，這使得分析工作面臨較多困難。本文針對典型的流行病學SIR模型，利用BST對其模型的有效性和穩定性進行了探討。首先進行動態特性分析，應用開源非線性解題器IPOPT對其進行仿真計算，獲得流行病學的各項指標的預測結果；利用BST方法將常規機理模型轉換為S系統，基于此模型分析速率常數對各項指標影響的靈敏度，以及三個關鍵參數對模型狀態的影響，依此評價模型的魯棒性，可以為今后模型的建立及評價等工作提出指導，同時，BST方法的易用性也大幅度降低了模型分析工作的難度。本研究拓寬了BST理論的應用范圍，但在不同的應用領域，其分析評價指標應給予針對性的考慮。
引言
在當下的信息時代，生物學家已經意識到傳統通過實驗的研究方式面臨成本高、效率低的問題。在計算機應用領域不斷拓展的背景下，計算生物學與系統生物學得到快速發展。依據長期積累的實驗數據，利用數學建模的方法去研究并預測不同生物過程的動態特性進而有效指導實際實驗研究方向，已成為生物研究領域的關注重點。隨著研究的不斷深入，不同類型生物過程的復雜程度也同時提高，與之對應的模型評價與分析工作也越來越困難，這勢必需要新的建模手段和分析方法去應對[1-4]。
流行病是在外部環境下，于短時間內在人群中廣泛傳播并極可能面臨大面積感染的傳染病。病原體的抗原變異在流行病傳播過程中是常見的現象，對流行病的防治是一個巨大的困難。疫苗可使機體產生特異性抗體，眾所周知，通過疫苗的接種可使人體獲得對某種病原體的免疫力，比如通過種痘的方式已在全世界范圍徹底抑制了天花。但目前，仍然存在很多沒有疫苗可應對的傳染病，這對流行病防治來說是極大的挑戰。通過建立流行病的數學模型，可較準確地評價和預測流行病的傳播、發展趨勢以及評估疫苗的有效性，有利于提前試行有效的防治措施，目前在這方面已有大量的工作報道[5-6]，但多集中在模型的建立和參數辨識上[7-8]，而模型的有效性與評價方法研究相對較少，由于流行病學的實際數據需要長期的積累，因此模型的可靠性的提高能夠保證數據應用的時效性。
針對流行病防治過程的“易感-感染-恢復”(Susceptible-Infective-Recovered，SIR)模型，在前期工作基礎上[9-10]，采用生化系統理論(Biochemical Systems Theory，BST)的方法對其進行評價，對模型特性進行分析并評價其魯棒性。
1 方法
不同類型生物過程的動態模型大多采用微分方程或微分代數方程來描述，針對難以分析的復雜模型，可利用BST理論方法建立與之對應的S系統來幫助分析。現有許多工具可用來求解，除MATLAB等常規軟件外，還可選用開源免費的軟件包，如微分/代數解題包SUNDIALS[11]和非線性規劃解題器IPOPT[9, 12]等。在使用計算機對生物數學模型進行分析時，為確保相關信息的正確描述，一般需要一種通用的數據描述、儲存規范。目前系統生物學標記語言(Systems Biology Markup Language，SBML)使用最為廣泛，當前已有多種公用模型數據庫[13]以SBML為主要的儲存格式。
1.1 動態模擬
不同規模的生物過程通常描述為復雜程度不同的非線性系統，其動態特性采用前期工作中開發的通用可擴展集成規劃工具GIEPT[9]求解。當生物系統的復雜性達到一定程度后，其數學模型不能采用常規的數值積分法求解，須將其轉化為非線性規劃問題(nonlinear programming，NLP)后采用動態優化的方法進行計算。通常NLP的標準形式如下：
| $ \begin{array}{l} min\int\limits_0^{{t_f}} {\Phi \left( {{z_d}\left( t \right),{z_a}\left( t \right),u\left( t \right),\theta ,\varepsilon } \right){\rm{d}}} t\\ {\rm{s}}{\rm{.t}}{\rm{.:}}\frac{{{\rm{d}}{z_d}}}{{{\rm{d}}t}} = {f_d}\left( {{z_d}\left( t \right),{z_a}\left( t \right),u\left( t \right),\theta ,\varepsilon } \right)\\ \;\;\;\;\;\;\;\;0 = {f_a}\left( {{z_d}\left( t \right),{z_a}\left( t \right),u\left( t \right),\theta ,\varepsilon } \right)\\ \;\;\;\;\;\;\;\;\;z_d^L \le {z_d}\left( t \right) \le z_d^U\\ \;\;\;\;\;\;\;\;\;z_a^L \le {z_a}\left( t \right) \le z_a^U\\ \;\;\;\;\;\;\;\;\;{u^L} \le u\left( t \right) \le {u^U}\\ \;\;\;\;\;\;\;\;\;z\left( 0 \right) = {z_0} \end{array} $ |
其中：zd(t)、za(t)和u(t)分別表示系統的微分狀態變量、代數狀態變量以及控制變量；θ表示系統內可變參數(為降低復雜度，只考慮經簡化的初始狀態)，ε描述了系統內具有固定值的參數。選取N個離散點；采用文獻[14-16]報道的方法將模型離散化，并構造非線性規劃的標準數學描述。
1.2 BST模型
BST[17-18]是一種可對系統進行模擬分析的數學計算方法。最初，BST主要是針對生化過程設計提出的，本文將其應用到相近的流行病傳播過程的分析中。
BST是基于動力學的方法，其最為常用的數學描述稱為S系統，它是由一種冪律函數描述的。取X代表目標系統中的某一狀態量，將對X有貢獻(流入)的因素合計為正向流量(記為V+)，使X有損失(流出)的因素合計為負向流量(記為V-)，則S系統可由下面的數學式表示：
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{i}}=V_{i}^{+}-V_{i}^{-}={{\alpha }_{i}}\underset{j=1}{\overset{n+m}{\mathop \prod }}\,X_{j}^{{{g}_{i,j}}}-{{\beta }_{i}}\underset{j=1}{\overset{n+m}{\mathop \prod }}\,X_{j}^{{{h}_{i,j}}},i=1,2,\cdots n $ |
其中Xi是系統的狀態變量(亦稱之為非獨立變量)；Xj是包括狀態變量在內的所有影響系統響應速率的因素；n為狀態變量的個數，m為獨立變量的個數；αi為流入(正向)速率常數，βi為流出(負向)速率常數；gi,j和hi,j分別為各因素的動力階。速率常數和動力階分別由下式給出：
| $ {g_{i,j}} = {\left. {\frac{{\partial V_i^ + }}{{\partial {X_j}}} \cdot \frac{{{X_j}}}{{V_i^ + }}} \right|_{op}},{h_{i,j}} = {\left. {\frac{{\partial V_i^ - }}{{\partial {X_j}}} \cdot \frac{{{X_j}}}{{V_i^j}}} \right|_{op}} $ |
| $ {\alpha _i} = {\left. {\frac{{V_i^ + }}{{_{\mathop \prod \limits_{j = 1}^{n + m} X_j^{{g_{i,j}}}}}}} \right|_{op}},{\beta _i} = {\left. {\frac{{\partial V_i^ - }}{{\mathop \prod \limits_{j = 1}^{n + m} X_j^{{h_{i,j}}}}}} \right|_{op}} $ |
其中OP為工作點，一般取自系統的穩態值。
由于S系統具有其他動態系統所不具備的特性，使得借助BST方法，能夠對生物系統進行數學分析。尤其體現在，當生物系統的規模逐步拓展，采用微分方程描述的機理模型的復雜性會無限提高，但與之對應的S系統的基本結構仍會保持不變。這就使得S系統能夠應對不同規模的系統[18]。
2 應用
當前研究較多的流行病模型大多以Kermack和McKendrick提出的經典倉室模型[5]為基礎，本文研究的模型即來源于典型的SIR流行病學模型[19]。該模型描述了疫苗接種防治過程中的免疫逃避動態過程：在某一群體中，通過接種疫苗使得原有土著致病菌株(百日咳)維持較低的發病率，引入突變菌株后，假定兩種菌株具有相同繁殖率，則這兩種菌株可達到共存的相對穩定狀態。文獻中原有數學模型狀態變量(或參數)與S系統因素的對應關系如表 1所示。
參數值來源于在線數據庫BioModels[13]。其中三個參數：疫苗覆蓋率(即新生兒接種率)p，接種后易感人群數量的下降率τ以及接種后傳染周期的縮減率η是接下來模型分析的考察目標。
表1
系統原參量與S系統因素的對應關系
Table1.
Correspondence between original variables and elements of S-system
2.1 動態過程仿真
使用GIEPT調用IPOPT進行求解。在利用IPOPT求解之前，如前所述，需將動態模型進行離散化，轉化成NLP。由于原有數學模型為單純的微分方程組，并未有顯式的目標函數，因此必須構造一個恰當的目標函數形式。這里，以離散點處斜率估算誤差總和最小為目標函數，如式(5)所示。
| $ \min :\sum\limits_{i}^{n}{\sum\limits_{j}^{N}{\left\| {{\overset{\centerdot }{\mathop{x}}\,}_{i,j}}-\frac{{{x}_{i,j+1}}-{{x}_{i,j}}}{{{t}_{i,j+1}}-{{t}_{i,j}}} \right\|}} $ |
式(5)中，x代表SIR系統數學模型中的狀態變量；n為狀態變量的個數；N為離散分段數。至此，形如式(1)的標準形式已經具備，可以采用非線性規劃的方法進行求解。為了避免數據點過多產生堆積，每隔6個時刻繪制一次數據點，如圖 1所示。
圖1
SIR系統的動態過程
子圖為
sub-figure describes the dynamic process of
從結果可以看出，該過程是一個非常復雜的非線性過程，伴隨有不同程度的振蕩，圖示中各指標經過振蕩后均能達到穩定，這種具有減幅振蕩的暫態過程表明系統的特征值具有非零的虛部。考察SIR系統的主要指標S(易感比率)、I(感染比率)以及R(恢復比率)，三種指標特別是土著菌株與突變菌株的感染率(I1、I2)的變化情況符合文獻[19]中的描述：兩種菌株達到共存的穩定狀態(圖 1中子圖)。
圖 1所示的SIR系統動態過程事實上是一種理想狀態，其描述了雙菌株在長期的致病及接種疫苗防治過程中達到穩定的情形，而在實際情況中，若再出現新的突變致病菌株，則有流行病大爆發的危險；另一方面，多菌株病原體的出現，也給疫苗研究工作帶來很大困難。因此，通過模型預測流行病的防治效果及潛在威脅就顯得尤為重要，而這其中的關鍵問題是模型是否可靠。
2.2 模型穩定性與靈敏度分析
利用前述介紹的BST理論方法，將原有的SIR模型微分方程中的因式分為正向和負向部分，依照表 1的參量對應關系，將其轉化為BST模型(即S系統)，如式(6)所示。
| $ \begin{align} & {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{1}}=0.0775X_{1}^{-0.0091}X_{4}^{0.0305}X_{5}^{0.0525}X_{6}^{0.0459}X_{10}^{0.8803}X_{12}^{-0.2459}- \\ & 1298.2{{X}_{1}}X_{2}^{0.2424}X_{3}^{0.3898}X_{7}^{0.0985}X_{8}^{0.1309}X_{9}^{0.1385} \\ \end{align} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{2}}=539.6934{{X}_{1}}X_{2}^{0.711}X_{7}^{0.289}-17.4{{X}_{2}} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{3}}=772.1643{{X}_{1}}X_{3}^{0.5914}X_{8}^{0.1986}X_{9}^{0.2101}-17.4{{X}_{3}} $ |
| ${{\overset{\centerdot }{\mathop{\text{X}}}\,}_{4}}=17.38{{X}_{2}}-243.0366X_{3}^{0.5391}X_{4}^{0.181}X_{9}^{0.1915} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{5}}=17.38{{X}_{3}}-113.3209X_{2}^{0.5987}{{X}_{5}}X_{7}^{0.2434} $ |
| $ \begin{align} & {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{6}}=0.0111X_{6}^{-0.2296}X_{12}^{-0.2459}- \\ & 17.7973X_{3}^{0.5194}{{X}_{6}}X_{8}^{0.1744}X_{9}^{0.1845}X_{11}^{-7.9048} \\ \end{align} $ |
| ${{\overset{\centerdot }{\mathop{\text{X}}}\,}_{7}}=269.8467X_{2}^{0.711}{{X}_{5}}X_{7}^{0.289}-34.78{{X}_{7}} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{8}}=386.0822X_{3}^{0.5914}{{X}_{4}}X_{8}^{0.1986}X_{9}^{0.2101}-34.78{{X}_{8}} $ |
| $ \begin{align} & {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{9}}=29.9152X_{3}^{0.5914}{{X}_{6}}X_{8}^{0.1986}X_{9}^{0.2101}X_{11}^{-9}- \\ & 69.5323{{X}_{9}}X_{13}^{0.9994} \\ \end{align} $ |
| $ {{\overset{\centerdot }{\mathop{\text{X}}}\,}_{10}}=133.9756X_{7}^{0.2678}X_{8}^{0.3558}X_{9}^{0.3764}X_{13}^{0.3764}-0.07{{X}_{10}} $ |
首先求解模型的特征值：－29.75，－34.78，－27.70，－0.44+3.51i，－0.44－3.51i，－0.14+1.69i，－0.14－1.69i，－0.39，－0.02，－0.09。由于所有特征值的實部均為負值，這表明該系統具有局部穩定的特點；而有4個特征值的虛部非零，這說明系統可能存在振蕩的特性，這點在圖 1所示的動態過程中已有體現。
模型的穩定性分析揭示了系統在非獨立變量(狀態變量)受到擾動時能否回到先前穩態的特性；而靈敏度分析則能夠預測獨立變量(參數)在發生微小改變后，系統所達到的新穩態與原穩態的差異，而這種相對擾動量的差異程度即可表征模型的魯棒性，通常情況下，魯棒的模型具有較低的參數靈敏度。借助S系統所具有的性質，可以很方便地分析系統參量的穩態靈敏度[20-21]。
2.2.1 速率常數引起的狀態變量的相對變化
靈敏度定義為：
| $ S\left( {{X}_{i,{{\alpha }_{j}}}} \right)=\partial \left( \ln {{X}_{i}} \right)/\partial \left( {{\alpha }_{j}} \right) $ |
負向速率常數的靈敏度S(Xi,βj)具有相同的定義，可以證明S(Xi,αj)與S(Xi,βj)絕對值相同，但符號相反[20]。因此只討論正向速率常數的靈敏度即可，計算結果如圖 2所示。
圖2
速率常數的靈敏度對比
Figure2.
Sensitivities of rate constant
從圖 2可以看出，在總共100個靈敏度中，絕對值均未超過5。其中對X1、X5和X9(對應的物理量分別為S、R2和IV2)的影響最小，靈敏度小于1；影響最大的為X2和X8(I1和J2)，靈敏度在4~5之間。參考常規的生化反應過程分析經驗[18]，靈敏度絕對值在0~5范圍內，表明其影響并不顯著，系統比較穩定，但對擾動仍會響應。
2.2.2 關鍵參數引起的狀態變量的相對變化
除了速率常數，模型的可靠性與參數的選取有密不可分的關系。為了考察參數的靈敏度，取原SIR模型中的三個關鍵參數τ、p和η作為S系統中的獨立變量X11、X12和X13來分析它們的可靠性。獨立變量的靈敏度一般稱之為對數增益(logarithmic gain)，其定義如下：
| $ L\left( {{X}_{i}},{{X}_{j}} \right)=\partial \left( \ln {{X}_{i}} \right)/\partial \left( {{X}_{j}} \right) $ |
可以看出，式(8)與式(7)有類似的數學描述。計算結果以柱狀圖給出，如圖 3所示。
圖3
關鍵參數的靈敏度對比
Figure3.
Sensitivities of key parameters
從結果中可以看出p和η對狀態變量的影響有限(靈敏度均不超過2)，相對地，參數τ對系統的影響作用較為顯著，尤其是對X6(V，已接種并對病原體1免疫)的靈敏度接近8，未超過10。文獻[18]指出，若對數增益在10以上，則模型對參數的擾動會有顯著的響應，文獻[19]中對τ取兩個不同值(0.7和0.9)后的結果顯示，兩者收斂狀態未見顯著差異。雖然參數τ并未導致系統穩態的顯著變化，但相對于另外兩個參數p和η來說，其對數增益明顯過高，這在實際模型評價及預測分析中應引起足夠的重視。通過對數增益分析可知，參數的設定對系統的可靠性至關重要，這對模型的建立與應用可起到有價值的指導作用。
通過分析系統特征值、速率常數的靈敏度以及關鍵參數對系統狀態影響的對數增益，該系統是局部穩定的，參量靈敏度也在可接受范圍，系統是魯棒的；而通過分析系統存在的振蕩特性可為疫苗的供給與接種以及流行病的防治工作給予指導。需要強調的是，前文的靈敏度分析參考了常規生化反應過程的經驗，而在流行病SIR模型中的分析評價指標需要在實踐中進一步驗證。
3 結語
通過對流行病的傳播過程進行數學建模分析，可對流行病的防治工作提供科學的指導，而流行病模型的常規分析方法需要應用復雜的數學理論，這使得分析工作面臨較多困難。本文針對典型的流行病學SIR模型，利用BST對其模型的有效性和穩定性進行了探討。首先進行動態特性分析，應用開源非線性解題器IPOPT對其進行仿真計算，獲得流行病學的各項指標的預測結果；利用BST方法將常規機理模型轉換為S系統，基于此模型分析速率常數對各項指標影響的靈敏度，以及三個關鍵參數對模型狀態的影響，依此評價模型的魯棒性，可以為今后模型的建立及評價等工作提出指導，同時，BST方法的易用性也大幅度降低了模型分析工作的難度。本研究拓寬了BST理論的應用范圍，但在不同的應用領域，其分析評價指標應給予針對性的考慮。
