Basic Experiments Workflow for Simple Epidemiological Models
Basic Experiments Workflow for Simple Epidemiological Models
Anton Antonov
March 2020
Introduction
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].
1
.Two separate infected populations: one is "severely symptomatic"; the other is "normally symptomatic."
2
.The monetary equivalent of lost productivity due to infected or deceased people is tracked.
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
Workflow
1
.Get one of the classical epidemiology models.
2
.Extend the equations of the model if needed or desired.
3
.Set relevant initial conditions for the populations.
4
.Pick model parameters to be adjusted and “played with.”
5
.Derive parametrized solutions of model’s system of equations (ODEs or DAEs).
5
.1
.Using the parameters of the previous step.
6
.Using an interactive interface, experiment with different values of the parameters.
6
.1
.In order to form “qualitative understanding.”
7
.Get real-life data.
7
.1
.Say, for the 2019–2020 coronavirus outbreak.
8
. Attempt manual or automatic calibration of the model.
8
.1
.This step will most likely require additional data transformations and programming.
8
.2
.Only manual calibration is shown in this notebook.
Load Packages of the Framework
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].
In[]:=
Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModels.m"]Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelModifications.m"]Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/WL/SystemDynamicsInteractiveInterfacesFunctions.m"]
Getting the Model Code
Getting the Model Code
In[]:=
modelSI2R=SEI2RModel[t,"InitialConditions"True,"RateRules"True];
We can show a tabulated visualization of the model using the function ModelGridTableForm from [AAp1]:
In[]:=
ModelGridTableForm[modelSI2R]
Out[]=
Stocks
,Rates
,Equations
,RateRules
,InitialConditions
# | 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 |
# | 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) |
# | Equation |
1 | ′ SP SP[t] μ[TP] INSP[t] SP[t] β[INSP] TP[0] ISSP[t] SP[t] β[ISSP] TP[0] |
2 | ′ EP EP[t] 1 aincp μ[TP] INSP[t] SP[t] β[INSP] TP[0] ISSP[t] SP[t] β[ISSP] TP[0] |
3 | ′ INSP INSP[t] aip EP[t] sspf[SP] aincp INSP[t] μ[INSP] |
4 | ′ ISSP ISSP[t] aip EP[t] sspf[SP] aincp ISSP[t] μ[ISSP] |
5 | ′ RP INSP[t] ISSP[t] aip RP[t] μ[TP] |
6 | ′ MLP lpcr[ISSP,INSP] RP[t] SP[t] |
# | 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 |
# | 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 |
Model Extensions and New Models
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:
1
.By adding a birthrate to the Susceptible Population equation (the birthrate is not included by default).
2
.By adding a new equation for the infected deceased population.
Adding Births Term
Adding Births Term
Here are the equations of SEI2R (from [AAp1]):
In[]:=
ModelGridTableForm[modelSI2R]["Equations"]
Out[]=
# | Equation |
1 | ′ SP SP[t] μ[TP] INSP[t] SP[t] β[INSP] TP[0] ISSP[t] SP[t] β[ISSP] TP[0] |
2 | ′ EP EP[t] 1 aincp μ[TP] INSP[t] SP[t] β[INSP] TP[0] ISSP[t] SP[t] β[ISSP] TP[0] |
3 | ′ INSP INSP[t] aip EP[t] sspf[SP] aincp INSP[t] μ[INSP] |
4 | ′ ISSP ISSP[t] aip EP[t] sspf[SP] aincp ISSP[t] μ[ISSP] |
5 | ′ RP INSP[t] ISSP[t] aip RP[t] μ[TP] |
6 | ′ MLP lpcr[ISSP,INSP] RP[t] SP[t] |
Here we find the position of the equation that corresponds to “Susceptible Population”:
In[]:=
pos=EquationPosition[modelSI2R["Equations"],First@GetPopulationSymbols[modelSI2R,"Susceptible Population"]]
Out[]=
1
Here we make the births term using a birthrate that is the same as the death rate:
In[]:=
birthTerm=GetPopulations[modelSI2R,"Total Population"]〚1〛*GetRates[modelSI2R,"Population death rate"]〚1〛
Out[]=
TP[t]μ[TP]
Here we add the births term to the equations of the new model:
In[]:=
modelSI2RNew=modelSI2R;modelSI2RNew["Equations"]=ReplacePart[modelSI2R["Equations"],{1,2}modelSI2R["Equations"]〚1,2〛+birthTerm];
Here we display the equations of the new model:
Adding Infected Deceased Population Equation
Adding Infected Deceased Population Equation
Here we add new population, equation and initial condition, which allow for tracking the deaths because of infection:
Here is how the model looks:
Parameters and Actual Simulation Equations Code
Parameters and Actual Simulation Equations Code
Here are the parameters we want to experiment with (or do calibration with):
Here we set custom rates and initial conditions:
Here is the system of ODEs we use to do parametrized simulations:
Simulation
Simulation
(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
Interactive Interface
Calibration over Real Data
Calibration over Real Data
It is important to calibrate these kinds of models with real data, or at least to give a serious attempt to do such a calibration. If the calibration is “too hard” or “impossible,” that would indicate that the model is not that adequate (if adequate at all).
The calibration efforts can be (semi-)automated using special model-to-data goodness-of-fit functions and a minimization procedure. (See, for example [AA2].)
In this section we just attempt to calibrate SEI2R over real data taken from [WRI1] using a specialized interactive interface.
Real Data
Real Data
Here is COVID-19 data taken from [WRI1] for the Chinese province Hubei:
The total population in Hubei is:
But we have to use a fraction of that population in order to produce good fits. That can be justified with the conjecture that the citizens of Hubei are spread out and it is mostly one city (Wuhan) where the outbreak is.
The real data has to be padded with a certain number of 0s related to the infectious and incubation periods in order to make good fits. Such padding is easy to justify: if we observe recovered people, that means that they passed through the incubation and infectious periods.
Calibration Interactive Interface
Calibration Interactive Interface
In this interface, we put the Infected Severely Symptomatic Population (ISSP) to zero. That way it is easier to compare the real data with the simulated results (and pick parameter values that give close fits). Also note that since SEI2R is simple in this interface, the system is always solved:
Maybe Good Enough Parameters
Maybe Good Enough Parameters
Basic reproduction number:
Basic reproduction number:
References
References
Articles
Articles
[Wk1] Wikipedia, "Compartmental Models in Epidemiology."
[HH1] Hethcote, Herbert W. (2000). "The Mathematics of Infectious Diseases." SIAM Review. 42(4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.
[AA1] Antonov, Anton. (2020). "Coronavirus Propagation Modeling Considerations." SystemModeling at GitHub.
[AA2] Antonov, Anton (2019). Answer of "Model Calibration with Phase Space Data." Mathematica StackExchange.
Repositories & Packages
Repositories & Packages
[AAr1] Antonov, Anton. (2020). “Coronavirus Propagation Dynamics Project.” SystemModeling at GitHub.
[AAp1] AAntonov, Anton. (2020). "Epidemiology Models Mathematica Package." SystemModeling at GitHub.
[AAp2] Antonov, Anton. (2020). "Epidemiology Models Modifications Mathematica Package." SystemModeling at GitHub.
[AAp3] Antonov, Anton. (2020). "System Dynamics Interactive Interfaces Functions Mathematica Package." SystemModeling at GitHub.