Basic Experiments Workflow for Simple Epidemiological Models

Anton Antonov
MathematicaForPrediction at WordPress​
​SystemModeling at GitHub
March 2020

Introduction

The primary purpose of this notebook is to give a “stencil workflow” for simulations using the packages in the project Coronavirus simulation dynamics [AAr1].
The model in this notebook—SEI2R—differs from the classical SEIR model with the following elements:
Two separate infected populations: one is "severely symptomatic"; the other is "normally symptomatic."
1.
The monetary equivalent of lost productivity due to infected or deceased people is tracked.
2.
Remark: We consider the coronavirus propagation models as instances of the more general System Dynamics (SD) models.
Remark: The SEI2R model is a modification of the classic epidemic model SEIR [Wk1].
Remark: The interactive interfaces in the notebook can be used for attempts to calibrate SEI2R with real data. (For example, data for the 2019–2020 coronavirus outbreak [WRI1].)

Workflow

Get one of the classical epidemiology models.
1.
Extend the equations of the model if needed or desired.
2.
Set relevant initial conditions for the populations.
3.
Pick model parameters to be adjusted and “played with.”
4.
Derive parametrized solutions of model’s system of equations (ODEs or DAEs).
5.
Using the parameters of the previous step.
5.1.
Using an interactive interface, experiment with different values of the parameters.
6.
In order to form “qualitative understanding.”
6.1.
Get real-life data.
7.
Say, for the 2019–2020 coronavirus outbreak.
7.1.
Attempt manual or automatic calibration of the model.
8.
This step will most likely require additional data transformations and programming.
8.1.
Only manual calibration is shown in this notebook.
8.2.

Load Packages of the Framework

The epidemiological models framework used in this notebook is implemented in the packages [AAp1, AAp2]; the interactive plots functions are from the package [AAp3].
Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModels.m"]​​Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelModifications.m"]​​Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/WL/SystemDynamicsInteractiveInterfacesFunctions.m"]
In[]:=

Getting the Model Code

Here we take the SEI2R model implemented in the package EpidemiologyModels.m [AAp1]:
modelSI2R=SEI2RModel[t,"InitialConditions"True,"RateRules"True];
In[]:=
We can show a tabulated visualization of the model using the function ModelGridTableForm from [AAp1]:
ModelGridTableForm[modelSI2R]
In[]:=
Stocks
#
Symbol
Description
1
TP[t]
Total Population
2
SP[t]
Susceptible Population
3
EP[t]
Exposed Population
4
INSP[t]
Infected Normally Symptomatic Population
5
ISSP[t]
Infected Severely Symptomatic Population
6
RP[t]
Recovered Population
7
MLP[t]
Money of Lost Productivity
,Rates
#
Symbol
Description
1
μ[TP]
Population death rate
2
μ[INSP]
Infected Normally Symptomatic Population death rate
3
μ[ISSP]
Infected Severely Symptomatic Population death rate
4
sspf[SP]
Severely Symptomatic Population Fraction
5
β[INSP]
Contact rate for the normally symptomatic population
6
β[ISSP]
Contact rate for the severely symptomatic population
7
aip
Average infectious period
8
aincp
Average incubation period
9
lpcr[ISSP,INSP]
Lost productivity cost rate (per person per day)
,Equations
#
Equation
1
′
SP
t-
SPt
μTP
-
INSP[t]
SP[t]
β[INSP]
TP[0]
-
ISSP[t]
SP[t]
β[ISSP]
TP[0]
2
′
EP
t-
EPt
1
aincp
+
μTP
+
INSP[t]
SP[t]
β[INSP]
TP[0]
+
ISSP[t]
SP[t]
β[ISSP]
TP[0]
3
′
INSP
t-
INSP[t]
aip
+
EP[t]
(1-
sspf[SP]
)
aincp
-
INSPt
μINSP
4
′
ISSP
t-
ISSP[t]
aip
+
EP[t]
sspf[SP]
aincp
-
ISSPt
μISSP
5
′
RP
t
INSP[t]
+
ISSP[t]
aip
-
RPt
μTP
6
′
MLP
[t]
lpcr[ISSP,INSP]
(-
RP[t]
-
SP[t]
+TP[0])
,RateRules
#
Symbol
Value
1
TP[0]
100000
2
μTP
1
45625
3
μISSP
0.035
aip
4
μINSP
0.01
aip
5
β[ISSP]
0.15
6
β[INSP]
0.15
7
aip
26
8
aincp
6
9
sspf[SP]
0.2
10
lpcr[ISSP,INSP]
1
,InitialConditions
#
Equation
1
SP[0]99998
2
EP[0]0
3
ISSP[0]1
4
INSP[0]1
5
RP[0]0
6
MLP[0]0

Out[]=

Model Extensions and New Models

