Python implementation of SIS model with a constant contact rate.

Introduction

Previously, a Python Code solved the SI Epidemiological Differential equation model. Here, we give an implementation of the SIS model. Then, I compare the computer generated solution to the analytic solution. In addition, the new code adds the ability to store the results in a comma delimited file that can be read by Excel.

Since we need a test case: the Hubei data will continue to serve that purpose.

Implementation of the SIS model.

We added the term $-\gamma_r I(t)$ to the right hand side of the differential equation for the number of people that are infected. This is in the function definition of f in the value of derivative. The other major addition is the code at the end that allows the user to store the information in a file.

Parameter selection

The Hubei data determined the parameter $\lambda$ . The recovery rate for COVID-19 determines the parameter $\gamma$ . After running the software, we used the comma delimited code as input to Excel to generate the graph. Below is the graphical output of the results. The population units are in thousands of people, so 40,000 on the vertical axis represents 40,000,000 infected people. The agreement between the formula and the computer solution (circles) is excellent.

Eventually, the infected population rises to just under 83% of the population. This does not mean that 17% of the population is unaffected. Eventually everyone gets the disease. What is true in this simulation is that at the end for every person that recovers another gets sick, so that at any time 83% of the population is sick. So this curve is not representing a herd immunity, it represents massive sickness.

This shows the excellent agreement between the implementation of the SIS model in Python and the analytic solution.

This represents one thousand more infections than the Hubei data indicates. No reasonable person would stand by and watch the entire population get infected without trying some defensive measures. In the absence of a vaccine the only way to reduce the infections is changing the infection rate. In the case of the model this means adopting measures that reduce the parameter plam in the code or $\lambda$ in the theory. We do that in the next code segment.

Links

To return to the main page for the SIS model go here. https://davidalarrabee.com/?page_id=754

If you want to see the modification for a varying $\lambda$ follow this link. (TBD)

The Python Code

Copy Code

                    
# ************** SIS Model differential equation Solver vs. 1.01 May 20, 2020
# miscellaneious constants and stuff 
Title  = "SIS Model Differential Equation Solver"
Vers   = "vs. 1.01 May 20, 2020  David A. Larrabee"
EOL = "\n"
Header = " day, Susceptible, infectious,   formula,     error,    % infected "
e=2.7182818284  # Euler's number (approx)

# ***************************************************************************
#   Imports and definitions
# **************************************************************************       
import numpy as np
from scipy.integrate import odeint

def f(y, t, parm):             # returns the derrivative of f(t) (y)
    I = y[0]                   # get the current value of y
    plam,gammar, Npop = parm             # get the parameter
    derivative = [plam*I*(1-I/Npop)] -gammar*I     # calculate the derrivative
    return derivative

#****************************************************************************
#   These are parameters the User can Input
#****************************************************************************
plam = 0.2587      # growth rate at start (infections per day per infector)
gammar = 0.045     # the daily probability of recovery from the virus
Npop = 58500       # Total "population" (Hubei in thousands)
I0 = 1.619         # Initial number of infected (hubei in thousdands)
Finaltime = 90.    # The sumulation goes from time = 0 to here (in days)
Increment = 1.0    # How often you want the result retured (in days)
#****************************************************************************
#     Set up for the simulation 
#****************************************************************************
print(Title)
print(Vers)
print(" ")

parm = [plam,gammar, Npop]   #ODE solver wants a list of Parameters
y0 = [I0]       #ODE also wants a list of initial values

# Make an array of the times at which we want a solution
times = np.arange(0.0, Finaltime, Increment)

#*************************************************************
#             Call the ODE solver
#*************************************************************
Solutions = odeint(f, y0, times, args=(parm,))

#*******************************************************************
# Compare the results to theory and print out the results to screen
#*******************************************************************

print(Header)
lambdastar = plam - gammar
Nstar=Npop*lambdastar/plam

for i in range(len(times)): # printout and compare with actual solution
    t = times[i]
    I = Solutions[i,0]
    S = Npop - I
    val = Nstar/(1+(Nstar-I0)/I0 * e**(-1.0*lambdastar*t) ) # analytic solution
    error = (val-I)/val
    pinfect = I/Npop*100.0
    print("{:4.1f}".format(t), "{:11.1f}".format(S),"{:11.1f}".format(I),
          "{:11.1f}".format(val),"  ","{:+8.4E}".format(error),
          "  ","{:8.4f}".format(pinfect) )
          
#*******************************************************************
# drop the results to a file?
#*******************************************************************
ans = input("insert filename to store results in a file. type N to finish: ")

