Calibrate multidose national model (#285)

This PR contains several bugfixes of the multidose national model and contains a calibration using data until 31/08

data/raw/sciensano/COVID19BE_HOSP.csv

@@ -5829,3 +5829,14 @@ DATE,PROVINCE,REGION,NR_REPORTING,TOTAL_IN,TOTAL_IN_ICU,TOTAL_IN_RESP,TOTAL_IN_E
 2021-08-26,VlaamsBrabant,Flanders,6,22,4,1,0,2,4
 2021-08-26,BrabantWallon,Wallonia,2,6,1,1,0,2,0
 2021-08-26,WestVlaanderen,Flanders,11,40,12,5,2,3,4
+2021-08-27,Antwerpen,Flanders,14,81,23,12,3,10,9
+2021-08-27,Brussels,Brussels,15,188,69,37,10,17,23
+2021-08-27,Hainaut,Wallonia,14,71,21,11,2,4,9
+2021-08-27,Limburg,Flanders,7,25,8,4,0,5,3
+2021-08-27,Liège,Wallonia,12,80,19,11,1,6,5
+2021-08-27,Luxembourg,Wallonia,3,7,2,0,0,0,0
+2021-08-27,Namur,Wallonia,6,11,5,1,1,1,1
+2021-08-27,OostVlaanderen,Flanders,14,68,17,10,0,9,5
+2021-08-27,VlaamsBrabant,Flanders,6,22,5,0,0,1,1
+2021-08-27,BrabantWallon,Wallonia,2,5,1,1,0,0,2
+2021-08-27,WestVlaanderen,Flanders,11,39,12,5,2,5,5

data/raw/sciensano/COVID19BE_MORT.csv

@@ -6268,12 +6268,18 @@ DATE,REGION,AGEGROUP,SEX,DEATHS
 2021-08-23,Brussels,85+,F,1
 2021-08-23,Wallonia,85+,F,1
 2021-08-24,Flanders,75-84,M,1
+2021-08-24,Flanders,85+,M,1
 2021-08-24,Brussels,85+,F,1
 2021-08-24,Wallonia,65-74,M,1
 2021-08-24,Wallonia,85+,M,1
 2021-08-24,Wallonia,85+,F,1
 2021-08-25,Flanders,75-84,M,1
 2021-08-25,Flanders,75-84,F,3
 2021-08-25,Wallonia,75-84,F,1
+2021-08-26,Flanders,45-64,M,1
+2021-08-26,Flanders,85+,M,1
+2021-08-26,Brussels,75-84,F,1
+2021-08-26,Brussels,85+,F,1
 2021-08-26,Wallonia,85+,M,1
 2021-08-26,Wallonia,85+,F,1
+2021-08-27,Wallonia,85+,F,1

notebooks/calibration/calibrate-COVID19_SEIRD_stratified_vacc-WAVE2.py

@@ -89,7 +89,7 @@
 # Serological data
 df_sero_herzog, df_sero_sciensano = sciensano.get_serological_data()
 # VOC data
-df_VOC_501Y = VOC.get_501Y_data()
+df_VOC_abc = VOC.get_abc_data()
 # Model initial condition on September 1st
 with open('../../data/interim/model_parameters/COVID19_SEIRD/calibrations/national/initial_states_2020-09-01.json', 'r') as fp:
     initial_states = json.load(fp)    
@@ -116,7 +116,7 @@
 # ---------------------------
 
 from covid19model.models.time_dependant_parameter_fncs import make_VOC_function
-VOC_function = make_VOC_function(df_VOC_501Y)
+VOC_function = make_VOC_function() # no argument = use of logistic function for alpha-gamma variants instead of data
 
 # -----------------------------------
 # Time-dependant vaccination function
@@ -134,6 +134,15 @@
 contact_matrix_4prev = make_contact_matrix_function(df_google, Nc_all)
 policies_WAVE2 = make_contact_matrix_function(df_google, Nc_all).policies_WAVE2_no_relaxation
 
+# -----------------------------------
+# Time-dependant seasonality function
+# -----------------------------------
+
+def seasonality(t,states,param,amplitude, peak_shift):
+    maxdate = pd.Timedelta(days=peak_shift) + pd.to_datetime('2021-01-01')
+    t = (t - pd.to_datetime(maxdate))/pd.Timedelta(days=1)/365
+    return param*(1+amplitude*np.cos( 2*np.pi*(t)))
+
 # --------------------
 # Initialize the model
 # --------------------
@@ -154,16 +163,16 @@
 # Add the time-dependant parameter function arguments
 # Social policies
 params.update({'l': 21, 'prev_schools': 0, 'prev_work': 0.5, 'prev_rest_lockdown': 0.5, 'prev_rest_relaxation':0.5, 'prev_home': 0.5})
-# VOC
-params.update({'t_sig': '2021-06-21', 'k': 0.06}) # Infliction point from genomic surveillance data
 # Vaccination
 params.update(
     {'initN': initN, 'vacc_order': np.array(range(9))[::-1], 'daily_first_dose': 55000,
      'refusal': 0.2*np.ones([9,2]), 'delay_immunity': 14, 'delay_doses': 3*7, 'stop_idx': 8}
 )
+# Seasonality
+params.update({'amplitude': 0.1, 'peak_shift': 0})
 # Initialize model
 model = models.COVID19_SEIRD_stratified_vacc(initial_states, params,
-                        time_dependent_parameters={'Nc': policies_WAVE2, 'N_vacc': vacc_strategy, 'alpha': VOC_function})
+                        time_dependent_parameters={'beta': seasonality, 'Nc': policies_WAVE2, 'N_vacc': vacc_strategy, 'alpha': VOC_function})
 
 # -----------------------------------
 # Define calibration helper functions
