<jupyter_start><jupyter_text>Modeling Covid-19 Pandemic Using an SIR Model<NAME>, <NAME> and <NAME>  Biological and Mathematical BackgroundMathematical models are important tools for analyzing the spread and control of infectious diseases. These models incorporate the major components of disease transmission in simple mathematical equations. Modeling can be used to identify important epidemiological parameters and to evaluate prevention and control strategies. The use of models have led to discovery of mechanisms for controlling the spread of infections: vaccination for smallpox, magamenent of mosquito populations for malaria and prophilactic treatment for HIV. Understanding the transmission characteristics of infectious diseases can lead to better approaches to decreasing the transmission of these diseases.Comparmental epidemiological models divide the population into compartments based on their health status: e.g. susceptible, infectious and recovered. If a small group of infected individuals is introduced into a large population, a basic problem is to describe the spread of the infection within the population as a function of time. Differential equations describe how people 'move' between compartments: e.g. how often a susceptible person becomes infectious, how long does it take for an infectious persons to recovered. Here we present the Susceptible-Infectious-Removed (SIR) model and use actual epidemiological data to estimate the parameters of the model. In the SIR model, $S(t)$ denotes the number of susceptible individuals, who are not infected but are susceptible to the disease. $I(t)$ denotes the number of infected individuals, who have been infected and are able to spread the disease by contact with susceptible individuals. $R(t)$ denotes the number of removed individuals, no longer susceptible or infectious. Removal is carried out either by recovery with full immunity or by death from the disease. A flow diagram of the SIR model looks like follows:The SIR model has three parameters: the average number of contacts between infectious and susceptible individuals per unit time $(N)$, the probability of infection per contact $(\gamma)$ and the average duration of the disease $\big(\frac{1}{\mu}\big)$. The equations of the SIR model are as follow: Equations:$\frac{dS}{dt} = - N \gamma S\frac{I}{S_0}$&nbsp;$\frac{dI}{dt} =  N \gamma S\frac{I}{S_0} - \mu I$&nbsp;$\frac{dR}{dt} = \mu I$ Parameters:$N \rightarrow$ number of encounter per unit time&nbsp;$\gamma \rightarrow$ infection probability per encounter&nbsp;$S_0 \rightarrow$ initial number of susceptible individuals&nbsp;$\mu \rightarrow$ removal rate from the infectious populationNote that $\frac{I}{S_0}$ denotes the fraction of infected individuals with respect to the total population hence the term $N\gamma \frac{I}{S_0}$ denotes the probability of having an infectious encounter per unit time. This is a basic model but we can make some highly relevant general comments about epidemics and adequately describe some specific epidemics with such a model. More complicated compartmental structures are also possible: models in which infected individuals go through a latent phase before becoming infectious (SEIR), models in which recovered individuals lose immunity and become susceptible again (SIRS/SEIRS) and models where certain fraction of the population is vaccinated against the disease (SVEIR). What level of complexity is appropriate for the model is determined by the questions the model is intended to address.Given the values of the parameters and the number of of susceptible and infectious individuals the model can be used to answer important questions regarding an epidemic: 1) wether the infections will spread or not, 2) how will it develop in time, 3) when it will start to decline, etc. The basic reproduction number, $R_0$ is defined as the expected number of disease cases (secondary infections) produced by an infected individual and can be determined from the parameters of the model: $R_0 = \frac{N \gamma}{\mu S_0}$$R_0$ is an indicator of how fast the infection is spreading, and can be used to determined wether the infection will spread or die out.  If every single infected individual produces more that one infection, $R_0 > 1$ and the epidemic ensues. The whole question of thresholds in epidemics is obviously important: the threshold $R_0 = 1$ is the dividing line between the infection dying out and the onset of an epidemic. $R_0$ also determines the number of individuals that will get infected over the course of the epidemic. One way to reduce the reproduction rate is to reduce the number of suceptibles ($S_0$). Vaccination is the common method of doing this and it has been successful in eradicating smallpox.  Instructions to Use the Application1. Load 'Active Cases' data by choosing region/country2. Estimate $\mu$ by zooming in to the downside of the 'Active Cases' curve 3. Estimate $N\gamma$ by zooming in to the upside of the 'Active Cases' curve4. Estimate $R_0$ and $S_0$5. Simulate the model with the estimated parameter values6. Fit the model to the 'Active Cases' data  Estimating Parameters from Experimental DataThe Covid-19 pandemic is one of the greatest public health challenges of our lifetimes. Along with publicly available epidemiological data, modeling can be used to determine how fast the disease is spreading, what would be the expected size of the pandemic (e.g. how many people will become infected) and what level of vaccionation rates and efficacy are required to control it. This application allows users to estimate the basic reproductive ratio $R_0$ along with other parameters such as the infection rate $\frac{N \gamma}{S_0}$ and the duration of the disease $\big(\frac{1}{\mu}\big)$. To estimate these parameters from experimental data, we assume that the SIR model is a good approximation at the early stages of the pandemic and include only data corrresponding to the first 'wave' of infections. Load Experimental DataTo estimate the parameters we need data for 'active cases' over time. Load the active cases data from the country/region to estimate the parameters from by choosing from the pulldown menu.<jupyter_code>Country_W = widgets.Dropdown(
    options=['Spain','Italy','New Zealand'],
    value='Italy',
    description='Country: ',
    disabled=False,
)

