基于復雜網絡的新冠病毒群體免疫數值仿真

王佳亮 李海濱 李海燕

摘要：鑒于構建流行病動力學模型、探索流行病傳播規律對疫情防控具有十分重要的理論意義和實際應用價值，在已有的均勻混合模型基礎上，針對個體接觸關系異質化越發明顯，且每個個體都處在不同的接觸關系中，建立了兼顧個體狀態與接觸追蹤的動態小世界網絡模型。模擬了新冠病毒在社會中的傳播過程。通過對比仿真結果，說明了所建模型的合理性。在此基礎上，仿真計算了網絡拓撲結構與接種免疫人數占比共同作用下對新冠病毒傳播的影響，分析得到群體免疫臨界值。說明所建傳播模型合理，接種疫苗實現群體免疫可行。

關鍵詞：小世界網絡；動力學建模；新冠病毒；群體免疫

中圖分類號： TP391.9文獻標識碼： A

收稿日期：2021-09-06；修回日期：2021-12-31

基金項目：國家自然科學基金（11962021）

第一作者：王佳亮（1996-），男，內蒙古赤峰人，碩士研究生，主要研究方向為系統動力學建模與仿真。

通信作者：李海濱（1973-），男，內蒙古呼和浩特人，博士，教授，主要研究方向為結構不確定性分析與量化、神經網絡計算、六維力傳感器設計。

李海燕（1980-），女，內蒙古呼和浩特人，本科，副主任護師，主要研究方向為內分泌疾病與糖足的臨床護理。

Numerical Simulation of the COVID-19 Herd Immunity Based on Complex Network Modeling

WANG Jialiang^1a， LI Haibin¹， LI Haiyan²

（1. a. College of Sciences; b. Engineering Training Center of Inner Mongolia University of Technology， Hohhot 010051， China; 2. Department of Endocrinology the First Affiliated Hospital of Inner Mongolia Medical University， Hohhot 010010， China）

Abstract：Constructing an epidemic dynamic model and exploring the spreading law of epidemic have very important theoretical significance for epidemic prevention and control. Based on the existing homogeneous mixing model， in view of the increasingly obvious heterogeneity of individual contact relationships， and each individual is in a different contact relationship， a dynamic small-world network model that takes into account individual status. Contact tracking has been established to simulate the spread of the COVID-19 in society. By comparing the simulation results， the rationality of the built model is explained. On this basis， the simulation calculated the impact of the network topology and the proportion of vaccinated people on the spread of the COVID-19， analyzed the critical value of herd immunity. The established propagation model is reasonable， and feasible to achieve herd immunization by vaccination.

Key words： small world network; dynamic modeling; COVID-19; herd immunity

0 引言

自從人類社會誕生以后，傳染病一直是困擾人類的難題之一，從早期的黑死病、天花、麻疹，到近期的埃博拉（EBOV）、寨卡病毒（Zika Virus）和新型冠狀病毒（COVID-19）^[1]。隨著航空、高鐵等交通網絡的飛速發展，便捷的交通為人類出行帶來方便，同時也為流行病迅速傳播提供了條件。建立流行病動力學模型、探索病毒傳播規律、提供有效預防措施是疫情防控的重要研究內容。