@@ -312,13 +321,13 @@
 start_calibration = '2020-09-01'
 # Last datapoint used to calibrate compliance and prevention
 if not args.enddate:
-    end_calibration = '2021-08-22'
+    end_calibration = datetime.datetime.strftime(df_sciensano.index[-1], '%Y-%m-%d')
 else:
     end_calibration = str(args.enddate)
 # PSO settings
 processes = int(mp.cpu_count()/2)
-multiplier = 3
-maxiter = 30
+multiplier = 5
+maxiter = 1000
 popsize = multiplier*processes
 # MCMC settings
 max_n = 500000
@@ -335,22 +344,24 @@
 # Define dataset
 # --------------
 
-data=[df_sciensano['H_in'][start_calibration:end_calibration]]
-states = ["H_in"]
-weights = [1]
+data=[df_sciensano['H_in'][start_calibration:end_calibration],df_sciensano['H_in']['2021-06-15':end_calibration]]
+states = ["H_in", "H_in"]
+weights = [1, 1]
 
 # -----------
 # Perform PSO
 # -----------
 
 # optimisation settings
-pars = ['beta','da','l', 'prev_schools', 'prev_work', 'prev_rest_lockdown', 'prev_rest_relaxation', 'prev_home', 'K_inf1', 'K_inf2']
-bounds=((0.010,0.030),(2,14),(2,12),(0.02,0.98),(0.02,0.98),(0.02,0.98),(0.02,0.98),(0.02,0.98),(1.4,1.7),(2,2.35))
+pars = ['beta','da','l', 'prev_schools', 'prev_work', 'prev_rest_lockdown', 'prev_rest_relaxation', 'prev_home', 'K_inf1', 'K_inf2', 'amplitude', 'peak_shift']
+bounds=((0.005,0.030),(2,14),(2,12),(0.05,0.95),(0.05,0.95),(0.05,0.95),(0.05,0.95),(0.05,0.95),(1.4,1.7),(2.1,2.6),(0,0.25), (-31, 31))
 # run optimization
-theta = pso.fit_pso(model, data, pars, states, bounds, weights, maxiter=maxiter, popsize=popsize,
-                    start_date=start_calibration, warmup=warmup, processes=processes)
-#theta = np.array([0.0134, 8.32, 4.03, 0.687, 0.118, 0.105, 0.50, 0.70, 1.52, 2.20])
-#theta = np.array([0.01489179, 6.52556664, 3.32749332, 0.75299559, 0.05099117, 0.2546443, 0.72560745, 0.63643327, 1.53328335, 2.32212406]) #-253281.68302163907
+#theta = pso.fit_pso(model, data, pars, states, bounds, weights, maxiter=maxiter, popsize=popsize,
+#                    start_date=start_calibration, warmup=warmup, processes=processes)
+
+theta = np.array([1.21910243e-02, 1.27648358e+01, 7.43853204e+00, 6.46597486e-01,
+                    5.02862404e-02, 5.08402716e-02, 9.50000000e-01, 3.59405878e-01,
+                    1.64860956e+00, 2.28768733e+00, 1.41699867e-01, 1.93592372e+01]) #-260517.1618694332
 
 # Assign estimate
 model.parameters = assign_PSO(model.parameters, pars, theta)
