<jupyter_start><jupyter_text>Test runs of the SEIR model for Belgium - coronavirus**Author:** <NAME>  **Date:** 14 march 2020  **Licence:** [CC-BY](https://creativecommons.org/licenses/by/4.0/)*Warning:* I am not an epidemiologist. I wrote this notebook to have a startingpoint for comparing estimates that I found online, for instance the one by NicolasVandewalle show below.I used a standard epidemiological model (SEIR) along with parameters from a recentarticle by a team of chinese scientists for the case of Wuhan.The absolute values are likely incorrect. The trends are good to know.![](figure_vdwnico_20200308.jpeg)Source: [@vdwnico](https://twitter.com/vdwnico/status/1236542044685959168) on Twitter, 8 march 2020.  Data from Belgian governmentBelow, I write down the number of new cases in Belgium since march 2.`data_I` is the total of infected cases (there are no recoveries recorded yet).<jupyter_code>new_I = [2, 6, 5, 10, 27, 59, 60, 31, 39, 28, 47, 85, 160, 133, 197, 172] # starting march 1
data_I = np.cumsum(new_I)
data_dates = np.arange('2020-03-01', len(new_I), dtype=np.datetime64)<jupyter_output><empty_output><jupyter_text>For estimating the healthcare network saturation, I use the estimate of Philippe Devos regarding theavailability of intensive care beds and the estimate of 6.1% of infected people needing intensive care(no other information is taken from that article).Article by <NAME> in La Libre Belgique, 2 march 2020:https://www.lalibre.be/debats/opinions/coronavirus-sans-mesures-de-precaution-drastiques-on-risque-d-avoir-850-000-personnes-infectees-et-50-000-morts-en-belgique-5e5cf60f9978e23106a0bfd9<jupyter_code>available_icu_beds = 140
fraction_of_icu_need = 0.061
saturation_of_icu_beds = 140/fraction_of_icu_need
print("Number of cases to saturate the ICU beds", saturation_of_icu_beds)<jupyter_output><empty_output><jupyter_text>Basic logarithmic plot and trendBelow, I plot the data on a lin-log plot to assess the exponentialcharacter of the number of confirmed infected cases.I superimpose an exponential trend for the coming month.<jupyter_code>plt.figure()

trend_start, trend_stop = 5, len(data_I)-1

plt.plot(data_dates, data_I, marker='o', label='reported I')

more_dates = np.arange(data_dates[0], data_dates[0]+30)
delta = np.array(more_dates - data_dates[trend_start], dtype=float)

current_exponential_rate = np.log(data_I[trend_stop]/data_I[trend_start]) / np.array(data_dates[trend_stop]-data_dates[trend_start], dtype=float)

plt.plot(more_dates, data_I[trend_start]*np.exp(delta*current_exponential_rate), marker='v')

plt.axhline(saturation_of_icu_beds)
plt.yscale('log')
plt.grid()

print(f"Current exponential rate {current_exponential_rate}")
print(f"Multiplication per day {np.exp(current_exponential_rate)}")
ax = plt.gca()
plt.setp(ax.get_xticklabels(), rotation=30, ha="right");<jupyter_output><empty_output><jupyter_text>Definition of the model and parametersI take the SEIR model defined in Wang et al. Cell Discovery (2020) 6:10 https://doi.org/10.1038/s41421-020-0148-0<jupyter_code># parameters for the model, page 6

σ = 1/5.2 # inverse of mean period
γ = 1/18  # time for which a case should be recovered

N = 11e6 / 3  # population size

# definition of model page 6
def SEIR(x, t, R0):
    S, E, I, R = x
    β = R0*γ
    return [
            - β*I*S/N,
            β*I*S/N - σ*E,
            σ*E - γ*I,
            γ*I,
           ]<jupyter_output><empty_output><jupyter_text>Estimate for BelgiumI start with the record of 2 infected cases on march 2 at https://www.info-coronavirus.be/fr/2020/03/01/deuxieme-infection-au-covid-19-en-belgique/ and list the new cases day per day (see above).The "day 0" is thus march 2 and the initial number of infected people is 2.To circumvent the low statistics of the first few days, I allow to delay the initial conditionof the SEIR model by an offset of `start_day`, taking the corresponding cumulated number of casesfor I0.<jupyter_code>start_day = 1
I0 = data_I[start_day]
E0 = 20*I0 # number of initially exposed people, see Wang et al

t = np.arange(0, 200, dtype=float) # running for 100 days
t_display = np.arange(data_dates[0], len(t), dtype=np.datetime64) + start_day

plt.figure()
plt.plot(data_dates, data_I, marker='o', label='reported I')

for R0 in [1.9, 3.1, 4.5, 7]:
    model = odeint(SEIR, [N-I0-E0, E0, I0, 0], t, args=(R0,))
    plt.plot(t_display, model[:,2], label=f'model I. R0={R0}')

plt.yscale('log')

plt.xlabel('time')

plt.axhline(saturation_of_icu_beds, label='saturation', c='k', ls='--')

ylo, yhi = plt.ylim()
plt.ylim(data_I[0], yhi)

plt.legend()
plt.grid()
ax = plt.gca()
plt.setp(ax.get_xticklabels(), rotation=30, ha="right");
start_day = 1
I0 = data_I[start_day]
E0 = 20*I0 # number of initially exposed people, see Wang et al

t = np.arange(0, 30, dtype=float) # running for 100 days
t_display = np.arange(data_dates[0], len(t), dtype=np.datetime64) + start_day

plt.figure()
plt.plot(data_dates, data_I, marker='o', label='reported I')

for R0 in [1.9, 3.1, 4.5, 7]:
    model = odeint(SEIR, [N-I0-E0, E0, I0, 0], t, args=(R0,))
    plt.plot(t_display, model[:,2], label=f'model I. R0={R0}')

plt.yscale('log')

plt.xlabel('time')

plt.axhline(saturation_of_icu_beds, label='saturation', c='k', ls='--')

ylo, yhi = plt.ylim()
plt.ylim(data_I[0], yhi)

plt.legend()
plt.grid()
ax = plt.gca()
plt.setp(ax.get_xticklabels(), rotation=30, ha="right");<jupyter_output><empty_output><jupyter_text>Risk of under-reportingAgain, as a test, I use below twice the number of reported cases tounderstand the risk related to under-reporting.If the actual number of cases is twice as high, the peak infection arisesearlier (mid-end july instead of mid-august) and the saturation of ICU bedsoccurs late march instead of mid-april.<jupyter_code>I0 = 2*data_I[start_day]
E0 = 20*I0

t = np.arange(0, 200, dtype=float) # running for 100 days
t_display = np.arange(data_dates[0], len(t), dtype=np.datetime64) + start_day

plt.figure()
plt.plot(data_dates, 2*data_I, marker='o', label='reported I')

for R0 in [1.9, 3.1, 4.5]:
    model = odeint(SEIR, [N-I0-E0, E0, I0, 0], t, args=(R0,))
    plt.plot(t_display, model[:,2], label=f'model I. R0={R0}')

plt.yscale('log')

plt.xlabel('time')

plt.axhline(saturation_of_icu_beds, label='saturation', c='k', ls='--')

plt.legend()
plt.grid()

ax = plt.gca()
plt.setp(ax.get_xticklabels(), rotation=30, ha="right");
start_day = 3
I0 = 2*data_I[start_day]
E0 = 20*I0

t = np.arange(0, 30, dtype=float) # running for 100 days
t_display = np.arange(data_dates[0], len(t), dtype=np.datetime64) + start_day

plt.figure()
plt.plot(data_dates, 2*data_I, marker='o', label='reported I')

for R0 in [1.9, 3.1, 4.5]:
    model = odeint(SEIR, [N-I0-E0, E0, I0, 0], t, args=(R0,))
    plt.plot(t_display, model[:,2], label=f'model I. R0={R0}')

plt.yscale('log')

plt.xlabel('time')

plt.axhline(saturation_of_icu_beds, label='saturation', c='k', ls='--')

plt.legend()
plt.grid()

ax = plt.gca()
plt.setp(ax.get_xticklabels(), rotation=30, ha="right");<jupyter_output><empty_output>