def load_data(Country):
    Active_Cases = np.genfromtxt('PP_%s.txt' % Country, delimiter=',')
    T = np.arange(len(Active_Cases))
    return(T,Active_Cases)

def plot_experimental_data(Country):
    global T
    global Active_Cases
    global CurrentCountry
    CurrentCountry = Country
    [T,Active_Cases] = load_data(Country.replace(' ','_'))
    plt.plot(T, Active_Cases,label='Data')
    plt.ylabel('Number of Cases')
    plt.xlabel('Time (days)')
    plt.title('Active Cases in %s' % Country)
    plt.grid()
    plt.legend()
    plt.show()

def on_value_change_data(change):
    for i in range(3):
        LW_beta.value = [0, max(T)]
        LW_beta.max = max(T)
        LW_mu.value = [0, max(T)]
        LW_mu.max = max(T)
        BU_mu.value = False
        BU_beta.disabled = True
        BU_beta.value = False
        BU_R0.disabled = True
        BU_R0.value = False
        BU_initial_guess.disabled = True
        BU_initial_guess.value = False
        BU_fit.disabled = True
        BU_fit.value = False

interact(plot_experimental_data, Country=Country_W)
Country_W.observe(on_value_change_data,'value')
print('Data retreived from John Hopkins University COVID-19 tracker repository: https://github.com/CSSEGISandData/COVID-19')<jupyter_output><empty_output><jupyter_text>Estimating $\mu$ You will learn: how to estimate $\mu$ by looking at the downside of the active cases curveTo estimate $\mu$ we assume that in late stages of the epidemic, the rate of new infections is small $\frac{N \gamma}{S_0}SI \approx 0$. This can the result of containment measusers such as lockdowns and mask-wearing, which reduce the number of contacts $N$ and the probability of infection for contact $\gamma$, but also because the number of remaining susceptible individuals is small $S(t) \approx 0$ (e.g. a large number of people have been infected). $\frac{dI}{dt} = \frac{N \gamma}{S_0}SI - \mu I \rightarrow \frac{dI}{dt} = - \mu I $The value of $\mu$ is given by the slope of the downside of the curve of active cases:$I_{tmax} = I_{tmin}e^{- \mu (tmax-tmin)} \rightarrow \mu = \frac{ln(I_{tmin})-ln(I_{tmax})}{tmax-tmin}$ [$days^{-1}$] If this approximation is reasonable, when we plot the data on a semi-log scale, it should look like a straight line. If the slope is not linear, then our approximation is not so good and we have to define the time windown in which we want to measure the slope.  Instructions1. In the following graph choose a time window that captures the downside of the curve of active cases by changing the time range using the slider2. Make sure that the data looks like a straight line3. Press the button to calculate $\mu$. The red line is the linear regression of the data from which the slope is estimated to calculate $\mu$.<jupyter_code>LW_mu = widgets.IntRangeSlider(
    value=[0, max(T)],
    min=0,
    max=max(T),
    description='Times:',
    disabled=False,
    continuous_update=False,
    orientation='horizontal',
    readout=True,
    readout_format='d',
)