@@ -382,28 +393,27 @@
 #density_da_norm = density_da/np.sum(density_da)
 
 # Setup uniform priors
-pars = ['beta', 'da', 'l', 'prev_schools', 'prev_work', 'prev_rest_lockdown', 'prev_rest_relaxation', 'prev_home','K_inf1', 'K_inf2']
-log_prior_fcn = [prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform]
-log_prior_fcn_args = [(0.001, 0.12), (0.01, 14), (0.1,14), (0.05,1), (0.05,1), (0.05,1), (0.05,1), (0.05,1),(1.3,1.8),(2.0,2.8)]
+pars = ['beta', 'da', 'l', 'prev_schools', 'prev_work', 'prev_rest_lockdown', 'prev_rest_relaxation', 'prev_home','K_inf1', 'K_inf2', 'amplitude', 'peak_shift']
+log_prior_fcn = [prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform, prior_uniform,  prior_uniform, prior_uniform]
+log_prior_fcn_args = [(0.001, 0.12), (0.01, 14), (0.1,14), (0.01,1), (0.01,1), (0.01,1), (0.01,1), (0.01,1),(1.3,1.65),(2.0,2.45),(0,0.30),(-31,31)]
 # Perturbate PSO Estimate
-pert = [5e-2, 5e-2, 5e-2, 100e-2, 100e-2, 100e-2, 100e-2, 100e-2, 10e-2, 10e-2]
+pert = [5e-2, 5e-2, 5e-2, 5e-2, 5e-2, 5e-2, 5e-2, 5e-2, 10e-2, 10e-2, 10e-2, 10e-2]
 ndim, nwalkers, pos = perturbate_PSO(theta, pert, 2)
 
-pos[:,3:8] = np.where(pos[:,3:8]<=0, 0, pos[:,3:8]) 
-pos[:,3:8] = np.where(pos[:,3:8]>=1, 1, pos[:,3:8]) 
+#pos[:,3:8] = np.where(pos[:,3:8]<=0, 0, pos[:,3:8]) 
+#pos[:,3:8] = np.where(pos[:,3:8]>=1, 1, pos[:,3:8]) 
 
 # Set up the sampler backend if needed
 if backend:
     filename = spatial_unit+'_R0_COMP_EFF_'+run_date
     backend = emcee.backends.HDFBackend(results_folder+filename)
     backend.reset(nwalkers, ndim)
 # Labels for traceplots
-labels = ['$\\beta$','$d_{a}$','$l$', '$\Omega_{schools}$', '$\Omega_{work}$', '$\Omega_{rest, lockdown}$', '$\Omega_{rest, relaxation}$', '$\Omega_{home}$', '$K_{inf,alpha}$', '$K_{inf,delta}$']
+labels = ['$\\beta$','$d_{a}$','$l$', '$\Omega_{schools}$', '$\Omega_{work}$', '$\Omega_{rest, lockdown}$', '$\Omega_{rest, relaxation}$', '$\Omega_{home}$', '$K_{inf,alpha}$', '$K_{inf,delta}$', '$\phi$', '$\\tau$']
 # Arguments of chosen objective function
 objective_fcn = objective_fcns.log_probability
 objective_fcn_args = (model, log_prior_fcn, log_prior_fcn_args, data, states, pars)
-objective_fcn_kwargs = {'weights': weights, 'start_date': start_calibration, 'warmup': warmup, 'poisson_offset': 0.1}
-
+objective_fcn_kwargs = {'weights': weights, 'start_date': start_calibration, 'warmup': warmup}
 
 # ----------------
 # Run MCMC sampler

notebooks/calibration/plot_fit_R0_COMP_EFF_WAVE2.py

@@ -91,7 +91,7 @@
 # Serological data
 df_sero_herzog, df_sero_sciensano = sciensano.get_serological_data()
 # VOC data
-df_VOC_501Y = VOC.get_501Y_data()
+df_VOC_abc = VOC.get_abc_data()
 # Start of data collection
 start_data = df_sciensano.idxmin()
 # Start of calibration warmup and beta
@@ -119,7 +119,7 @@
 # ---------------------------
 
 from covid19model.models.time_dependant_parameter_fncs import make_VOC_function
-VOC_function = make_VOC_function(df_VOC_501Y)
+VOC_function = make_VOC_function()
 
 # ---------------------------------------------------------
 # Time-dependant vaccination function and sampling function
