<a href="https://colab.research.google.com/github/JesperDramsch/corona/blob/master/Flattening_the_Curve_Interactively.ipynb" target="_parent"></a>

# Infection Modeling

The [Numberphile video](https://www.youtube.com/watch?v=k6nLfCbAzgo) on the corona curve inspired me to try solve these ODEs in Python and make it interactive.

The original video shows the SIR model, a ["compartmental model"](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology) for epidemology, deterministically modeling:

- **S**usceptible to the disease
- **I**nfectious with the disease
- **R**emoved from the disease and immune

This is done by solving these three ordinary differential equations (ODEs) and luckily SciPy has just the tools for this. The change of susceptible population is dependent on the transmission coefficient $\beta$, the number of infected and suscptible people $(I, S)$, as well as the population size $N$:

$
\frac{dS}{dt} = - \frac{\beta I S}{N},
$

the change of infected people is given by the influx of the infected people from the susceptible group from before and the transfer of infected people multiplied by the recovery factor $\gamma$:

$
\frac{dI}{dt} = \frac{\beta I S}{N}- \gamma I, 
$

the change of recovered people is simply given by the influx of infected people modified by the recovery factor from before:

$
\frac{dR}{dt} = \gamma I,
$

We can see that eventually all people will have gone through the infected group to the recovered group. This model simply depends on the assumption that the population N never changes (steady state).


```
import numpy as np
from scipy.integrate import solve_ivp

import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, FloatSlider
import ipywidgets as widgets

%matplotlib inline
```

These following start values show a simple percentage ($N=100$) and the values for the novel Corona virus that causes COVID-19.

Additionally, we'll define some time we want to analyze and the sampling of time for our simulation. Coarser is faster but finer is more stable.


```
N = 100

transmisision_rate = 3.2
recovery_rate = 0.23
waning_rate = 0.00

max_time = 5e1
dt = 1e-2
```

Then we need some initially infected people from our population.


```
# Population stats

def SIR_init(N=1, infected_ratio=0.01):
    """ Initialize S, I, R
    Initial Susceptible, Infected and Removed Population

    N (int): number of people in population
    infected_ratio (float): ratio of initially infected people
    """
    I = (infected_ratio * N)
    S = ((1 - infected_ratio) * N)
    R = 0
    return S, I, R
    
S_ini, I_ini, R_ini = SIR_init(N)
print(SIR_init(N))
```

    (99.0, 1.0, 0)


# Solve Differential Equations

We use a slightly modified version of the SIR formulation. Influenza and influenza-like viruses (possibly the corona virus?) can cause a loss of immunity. That means the SIRS model, which has a rate that changes from recovered to susceptible, is possibly better suited to model this virus. That changes the aforementioned ODEs to: 

$
\begin{align}
& \frac{dS}{dt} = - \frac{\beta I S}{N} + \xi R \\[6pt]
& \frac{dI}{dt} = \frac{\beta I S}{N}- \gamma I \\[6pt]
& \frac{dR}{dt} = \gamma I - \xi R
\end{align}
$

adding a change to $\frac{dR}{dt}$ of $R$ modified by the waning ratio $\xi$ that is added to the susceptible population $S$. As long as we keep $\xi=0$ the SIRS model is equivalent to the SIR model.


```

def dSIR_dt(t, SIR, N, transmisision_rate, recovery_rate, waning_rate=0): 
#def dSIR_dt(t, S, I, R, N, transmisision_rate, recovery_rate, waning_rate=0): 
    S, I, R = SIR

    infected = transmisision_rate * S * I / N
    removed = recovery_rate * I
    waned = waning_rate * R

    S_new = - infected + waned
    I_new = infected - removed
    R_new = removed - waned
    return (S_new, I_new, R_new)

```

The Scipy package provides us with several ODE solvers for "Initial Value Provlems", which is our kind of problem. These solvers are neatly packed in [solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp).

This function can be used to portray all the data we want in the entire time we simulate. However, it can solve specific points as well without providing that dense solution. Pretty neat:


```
def solve_SIR(N, max_time, dt, transmisision_rate, recovery_rate, waning_rate):
    t = np.arange(0, max_time, dt)
    SIR = solve_ivp(dSIR_dt, 
                    [0, max_time], 
                    SIR_init(N), 
                    args=(N, transmisision_rate, recovery_rate, waning_rate), 
                    dense_output=True)
    return t, SIR.sol(t)
```

# Interactive Plot

A nice big chunk of code to plot three lines. But now we can interactively explore the data and how we can flatten the curve.

Social Distancing reduces the transmission coefficient, try it out!


```
t, SIR = solve_SIR(N, max_time, dt, transmisision_rate, recovery_rate, waning_rate)

fig, ax = plt.subplots(figsize=(10,5))
S_plot, = ax.plot(t, SIR[0],  label='Susceptible')
I_plot, = ax.plot(t, SIR[1],  label='Infectious')
R_plot, = ax.plot(t, SIR[2],  label='Removed')
_ = ax.legend(loc='best')
_ = ax.set_xlabel('Time')
_ = ax.set_ylabel('Population')
_ = ax.set_title('Deterministic SIR(S) model')
r_o = plt.text(40, 0, f'$R_O$={transmisision_rate/recovery_rate:.2f}', fontsize=10)
plt.close()

def plot_SIR(transmission=3.2, recovery=.23, wane=.05):
    t, SIR = solve_SIR(N, max_time, dt, transmission, recovery, wane)
    S_plot.set_ydata(SIR[0])
    I_plot.set_ydata(SIR[1])
    R_plot.set_ydata(SIR[2])
    r_o.set_text(f'$R_O$={transmission/recovery:.2f}')
    fig.canvas.draw()
    display(fig)

style = {'description_width': 'initial'}
interactive_plot = interact(plot_SIR,
         transmission=FloatSlider(value=3.2, min=0, max=5, step=1e-2, continuous_update=False, description="Transmission Rate", style=style),
         recovery=FloatSlider(value=.23, min=0, max=1, step=1e-2, continuous_update=False, description="Recovery Rate", style=style),
         wane=FloatSlider(value=0, min=0, max=1, step=1e-2, continuous_update=False, description="Immunity Loss", style=style))

```


    interactive(children=(FloatSlider(value=3.2, continuous_update=False, description='Transmission Rate', max=5.0…



```

```

9438328e8e7f26946dd3f46ac1c0001fc08b5df4
ipynb
Jupyter Notebook
Flattening_the_Curve_Interactively.ipynb
JesperDramsch/corona
24e9792dad7313844c7f7d80ed8a80cc4722a055
MIT
Flattening_the_Curve_Interactively.ipynb
JesperDramsch/corona
24e9792dad7313844c7f7d80ed8a80cc4722a055
MIT
Flattening_the_Curve_Interactively.ipynb
JesperDramsch/corona
24e9792dad7313844c7f7d80ed8a80cc4722a055
MIT
Qwen/Qwen-72B
1. YES
2. YES