The framework implemented with the packages [AAp1, AAp2, AAp3] can be utilized using custom-made data structures that follow the structure of the models in [AAp1].
Of course, we can also just extend the models from [AAp1]. In this section, we show how SEI2R can be extended in two ways:
By adding a birthrate to the Susceptible Population equation (the birthrate is not included by default).
1.
By adding a new equation for the infected deceased population.
2.

Adding Births Term

Here are the equations of SEI2R (from [AAp1]):
ModelGridTableForm[modelSI2R]["Equations"]
In[]:=
#
Equation
1
′
SP
t-
SPt
μTP
-
INSP[t]
SP[t]
β[INSP]
TP[0]
-
ISSP[t]
SP[t]
β[ISSP]
TP[0]
2
′
EP
t-
EPt
1
aincp
+
μTP
+
INSP[t]
SP[t]
β[INSP]
TP[0]
+
ISSP[t]
SP[t]
β[ISSP]
TP[0]
3
′
INSP
t-
INSP[t]
aip
+
EP[t]
(1-
sspf[SP]
)
aincp
-
INSPt
μINSP
4
′
ISSP
t-
ISSP[t]
aip
+
EP[t]
sspf[SP]
aincp
-
ISSPt
μISSP
5
′
RP
t
INSP[t]
+
ISSP[t]
aip
-
RPt
μTP
6
′
MLP
[t]
lpcr[ISSP,INSP]
(-
RP[t]
-
SP[t]
+TP[0])
Out[]=
Here we find the position of the equation that corresponds to “Susceptible Population”:
pos=EquationPosition[modelSI2R["Equations"],First@GetPopulationSymbols[modelSI2R,"Susceptible Population"]]
In[]:=
1
Out[]=
Here we make the births term using a birthrate that is the same as the death rate:
birthTerm=GetPopulations[modelSI2R,"Total Population"]〚1〛*GetRates[modelSI2R,"Population death rate"]〚1〛
In[]:=
TP[t]μ[TP]
Out[]=
Here we add the births term to the equations of the new model:
modelSI2RNew=modelSI2R;​​modelSI2RNew["Equations"]=ReplacePart[modelSI2R["Equations"],{1,2}modelSI2R["Equations"]〚1,2〛+birthTerm];
In[]:=
Here we display the equations of the new model:
ModelGridTableForm[modelSI2RNew]["Equations"]
In[]:=
#
Equation
1
′
SP
t-
SPt
μTP
+
TPt
μTP
-
INSP[t]
SP[t]
β[INSP]
TP[0]
-
ISSP[t]
SP[t]
β[ISSP]
TP[0]
2
′
EP
t-
EPt
1
aincp
+
μTP
+
INSP[t]
SP[t]
β[INSP]
TP[0]
+
ISSP[t]
SP[t]
β[ISSP]
TP[0]
3
′
INSP
t-
INSP[t]
aip
+
EP[t]
(1-
sspf[SP]
)
aincp
-
INSPt
μINSP
4
′
ISSP
t-
ISSP[t]
aip
+
EP[t]
sspf[SP]
aincp
-
ISSPt
μISSP
5
′
RP
t
INSP[t]
+
ISSP[t]
aip
-
RPt
μTP
6
′
MLP
[t]
lpcr[ISSP,INSP]
(-
RP[t]
-
SP[t]
+TP[0])
Out[]=

Adding Infected Deceased Population Equation

