notebook.community

Edit and run



In [1]:

    
import sympy as sym
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt



In [2]:

    
S, I, R = sym.Function("S"), sym.Function("I"), sym.Function("V")
N, mu, alpha, beta, t = sym.symbols("N, mu, alpha, beta, t")



In [3]:

    
eq1 = sym.Derivative(S(t), t) - (- alpha * S(t) * I(t) - mu * R(t))
eq2 = sym.Derivative(I(t), t) - (alpha * I(t) * S(t) / N - beta * I(t))
eq3 = sym.Derivative(R(t), t) - (beta * I(t) + mu * R(t))



In [4]:

    
sym.dsolve((eq1, eq2, eq3))









    



---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-4-f5e5bded6978> in <module>()
----> 1 sym.dsolve((eq1, eq2, eq3))

~/anaconda3/envs/cfm/lib/python3.6/site-packages/sympy/solvers/ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
    606             "number of functions being equal to number of equations")
    607         if match['type_of_equation'] is None:
--> 608             raise NotImplementedError
    609         else:
    610             if match['is_linear'] == True:

NotImplementedError:

Further investigation shows that an exact solution to this system of differential equations is difficult. Let us do this numerically:



In [5]:

    
def dx(x, t, alpha, beta, mu):

    return (- alpha * x[1] * x[0] - mu * x[0],
            alpha * x[1] * x[0]  - beta * x[1],
            beta * x[1] + mu * x[0])



In [6]:

    
alpha = 1 / 1000  # Every 1000 interactions leads to infection
beta = 1 / 5  # take 5 time units to recover
N = 10 ** 4  # Population of 10 thousand people
mu = 0  # 0 vaccination percentage

ts = np.linspace(0, 10, 5000)
xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
S, I, R = xs.T
plt.figure()
plt.plot(ts, S, label="Susceptibles")
plt.plot(ts, I, label="Infected")
plt.plot(ts, R, label="Recovered")
plt.legend()
plt.title(f"$\max(I)={round(max(I))}$ ($\\alpha={alpha}$, $\\beta={beta}$, $\mu={mu}$)")
plt.savefig("base_scenario.pdf");



In [7]:

    
mu = 1 / 2  # Vaccinate half the population
ts = np.linspace(0, 10, 5000)
xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
S, I, R = xs.T
plt.figure()
plt.plot(ts, S, label="Susceptibles")
plt.plot(ts, I, label="Infected")
plt.plot(ts, R, label="Recovered")
plt.legend()
plt.title(f"$\max(I)={round(max(I))}$ ($\\alpha={alpha}$, $\\beta={beta}$, $\mu={mu}$)")
plt.savefig("moderate_vaccination_rate.pdf");



In [8]:

    
mu = 99 / 100  # Vaccinate 99% of the population
ts = np.linspace(0, 10, 5000)
xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
S, I, R = xs.T
plt.figure()
plt.plot(ts, S, label="Susceptibles")
plt.plot(ts, I, label="Infected")
plt.plot(ts, R, label="Recovered")
plt.legend()
plt.title(f"$\max(I)={round(max(I))}$ ($\\alpha={alpha}$, $\\beta={beta}$, $\mu={mu}$)")
plt.savefig("high_vaccination_rate.pdf");



In [9]:

    
vaccination_rates = np.linspace(0, 1, 500)
max_percent_of_infected = []
for mu in vaccination_rates:
    xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
    S, I, R = xs.T
    max_percent_of_infected.append(max(I) / N)
plt.figure()
plt.plot(vaccination_rates, max_percent_of_infected)
plt.xlabel("$\mu$")
plt.ylabel("% of population infected")
plt.savefig("effect_of_vaccination_rate.pdf");