經典的均勻混合動力學模型由Kermack和McKendrick于1926年提出，隨后基于此工作，出現了SEIR（Susceptible Exposed Infected Removed）、SEIQR（Susceptible Exposed Infected Quarantined Removed）等一系列擴展模型用來研究流行病傳播，并取得了諸多重要結論和成果。針對新冠疫情，Fang等^[2]應用SEIR模型進行仿真模擬和數據擬合，預測了武漢疫情確診病例數的峰值，仿真結果表明所建模型具有較高的擬合精度。范如國等^[3]利用SEIR模型，模擬了5，7，10天3種不同潛伏期下的武漢疫情，預計此次疫情拐點在第71，74，78天出現。鐘南山院士團隊^[4]利用優化的SEIR模型和人工智能方法預測國內疫情將在2月底達到高峰，4月底趨于平緩。如果封城等管控措施推遲5天實施，中國國內的疫情規模將擴大3倍。Gu等^[5]利用SEIR模型研究了COVID-19的傳播速度、空間范圍和動力學機制，并對國內疫情的擴散趨勢進行了預測。Mwalili等^[6]利用改進的SEIR模型，探索了新冠疫情在不采取防控措施下的蔓延情況，結果表明，若不加以防控，新冠疫情將在世界范圍內持續擴散，從而產生毀滅性影響。的確經典的均勻混合動力學模型已經為流行病的防控做出了應有的貢獻，但是在個體間接觸關系異質化越來越明顯的情況下，均勻混合動力學模型已經無法描述目標人群內部復雜的動力學性質。基于網絡建立的流行病動力學模型能夠很好反映個體間接觸關系，從而出現了利用復雜網絡研究流行病傳播的動力學模型。楊洪勇等^[7]利用小世界網絡模型模擬了禽流感病毒的傳播過程，通過與實際數據對比來說明所構建模型的有效性。程靜等^[8]用小世界網絡模型對埃博拉病毒的傳播情況進行了仿真，仿真結果與實際報道的傳播情況相符。劉漢卿等^[9]利用小世界網絡模型研究了COVID-19在其網絡中的傳播規律以及通過控制節點行為對其傳播的影響，數據仿真結果驗證了模型的適用性。上述工作表明基于小世界網絡的流行病動力學模型不僅能描述個體狀態轉化關系，還能夠描述個體間的有效接觸關系。因此，本文嘗試采用小世界網絡構建流行病傳播模型，以更好模擬新冠病毒在現實中的接觸傳播過程。

目前，中國的防控措施已從積極發現并隔離病例、追蹤隔離密切接觸者、控制人群流動，發展到針對局部疫情暴發迅速采取多輪大規模人群病毒核酸檢測，形成了行之有效的防控策略^[10]。然而，非藥物干預措施下的疫情防控仍然面臨著外來輸入的風險，經過科研人員的不懈努力，中國已經成功研制出針對新型冠狀病毒的疫苗，中國國家藥品監督管理局于2020年12月30日批準了國藥集團中國生物北京生物制品研究所新冠病毒滅活疫苗的注冊申請^[11]。在中國新冠疫情防控進入疫苗時代，通過構建群體免疫屏障來戰勝新冠病毒成為可能。所謂群體免疫，就是在人群中讓大多數人對某種傳染病產生免疫力，在社會層面形成一道保護屏障，當一個或者多個傳染病患者進入這樣的人群時，也不會出現傳染病暴發的情況^[12]。那么，接種疫苗人群占比多少能夠達到群體免疫效果是一個很有研究意義的問題。針對這一問題，吳丹等^[13]根據新冠病毒的基本再生數R0=1.6～6.5，通過基本再生數與群體免疫臨界值的關系式計算得到群體免疫臨界值為38%～85%。目前給出的群體免疫臨界值都是根據基本再生數通過關系式計算得到，屬于經驗估計，缺乏必要的實驗或仿真驗證。因此，本文將在動態小世界網絡中仿真計算免疫個體占比對新冠病毒傳播的影響，根據仿真結果給出群體免疫臨界值。

1 動態小世界網絡上的流行病動力學建模與仿真

在動態小世界網絡中研究新冠病毒擴散動力學過程，包括：在具有復雜拓撲結構的網絡中研究新冠病毒的傳播，以及由于新冠病毒傳播對網絡中不同狀態節點的影響。不僅需要考慮網絡的拓撲結構，還需要確定新冠病毒的傳播模式。本節將用小世界網絡模型模擬新冠病毒傳播過程，并與已有均勻混合模型做對比，來說明所建模型的合理性。

1.1 小世界網絡上的新冠病毒傳播模型

小世界網絡模型最早由Watts和Strogatz提出，具有N個節點的小世界網絡模型構造規則為^[14]：1）首先將N個節點排列成一個圓并順時針按序編號，每個節點與其左右最近鄰的各K/2個節點連線，生成一個具有N·K/2條邊的最近鄰規則網絡；2）將步驟1）生成的規則網絡中的每一條邊以概率p進行斷線重連處理，重連過程中，保證網絡沒有重復邊和自連接的情況。當p=0時，對應規則網絡；當p=1時，對應隨機網絡；當0