BU_mu = widgets.ToggleButton(
    value=False,
    description='Calculate mu',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Description',
    icon='check' # (FontAwesome names without the `fa-` prefix)
)

def do_linear_fit_mu(tmin,tmax):
    global mu
    x = T[tmin:tmax]
    y = np.log(Active_Cases[tmin:tmax])
    mu = np.abs(np.polyfit(x, y, 1)[0])
    line = x*np.polyfit(x, y, 1)[0] + np.polyfit(x, y, 1)[1]
    plt.plot(x,np.exp(line),'r',label='Linear Fit')
    print('mu = %.2f per day' % mu)
#     print('mu2 =',(np.log(Active_Cases[tmin]/Active_Cases[tmax]))/(tmax-tmin))
    print('Duration of the infection = %.2f days' % (1.0/mu))

def on_value_change_mu(change):
    BU_mu.value = False
    LW_beta.value = [0, max(T)]
    LW_beta.max = max(T)
    BU_beta.disabled = True
    BU_beta.value = False
    BU_R0.disabled = True
    BU_R0.value = False
    BU_initial_guess.disabled = True
    BU_initial_guess.value = False
    BU_fit.disabled = True
    BU_fit.value = False
    
def plot_data_mu(L,B):
    MinAct = min(Active_Cases[L[0]:L[1]])
    MaxAct = max(Active_Cases[L[0]:L[1]])
    plt.subplot(2,1,1)
    plt.title('Active Cases in %s' % CurrentCountry)
    plt.plot(T, Active_Cases,label='Data')
    plt.gca().add_patch(Rectangle((L[0], MinAct), 
                        L[1]-L[0], MaxAct-MinAct, 
                        fc ='none',  
                        ec ='red', 
                        lw = 2) ) 
    plt.yscale('log')
    plt.grid()
    plt.ylabel('Log Number of Cases')
    plt.subplot(2,1,2)
    plt.plot(T[L[0]:L[1]], Active_Cases[L[0]:L[1]],label='Data')
    plt.yscale('log')
    plt.ylabel('Log Number of Cases')
    plt.xlabel('Time (days)')
    plt.gca().spines['top'].set_color('red')
    plt.gca().spines['top'].set_linewidth(2)
    plt.gca().spines['bottom'].set_color('red')
    plt.gca().spines['bottom'].set_linewidth(2)
    plt.gca().spines['left'].set_color('red')
    plt.gca().spines['left'].set_linewidth(2)
    plt.gca().spines['right'].set_color('red')
    plt.gca().spines['right'].set_linewidth(2)
    plt.grid()
    if B:
        BU_beta.disabled = False
        do_linear_fit_mu(L[0],L[1])
    plt.legend()
    plt.show()

