Message Boards

GROUPS:

Staff Picks Biological Sciences Medical Sciences Operations Research Social Science Dynamic Interactivity Equation Solving Interface Construction Modeling

Anton Antonov

[NB] Basic experiments workflow for simple epidemiological models

Anton Antonov, Accendo Data LLC

Posted 2 years ago

13682 Views

10 Replies

12 Total Likes

Follow this post

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus

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:

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 died 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–20 coronavirus outbreak, [WRI1].)

Workflow

Get one of the classical epidemiology models.

Extend the equations of model if needed or desired.

Set relevant initial conditions for the populations.

Pick model parameters to be adjust and “play with.”

Derive parametrized solutions of model’s system of equations (ODE’s or DAE’s.)

Using the parameters of the previous step.

Using an interactive interface experiment with different values of the parameters.

In order to form “qualitative understanding.”

Get real life data.

Say, for the 2019-20 coronavirus outbreak.

Attempt manual or automatic calibration of the model.

This step will most likely require additional data transformations and programming.

Only manual calibration is shown in this notebook.

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-dynamics/WL/EpidemiologyModels.m"]Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelModifications.m"]Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/WL/SystemDynamicsInteractiveInterfacesFunctions.m"]

Getting the model code

Here we take the SEI2R model implemented in the package "EpidemiologyModels.m", [AAp1]:

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

#	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]	6
6	β[INSP]	3
7	aip	28
8	aincp	6
9	sspf[SP]	0.2
10	lpcr[ISSP,INSP]	600

,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



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 birth rate added to the Susceptible Population equation (the birth rate is not included by default)

By adding a new equation for the infected deceased population.

Adding births term

Here are the equations of SEI2R (from [AAp1]):

In[]:=

ModelGridTableForm[modelSI2R]["Equations"]

Out[]=

#	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])

Here we find the position of the equation that corresponds to “Susceptible Population”:

In[]:=

pos=EquationPosition[modelSI2R["Equations"],First@GetPopulationSymbols[modelSI2R,"Susceptible Population"]]

Out[]=

Here we make the births term using a birth rate 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 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:

In[]:=

ModelGridTableForm[modelSI2RNew]["Equations"]

Out[]=

#	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])

Adding infected deceased population equation

Here we add new population, equation, and initial condition that allow for tracking the deaths because of infection:

In[]:=

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];

Here is how the model looks like:

In[]:=

ModelGridTableForm[modelSI2R]["Equations"]

Out[]=

#	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]

Parameters and actual simulation equations code

Here are the parameters we want to experiment with (or do calibration with):

In[]:=

lsFocusParams={aincp,aip,sspf[SP],β[ISSP],β[INSP]};

Here we set custom rates and initial conditions:

In[]:=

population=58160000/400;modelSI2R=SetRateRules[modelSI2R,<|TP[0]population|>];modelSI2R=SetInitialConditions[modelSI2R,<|SP[0]population-1,ISSP[0]0,INSP[0]1|>];

Here is the system of ODE’s we use with to do parametrized simulations:

In[]:=

lsActualEquations=Join[modelSI2R["Equations"]//.KeyDrop[modelSI2R["RateRules"],lsFocusParams],modelSI2R["InitialConditions"]];ResourceFunction["GridTableForm"][List/@lsActualEquations]

Out[]=

#	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]600(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

Simulation

Straightforward simulation for one year with using ParametricNDSolve :

In[]:=

aSol=Association@Flatten@ParametricNDSolve[lsActualEquations,Head/@Keys[modelSI2R["Stocks"]],{t,0,365},lsFocusParams]

Out[]=

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]}





(The advantage having parametrized solutions is that we can quickly compute simulation results with new parameter values without solving model’s system of ODE’s; see the interfaces below.)

Interactive interface

In[]:=

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 symptomatic population fraction"},0,1,0.025,Appearance{"Open"}},{{crisp,0.56,"Contact rate of the infected severely symptomatic population"},0,20,0.01,Appearance{"Open"}},{{criap,0.56,"Contact rate of the infected normally symptomatic population"},0,20,0.01,Appearance{"Open"}},{{ndays,160,"Number of days"},1,365,1,Appearance{"Open"}},{{popTogetherQ,True,"Plot populations together"},{False,True}},{{popDerivativesQ,False,"Plot populations derivatives"},{False,True}},{{popLogPlotQ,False,"LogPlot populations"},{False,True}},{{econTogetherQ,False,"Plot economics functions together"},{False,True}},{{econDerivativesQ,False,"Plot economics functions derivatives"},{False,True}},{{econLogPlotQ,False,"LogPlot economics functions"},{False,True}},{{nPlotColumns,1,"Number of plot columns"},Range[5]},ControlPlacementLeft,ContinuousActionFalse]

Out[]=

Average incubation period (days)

Average infectious period (days)