流行病在傳播過程中，個體的狀態可分為易感態S、潛伏態E、有癥狀的患病態I、無癥狀的患病態A、隔離態Q、治愈態R和死亡態D等幾種基本狀態。其中，除了易感態S和有癥狀的患病態I是必選之外，其它狀態要根據實際需要加以選擇^[15]。根據新冠病毒傳播特性，本文將目標人群中的個體狀態分為5類：易感態S（未患病但可以被感染）；非傳染潛伏態E1（無癥狀患病態但不具有傳染性）；傳染潛伏態E2（無癥狀患病態且具有傳染性）；患病態I（有癥狀患病態且具有傳染性）；隔離態Q（被隔離的有癥狀患病態）。根據新冠病毒傳播規律，易感態個體被接觸感染后變為潛伏態，文獻[16]表明，潛伏態個體在潛伏期初期不具有傳染性，到了潛伏期后期開始具有傳染性。因此，將潛伏態個體分為非傳染潛伏態和傳染潛伏態2個相互獨立的狀態。潛伏態個體經過潛伏期后變為患病態，之后以一定概率被隔離。至此，完成了對任一易感態個體從感染到隔離全過程的運動軌跡跟蹤描述。所以，本文確定網絡中節點狀態轉化過程如下：當網絡中的節點處于易感態時，與其相連的每一個傳染潛伏態節點以概率β_E對其進行感染，與其相連的每一個患病態節點以概率β_I對其進行感染，同時需要判斷該節點是否被感染為非傳染潛伏態；當網絡中的節點處于非傳染潛伏態時，該節點經過非傳染潛伏期τ_E1后轉化為傳染潛伏態；當網絡中的節點處于傳染潛伏態時，該節點經過傳染潛伏期τ_E2后轉化為患病態；當網絡中的節點處于患病態時，該節點將以概率γ_Q被隔離轉化為隔離態。上述節點狀態轉化過程如圖1所示，對應動力學方程如式（1）所示。

綜合上述分析，將小世界網絡上的新冠病毒傳播過程總結為：1）初始化，生成一個具有N·K/2條邊的最近鄰規則網絡。2）小世界網絡，將規則網絡中的每一條邊以概率p進行斷線重連處理，重連過程中，保證網絡沒有重復邊和自連接的情況。3）引入傳染源，在網絡中隨機選取比例為E0的非傳染潛伏態個體作為初始傳染源，其余為易感態個體。4）開始傳染，網絡拓撲結構變化與節點狀態轉化同步進行，每同步進行一次為一個迭代步，周期為1 d。5）記錄數據，統計每一個迭代步所有節點的狀態，在觀測時間點上統計數據并作圖。

1.2 小世界網絡上新冠病毒的傳播仿真

本節將通過與文獻[17]中的時滯均勻混合模型仿真對比來驗證所構建模型的合理性。在仿真之前首先要確定模型參數。

待定的模型參數有β_E，β_I，τ_E1，τ_E2，γ_Q。其中，接觸感染率β_E，β_I和潛伏期時間常數τ_E1，τ_E2等參數反映的是病毒傳播特性。文獻[15]根據武漢某醫院新冠患者接觸治療情況給出了一次接觸感染率β_E，β_I；文獻[18]給出了新冠病毒潛伏期參數估計的研究報道，潛伏期中位數為5.5 d，95%患者的潛伏期不超過13 d。據此本文設定新冠病毒的潛伏期為12～13 d，大致分為兩個階段，即非傳染潛伏期τ_E1為5～6 d、傳染潛伏期τ_E2為7～8 d。將上述各參數列于表1。隔離率γ_Q的設定取決于對疫情防控的響應速度。