@@ -147,6 +147,16 @@
 contact_matrix_4prev = make_contact_matrix_function(df_google, Nc_all)
 policies_WAVE2 = make_contact_matrix_function(df_google, Nc_all).policies_WAVE2_no_relaxation
 
+# -----------------------------------
+# Time-dependant seasonality function
+# -----------------------------------
+
+def seasonality(t,states,param,amplitude, peak_shift):
+    maxdate = pd.Timedelta(days=peak_shift) + pd.to_datetime('2021-01-01')
+    t = (t - pd.to_datetime(maxdate))/pd.Timedelta(days=1)/365
+    return param*(1+amplitude*np.cos( 2*np.pi*(t)))
+
+
 # ---------------------------------------------------
 # Function to add poisson draws and sampling function
 # ---------------------------------------------------
@@ -172,13 +182,13 @@
     params.update({'N_vacc': np.zeros([9,3])})
     # Social policies
     params.update({'l': 21, 'prev_schools': 0, 'prev_work': 0.5, 'prev_rest_lockdown': 0.5, 'prev_rest_relaxation':0.5, 'prev_home': 0.5})
+    # Seasonality
+    params.update({'amplitude': 0.1, 'peak_shift': 0})
 else:
     # Social policies
     params.update({'l': 21, 'prev_schools': 0, 'prev_work': 0.5, 'prev_rest': 0.5, 'prev_home': 0.5})
 
 # Add the remaining time-dependant parameter function arguments