if ans != "N":
     fhandler=open(ans,"a")
     strng=Title + EOL + Vers+EOL
     strng=strng+'Total number of people =, '+str(Npop)+ EOL
     strng=strng+"Total number of people initially infected =, "
     strng=strng+str(I0)+EOL
     strng=strng+" lanmbda =,"+str(plam)+EOL
     strng=strng+" gammar =,"+str(gammar)+EOL     
     strng=strng+" Ending time, "+str(Finaltime)+EOL
     strng=strng+" fraction of a day in a timestep, "+str(Increment)+EOL 
     strng = strng+Header+EOL
     fhandler.write(strng)
     for i in range(len(times)): # printout and compare with actual solution
           t = times[i]
           I = Solutions[i,0]
           S = Npop - I
           val = Nstar/(1+(Nstar-I0)/I0 * e**(-1.0*lambdastar*t) ) # analytic solution
           error = (val-I)/val
           pinfect = I/Npop*100.0
           
           strng = str(times[i])+","
           strng = strng+str(S)+","
           strng = strng+str(I)+","
           strng = strng+str(val)+","
           strng = strng+str(error)+","
           strng = strng+str(pinfect)+EOL
           
           fhandler.write(strng)
           
     fhandler.close()
    
print("All Done")

Code output

Copy Code