在上述參數條件下，分別利用小世界網絡模型和文獻[17]中的時滯均勻混合模型進行疫情動力學仿真。取網絡節點總數N=5 000，每個節點初始擁有的鄰居節點個數K=28，斷線重連率p=0.05，非傳染潛伏態個體占比E0=0.01，其他參數見表1，以180 d為觀測終點進行仿真。在上述各態不考慮隔離措施情況下，圖2中的“—”模擬了新冠疫情在動態小世界網絡中傳播的全過程；圖2中的“○”為同參數條件下文獻[17]中的時滯均勻混合動力學模型仿真所得結果。

對比圖2所得結果可以看到，兩種方法在疫情傳播過程中各態節點變化趨勢高度吻合，可以說明本文所建模型的合理性。由于文獻[17]中的均勻混合傳播模型是以時滯微分方程形式描述的，模型給定初始條件后，應用MATLAB中的“dde23”可以求出確定的解；且均勻混合傳播模型對應“全連通網絡”，表示為任意一個患病態個體可以把病毒傳染給網絡中任意一個易感態個體。本文所構建的動態小世界網絡傳播模型，既考慮了新冠病毒傳播過程中的隨機性，又根據網絡拓撲結構來說明實際的患病態個體只把病毒傳染給了所接觸的有限個體，從而造成了仿真結果的差異。基于網絡建立的傳播模型，既涵蓋了病毒的傳播特性，又包含了自身結構的特性，可以使人們對新冠病毒接觸傳播的過程認識更加清晰。

2 基于網絡的群體免疫數值仿真

以最小的代價去控制流行病的傳播，將流行病擴散造成的影響和損失降到最低，可以通過對人群中大部分人進行免疫來實現，這樣不僅接種疫苗的免疫者能夠得到保護，全人群也可能獲得保護。因此，本節將在動態小世界網絡中討論免疫個體占比對新冠病毒傳播的影響。主要包括兩方面工作：1）仿真對比有、無免疫個體兩種情況下，新冠病毒在動態小世界網絡中的傳播過程，通過仿真結果來說明，由于免疫個體的存在對病毒傳播造成的影響。2）在網絡拓撲結構與初始免疫人數占比共同作用下，在不同觀測節點上仿真計算了新冠病毒傳播的最終感染人數占比；根據感染人數占比計算有效再生系數，分析給出了群體免疫臨界值。

首先通過仿真來說明免疫個體的存在對新冠病毒傳播的影響。取網絡節點總數N=5 000，每個節點初始擁有的鄰居節點個數K=24，斷線重連率p=0.05，非傳染潛伏態個體占比E0=0.01，初始免疫個體占比R=0.2，其他參數見表1，以90 d為觀測節點進行仿真。圖3、圖4分別模擬了有、無免疫個體情況下新冠病毒的傳播過程。

分析圖3與圖4所得結果，在90 d觀測節點上，無免疫個體時最終感染人數占比約為51.68%，有免疫個體時最終感染人數占比約為17.66%，其中感染人數占比為非傳染潛伏態、傳染潛伏態、患病態個體占比的總和。通過仿真結果可以看出，免疫個體的存在能夠大幅減少病毒擴散。

在網絡中討論接種免疫人數占比對病毒傳播的影響，需要考慮網絡的拓撲結構，以及初始免疫人數的占比。假設表征網絡拓撲結構的參數用平均路徑長度L來描述，平均路徑長度是指網絡中連接任意兩個節點所要經過的最少連邊個數的平均值^[19]。由于圖3所得結果僅選擇了一組參數，下面將對模型參數取不同的值分別進行仿真計算，以探索表征網絡拓撲結構的參數L與初始免疫人數占比R共同作用下，對網絡中新冠病毒傳播最終感染人數G占比的影響。取網絡節點總數N=5 000，非傳染潛伏態個體占比E0=0.01保持不變，為了得到不同觀測節點上感染人數G的占比，將在隨機變量L和R組成的各離散點上，分別以30 d，60 d，90 d為觀測節點進行仿真計算，所得結果如圖5所示。

分析圖5所得結果，當網絡的平均路徑長度L較小時，病毒在網絡中的傳播速度較快；隨著免疫人數占比的逐漸增加，病毒在網絡中的傳播速度減慢，最終感染人數占比減少。下面將根據感染人數占比計算有效再生系數，分析給出不同觀測節點上的群體免疫臨界值。

