基本的傳染病模型：SI、SIS、SIR及其Python代碼實現

本文主要參考博客：http://chengjunwang.com/en/2013/08/learn-basic-epidemic-models-with-python/。該博客有一些筆誤，並且有些地方表述不准確，推薦大家閱讀Albert-Laszlo Barabasi寫得書Network Science，大家可以在如下網站直接閱讀傳染病模型這一章：http://barabasi.com/networksciencebook/chapter/10#contact-networks。Barabasi是一位復雜網絡科學領域非常厲害的學者，大家也可以在他的官網上查看作者的一些相關工作。

　　下面我就直接把SIS模型和SIR模型的代碼放上來一起學習一下。我的Python版本是3.6.1，使用的IDE是Anaconda3。Anaconda3這個IDE我個人推薦使用，用起來很方便，而且提供了Jupyther Notebook這個很好的交互工具，大家可以嘗試一下，可在官網下載：https://www.continuum.io/downloads/。

     在Barabasi寫得書中，有兩個Hypothesis：1,Compartmentalization; 2, Homogenous Mixing。具體看教材。

默認條件：1, closed population; 2, no births; 3, no deaths; 4, no migrations.

　　1.    SI model 

# -*- coding: utf-8 -*-

import scipy.integrate as spi
import numpy as np
import pylab as pl

beta=1.4247
"""the likelihood that the disease will be transmitted from an infected to a susceptible
individual in a unit time is β"""
gamma=0
#gamma is the recovery rate and in SI model, gamma equals zero
I0=1e-6
#I0 is the initial fraction of infected individuals
ND=70
#ND is the total time step
TS=1.0
INPUT = (1.0-I0, I0)

def diff_eqs(INP,t):
    '''The main set of equations'''
    Y=np.zeros((2))
    V = INP
    Y[0] = - beta * V[0] * V[1] + gamma * V[1]
    Y[1] = beta * V[0] * V[1] - gamma * V[1]
    return Y   # For odeint

t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)
"""RES is the result of fraction of susceptibles and infectious individuals at each time step respectively"""
print(RES)

#Ploting
pl.plot(RES[:,0], '-bs', label='Susceptibles')
pl.plot(RES[:,1], '-ro', label='Infectious')
pl.legend(loc=0)
pl.title('SI epidemic without births or deaths')
pl.xlabel('Time')
pl.ylabel('Susceptibles and Infectious')
pl.savefig('2.5-SI-high.png', dpi=900) # This does increase the resolution.
pl.show()

結果如下圖所示

　　在早期，受感染個體的比例呈指數增長， 最終這個封閉群體中的每個人都會被感染，大概在t=16時，群體中所有個體都被感染了。

　　2.    SIS model 

# -*- coding: utf-8 -*-

import scipy.integrate as spi
import numpy as np
import pylab as pl

beta=1.4247
gamma=0.14286
I0=1e-6
ND=70
TS=1.0
INPUT = (1.0-I0, I0)

def diff_eqs(INP,t):
    '''The main set of equations'''
    Y=np.zeros((2))
    V = INP
    Y[0] = - beta * V[0] * V[1] + gamma * V[1]
    Y[1] = beta * V[0] * V[1] - gamma * V[1]
    return Y   # For odeint

t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)

print(RES)

#Ploting
pl.plot(RES[:,0], '-bs', label='Susceptibles')
pl.plot(RES[:,1], '-ro', label='Infectious')
pl.legend(loc=0)
pl.title('SIS epidemic without births or deaths')
pl.xlabel('Time')
pl.ylabel('Susceptibles and Infectious')
pl.savefig('2.5-SIS-high.png', dpi=900) # This does increase the resolution.
pl.show()

運行之后得到結果如下圖：

 

 

 

　　由於個體被感染后可以恢復，所以在一個大的時間步，上圖是t=17，系統達到一個穩態，其中感染個體的比例是恆定的。因此，在穩定狀態下，只有有限部分的個體被感染，此時並不意味着感染消失了，而是此時在任意一個時間點，被感染的個體數量和恢復的個體數量達到一個動態平衡，雙方比例保持不變。請注意，對於較大的恢復率gamma，感染個體的數量呈指數下降，最終疾病消失，即此時康復的速度高於感染的速度，故根據恢復率gamma的大小，最終可能有兩種可能的結果。

　　3.    SIR model 

# -*- coding: utf-8 -*-

import scipy.integrate as spi
import numpy as np
import pylab as pl

beta=1.4247
gamma=0.14286
TS=1.0
ND=70.0
S0=1-1e-6
I0=1e-6
INPUT = (S0, I0, 0.0)

def diff_eqs(INP,t):
	'''The main set of equations'''
	Y=np.zeros((3))
	V = INP
	Y[0] = - beta * V[0] * V[1]
	Y[1] = beta * V[0] * V[1] - gamma * V[1]
	Y[2] = gamma * V[1]
	return Y   # For odeint

t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)

print(RES)

#Ploting
pl.plot(RES[:,0], '-bs', label='Susceptibles')  # I change -g to g--  # RES[:,0], '-g',
pl.plot(RES[:,2], '-g^', label='Recovereds')  # RES[:,2], '-k',
pl.plot(RES[:,1], '-ro', label='Infectious')
pl.legend(loc=0)
pl.title('SIR epidemic without births or deaths')
pl.xlabel('Time')
pl.ylabel('Susceptibles, Recovereds, and Infectious')
pl.savefig('2.1-SIR-high.png', dpi=900) # This does, too
pl.show()

所得結果如下圖：