Epidemic modeling  Part 2
A deterministic numerical SEIR model
 Motivation for writeup
 Recall SEIR model equations
 Numerical solution to this deteministic population level model
 Qualitative analysis of $\beta$, $\sigma$, and $\gamma$
 Discussion
This is the 2nd part of a multipart series blog post on modeling in epidemiology.
The COVID19 pandemic has brought a lot of attention to the study of epidemiology and more specifically to the various mathematical models that are used to inform public health policies. Everyone has been trying to understand the growth or slowing of new cases and trying to predict the necessary sanitary resources. This blog post attempts to explain the foundations for some of the most used models and enlighten the reader on two key points.
After introducing the concepts of compartmentalization and disease dynamics in the first blog post, this second part is focused on developing a deterministic numerical solution for the SEIR model discussed there.
While normally the goal is to use realworld data to infer characteristics of the underlying disease (as will be done in later blog posts), here we want to use simulate the spread of a COVID19 like disease in a population of 10000, and look at the effects of the different parameters on the spread.
See the first blog post for derivation.

Continuoustime:
 $\frac{ds(t)}{dt}=\beta i(t) s(t)$
 $\frac{de(t)}{dt}=\beta i(t) s(t)  \sigma e(t)$
 $\frac{di(t)}{dt}=\sigma e(t)  \gamma i(t)$
 $\frac{dr(t)}{dt}=\gamma i(t)$

Discretetime:
 $\Delta S = \beta I S \Delta T$
 $\Delta E = (\beta I S \sigma E) \Delta T$
 $\Delta I = (\sigma E  \gamma I) \Delta T$
 $\Delta R = \gamma I \Delta T$
Coding the SEIR model
To build the SEIR model we simply use the discretetime set of equations above.
The model will thus take as input the following:
 Initial proportion of S, E, I, and R in the population
 $\beta$ parameter pertaining to the population in question
 $\sigma$ and $\gamma$ parameters pertaining to the disease
 Numbers of days to run the simulation
#collapse_hide
# Import required libraries
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
# Let's build a numerical solution
def seir_model(init, parms, days):
S_0, E_0, I_0, R_0 = init
Epd, Ipd, Rpd = [0], [0], [0]
S, E, I, R = [S_0], [E_0], [I_0], [R_0]
dt=0.1
t = np.linspace(0,days,int(days/dt))
sigma, beta, gam = parms
for _ in t[1:]:
next_S = S[1]  beta*S[1]*I[1]*dt
Epd.append(beta*S[1]*I[1]*dt)
next_E = E[1] + (beta*S[1]*I[1]  sigma*E[1])*dt
Ipd.append(sigma*E[1]*dt)
next_I = I[1] + (sigma*E[1]  gam*I[1])*dt
Rpd.append(gam*I[1]*dt)
next_R = R[1] + (gam*I[1])*dt
S.append(next_S)
E.append(next_E)
I.append(next_I)
R.append(next_R)
return np.stack([S, E, I, R, Epd, Ipd, Rpd]).T
Simulation parameters used for plot below:
 Days = 100
 Population = 10000
 Number of susceptible people on day 0 = 9999
 Number of exposed people on day 0 = 1
 No infected or recovered people on day 0
A lot of research is ongoing into the COVID19 characteristics of $\beta$, $\sigma$, and $\gamma$.
However, these are complex studies that require a lot of data and so far we have little information to go on.
The literature suggests the following:
 $\underline{T_{Incubation}}$:
The mean is 56 days but it can range anywhere from 214 days ^{1} ^{2}
Another paper reports a mean incubation period of 5.2 days and the 95th percentile at 12.5 days ^{3}.
There are reports of presymptomatic infections[^2], but these are reportedly rare [^1] so in the following models we will assume: $$T_{Incubation} = T_{Latent}$$ And so: $$\sigma = \frac{1}{5.2} days^{1}$$
 $\underline{T_{Infectious}}$:
Again it is very difficult to say for sure and the period of communicability is very uncertain for COVID19.
Research suggests a median of 20 days of viral shedding after onset of symptoms ^{4}.
Ranging from 8 to 37 days in survivors.
While it is noted PCR positivity does not necessarily reflect the infectious period (virus may not be viable but the PCR amplification will result in a positive), for the purpose of this blog post we will assume the following: $$T_{Infectious} = T_{Clinical}$$ To obtain an exponential distribution with median M, the scale A is calculated as follows: $$A = \frac{M}{\ln2} = \frac{20}{\ln2}$$ This results in $$\gamma = \frac{\ln2}{20} = \frac{1}{28.85}\ days^{1}$$
 $\underline{Beta= \beta}$:
While difficult to estimate this parameter as there is a lot of variation between countries, cultures, societal norms, etc.. a little thought experiment can help us evaluate the value for $\beta = r\rho$ in Switerland or France for example.
If no control measures are put in place and people do not change habits (as is the case in this blog post), we can expect the following:
 Average number of contacts per day:
$$r = 10\ contacts\ per\ day$$
 Average probability of transmission from contact:
$$\rho = 5\%$$
And so: $$\beta = r\rho = 0.5$$
#collapse_hide
#Define parameters
days = 200
N = 10000
init = 1  1/N, 1/N, 0, 0
sigma = 1/5.2
beta = 0.5
gam = 1/28.85
parms = sigma, beta, gam
# Run simulation
results_avg = seir_model(init, parms, days)
#collapse_hide
fig = go.Figure(data=[
go.Scatter(name='S', x=np.linspace(0,days,days*10), y=results_avg.T[0], line={'dash':'solid', 'color':'blue'}),
go.Scatter(name='E', x=np.linspace(0,days,days*10), y=results_avg.T[1], line={'dash':'solid', 'color':'yellow'}),
go.Scatter(name='I', x=np.linspace(0,days,days*10), y=results_avg.T[2], line={'dash':'solid', 'color':'red'}),
go.Scatter(name='R', x=np.linspace(0,days,days*10), y=results_avg.T[3], line={'dash':'solid', 'color':'green'}),
])
fig.update_layout(
xaxis_title = 'Day',
yaxis_title = 'Proportion of population',
title={
'text':'Deterministic SEIR model  COVID19 parameters',
'x':0.5,
'xanchor':'center'
}
)
fig.show()
#collapse_hide
## Let's try to see how the model changes
days = 1000
N = 10000
init = 1  1/N, 1/N, 0, 0
sigma_high = 1 # 1 > Average 1 day from E > I (ressembles SIR model)
sigma_low = 1/100 #10 days on average, twice as long as COVID19
sigma_covid = 1/5.2
beta = 0.5
gam = 1/28.85
parms_fastEI = sigma_high, beta, gam
parms_slowEI = sigma_low, beta, gam
parms_avg = sigma_covid, beta, gam
# Run simulation
results_fastEtoI = seir_model(init, parms_fastEI, days)
results_slowEtoI = seir_model(init, parms_slowEI, days)
results_avg = seir_model(init, parms_avg, days)
#collapse_hide
fig = go.Figure(data=[
go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[0], line={'dash':'solid', 'color':'blue'}, legendgroup="COVID", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[1], line={'dash':'solid', 'color':'yellow'}, legendgroup="COVID", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[2], line={'dash':'solid', 'color':'red'}, legendgroup="COVID", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[3], line={'dash':'solid', 'color':'green'}, legendgroup="COVID", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{high}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[0], line={'dash':'dash','color':'blue'}, legendgroup="high", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{high}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[1], line={'dash':'dash', 'color':'yellow'}, legendgroup="high", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{high}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[2], line={'dash':'dash', 'color':'red'}, legendgroup="high", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{high}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[3], line={'dash':'dash', 'color':'green'}, legendgroup="high", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{low}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[0], line={'dash':'dot', 'color':'blue'}, legendgroup="slow", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{low}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[1], line={'dash':'dot', 'color':'yellow'}, legendgroup="slow", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{low}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[2], line={'dash':'dot', 'color':'red'}, legendgroup="slow", hoverinfo='x+y'),
go.Scatter(name=r'$S:\sigma_{low}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[3], line={'dash':'dot', 'color':'green'}, legendgroup="slow", hoverinfo='x+y'),
])
fig.update_layout(
xaxis_title = 'Day',
yaxis_title = 'Proportion of population',
title={
'text':r'$\text{Effect of } \sigma \ \text{on Deterministic SEIR model}$',
'x':0.5,
'xanchor':'center'
}
)
fig.show()