通過式（2）^[20]：

計算有效再生系數，其中，total（t）為第t天的感染人數占比（見圖5），t為傳播天數，取值為30，60，90；n為初始的感染人數占比，取值為0.01；k為有效再生系數；將圖5中所得的每種情況的total（t）分別對有效再生系數k進行計算，所得結果如圖6所示。

綜合分析圖5、圖6所得結果，在30 d觀測節點上，當免疫人數占比達到35%左右時，感染人數占比最多為1.88%，有效再生系數k=0.468 08，小于0.5，可以認為達到群體免疫效果。在60 d觀測節點上，當免疫人數占比達到65%左右時，感染人數占比最多為1.81%，有效再生系數k=0.447 51，小于0.5，可以認為達到群體免疫效果。在90 d觀測節點上，當免疫人數占比達到70%左右時，感染人數占比最多為1.84%，有效再生系數k=0.456 52，小于0.5，可以認為達到群體免疫效果。

3 總結

針對流行病動力學建模，本文建立了兼顧個體狀態與接觸追蹤的動態小世界網絡上的新冠病毒傳播模型，與均勻混合動力學模型相比，復雜網絡上的新冠病毒傳播模型能夠更好地反映個體接觸的傳染情況，對病毒傳播過程的追蹤更容易。針對疫情防控策略的群體免疫，首先仿真計算了有、無免疫個體對網絡中新冠病毒傳播的影響。仿真結果表明，免疫個體的存在能夠大幅減少病毒擴散。之后對模型參數取不同的值進行仿真計算，定量探索了初始免疫人數占比與網絡拓撲結構共同作用對新冠病毒傳播的影響。根據感染人數占比通過式（2）計算有效再生系數，綜合分析得到30 d，60 d，90 d觀測節點上的群體免疫臨界值分別為35%，65%，70%，即當初始免疫人數占比分別達到35%，65%，70%時，疫情在網絡中傳播30 d，60 d，90 d也不會出現暴發的情況。隨著觀測時間向后推移，感染人數占比呈明顯上升趨勢，可知疫情長期防控比短期防控難度更大，但是隨著疫苗的優化與大規模接種，群體免疫屏障會越來越堅固，人類終將戰勝新冠病毒。

參考文獻：

[1]SARKAR A， LIU G， JIN Y， et al. Public health preparedness and responses to the coronavirus disease 2019（COVID-19） pandemic in South Asia： a situation and policy analysis[J]. Glob Health J， 2020， 4（4）： 121-132.

[2]FANG Y Q， NIE Y T， Penny M. Transmission dynamics of the COVID-19 outbreak and effectiveness of government interventions： a data-driven analysis[J]. Journal of medical virology， 2020， 92（6）： 645-659.

[3]范如國， 王奕博， 羅明， 等. 基于SEIR的新冠肺炎傳播模型及拐點預測分析[J]. 電子科技大學學報， 2020， 49（3）： 369-374.

FAN R G， WANG Y B， LUO M， et al. SEIR-based COVID-19 transmission model and inflection point prediction analysis[J]. Journal of University of Electronic Science and Technology of China， 2020， 49（3）： 369-374.

[4]YANG Z， ZENG Z， WANG K， et al. Modified SEIR and AI prediction of the epidemics trend of COVID-19 in China under public health interventions[J]. Journal of Thoracic Disease， 2020， 12（3）： 165-174.

[5]GU B R. Forecast and analysis of COVID-19 epidemic based on improved SEIR model[J]. Journal of Physics：Conference Series， 2021， 1802（4）： 042050.

[6]MWALILI S， KIMATHI M， OJIAMBO V， et al. SEIR model for COVID-19 dynamics incorporating the environment and social distancing[J]. BMC Research Notes， 2020， 13（1）： 1-5.

[7]楊洪勇， 張嗣瀛. 基于復雜網絡的禽流感病毒傳播[J]. 系統仿真學報， 2008， 20（18）： 51-55.