interact(plot_data_mu, L=LW_mu, B=BU_mu)
LW_mu.observe(on_value_change_mu)<jupyter_output><empty_output><jupyter_text>Estimating $N\gamma$ You will learn: how to estimate $N\gamma$ by looking at the upside of the active cases curveNow that we have estimated $\mu$, we can estimate $N\gamma$. Note that given the available data, we cannot independently estimate $N$ and $\gamma$, we can only estimate their product (often called $\beta$). To estimate $\beta$ we assume that at the beginning of the pandemic the whole population is susceptible because the number of infected individuals is small $\frac{S}{S_0} \approx 1$. $\beta \equiv N \gamma$$\frac{dI}{dt} = \frac{\beta}{S_0}SI - \mu I \rightarrow \frac{dI}{dt} = (\beta - \mu)I $The value of $\beta - \mu$ is given by the slope of the upside of the curve of active cases:$I(t) = I_0e^{(\beta - \mu )t} \rightarrow \beta = \frac{ln(I(t))-ln(I_{0})}{t-t_0} + \mu$ [$days^{-1}$]If this approximation is reasonable, when we plot the data on a semi-log scale, it should look like a straight line. If the slope is not linear, then our approximation is not so good and we have to define the time window in which we want to measure the slope.  Instructions1. In the following graph choose a time window that captures the upside of the curve of active cases by changing the time range using the slider2. Make sure that the data looks like a straight line. Once the appropiate time range has been selected3. Press the button to calculate $\beta$. The red line is the linear regression of the data from which the slope is estimated to calculate $\beta$.<jupyter_code>LW_beta = widgets.IntRangeSlider(
    value=[0, max(T)],
    min=0,
    max=max(T),
    description='Times:',
    disabled=False,
    continuous_update=False,
    orientation='horizontal',
    readout=True,
    readout_format='d',
)

BU_beta = widgets.ToggleButton(
    value=False,
    description='Calculate Beta',
    disabled=True,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Description',
    icon='check' # (FontAwesome names without the `fa-` prefix)
)

def do_linear_fit_beta(tmin,tmax):
    global beta
    x = T[tmin:tmax]
    y = np.log(Active_Cases[tmin:tmax])
    upside_slope = np.abs(np.polyfit(x, y, 1)[0])
    beta = upside_slope + mu
    line = x*np.polyfit(x, y, 1)[0] + np.polyfit(x, y, 1)[1]
    plt.plot(x,np.exp(line),'r',label='Linear Fit')
    print('Beta = %.2f per day' % beta)
#     print('Beta2 =',np.log(Active_Cases[tmax]/Active_Cases[tmin])/(tmax-tmin)+mu)

def on_value_change_beta(change):
    BU_beta.value = False
    BU_R0.disabled = True
    BU_R0.value = False
    BU_initial_guess.disabled = True
    BU_initial_guess.value = False
    BU_fit.disabled = True
    BU_fit.value = False

def plot_data_beta(L,B):
    MinAct = min(Active_Cases[L[0]:L[1]])
    MaxAct = max(Active_Cases[L[0]:L[1]])
    plt.subplot(2,1,1)
    plt.title('Active Cases in %s' % CurrentCountry)
    plt.plot(T, Active_Cases,label='Data')
    plt.gca().add_patch(Rectangle((L[0], MinAct), 
                        L[1]-L[0], MaxAct-MinAct, 
                        fc ='none',  
                        ec ='red', 
                        lw = 2) ) 
    plt.yscale('log')
    plt.ylabel('Log Number of Cases')
    plt.grid()
    plt.subplot(2,1,2)
    plt.plot(T[L[0]:L[1]], Active_Cases[L[0]:L[1]],label='Data')
    plt.yscale('log')
    plt.ylabel('Log Number of Cases')
    plt.xlabel('Time (days)')
    plt.gca().spines['top'].set_color('red')
    plt.gca().spines['top'].set_linewidth(2)
    plt.gca().spines['bottom'].set_color('red')
    plt.gca().spines['bottom'].set_linewidth(2)
    plt.gca().spines['left'].set_color('red')
    plt.gca().spines['left'].set_linewidth(2)
    plt.gca().spines['right'].set_color('red')
    plt.gca().spines['right'].set_linewidth(2)
    plt.grid()
    if B:
        BU_R0.disabled = False
        do_linear_fit_beta(L[0],L[1])
    plt.legend()
    plt.show()
    