Here we add new population, equation and initial condition, which allow for tracking the deaths because of infection:
AppendTo[modelSI2R["Equations"],IDP'[t]μ[INSP]*INSP[t]+μ[ISSP]*ISSP[t]];​​AppendTo[modelSI2R["Stocks"],IDP[t]"Infected Deceased Population"];​​AppendTo[modelSI2R["InitialConditions"],IDP[0]0];
In[]:=
Here is how the model looks:
ModelGridTableForm[modelSI2R]["Equations"]
In[]:=
#
Equation
1
′
SP
t-
SPt
μTP
-
INSP[t]
SP[t]
β[INSP]
TP[0]
-
ISSP[t]
SP[t]
β[ISSP]
TP[0]
2
′
EP
t-
EPt
1
aincp
+
μTP
+
INSP[t]
SP[t]
β[INSP]
TP[0]
+
ISSP[t]
SP[t]
β[ISSP]
TP[0]
3
′
INSP
t-
INSP[t]
aip
+
EP[t]
(1-
sspf[SP]
)
aincp
-
INSPt
μINSP
4
′
ISSP
t-
ISSP[t]
aip
+
EP[t]
sspf[SP]
aincp
-
ISSPt
μISSP
5
′
RP
t
INSP[t]
+
ISSP[t]
aip
-
RPt
μTP
6
′
MLP
[t]
lpcr[ISSP,INSP]
(-
RP[t]
-
SP[t]
+TP[0])
7
′
IDP
[t]
INSP[t]
μ[INSP]
+
ISSP[t]
μ[ISSP]
Out[]=

Parameters and Actual Simulation Equations Code

Here are the parameters we want to experiment with (or do calibration with):
lsFocusParams={aincp,aip,sspf[SP],β[ISSP],β[INSP]};
In[]:=
Here we set custom rates and initial conditions:
population=58160000/400;​​modelSI2R=SetRateRules[modelSI2R,<|TP[0]population|>];​​modelSI2R=SetInitialConditions[modelSI2R,<|SP[0]population-1,ISSP[0]0,INSP[0]1|>];
In[]:=
Here is the system of ODEs we use to do parametrized simulations:
lsActualEquations=Join[modelSI2R["Equations"]//.KeyDrop[modelSI2R["RateRules"],lsFocusParams],modelSI2R["InitialConditions"]];​​ResourceFunction["GridTableForm"][List/@lsActualEquations]
In[]:=
#
1
1
′
SP
t-
SP[t]
45625
-
INSP[t]SP[t]β[INSP]
145400
-
ISSP[t]SP[t]β[ISSP]
145400
2
′
EP
t-
1
45625
+
1
aincp
EPt+
INSP[t]SP[t]β[INSP]
145400
+
ISSP[t]SP[t]β[ISSP]
145400
3
′
INSP
t-
1.01INSP[t]
aip
+
EP[t](1-sspf[SP])
aincp
4
′
ISSP
t-
1.035ISSP[t]
aip
+
EP[t]sspf[SP]
aincp
5
′
RP
t
INSP[t]+ISSP[t]
aip
-
RP[t]
45625
6
′
MLP
[t]145400-RP[t]-SP[t]
7
′
IDP
t
0.01INSP[t]
aip
+
0.035ISSP[t]
aip
8
SP[0]145399
9
EP[0]0
10
ISSP[0]0
11
INSP[0]1
12
RP[0]0
13
MLP[0]0
14
IDP[0]0
Out[]=

Simulation

Straightforward simulation for one year using ParametricNDSolve:
aSol=​​Association@Flatten@​​ParametricNDSolve[lsActualEquations,Head/@Keys[modelSI2R["Stocks"]],{t,0,365},lsFocusParams]
In[]:=
TPParametricFunction
Expression: TP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}
,SPParametricFunction
Expression: SP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}
,EPParametricFunction
Expression: EP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}
,INSPParametricFunction
Expression: INSP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}
,ISSPParametricFunction
Expression: ISSP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}
,RPParametricFunction
Expression: RP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}
,MLPParametricFunction
Expression: MLP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}
,IDPParametricFunction
Expression: IDP
Parameters: {aincp,aip,sspf[SP],β[ISSP],β[INSP]}

Out[]=
(The advantage of having parametrized solutions is that we can quickly compute simulation results with new parameter values without solving the model’s system of ODEs; see the interfaces below.)

Interactive Interface

opts={PlotRangeAll,PlotLegendsNone,PlotTheme"Detailed",PerformanceGoal"Speed",ImageSize300};​​lsPopulationKeys=GetPopulationSymbols[modelSI2R,__~~"Population"];​​lsEconKeys={MLP};​​Manipulate[​​DynamicModule[{lsPopulationPlots,lsEconPlots,lsRestPlots},​​​​lsPopulationPlots=​​ParametricSolutionsPlots[​​modelSI2R["Stocks"],​​KeyTake[aSol,lsPopulationKeys],​​{aincp,aip,spf,crisp,criap},ndays,​​"LogPlot"popLogPlotQ,"Together"popTogetherQ,"Derivatives"popDerivativesQ,"DerivativePrefix""Δ",opts];​​​​lsEconPlots=​​ParametricSolutionsPlots[​​modelSI2R["Stocks"],​​KeyTake[aSol,lsEconKeys],​​{aincp,aip,spf,crisp,criap},ndays,​​"LogPlot"econLogPlotQ,"Together"econTogetherQ,"Derivatives"econDerivativesQ,"DerivativePrefix""Δ",opts];​​​​lsRestPlots=​​ParametricSolutionsPlots[​​modelSI2R["Stocks"],​​KeyDrop[aSol,Join[lsPopulationKeys,lsEconKeys]],​​{aincp,aip,spf,crisp,criap},ndays,​​"LogPlot"econLogPlotQ,"Together"econTogetherQ,"Derivatives"econDerivativesQ,"DerivativePrefix""Δ",opts];​​​​Multicolumn[Join[lsPopulationPlots,lsEconPlots,lsRestPlots],nPlotColumns,DividersAll,FrameStyleGrayLevel[0.8]]​​],​​{{aincp,6.,"Average incubation period (days)"},1,60.,1,Appearance{"Open"}},​​{{aip,32.,"Average infectious period (days)"},1,100.,1,Appearance{"Open"}},​​{{spf,0.2,"Severely symptom