YANG H Y， ZHANG S Y. Viruses epidemics of avian influenza based on complex networks[J]. Journal of System Simulation， 2008， 20（18）： 51-55.

[8]程靜， 黃青， 謝銘杰， 等. 小世界網絡中埃博拉病毒傳播的研究[J]. 生物醫學工程學進展， 2015， 36（2）： 91-94. CHENG J， HUANG Q， XIE M J， et al. Study on the transmission of ebola virus in small world network[J]. Advances in Biomedical Engineering， 2015， 36（2）： 91-94.

[9]劉漢卿， 康曉東， 高萬春， 等. 基于多模型的COVID-19傳播研究[J]. 計算機科學， 2021， 48（6）： 196-202.

LIU H Q， KANG X D， GAO W C， et al. Research on propagation of COVID-19 based on multiple models[J]. Computer Science， 2021， 48（6）： 196-202.

[10] LI Z， CHEN Q， FENG L， et al. Active case finding with case management： the key to tackling the COVID-19 pandemic[J]. The Lancet， 2020， 396（10243）： 63-70.

[11] 國藥集團. 國藥集團中國生物新冠滅活疫苗獲批附條件上市[EB/OL]. [2021-03-02]. http：//www.sinopharm.com/s/1223-4126-38840.html.

SINOPHARM. Sinopharm China biotech covid-19 inactivated vaccine approved for conditional listing[EB/OL]. [2021-03-02]. http：//www.sinopharm.com/s/1223-4126-38840.html.

[12] FINE P， EAMES K， HEYMANND D L. "Herd immunity" ： a rough guide[J]. Clinical Infectious Diseases， 2011， 52（7）： 911-916.

[13] 吳丹， 鄭徽， 李藝星， 等. 群體免疫及其對傳染病防控的意義[J]. 中國疫苗和免疫， 2020， 26（4）： 123-127.

WU D， ZHENG H， LI Y X， et al. Herd immunity and its importance in infectious disease prevention and control[J]. Chinese Journal of Vaccines and Immunization， 2020， 26（4）： 123-127.

[14] WATTS D J， STROGATZ S H. Collective dynamics of small-world networks[J]. Nature， 1998， 393（6684）： 440-442.

[15] 李海濱. 基于社會分工的流行病動力學建模與仿真研究[J]. 系統仿真學報， 2020， 32（5）： 745-758.

LI H B. Modeling and simulation on dynamics of epidemic disease based on social division of labor[J]. Journal of System Simulation， 2020， 32（5）： 745-758.

[16] XI H， ERIC H， LAU Y， et al. Temporal dynamics in viral shedding and transmissibility of COVID-19[J]. Nature Medicine， 2020， 26（5）： 672-675.

[17] 李海濱， 王佳亮， 李海燕. 新冠后疫情時代復學風險評估的不確定性量化分析[J]. 系統仿真學報， 2021， 33（1）： 13-23.

LI H B， WANG J L， LI H Y.Uncertainty quantitative analysis in risk assessment of returning to school in the post-COVID-19 era[J]. Journal of System Simulation， 2021， 33（1）： 13-23.

[18] 邱明， 悅胡濤， 崔恒建. 雙區間刪失下新冠病毒肺炎潛伏期分布的參數估計[J]. 應用數學學報， 2020， 43（2）： 200-210.

QIU M， YUE H T， CUI H J. Parameter estimation of incubation period distribution of new coronavirus pneumonia under double interval censoring[J]. Journal of Applied Mathematics， 2020， 43（2）： 200-210.

[19] 汪小帆. 復雜網絡理論及其應用[M]. 北京：清華大學出版社， 2006.

[20] 王志心， 劉治， 劉兆軍. 基于機器學習的新型冠狀病毒（COVID-19）疫情分析及預測[J]. 生物醫學工程研究， 2020， 39（1）： 1-5.

WANG Z X， LIU Z， LIU Z J. COVID-19 analysis and forecast based on machine learning[J]. Journal of Biomedical Engineering Research， 2020， 39（1）： 1-5.

（責任編輯 李 進）