interact(plot_data_beta, L=LW_beta, B=BU_beta)
LW_beta.observe(on_value_change_beta)<jupyter_output><empty_output><jupyter_text>Estimating $R_0$ and $S_0$$R_0$ is the number of individuals an infected individual will infect through the duration of the infection. $R_0$ is directly proportional to the number of contacts between susceptibles and infected individuals per unit time ($N$) and the propbability of infection per contact ($\gamma$). $R_0$ is also proportional to the duration of the infection ($\frac{1}{\mu}$). Now that we have $N \gamma$ and $\mu$ we can estimate $R_0$.$R_0 = \frac{N \gamma}{\mu}$Once we have estimated $R_0$, we can estimate the effective size  $S_0$ from the maximun of the active cases $I_{max}$.$S_0 = \frac{I_{max} R_0}{R_0 - 1 - ln(R_0)}$ Instructions1. Press the button to calculate $R_0$ and $S_0$.<jupyter_code>BU_R0 = widgets.ToggleButton(
    value=False,
    description='Calculate R0',
    disabled=True,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Description',
    icon='check' # (FontAwesome names without the `fa-` prefix)
)

def print_params(B):
    global R0
    global S0
    if B:
        BU_initial_guess.disabled = False
        BU_initial_guess.value = False
        
        #Estimating R0
        R0 = beta/mu
        print('R0 = %.2f' % R0)

        #Estimating S0
        Imax = max(Active_Cases)
        S0 = Imax / (1 -(1+ np.log(R0))/R0)
        print('S0 = %.2f' % S0)
        
interact(print_params, B=BU_R0)<jupyter_output><empty_output><jupyter_text>Simulating SIR Model with Parameter EstimatesOnce we estimate the basic model parameters, we can plug in those values into the model equations, simulate the model and compare it to the experimental data. The chi-squared error is the sum of the squared differences between the epidemiological data and the simulation results. Press the button to simulate the model with the estimated parameter value and compare simulation results to the epidemiological data.<jupyter_code>toFit = ['S1','I1','b','k1']
nParameters = len(toFit)
params = lmfit.Parameters()

BU_initial_guess = widgets.ToggleButton(
    value=False,
    description='Simulate Model',
    disabled=True,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Description',
    icon='check' # (FontAwesome names without the `fa-` prefix)
)

def plot_initial_guess(B):
    if B:
        BU_fit.disabled = False
        params.add('S1', value=S0, min=S0/10.0, max=S0*20.0)
        params.add('I1', value=1.0, min=0, max=S0/10.0)
        params.add('b', value=beta/S0, min=0, max=1)
        params.add('k1', value=mu, min=0, max=1)
        cost_fun(params)

interact(plot_initial_guess, B=BU_initial_guess)<jupyter_output><empty_output><jupyter_text>Fitting SIR Model to Epidemiological DataNow that we have an original guess of the parameters of the model, we can use an optimizer to fit the model to the epidemiological data and get a more accurate value of the parameters. The role of the optimizer is to minimize the Chi-squared error by exploring different combinations of parameter values starting from the initial guesses we estimated above. Here we use a local 'nelder-mead' optimization strategy because we have good initial guesses. More complicated strategies can be employed when initial guesses of the parameters cannot be estimated easily from available data. Press the button to fit the model and obtain a better estimate of the parameter values.<jupyter_code>minimizer = lmfit.Minimizer(cost_fun, params)

BU_fit = widgets.ToggleButton(
    value=False,
    description='Fit Model',
    disabled=True,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Description',
    icon='check' # (FontAwesome names without the `fa-` prefix)
)

def plot_fit(B):
    if B:
        result = minimizer.minimize(method='nelder')
        lmfit.report_fit(result.params, min_correl=0.5)
        k1 = result.params['k1'].value
        b = result.params['b'].value
        S1 = result.params['S1'].value
        print('Estimated R0 = %.2f, Fitted R0 = %.2f' % (R0,b*S1/k1))
        print('Duration of the infection = %.2f days' % (1.0/k1))

interact(plot_fit, B=BU_fit)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10673930589184827, 'fraction_numerical': 0.01104721549636804, 'mean_word_length': 3.9736578023080784, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 29, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '542180', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6102, 'n_tokens_neox': 5856, 'n_words': 2196, 'path': 'App_Prototype_7.ipynb', 'repo_name': 'joaponte/tellurium-nanohub-base-main', 'size': '19851'}