This article mainly refer to the blog: http://chengjunwang.com/en/2013/08/learn-basic-epidemic-models-with-python/. The blog has some clerical errors, and some places are not accurate, recommended that you read Albert-laszlo Barabasi write a book NetworkScience, you can read the following website directly to the infectious Disease Model chapter:/HTTP Barabasi.com/networksciencebook/chapter/10#contact-networks. Barabasi is a very powerful scholar in the field of complex network science, and you can also view some of the author's work on his website.
Below I will directly put the SIS model and the Code of the Sir model to study together. My Python version is 3.6.1, and the IDE used is ANACONDA3. Anaconda3 This IDE I personally recommend the use, it is very convenient, and provides jupyther notebook this very good interactive tool, you can try it, can be downloaded on the official website: https://www.continuum.io/downloads/.
There are two hypothesis:1,compartmentalization in the book written by Barabasi; 2, homogenous Mixing. Read the textbook in detail.
Default condition: 1, closed population; 2, no births; 3, no deaths; 4, no migrations.
1. SI model
1 #-*-coding:utf-8-*-2 3 ImportScipy.integrate as SPI4 ImportNumPy as NP5 ImportPylab as PL6 7beta=1.42478 """The likelihood that the disease would be transmitted from an infected to a susceptible9 Individual in a unit time isβ"""TenGamma=0 One #Gamma is the recovery rate and in SI model, gamma equals zero AI0=1e-6 - #I0 is the initial fraction of infected individuals -Nd=70 the #ND is the total time step -ts=1.0 -INPUT = (1.0-I0, I0) - + defDiff_eqs (inp,t): - " "The main set of equations" " +Y=np.zeros ((2)) AV =INP atY[0] =-beta * v[0] * v[1] + gamma * v[1] -Y[1] = Beta * v[0] * v[1]-gamma * v[1] - returnY#For Odeint - -T_start = 0.0; T_end = ND; T_inc =TS -T_range = Np.arange (T_start, t_end+t_inc, T_inc) inRES =spi.odeint (Diff_eqs,input,t_range) - """RES is the result of fraction of susceptibles and infectious individuals at each time step respectively""" to Print(RES) + - #ploting thePl.plot (res[:,0],'-bs', label='susceptibles') *Pl.plot (res[:,1],'-ro', label='Infectious') $Pl.legend (loc=0)Panax NotoginsengPl.title ('SI epidemic without births or deaths') -Pl.xlabel (' Time') thePl.ylabel ('susceptibles and Infectious') +Pl.savefig ('2.5-si-high.png', dpi=900)#This does increase the resolution. APl.show ()
The result is as shown
In the early days, the proportion of infected individuals grew exponentially, and eventually everyone in the closed group was infected, presumably at T=16, when all the individuals in the group were infected.
2. SIS model
1 #-*-coding:utf-8-*-2 3 ImportScipy.integrate as SPI4 ImportNumPy as NP5 ImportPylab as PL6 7beta=1.42478gamma=0.142869I0=1e-6TenNd=70 Onets=1.0 AINPUT = (1.0-I0, I0) - - defDiff_eqs (inp,t): the " "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] - returnY#For Odeint + AT_start = 0.0; T_end = ND; T_inc =TS atT_range = Np.arange (T_start, t_end+t_inc, T_inc) -RES =spi.odeint (Diff_eqs,input,t_range) - - Print(RES) - - #ploting inPl.plot (res[:,0],'-bs', label='susceptibles') -Pl.plot (res[:,1],'-ro', label='Infectious') toPl.legend (loc=0) +Pl.title ('SIS epidemic without births or deaths') -Pl.xlabel (' Time') thePl.ylabel ('susceptibles and Infectious') *Pl.savefig ('2.5-sis-high.png', dpi=900)#This does increase the resolution. $Pl.show ()
After running, the results are as follows:
Since individuals can recover after being infected, in a large time step, it is t=17 that the system reaches a steady state in which the proportion of infected individuals is constant. Therefore, in a stable state, only a limited number of individuals are infected, at this time does not mean that the infection disappeared, but at any point in time, infected individuals and the number of individuals recovered to a dynamic equilibrium, the proportion of the two sides remained unchanged. Please note that for larger recovery rate Gamma, the number of infected individuals declines exponentially, eventually the disease disappears, that is, at this point the rate of recovery is higher than the rate of infection, so depending on the size of the recovery rate gamma, there may eventually be two possible outcomes.
3. SIR Model
#-*-Coding:utf-8-*-import scipy.integrate as Spiimport numpy as Npimport Pylab as plbeta=1.4247gamma=0.14286ts=1.0nd=7 0.0s0=1-1e-6i0=1e-6input = (S0, I0, 0.0) def diff_eqs (inp,t): ' The main set of equations ' Y=np.zeros ((3)) V = inpy[0] =- Beta * v[0] * v[1]y[1] = beta * v[0] * v[1]-gamma * v[1]y[2] = gamma * V[1]return Y # for odeintt_start = 0.0; t_end = ND; T_inc = Tst_range = Np.arange (T_start, T_end+t_inc, t_inc) RES = Spi.odeint (diff_eqs,input,t_range) print (res) # Plotingpl.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, Toopl.show ()
The results are as follows:
Basic Infectious Disease model: SI, SIS, Sir and its Python code implementation