<jupyter_start><jupyter_text>SEIRD Approach in Modeling of Novel Coronavirus Disease 2019 Dynamics Part 1 **Overview.** This report serves as a high-level overview of the compartmental epidemiological deterministic model, more concretely, Susceptible-Exposed-Infected-Recovered-Deceased (SEIRD) model, detail mathematical description is to be found elsewhere. Additionaly, uncertainty quantification of the SEIRD output is achieved using Markov chain Monte Carlo (MCMC) method. The goal is to examine how the variability of input epidemiological parameters affects the overall output of the model by treating input epidemiological parameters as random variables and recast the original deterministic model as stochastic model without affecting the configuration of the model itself.**SEIRD model**\begin{aligned}\frac{S(t + \Delta t)}{\Delta t} &= -\frac{\beta \cdot S(t) \cdot I(t)}{N} - \langle \delta \rangle \cdot S(t) \cdot E(t) \\\frac{E(t + \Delta t)}{\Delta t} &= \frac{\beta \cdot S(t) \cdot I(t)}{N} - \langle \alpha \rangle \cdot E(t) + \langle \delta \rangle \cdot S(t) \cdot E(t) \\\frac{I(t + \Delta t)}{\Delta t} &= \langle \alpha \rangle \cdot E(t) - \langle \gamma \rangle \cdot I(t) - \mu \cdot I(t) \\\frac{R(t + \Delta t)}{\Delta t} &= \langle \gamma \rangle \cdot I(t) \\\frac{D(t + \Delta t)}{\Delta t} &= \langle \mu \rangle \cdot I(t) \\ \\\beta &= 0.5 \\N &= S(t) + E(t) + I(t) + R(t) + D(t)\end{aligned}where $\beta$ is the transition  or  infectious  rate, $\gamma$ is the recovery rate and is determined by the average duration of an infection, $\alpha$ is the incubation rate or the reciprocal value of the incubation period. Since, for  the  case  of  COVID-19  disease,  an  exposed individual is infectious even in the incubation period, $\delta$ represents the direct transition rate between susceptible and exposed compartment. Finally, $\mu$ is the mortality rate.Each compartment is a function of time, where the sum of all compartments yields to the population number.Unknown epidemiological input parameters, $\delta$, $\alpha$, $\gamma$ and $\mu$, are all treated as random variables with expected values obtained through inverse system identification using L-BFGS-B optimization of the loss function defined as the difference between *a priori* known confirmed numbers of infectious, recovered and deceased individuals and associated predicted curves obtained via the SEIRD model. The variation coefficient (standard deviation or spread) is defined arbitrarily and is fixed to $0.2$ for the further analysis.<jupyter_code>import os
import datetime

from coropy.compartmental_models import SEIRDModel
from matplotlib import dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp

import utils
np.random.seed(42)  # for the reproducibility of results
utils.configure(grid=False)  # configure paper-quality plots
b = utils.default_colors('blue')
o = utils.default_colors('orange')
g = utils.default_colors('green')
r = utils.default_colors('red')
p = utils.default_colors('purple')
# all data
DATA_DIR = os.path.join(os.pardir, 'data')
FILE_NAME = '1st_wave_data_cro.csv' 
FILE_PATH = os.path.join(DATA_DIR, FILE_NAME)
df = pd.read_csv(FILE_PATH)

# derived data
cum_positives = np.cumsum(df.confirmed_positives.values)
cum_recovered = df.estimate_recovered.values
cum_deceased = np.cumsum(df.confirmed_deceased.values)
tests = df.total_tests.values
active = cum_positives - cum_recovered - cum_deceased
# initial conditions
S0 = cum_positives[-1] + cum_deceased[-1]
E0 = 10
I0 = 1
R0 = 0
D0 = 0
IC = [S0, E0, I0, R0, D0]

# dates
t = np.arange(cum_positives.size)
start_dt = datetime.datetime(2020, 2, 25)
end_dt = start_dt + datetime.timedelta(int(t.max())+1)
dates = mdates.drange(start_dt, end_dt, datetime.timedelta(days=1))

# transmission rate
beta = 0.5

# system of odes
def SEIRD_ode(t, y, delta, alpha, gamma, mu):
    S, E, I, R, D = y
    N = S + E + I + R + D
    return [
        -beta*S*I/N - delta*S*E, 
        beta*S*I/N - alpha*E + delta*S*E, 
        alpha*E - gamma*I - mu*I, 
        gamma*I,
        mu*I,
    ]

model = SEIRDModel(ode_system=SEIRD_ode)
params, loss = model.fit(cum_positives, cum_recovered, cum_deceased, IC, guess=[0.1, 0.1, 0.1, 0.1])
S, E, I, R, D = model.simulate()
# metropolis-hastings

n_iters = 1000
I_sampled = np.zeros(shape=(n_iters, t.size))

for i in range(n_iters):
    noisy_S0 = np.random.normal(loc=S0, scale=0.2*S0)
    noisy_E0 = np.random.normal(loc=E0, scale=0.2*E0)
    noisy_params = np.random.normal(loc=params, scale=0.2*params, size=params.size)
    noisy_sol = solve_ivp(
        SEIRD_ode,
        t_span=(0, int(t[-1])),
        y0=[noisy_S0, noisy_E0, I0, R0, D0],
        args=noisy_params,
        method='RK45',
        t_eval=t,
        vectorized=True,
    )
    if np.mean(noisy_sol.y[2] / I):
        I = noisy_sol.y[2]
    
    I_sampled[i, :] = I  
I_exp = np.mean(I_sampled, axis=0)
I_err = np.std(I_sampled, axis=0)
pcr_sens = 0.85
pcr_spec = 0.95
active_lb_ci, active_lb, active_ub, active_ub_ci = SEIRDModel.calculate_ci(
    pcr_sens,
    pcr_spec,
    df.confirmed_positives.values,
    cum_recovered + cum_deceased,
    tests)
sre_fit = np.sqrt(np.mean((active - I_exp)**2))
fig = plt.figure()

# active
eb = plt.errorbar(dates, active, 
                  yerr=np.r_[
                      (active-active_lb).reshape(1, -1),
                      (active_ub-active).reshape(1, -1)],
                  markevery=5,
                  errorevery=5, capsize=3,
                  marker='s', linestyle='None', mew=1, mfc='none',
                  c=p,
                  label=f'$I[t]$')
exp = plt.plot(dates, I_exp, '-', c=b, alpha=0.5, label=f'$\\langle \\hat{{I}}(t) \\rangle$')
plt.plot(dates, I_sampled.T, c=b, alpha=0.01)
plt.ylabel('$N$')
plt.legend()
_ = fig.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y/%m/%d'))
_ = fig.gca().xaxis.set_major_locator(mdates.DayLocator(bymonthday=[10, 25]))
_ = plt.gcf().autofmt_xdate()
plt.yticks([0, 1000, 2000, 3000])
plt.ylim([-250, 3000])
plt.show()
fname = f'MCMC_SEIRD_rigid'
fig.savefig(fname=os.path.join(os.pardir, 'figures', fname + '.png'),
            format='png', bbox_inches='tight', dpi=1000)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13339700732250875, 'fraction_numerical': 0.02037567653613499, 'mean_word_length': 3.7311746987951806, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '222459', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2321, 'n_tokens_neox': 2159, 'n_words': 622, 'path': 'notebooks/02_SEIRD_mcmc.ipynb', 'repo_name': 'antelk/cmj-covid-19-paper', 'size': '6282'}