<jupyter_start><jupyter_text>Basic SIR model - with an extensionAfter presenting the basic model, we introduce a modification that allows for:- Recoveries and Deaths- Serious cases and Asymptomatic cases- differential rates of recruitment/infection Key result: _If infection travels faster among asymptomatic cases, "flattening the curve" may result in more deaths._Begin with a review and presentation of the basic SIR model.   A simple Susceptible-Infected-Recovered (SIR) modelThe SIR model is explained in many places (e.g. [here](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiologyThe_SIR_model)). Economists sometimes embed this in a larger economics model, for example this [NBER paper](https://www.nber.org/papers/w26902) working paper and references cited within. Python code below by <NAME> modified by <NAME>, guided by this ['Coronavirus Curve - Numberphile](https://www.youtube.com/watch?v=k6nLfCbAzgo) video.  If you're not doing so already you can run this as an interactive [web app](https://ricardian.herokuapp.com/).  The notebook with code is [here](https://github.com/jhconning/Econ-Teach/blob/master/notebooks/epidemic/SIRmodel.ipynb)  ModelThe proportion of the population are infected $I_t$, susceptible $S$, and recovered $R$ evolve over time according to these equation swhich depend on the transmission rate $\beta$ and the recovery rate $\gamma$:$$\begin{align}\frac{dI}{dt} &= \beta \cdot S \cdot I - \gamma \cdot I  \\\frac{dS}{dt} &=-\beta \cdot S \cdot I  \\\frac{dR}{dt} &= \gamma \cdot I\end{align}$$<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from ipywidgets import interact
N = 1          # Size of the population (so everything in proportions)
I0 = 0.01      # Initial proportion of the population infected
S0 = N - I0    # Initial proportion of the population susceptible
R0 = 0.0       # Initial proportion of the population recovered
maxT = 25      # max number of periods in simulation

beta = 0.5      # transmission rate
gamma = 0.1     # recovery rate
def SIR(y, t, beta, gamma):
    '''the SIR model'''
    S, I, R = y
    dSdt = -beta*S*I
    dIdt = beta*S*I - gamma*I
    dRdt = gamma*I
    return([dSdt, dIdt, dRdt])
def plotSIR(beta = beta, gamma = gamma, maxT = maxT):
    '''Solve differential equations in SIR and plot'''
    t = np.linspace(0, maxT, 1000)
    soln = odeint(SIR,[S0,I0,R0], t, args=(beta, gamma))
    soln = np.array(soln)
    plt.figure(figsize=[8,6])

    plt.plot(t, soln[:,0], linewidth=3, label = 'S(t)')
    plt.plot(t, soln[:,1], linewidth=3, label = 'I(t)')
    plt.plot(t, soln[:,2], linewidth=3, label = 'R(t)')
    plt.grid()
    plt.legend()
    plt.xlabel("Time")
    plt.ylabel("proportions")
    plt.title("SIR model")
    plt.show()<jupyter_output><empty_output><jupyter_text>Change the parameters with the sliders below (you'll only have interactivity with a running jupyter notebook server).<jupyter_code>interact(plotSIR, beta=(0,1,0.05), gamma=(0,1,0.05), maxT=(5,100,5));<jupyter_output><empty_output><jupyter_text>Below is a plot with the default parameters ($\beta=0.5$, $\gamma=0.1$) to have a graphic in case widget above does not display.<jupyter_code>plotSIR(beta, gamma, maxT)<jupyter_output><empty_output><jupyter_text>A More complex model In this model, social distancing for all can be counterproductive!One feature of the Coronavirus is that some cases are almost Asymptomatic. In general, these cases recover faster. Also, one might expect that these cases move faster to begin with. How would we do this? Let's make a model in which there are two types of cases and also deaths. The first type of cases are "asymptomatic" which we imagine are light cases, subscripted $a$. The second type are "bad" cases, which are subscripted $b$. Each type of case has its own infection rate and recovery rate, and one would expect in general that these rates, denoted by $\beta_a,\beta_b$ and $\gamma_a,\gamma_b$ would obey:$$\beta_a>\beta_b$$This is because people who are asymptomatic are more willing and able to move around; in a vacuum, they don't change their actions because they don't know they are sick. Also:$$\gamma_a>\gamma_b$$Recovery is easier from a lighter illness. Therefore we have the following equations: $$\begin{align}\frac{dI_a}{dt} &= \beta_a \cdot S \cdot I_a - \gamma_a \cdot I_a  \\\frac{dI_b}{dt} &= \beta_b \cdot S \cdot I_b - \gamma_b \cdot I_b - \delta \cdot I_b\\\frac{dS}{dt} &=-\beta_a \cdot S \cdot I_a - \beta_b \cdot S \cdot I_b  \\\frac{dR}{dt} &= \gamma_a \cdot I_a + \gamma_b \cdot I_b \\\frac{dD}{dt} &= \delta I_b\end{align}$$Note the further addition of an equation describing deaths, denoted $D$. It is assumed deaths occur at a rate $\delta$, and that deaths only occur among bad cases. Programs:<jupyter_code>N = 1                  # Size of the population (so everything in proportions)
Ia0 = 0.01             # Initial proportion of the population infected
Ib0 = 0.01
S0 = N - Ia0 - Ib0     # Initial proportion of the population susceptible
R0 = 0.0               # Initial proportion of the population recovered
D0 = 0.0               # Initial proportion of the population dead
maxT = 100              # max number of periods in simulation

betaa = 0.5           # transmission rate - a cases
betab = 0.5           # transmission rate - b cases
gamma = 0.15           # recovery rate - a cases
gammb = 0.15           # recovery rate - b cases

delta = 0.1            # death rate from b cases
def SIR_expanded(y, t, betaa, betab, gamma, gammb, delta):
    '''the SIR model with asymptomatic and severe (bad) cases'''
    S, Ia, Ib, R, D = y
    dSdt = -betaa*S*Ia - betab*S*Ib
    dIadt = betaa*S*Ia - gamma*Ia
    dIbdt = betab*S*Ib - gammb*Ib - delta*Ib
    dRdt = gamma*Ia + gammb*Ib
    dDdt = delta*Ib
    return([dSdt, dIadt, dIbdt, dRdt, dDdt])
def plotSIR_expanded(betaa = betaa, betab=betab, gamma = gamma, gammb=gammb, delta=delta, maxT = maxT):
    '''Solve differential equations in SIR and plot'''
    t = np.linspace(0, maxT, 1000)
    soln = odeint(SIR_expanded,[S0,Ia0,Ib0, R0, D0], t, args=(betaa, betab, gamma, gammb, delta))
    soln = np.array(soln)
    plt.figure(figsize=[8,6])

    plt.plot(t, soln[:,0], linewidth=3, label = 'S(t)')
    plt.plot(t, soln[:,1], linewidth=3, label = 'Ia(t)')
    plt.plot(t, soln[:,2], linewidth=3, label = 'Ib(t)')
    plt.plot(t, soln[:,3], linewidth=3, label = 'R(t)')
    plt.plot(t, soln[:,4], linewidth=3, label = 'D(t)')
    plt.grid()
    plt.legend()
    plt.xlabel("Time")
    plt.ylabel("proportions")
    plt.suptitle("SIR model with asymptomatic and bad cases")
    plt.title('Fraction of population dead: %0.4f - Fraction susceptible %0.4f' %(soln[-1, 4], soln[-1, 0]))
    plt.show()<jupyter_output><empty_output><jupyter_text>Note that the ONLY thing that differs between the graph below and the original one is that the b cases have a death rate. The take up rates and the recovery rates are the same for the two types of cases.<jupyter_code>plotSIR_expanded(betaa, betab, gamma, gammb, delta)<jupyter_output><empty_output><jupyter_text>Note how the inclusion of a death rate drastically alters the recruitment rates for the two types of affliction. Let's make an interactive plot:<jupyter_code>interact(plotSIR_expanded, betaa=(0,1,0.05), betab=(0,1,0.05), gamma=(0,1,0.05), gammb=(0,1,0.05), delta=(0,1,0.05), maxT=(5,100,5));<jupyter_output><empty_output><jupyter_text>One of the most interesting things about this model is the very complex role played by the deadliness of the disease. An increase in the death rate does not necessarily mean that the number of deaths in aggregate will increase! This is because deaths reduce recruitment into the population, especially in the presence of asymptomatic cases.   Flattening (or is that fattening) the curveSo, flattening the curve can mean a number of things, but let's suppose through measures like social distancing and staying home, we reduce the takeup rate in a particular fashion. In particular, suppose that we can reduce the rate $\beta$ by some amount $(1-f)$, where $f$ is our flattening parameter. But let's also introduce a "ham-handedness" parameter, which acknowledges that the flattening is likely to disproportionately impact asymptomatic cases. This is because if you have a bad case, you are staying home anyways. If you have a mild case, you wouldn't be staying home, but under the PAUSE policy, you now have to. So, the takeup rates for the two types of cases will converge. Suppose that the rate for bad cases is reduced by $f$, while the rate for asymptomatic cases is reduced by, so it becomes, effectively:$$\tilde{\beta}_b=(1-f)\beta_b$$For asymptomatic cases:$$\tilde{\beta}_a=\left(1-fh-(1-h)f\frac{\beta_b}{\beta_a}\right)$$Our parameter $h$ when 0 impacts both types proportionately, but when $h=1$, it now applies $f$ in such a way so as to force the types of cases to have the same take up rate.<jupyter_code>def plotSIR_expandedflat(betaa = betaa, betab=betab, gamma = gamma, gammb=gammb, delta=delta, maxT = maxT, f = 0, h = 0):
    '''Solve differential equations in SIR and plot, with flattening added in'''
    t = np.linspace(0, maxT, 1000)
    soln = odeint(SIR_expanded,[S0, Ia0, Ib0, R0, D0], t, args=(betaa*(1-f*(1-h)-f*h*betab/betaa), betab*(1-f), gamma, gammb, delta))
    soln = np.array(soln)
    plt.figure(figsize=[8,6])

    plt.plot(t, soln[:,0], linewidth=3, label = 'S(t)')
    plt.plot(t, soln[:,1], linewidth=3, label = 'Ia(t)')
    plt.plot(t, soln[:,2], linewidth=3, label = 'Ib(t)')
    plt.plot(t, soln[:,3], linewidth=3, label = 'R(t)')
    plt.plot(t, soln[:,4], linewidth=3, label = 'D(t)')
    plt.grid()
    plt.legend()
    plt.xlabel("Time")
    plt.ylabel("proportions")
    plt.suptitle("SIR model with asymptomatic and bad cases")
    plt.title('Fraction of population dead: %0.4f - Fraction susceptible %0.4f' %(soln[-1, 4], soln[-1, 0]))
    plt.show()
plotSIR_expandedflat()<jupyter_output><empty_output><jupyter_text>Interactive plot with "Flattening"Here, we allow one to adjust the "flattening parameter and the ham-handedness parameter.<jupyter_code>interact(plotSIR_expandedflat, betaa=(0,1,0.05), betab=(0,1,0.05), gamma=(0,1,0.05), gammb=(0,1,0.05), delta=(0,1,0.05), 
         maxT=(5,100,5), f=(0,1,0.05), h=(0,1,0.05));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1240833654959475, 'fraction_numerical': 0.022288691624855268, 'mean_word_length': 4.100885826771654, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 37, 'https://': 5, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '331405', 'license': 'MIT', 'n_tokens_mistral': 3595, 'n_tokens_neox': 3415, 'n_words': 1214, 'path': 'notebooks/SIRmodel-Extended.ipynb', 'repo_name': 'hunterecon2/epidemics', 'size': '10381'}