-# VOC
-params.update({'t_sig': '2021-06-21', 'k': 0.06})
 # Vaccination
 params.update(
     {'initN': initN, 'vacc_order': np.array(range(9))[::-1], 'daily_first_dose': 55000,
@@ -187,7 +197,7 @@
 # Initialize model
 if args.vaccination_model == 'stratified':
     model = models.COVID19_SEIRD_stratified_vacc(initial_states, params,
-                        time_dependent_parameters={'Nc': policies_WAVE2, 'N_vacc': vacc_strategy, 'alpha':VOC_function})
+                        time_dependent_parameters={'beta': seasonality, 'Nc': policies_WAVE2, 'N_vacc': vacc_strategy, 'alpha':VOC_function})
 else:
     model = models.COVID19_SEIRD_vacc(initial_states, params,
                         time_dependent_parameters={'Nc': policies_WAVE2, 'N_vacc': vacc_strategy, 'alpha':VOC_function})

src/covid19model/data/VOC.py

@@ -2,14 +2,19 @@
 import datetime
 import pandas as pd
 
-def get_501Y_data():
+def get_abc_data():
     # Load and format B1.1.7 VOC data 
     rel_dir = os.path.join(os.path.dirname(__file__), '../../../data/raw/VOCs/sequencing_501YV1_501YV2_501YV3.csv')
-    df_VOC_501Y = pd.read_csv(rel_dir, parse_dates=True).set_index('collection_date', drop=True).drop(columns=['sampling_week','year', 'week'])
+    df_VOC_abc = pd.read_csv(rel_dir, parse_dates=True).set_index('collection_date', drop=True).drop(columns=['sampling_week','year', 'week'])
     # Converting the index as date
-    df_VOC_501Y.index = pd.to_datetime(df_VOC_501Y.index)
+    df_VOC_abc.index = pd.to_datetime(df_VOC_abc.index)
     # Extrapolate missing dates to obtain a continuous index
-    df_VOC_501Y = df_VOC_501Y.resample('D').interpolate('linear')
-    df_VOC_501Y['baselinesurv_f_501Y.V1_501Y.V2_501Y.V3'] = (df_VOC_501Y['baselinesurv_n_501Y.V1']+df_VOC_501Y['baselinesurv_n_501Y.V2']+df_VOC_501Y['baselinesurv_n_501Y.V3'])/df_VOC_501Y['baselinesurv_total_sequenced']
+    df_VOC_abc['baselinesurv_f_501Y.V1_501Y.V2_501Y.V3'] = (df_VOC_abc['baselinesurv_n_501Y.V1']+df_VOC_abc['baselinesurv_n_501Y.V2']+df_VOC_abc['baselinesurv_n_501Y.V3'])/df_VOC_abc['baselinesurv_total_sequenced']
     # Assign data to class
-    return df_VOC_501Y
+    return df_VOC_abc
+
+def get_delta_data():
+    #Copied from the molecular surveillance reported in the weekly Sciensano bulletins: https://covid-19.sciensano.be/sites/default/files/Covid19/COVID-19_Weekly_report_NL.pdf
+    dates = pd.date_range(start = '2021-06-16', end = '2021-08-18', freq='7D')
+    df_VOC_delta = pd.DataFrame(data=[0.284, 0.456, 0.675, 0.790, 0.902, 0.950, 0.960, 0.986, 0.992, 1], index=dates, columns=['delta'])
+    return df_VOC_delta

src/covid19model/data/model_parameters.py

@@ -318,11 +318,11 @@ def get_COVID19_SEIRD_parameters(age_stratified=True, spatial=None, vaccination=
         # vaccination
         if vaccination == True:
             pars_dict['N_vacc'] = np.zeros(9) # Default: no vaccination at simulation start
-            pars_dict['e_s'] = 0.95 # Default: 95% lower susceptibility to SARS-CoV-2 on a per contact basis
-            pars_dict['e_h'] = 0.90 # Default: 100% protection against severe COVID-19
-            pars_dict['e_a'] = 1.00 # Default: vaccination works in 100% of people
-            pars_dict['e_i'] = 0.90 # Default: vaccinated infectious individual is equally infectious as non-vaccinated individual
-            pars_dict['d_vacc'] = 36*30 # Default: 12 month coverage of vaccine
+            pars_dict['e_s'] = np.array([0.80, 0.80, 0.75]) # Default: 95% lower susceptibility to SARS-CoV-2 on a per contact basis
+            pars_dict['e_h'] = np.array([0.95, 0.95, 0.95]) # Default: 100% protection against severe COVID-19
+            pars_dict['e_a'] = 1.00*np.ones(3) # Default: vaccination works in 100% of people
+            pars_dict['e_i'] = 0.5*np.ones(3)# Default: vaccinated infectious individual is equally infectious as non-vaccinated individual
+            pars_dict['d_vacc'] = 36*30 # Default: 36 months coverage of vaccine
 
     else:
         pars_dict['Nc'] = np.array([17.65]) # Average interactions assuming weighing by age, by week/weekend and the inclusion of supplemental professional contacts (SPC)
@@ -413,7 +413,7 @@ def get_COVID19_SEIRD_parameters(age_stratified=True, spatial=None, vaccination=
         pars_dict['alpha'] = [1, 0, 0] # Must be a list so we can check if "sum(alpha) == 1" 
         pars_dict['K_inf1'] = 1.45 # British variant infectivity gain
         pars_dict['K_inf2'] = 1.45*1.5 # Indian variant infectivity gain
-        pars_dict['K_hosp'] = [1, 1.4, 1.4]
+        pars_dict['K_hosp'] = np.ones(3)
 
     return pars_dict
 

src/covid19model/models/models.py

@@ -464,22 +464,58 @@ def integrate(t, S, E, I, A, M, C, C_icurec, ICU, R, D, H_in, H_out, H_tot,
         # Compute the vaccination transitionings and waning of immunity
         # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+        dS = np.zeros(S.shape)
+        dE = np.zeros(E.shape)
+        dI = np.zeros(I.shape)
+        dA = np.zeros(A.shape)
+        dR = np.zeros(R.shape)
+
         # 0 doses --> 1 dose and 0 doses --> 2 doses
-        VE = S[:,0] + R[:,0]
+        VE = S[:,0] + E[:,0] + I[:,0] + A[:,0] + R[:,0]
+        #VE = np.where(VE<=0, 1e-3, VE)
         f_S = S[:,0]/VE
+        f_E = E[:,0]/VE
+        f_I = I[:,0]/VE
+        f_A = A[:,0]/VE
         f_R = R[:,0]/VE
-        S[:,0] = S[:,0] - N_vacc[:,0]*f_S -  N_vacc[:,2]*f_S 
-        R[:,0] = R[:,0] - N_vacc[:,0]*f_R - N_vacc[:,2]*f_R
-        S[:,1] = S[:,1] + N_vacc[:,0] # 0 --> 1 dose
-        S[:,2] = S[:,2] + N_vacc[:,2] # 0 --> 2 doses
+        dS[:,0] =  - N_vacc[:,0]*f_S -  N_vacc[:,2]*f_S 
+        dE[:,0] =  - N_vacc[:,0]*f_E -  N_vacc[:,2]*f_E
+        dI[:,0] =  - N_vacc[:,0]*f_I -  N_vacc[:,2]*f_I
+        dA[:,0] =  - N_vacc[:,0]*f_A -  N_vacc[:,2]*f_A
+        dR[:,0] = - N_vacc[:,0]*f_R - N_vacc[:,2]*f_R
+
+        dS[:,1] =  N_vacc[:,0]*f_S # 0 --> 1 dose
+        dE[:,1] =  N_vacc[:,0]*f_E # 0 --> 1 dose
+        dI[:,1] =  N_vacc[:,0]*f_I # 0 --> 1 dose
+        dA[:,1] =  N_vacc[:,0]*f_A # 0 --> 1 dose
+        dR[:,1] =  N_vacc[:,0]*f_R # 0 --> 1 dose
+
+        dS[:,2] =  N_vacc[:,2]*f_S # 0 --> 2 doses
+        dE[:,2] =  N_vacc[:,2]*f_E # 0 --> 2 doses
+        dI[:,2] =  N_vacc[:,2]*f_I # 0 --> 2 doses
+        dA[:,2] =  N_vacc[:,2]*f_A # 0 --> 2 doses
+        dR[:,2] =  N_vacc[:,2]*f_R # 0 --> 2 doses
+
         # 1 dose --> 2 doses
-        VE = S[:,1] + R[:,1]
-        VE = np.where(VE==0, 1, VE) 
+        VE = S[:,1] + E[:,1] + I[:,1] + A[:,1] + R[:,1]
+        #VE = np.where(VE<=0, 1e-3, VE)
         f_S = S[:,1]/VE
+        f_E = E[:,1]/VE
+        f_I = I[:,1]/VE
+        f_A = A[:,1]/VE
         f_R = R[:,1]/VE
-        S[:,1] = S[:,1] - N_vacc[:,1]*f_S
-        R[:,1] = R[:,1] - N_vacc[:,1]*f_R
-        S[:,2] = S[:,2] + N_vacc[:,1]
+
+        dS[:,1] = dS[:,1] - N_vacc[:,1]*f_S
+        dE[:,1] = dE[:,1] - N_vacc[:,1]*f_E
+        dI[:,1] = dI[:,1] - N_vacc[:,1]*f_I
+        dA[:,1] = dA[:,1] - N_vacc[:,1]*f_A
+        dR[:,1] = dR[:,1] - N_vacc[:,1]*f_R
+
+        dS[:,2] = dS[:,2] + N_vacc[:,1]*f_S
+        dE[:,2] = dE[:,2] + N_vacc[:,1]*f_E
+        dI[:,2] = dI[:,2] + N_vacc[:,1]*f_I
+        dA[:,2] = dA[:,2] + N_vacc[:,1]*f_A
+        dR[:,2] = dR[:,2] + N_vacc[:,1]*f_R
 
         # Compute infection pressure (IP) of all variants
         # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -490,16 +526,16 @@ def integrate(t, S, E, I, A, M, C, C_icurec, ICU, R, D, H_in, H_out, H_tot,
         # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         h_acc = (1-e_h)*h
 
-        dS  = - IP*S*(1-e_s)
-        dE  = IP*S*(1-e_s) - E/sigma 
-        dI = (1/sigma)*E - (1/omega)*I 
-        dA = (a/omega)*I - A/da
+        dS  = dS - IP*(S+dS)*(1-e_s)
+        dE  = dE + IP*(S+dS)*(1-e_s) - E/sigma 
+        dI = dI + (1/sigma)*E - (1/omega)*I 
+        dA = dA + (a/omega)*I - A/da
         dM = ((1-a)/omega)*I - M*((1-h_acc)/dm) - M*h_acc/dhospital
         dC = M*(h_acc/dhospital)*c - (1-m_C)*C*(1/(dc_R)) - m_C*C*(1/(dc_D))
         dICUstar = M*(h_acc/dhospital)*(1-c) - (1-m_ICU)*ICU/(dICU_R) - m_ICU*ICU/(dICU_D)
 
         dC_icurec = (1-m_ICU)*ICU/(dICU_R) - C_icurec*(1/dICUrec)
-        dR  = A/da + ((1-h_acc)/dm)*M + (1-m_C)*C*(1/(dc_R)) + C_icurec*(1/dICUrec)
+        dR  = dR + A/da + ((1-h_acc)/dm)*M + (1-m_C)*C*(1/(dc_R)) + C_icurec*(1/dICUrec)
         dD  = (m_ICU/(dICU_D))*ICU + (m_C/(dc_D))*C 
         dH_in = M*(h_acc/dhospital) - H_in
         dH_out =  (1-m_C)*C*(1/(dc_R)) +  m_C*C*(1/(dc_D)) + m_ICU/(dICU_D)*ICU + C_icurec*(1/dICUrec) - H_out