SIS Model Differential Equation Solver vs. 1.01 May 20, 2020 David A. Larrabee day, Susceptible, infectious, formula, error, % infected 0.0 58498.4 1.6 1.6 +0.0000E+00 0.0028 1.0 58498.0 2.0 2.0 +8.5523E-09 0.0034 2.0 58497.5 2.5 2.5 +4.4555E-09 0.0042 3.0 58496.9 3.1 3.1 -1.5477E-09 0.0053 4.0 58496.2 3.8 3.8 -1.6275E-08 0.0065 5.0 58495.3 4.7 4.7 -2.1809E-08 0.0081 6.0 58494.2 5.8 5.8 -2.4312E-08 0.0100 7.0 58492.8 7.2 7.2 -2.7152E-08 0.0124 8.0 58491.1 8.9 8.9 -3.2746E-08 0.0153 9.0 58488.9 11.1 11.1 -3.0027E-08 0.0189 10.0 58486.3 13.7 13.7 -3.1853E-08 0.0234 11.0 58483.0 17.0 17.0 -3.3154E-08 0.0290 12.0 58479.0 21.0 21.0 -3.4854E-08 0.0359 13.0 58474.0 26.0 26.0 -3.8679E-08 0.0445 14.0 58467.8 32.2 32.2 -4.2226E-08 0.0551 15.0 58460.1 39.9 39.9 -4.6448E-08 0.0682 16.0 58450.6 49.4 49.4 -4.9390E-08 0.0845 17.0 58438.8 61.2 61.2 -5.2512E-08 0.1045 18.0 58424.3 75.7 75.7 -5.6076E-08 0.1294 19.0 58406.3 93.7 93.7 -5.8664E-08 0.1602 20.0 58384.0 116.0 116.0 -6.0944E-08 0.1983 21.0 58356.5 143.5 143.5 -6.3416E-08 0.2454 22.0 58322.4 177.6 177.6 -6.5065E-08 0.3036 23.0 58280.3 219.7 219.7 -6.6262E-08 0.3756 24.0 58228.2 271.8 271.8 -6.8173E-08 0.4646 25.0 58163.9 336.1 336.1 -6.7771E-08 0.5745 26.0 58084.5 415.5 415.5 -6.5615E-08 0.7102 27.0 57986.6 513.4 513.4 -6.1783E-08 0.8776 28.0 57865.9 634.1 634.1 -5.6112E-08 1.0840 29.0 57717.3 782.7 782.7 -4.7385E-08 1.3380 30.0 57534.5 965.5 965.5 -3.4910E-08 1.6505 31.0 57310.1 1189.9 1189.9 -1.8972E-08 2.0340 32.0 57035.2 1464.8 1464.8 -1.9531E-09 2.5039 33.0 56699.2 1800.8 1800.8 -5.0459E-09 3.0782 34.0 56289.8 2210.2 2210.2 +3.9211E-09 3.7781 35.0 55792.7 2707.3 2707.3 +7.6439E-09 4.6278 36.0 55191.9 3308.1 3308.1 +1.2609E-08 5.6549 37.0 54469.5 4030.5 4030.5 +1.4466E-08 6.8898 38.0 53606.4 4893.6 4893.6 +1.3057E-08 8.3651 39.0 52583.3 5916.7 5916.7 +3.8318E-09 10.1140 40.0 51381.3 7118.7 7118.7 -1.0640E-08 12.1687 41.0 49984.1 8515.9 8515.9 -2.7129E-08 14.5570 42.0 48380.1 10119.9 10119.9 -2.5799E-08 17.2990 43.0 46564.5 11935.5 11935.5 -2.6669E-08 20.4025 44.0 44542.3 13957.7 13957.7 -3.2034E-08 23.8594 45.0 42329.6 16170.4 16170.4 -2.8104E-08 27.6417 46.0 39955.4 18544.6 18544.6 -2.6947E-08 31.7001 47.0 37460.8 21039.2 21039.2 -8.7815E-09 35.9645 48.0 34896.5 23603.5 23603.5 +1.0488E-08 40.3478 49.0 32319.6 26180.4 26180.4 +9.4684E-09 44.7528 50.0 29788.1 28711.9 28711.9 +8.4971E-09 49.0801 51.0 27356.1 31143.9 31143.9 +8.6682E-09 53.2374 52.0 25069.3 33430.7 33430.7 +2.8838E-09 57.1465 53.0 22961.9 35538.1 35538.1 +4.5598E-10 60.7489 54.0 21055.6 37444.4 37444.4 -7.8452E-10 64.0075 55.0 19360.1 39139.9 39139.9 -2.0320E-09 66.9058 56.0 17874.5 40625.5 40625.5 -6.1598E-10 69.4453 57.0 16589.8 41910.2 41910.2 +2.3381E-09 71.6413 58.0 15491.5 43008.5 43008.5 +4.1657E-09 73.5189 59.0 14561.5 43938.5 43938.5 +5.7158E-09 75.1085 60.0 13780.6 44719.4 44719.4 +6.0279E-09 76.4434 61.0 13129.4 45370.6 45370.6 +1.1957E-08 77.5565 62.0 12589.5 45910.5 45910.5 +1.7243E-08 78.4794 63.0 12144.0 46356.0 46356.0 +2.1152E-08 79.2410 64.0 11777.9 46722.1 46722.1 +2.2026E-08 79.8668 65.0 11478.0 47022.0 47022.0 +2.1090E-08 80.3796 66.0 11232.9 47267.1 47267.1 +1.9782E-08 80.7985 67.0 11033.1 47466.9 47466.9 +1.8050E-08 81.1399 68.0 10870.6 47629.4 47629.4 +1.5862E-08 81.4178 69.0 10738.5 47761.5 47761.5 +1.7373E-08 81.6437 70.0 10631.2 47868.8 47868.8 +1.7503E-08 81.8270 71.0 10544.3 47955.7 47955.7 +1.7213E-08 81.9756 72.0 10473.8 48026.2 48026.2 +1.7399E-08 82.0960 73.0 10416.8 48083.2 48083.2 +1.7294E-08 82.1935 74.0 10370.6 48129.4 48129.4 +1.6992E-08 82.2724 75.0 10333.3 48166.7 48166.7 +1.7758E-08 82.3363 76.0 10303.1 48196.9 48196.9 -6.0195E-10 82.3879 77.0 10278.7 48221.3 48221.3 -2.8438E-09 82.4297 78.0 10258.9 48241.1 48241.1 +9.7570E-10 82.4634 79.0 10243.0 48257.0 48257.0 +8.5952E-10 82.4907 80.0 10230.1 48269.9 48269.9 -6.7638E-10 82.5127 81.0 10219.6 48280.4 48280.4 +6.5238E-10 82.5305 82.0 10211.2 48288.8 48288.8 -6.4362E-11 82.5449 83.0 10204.4 48295.6 48295.6 -9.1762E-10 82.5565 84.0 10198.9 48301.1 48301.1 -3.4805E-10 82.5659 85.0 10194.5 48305.5 48305.5 -3.5870E-10 82.5735 86.0 10190.9 48309.1 48309.1 -5.8520E-10 82.5796 87.0 10188.0 48312.0 48312.0 -4.5463E-10 82.5846 88.0 10185.7 48314.3 48314.3 -3.5481E-10 82.5886 89.0 10183.8 48316.2 48316.2 -5.1384E-10 82.5918 insert filename to store results in a file. type N to finish: SISfl01.csv All Done

Listening for the Cry of the Earth

Part 3.1.1 Python Implementation of SIS Model with constant Lambda

Introduction

Implementation of the SIS model.

Parameter selection

Links

The Python Code

Code output

A dialog between Science, economics and Religion