Investigating COVID-19 with R: data analysis and simulations

Posted on May 29, 2020 by Anton Antonov Antonov

Methodological presentation
R-Ladies Miami Meetup, May 28th 2020

The extended abstract of the presentation was loosely followed. Here is the presentation mind-map:

(Note that mind-map’s PDF has hyperlinks. Also, see the folder Presentation-aids. )

The organizers and I did a poll for what people want to hear. After discussing the results of the 15 votes from that poll we decided the presentation to be a methodological one instead of a know-how one.

Approximately 30% of the presentation was based on the R-project “COVID-19-modeling-in-R”, [AA1].

Approximately 30% of the presentation was based on an R-programmed software monad for epidemiology compartmental models, ECMMon-R, [AAr2].

For the rest were used frameworks, simulations, and graphics made with Mathematica, [AAr1].

The presentation was given online (because of COVID-19) using Zoom. 90 people registered. Nearly 40 showed up (and maybe 20 stayed throughout.)

Here is a link to the video recording.

Screenshots

Here are screenshots of statistics used in the introduction:

References

Coronavirus

[Wk1] Wikipedia entry, Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2).

[Wk2] Wikipedia entry, Coronavirus disease 2019.

Modeling

[Wk3] Wikipedia entry, Compartmental models in epidemiology.

[Wk4] Wikipedia entry, System dynamics.

R code/software

[KS1] Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer, “deSolve: Solvers for Initial Value Problems of Differential Equations (‘ODE’, ‘DAE’, ‘DDE’)”, CRAN.

[AA1] Anton Antonov, “COVID-19-modeling-in-R”, 2020, SystemModeling at GitHub.

[AAr1] Anton Antonov, Coronavirus-propagation-dynamics, 2020, SystemModeling at GitHub.

[AAr2] Anton Antonov, Epidemiology Compartmental Modeling Monad in R, 2020, ECMMon-R at GitHub.

Coronavirus propagation modeling, useR! Boston April 2020

Posted on April 15, 2020 by Anton Antonov Antonov

Tutorial presentation

The extended abstract of the presentation was loosely followed. Here is the (main) presentation mind-map:

(Note that mind-map’s PDF has hyperlinks. Also, see the folder Presentation-aids. )

Approximately 70% of the presentation was based on an R-programmed software monad for epidemiology compartmental models, ECMMon-R, [AAr2]. For the rest were used frameworks, simulations, and graphics made with Mathematica, [AAr1], and Wolfram System Modeler .

The presentation was given online (because of COVID-19) using Zoom. 190 people registered. Nearly 70 showed up (and maybe 60 stayed throughout.)

Here is a link to the video recording.

A similar presentation was given two weeks prior for OMLDS.

References

Coronavirus

[Wk1] Wikipedia entry, Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2).

[Wk2] Wikipedia entry, Coronavirus disease 2019.

Modeling

[Wk3] Wikipedia entry, Compartmental models in epidemiology.

[Wk4] Wikipedia entry, System dynamics.

[JD1] Jim Duggan, System Dynamics Modeling with R, 2016, Springer.

R code/software

[KS1] Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer, “deSolve: Solvers for Initial Value Problems of Differential Equations (‘ODE’, ‘DAE’, ‘DDE’)”, CRAN.

[JD2] Jim Duggan, “SDMR”, 2016, GitHub.
(Resources for text book “System Dynamics Modeling with R”.)

[AA1] Anton Antonov, “COVID-19-modeling-in-R”, 2020, SystemModeling at GitHub.

[AAr1] Anton Antonov, Coronavirus-propagation-dynamics, 2020, SystemModeling at GitHub.

[AAr2] Anton Antonov, Epidemiology Compartmental Modeling Monad in R, 2020, ECMMon-R at GitHub.

SEI2HR-Econ model with quarantine and supplies scenarios

Posted on April 11, 2020 by Anton Antonov Antonov

Introduction

The epidemiology compartmental model, [Wk1], presented in this notebook – SEI2HR-Econ – deals with all three rectangles in this diagram:

ImageResize[Import["https://github.com/antononcube/SystemModeling/raw/master/Projects/Coronavirus-propagation-dynamics/Diagrams/Coronavirus-propagation-simple-dynamics.jpeg"], 900]

“SEI2HR” stands for “Susceptible, Exposed, Infected two, Hospitalized, Recovered” (populations.) “Econ” stands for “Economic”.

In this notebook we also deal with both quarantine scenarios and medical supplies scenarios. In the notebook [AA4] we deal with quarantine scenarios over a simpler model, SEI2HR.

Remark: We consider the contagious disease propagation models as instances of the more general System Dynamics (SD) models. We use SD terminology in this notebook.

The models

SEI2R

The model SEI2R is introduced and explained in the notebook [AA2]. SEI2R differs from the classical SEIR model, [Wk1, HH1], with the following elements:

Two separate infected populations: one is “severely symptomatic”, the other is “normally symptomatic”
The monetary equivalent of lost productivity due to infected or died people is tracked

SEI2HR

For the formulation of SEI2HR we use a system of Differential Algebraic Equations (DAE’s). The package [AAp1] allows the use of a formulation that has just Ordinary Differential Equations (ODE’s).

Here are the unique features of SEI2HR:

People stocks
- There are two types of infected populations: normally symptomatic and severely symptomatic.
- There is a hospitalized population.
- There is a deceased from infection population.
Hospital beds
- Hospital beds are a limited resource that determines the number of hospitalized people.
- Only severely symptomatic people are hospitalized according to the available hospital beds.
- The hospital beds stock is not assumed constant, it has its own change rate.
Money stocks
- The money from lost productivity is tracked.
- The money for hospital services is tracked.

SEI2HR-Econ

SEI2HR-Econ adds the following features to SEI2HR:

Medical supplies
- Medical supplies production is part of the model.
- Medical supplies delivery is part of the model..
- Medical supplies accumulation at hospitals is taken into account.
- Medical supplies demand tracking.
Hospitalization
- Severely symptomatic people are hospitalized according to two limited resources: hospital beds and medical supplies.
Money stocks
- Money for medical supplies production is tracked.

SEI2HR-Econ’s place a development plan

This graph shows the “big picture” of the model development plan undertaken in [AAr1] and SEI2HR (discussed in this notebook) is in that graph:

Notebook structure

The rest of notebook has the following sequence of sections:

Package load section
SEI2HR-Econ structure in comparison of SEI2HR
Explanations of the equations of SEI2HR-Econ
Quarantine scenario modeling preparation
Medical supplies production and delivery scenario modeling preparation
Parameters and initial conditions setup
- Populations, hospital beds, quarantine scenarios, medical supplies scenarios
Simulation solutions
Interactive interface
Sensitivity analysis

Load packages

The epidemiological models framework used in this notebook is implemented with the packages [AAp1-AAp4, AA3]; many of the plot functions are from the package [AAp5].

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/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelingVisualizationFunctions.m"];
Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelingSimulationFunctions.m"];
Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/WL/SystemDynamicsInteractiveInterfacesFunctions.m"];

SEI2HR-Econ extends SEI2HR

The model SEI2HR-Econ is an extension of the model SEI2HR, [AA4].

Here is SEI2HR:

reprTP = "AlgebraicEquation";
lsModelOpts = {"Tooltips" -> True, 
   TooltipStyle -> {Background -> Yellow, CellFrameColor -> Gray, 
     FontSize -> 20}};
modelReference = 
  SEI2HRModel[t, "InitialConditions" -> True, "RateRules" -> True, 
   "TotalPopulationRepresentation" -> reprTP];
ModelGridTableForm[modelReference, lsModelOpts]

Here is SEI2HR-Econ:

modelSEI2HREcon = 
  SEI2HREconModel[t, "InitialConditions" -> True, "RateRules" -> True,
    "TotalPopulationRepresentation" -> reprTP];
ModelGridTableForm[modelSEI2HREcon, lsModelOpts]

Here are the “differences” between the two models:

ModelGridTableForm@
 Merge[{modelSEI2HREcon, modelReference}, 
  If[AssociationQ[#[[1]]], KeyComplement[#], Complement @@ #] &]

Equations explanations

In this section we provide rationale for the equations of SEI2HR-Econ.

The equations for Susceptible, Exposed, Infected, Recovered populations of SEI2R are “standard” and explanations about them are found in [WK1, HH1]. For SEI2HR those equations change because of the stocks Hospitalized Population and Hospital Beds. For SEI2HR-Econ the SEI2HR equations change because of the stocks Medical Supplies, Medical Supplies Demand, and Hospital Medical Supplies.

The equations time unit is one day. The time horizon is one year. Since we target COVID-19, [Wk2, AA1], we do not consider births.

Remark: For convenient reading the equations in this section have tooltips for the involved stocks and rates.

Verbalization description of the model

We start with one infected (normally symptomatic) person, the rest of the people are susceptible. The infected people meet other people directly or get in contact with them indirectly. (Say, susceptible people touch things touched by infected.) For each susceptible person there is a probability to get the decease. The decease has an incubation period: before becoming infected the susceptible are (merely) exposed. The infected recover after a certain average infection period or die. A certain fraction of the infected become severely symptomatic. The severely symptomatic infected are hospitalized if there are enough hospital beds and enough medical supplies. The hospitalized severely infected have different death rate than the non-hospitalized ones. The number of hospital beds might change: hospitals are extended, new hospitals are build, or there are not enough medical personnel or supplies.

The different types of populations (infected, hospitalized, recovered, etc.) have their own consumption rates of medical supplies. The medical supplies are produced with a certain rate (units per day) and delivered after a certain delay period. The hospitals have their own storage for medical supplies. Medical supplies are delivered to the hospitals only, non-hospitalized people go to the medical supplies producer to buy supplies. The hospitals have precedence for the medical supplies: if the medical supplies are not enough for everyone, the hospital needs are covered first (as much as possible.)

The medical supplies producer has a certain storage capacity (for supplies.) The medical supplies delivery vehicles have a certain – generally speaking, smaller – capacity. The hospitals have a certain capacity to store medical supplies. It is assumed that both producer and hospitals have initial stocks of medical supplies. (Following a certain normal, general preparedness protocol.)

The combined demand from all populations for medical supplies is tracked (accumulated.) The deaths from infection are tracked (accumulated.) Money for medical supplies production, money for hospital services, and money from lost productivity are tracked (accumulated.)

The equations below give mathematical interpretation of the model description above.

Code for the equations

Each equation in this section are derived with code like this:

ModelGridTableForm[modelSEI2HREcon, lsModelOpts]["Equations"][[1, 
 EquationPosition[modelSEI2HREcon, RP] + 1, 2]]

and then the output cell is edited to be “DisplayFormula” and have CellLabel value corresponding to the stock of interest.

The infected and hospitalized populations

SEI2HR has two types of infected populations: a normally symptomatic one and a severely symptomatic one. A major assumption for SEI2HR is that only the severely symptomatic people are hospitalized. (That assumption is also reflected in the diagram in the introduction.)

Each of those three populations have their own contact rates and mortality rates.

Here are the contact rates from the SEI2HR-Econ dictionary

ColumnForm@
 Cases[Normal@modelSEI2HREcon["Rates"], 
  HoldPattern[\[Beta][_] -> _], \[Infinity]]

Here are the mortality rates from the SEI2HR-Econ dictionary

ColumnForm@
 Cases[Normal@modelSEI2HREcon["Rates"], 
  HoldPattern[\[Mu][_] -> _], \[Infinity]]

Remark: Below with “Infected Population” we mean both stocks Infected Normally Symptomatic Population (INSP) and Infected Severely Symptomatic Population (ISSP).

Total Population

In this notebook we consider a DAE’s formulation of SEI2HR-Econ. The stock Total Population has the following (obvious) algebraic equation:

Note that with Max we specified that the total population cannot be less than $0$ .

Remark: As mentioned in the introduction, the package [AAp1] allows for the use of non-algebraic formulation, without an equation for TP.

Susceptible Population

The stock Susceptible Population (SP) is decreased by (1) infections derived from stocks Infected Populations and Hospitalized Population (HP), and (2) morality cases derived with the typical mortality rate.

Because we hospitalize the severely infected people only instead of the term

we have the terms

The first term is for the infections derived from the hospitalized population. The second term for the infections derived from people who are infected severely symptomatic and not hospitalized.

Births term

Note that we do not consider in this notebook births, but the births term can be included in SP’s equation:

Block[{m = SEI2HREconModel[t, "BirthsTerm" -> True]},
 ModelGridTableForm[m]["Equations"][[1, EquationPosition[m, SP] + 1, 
  2]]
 ]

The births rate is the same as the death rate, but it can be programmatically changed. (See [AAp2].)

Exposed Population

The stock Exposed Population (EP) is increased by (1) infections derived from the stocks Infected Populations and Hospitalized Population, and (2) mortality cases derived with the typical mortality rate. EP is decreased by (1) the people who after a certain average incubation period (aincp) become ill, and (2) mortality cases derived with the typical mortality rate.

Infected Normally Symptomatic Population

INSP is increased by a fraction of the people who have been exposed. That fraction is derived with the parameter severely symptomatic population fraction (sspf). INSP is decreased by (1) the people who recover after a certain average infection period (aip), and (2) the normally symptomatic people who die from the disease.

Infected Severely Symptomatic Population

ISSP is increased by a fraction of the people who have been exposed. That fraction is corresponds to the parameter severely symptomatic population fraction (sspf). ISSP is decreased by (1) the people who recover after a certain average infection period (aip), (2) the hospitalized severely symptomatic people who die from the disease, and (3) the non-hospitalized severely symptomatic people who die from the disease.

Note that we do not assume that severely symptomatic people recover faster if they are hospitalized, only that they have a different death rate.

Hospitalized Population

The amount of people that can be hospitalized is determined by the available Hospital Beds (HB) – the stock Hospitalized Population (HP) is subject to a resource limitation by the stock HB.

The equation of the stock HP can be easily understood from the following dynamics description points:

If the number of hospitalized people is less that the number of hospital beds we hospitalize the new ISSP people.
The Available Hospital Beds (AHB) are determined by the minimum of (i) the non-occupied hospital beds, and (ii) the hospital medical supplies divided by the ISSP consumption rate.
If the new ISSP people are more than AHB the hospital takes as many as AHB.
Hospitalized people have the same average infection period (aip).
Hospitalized (severely symptomatic) people have their own mortality rate.

Here is the HP equation:

Note that although we know that in a given day some hospital beds are going to be freed they are not considered in the hospitalization plans for that day. Similarly, we know that new medical supplies are coming but we do not include them into AHB.

Recovered Population

The stock Recovered Population (RP) is increased by the recovered infected people and decreased by mortality cases derived with the typical mortality rate.

Deceased Infected Population

The stock Deceased Infected Population (DIP) accumulates the deaths of the people who are infected. Note that we utilize the different death rates for HP and ISSP.

Hospital Beds

The stock Hospital Beds (HB) can change with a rate that reflects the number of hospital beds change rate (nhbcr) per day. Generally speaking, using nhbcr we can capture scenarios, like, extending hospitals, building new hospitals, recruitment of new medical personnel, loss of medical personnel (due to infections.)

Hospital Medical Supplies

The Hospital Medical Supplies (HMS) are decreased according to the medical supplies consumption rate (mscr) of HP and increased by a Medical Supplies (MS) delivery term (to be described next.)

The MS delivery term is build with the following assumptions / postulates:

Every day the hospital attempts to order MS that correspond to HB multiplied by mscr.
The hospital has limited capacity of MS storage, $\kappa [\text{HMS}]$ .
The MS producer has limited capacity for delivery, $\kappa [\text{MDS}]$ .
The hospital demand for MS has precedence over the demands for the non-hospitalized populations.
Hence, if the MS producer has less stock of MS than the demand of the hospital then MS producer’s whole amount of MS goes to the hospital.
The supplies are delivered with some delay period: the medical supplies delivery period (msdp).

Here is the MS delivery term:

Here is the corresponding HMS equation:

Medical Supplies

The equation of the Medical Supplies (MS) stock is based on the following assumptions / postulates:

The non-hospitalized people go to the MS producer to buy supplies. (I.e. MS delivery is to the hospital only.)
The MS producer vehicles have certain capacity, $\kappa [\text{MSD}]$ .
The MS producer has a certain storage capacity (for MS stock.)
Each of the populations INSP, ISSP, and HP has its own specific medical supplies consumption rate (mscr). EP, RP, and TP have the same mscr.
The hospital has precedence in its MS order. I.e. the demand from the hospital is satisfied first, and then the demand of the rest of the populations.

Here is the MS delivery term described in the previous section:

Here is the MS formula with the MS delivery term replaced with “Dlvr”:

ModelGridTableForm[modelSEI2HREcon, "Tooltips" -> False][
   "Equations"][[1, EquationPosition[modelSEI2HREcon, MS] + 1, 2]] /. 
 dlvr -> Dlvr

We can see from that equation that MS is increased by medical supplies production rate (mspr) with measuring dimension number of units per day. The production is restricted by the storage capacity, $\kappa [\text{MS}]$ :

(*Min[mspr[HB], -MS[t] + \[Kappa][MS]]*)

MS is decreased by the MS delivery term and the demand from the non-hospitalized populations. Because the hospital has precedence, we use this term form in the equation:

(*Min[-Dlvr + MS[t], "non-hospital demand"]*)

Here is the full MS equation:

Medical Supplies Demand

The stock Medical Supplies Demand (MSD) simply accumulates the MS demand derived from population stocks and their corresponding mscr:

Money for Hospital Services

The stock Money for Hospital Services (MHS) simply tracks expenses for hospitalized people. The parameter hospital services cost rate (hscr) with unit money per bed per day simply multiplies HP.

Money from Lost Productivity

The stock Money from Lost Productivity (MLP) simply tracks the work non-availability of the infected and died from infection people. The parameter lost productivity cost rate (lpcr) with unit money per person per day multiplies the total count of the infected and dead from infection.

Quarantine scenarios

In order to model quarantine scenarios we use piecewise constant functions for the contact rates $\beta [\text{ISSP}]$ and $\beta [\text{INSP}]$ .

Remark: Other functions can be used, like, functions derived through some statistical fitting.

Here is an example plot :

Block[{func = \[Beta]*
    Piecewise[{{1, t < qsd}, {qcrf, qsd <= t <= qsd + ql}}, 1]},
 Legended[
  Block[{\[Beta] = 0.56, qsd = 60, ql = 8*7, qcrf = 0.25},
   ListLinePlot[Table[func, {t, 0, 365}], PlotStyle -> "Detailed"]
   ], func]]

To model quarantine with a piecewise constant function we use the following parameters:

$\text{qsd}$ for quarantine’s start
$\text{ql}$ for quarantines duration
$\text{qcrf}$ for the effect on the quarantine on the contact rate

Medical supplies scenarios

We consider three main scenarios for the medical supplies:

Constant production rate and consistent delivery (constant delivery period)
Constant production rate and disrupted delivery
Disrupted production and disrupted delivery

The disruptions are simulated with simple pulse functions – we want to see how the system being modeled reacts to singular, elementary disruption.

Here is an example plot of a disruption of delivery period plot :

Block[{func = 
   dbp*Piecewise[{{1, t < dds}, {dsf, dds <= t <= dds + ddl}}, 1]},
 Legended[
  Block[{dbp = 1, dds = 70, ddl = 7, dsf = 1.8},
   ListLinePlot[Table[func, {t, 0, 365}], PlotStyle -> "Detailed"]
   ], func]]

To model disruption of delivery with a piecewise constant function we use the following parameters:

$\text{dbp}$ for the delivery base period
$\text{dds}$ for delivery disruption start
$\text{ddl}$ for delivery disruption duration
$\text{dsf}$ for how much slower or faster the delivery period becomes

Parameters and actual simulation equations code

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

lsFocusParams = {aincp, aip, sspf[SP], \[Beta][HP], qsd, ql, qcrf, 
   nhbcr[ISSP, INSP], nhbr[TP], mspr[HB]};

Here we set custom rates and initial conditions:

population = 10^6;
modelSEI2HREcon =
  SetRateRules[
   modelSEI2HREcon,
   <|
    TP[0] -> population,
    \[Beta][ISSP] -> 
     0.5*Piecewise[{{1, t < qsd}, {qcrf, qsd <= t <= qsd + ql}}, 1],
    \[Beta][INSP] -> 
     0.5*Piecewise[{{1, t < qsd}, {qcrf, qsd <= t <= qsd + ql}}, 1],
    qsd -> 60,
    ql -> 8*7,
    qcrf -> 0.25,
    \[Beta][HP] -> 0.01,
    \[Mu][ISSP] -> 0.035/aip,
    \[Mu][INSP] -> 0.01/aip,
    nhbr[TP] -> 3/1000,
    lpcr[ISSP, INSP] -> 1,
    hscr[ISSP, INSP] -> 1,
    msdp[HB] -> 
     dbp*Piecewise[{{1, t < dds}, {dsf, dds <= t <= dds + ddl}}, 1],
    dbp -> 1,
    dds -> 70,
    ddl -> 7,
    dsf -> 2
    |>
   ];

Remark: Note the piecewise functions for $\beta [\text{ISSP}]$ , $\beta [\text{INSP}]$ , and $\text{msdp}[\text{HB}]$ .

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

lsActualEquations = 
  ModelNDSolveEquations[modelSEI2HREcon, 
   KeyDrop[modelSEI2HREcon["RateRules"], lsFocusParams]];
ResourceFunction["GridTableForm"][List /@ lsActualEquations]

lsActualEquations = 
  ModelNDSolveEquations[modelSEI2HREcon, modelSEI2HREcon["RateRules"]];
ResourceFunction["GridTableForm"][List /@ lsActualEquations]

Simulation

Instead of using ParametricNDSolve as in [AA4] in this notebook we use ModelNDSolve and SetRateRules from the package [AAp4].

Different quarantine starts

Here we compute simulation solutions for a set of quarantine starts:

AbsoluteTiming[
 aVarSolutions =
   Association@
    Map[
     Function[{qsdVar},
      qsdVar -> 
       Association[
        ModelNDSolve[
          SetRateRules[
           modelSEI2HREcon, <|ql -> 56, qsd -> qsdVar|>], {t, 365}][[
         1]]]
      ],
     Range[40, 120, 20]
     ];
 ]

(*{0.366168, Null}*)

Here we plot the results for ISSP only:

SeedRandom[2532]
aVals = #[ISSP][Range[0, 365]] & /@ aVarSolutions;
ListLinePlot[
 KeyValueMap[
  Callout[Tooltip[#2, #1], #1, {If[#1 <= 70, 
      RandomInteger[{120, 200}], RandomInteger[{80, 110}]], Above}] &,
   aVals], PlotLegends -> 
  SwatchLegend[Keys[aVals], LegendLabel -> "Quarantine start"], 
 PlotRange -> All, ImageSize -> Large, 
 PlotLabel -> ISSP[t] /. modelSEI2HREcon["Stocks"]]

Remark: We use the code in this section to do the computations in the section “Sensitivity Analysis”.

Interactive interface

Using the interface in this section we can interactively see the effects of changing parameters. (This interface is programmed without using parametrized NDSolve solutions in order to be have code that corresponds to the interface implementations in [AAr2].)

opts = {PlotRange -> All, PlotLegends -> None, 
   PlotTheme -> "Detailed", PerformanceGoal -> "Speed", 
   ImageSize -> 400};
lsPopulationKeys = {TP, SP, EP, ISSP, INSP, HP, RP, DIP, HB};
lsSuppliesKeys = {MS, MSD, HMS};
lsMoneyKeys = {MHS, MLP, MMSP};
Manipulate[
 DynamicModule[{modelLocal = modelSEI2HREcon, 
   aStocks = modelSEI2HREcon["Stocks"], aSolLocal = aParSol, 
   lsPopulationPlots, lsMoneyPlots, lsSuppliesPlots},
  
  modelLocal = 
   SetRateRules[
    modelLocal, <|aincp -> aincpM, aip -> aipM, 
     sspf[SP] -> sspfM, \[Beta][HP] -> crhpM, qsd -> qsdM, ql -> qlM, 
     qcrf -> qcrfM, nhbr[TP] -> nhbrM/1000, 
     nhbcr[ISSP, ISNP] -> nhbcrM, mspr[HB] -> msprM, 
     msdp[HB] -> msdpM|>];
  aSolLocal = Association[ModelNDSolve[modelLocal, {t, ndays}][[1]]];
  
  lsPopulationPlots =
   Quiet@ParametricSolutionsPlots[
     aStocks,
     KeyTake[aSolLocal, 
      Intersection[lsPopulationKeys, displayPopulationStocks]],
     None, ndays,
     "LogPlot" -> popLogPlotQ, "Together" -> popTogetherQ, 
     "Derivatives" -> popDerivativesQ, 
     "DerivativePrefix" -> "\[CapitalDelta]", opts, 
     Epilog -> {Gray, Dashed, 
       Line[{{qsdM, 0}, {qsdM, 1.5*population}}], 
       Line[{{qsdM + qlM, 0}, {qsdM + qlM, 1.5*population}}]}];
  
  lsSuppliesPlots =
   If[Length[
      KeyDrop[aSolLocal, Join[lsPopulationKeys, lsMoneyKeys]]] == 
     0, {},
    (*ELSE*)
    Quiet@ParametricSolutionsPlots[
      aStocks,
      KeyTake[KeyDrop[aSolLocal, Join[lsPopulationKeys, lsMoneyKeys]],
        displaySupplyStocks],
      None, ndays,
      "LogPlot" -> supplLogPlotQ, "Together" -> supplTogetherQ, 
      "Derivatives" -> supplDerivativesQ, 
      "DerivativePrefix" -> "\[CapitalDelta]", opts]
    ];
  
  lsMoneyPlots =
   Quiet@ParametricSolutionsPlots[
     aStocks,
     KeyTake[aSolLocal, Intersection[lsMoneyKeys, displayMoneyStocks]],
     None, ndays,
     "LogPlot" -> moneyLogPlotQ, "Together" -> moneyTogetherQ, 
     "Derivatives" -> moneyDerivativesQ, 
     "DerivativePrefix" -> "\[CapitalDelta]", opts];
  
  Multicolumn[Join[lsPopulationPlots, lsSuppliesPlots, lsMoneyPlots], 
   nPlotColumns, Dividers -> All, FrameStyle -> GrayLevel[0.8]],
  SaveDefinitions -> True
  ],
 {{ndays, 365, "Number of days"}, 1, 365, 1, Appearance -> {"Open"}},
 Delimiter,
 {{aincpM, 6., "Average incubation period (days)"}, 1, 60., 1, 
  Appearance -> {"Open"}},
 {{aipM, 21., "Average infectious period (days)"}, 1, 60., 1, 
  Appearance -> {"Open"}},
 {{sspfM, 0.2, "Severely symptomatic population fraction"}, 0, 1, 
  0.025, Appearance -> {"Open"}},
 {{crhpM, 0.1, "Contact rate of the hospitalized population"}, 0, 30, 
  0.1, Appearance -> {"Open"}},
 Delimiter,
 {{qsdM, 55, "Quarantine start days"}, 0, 365, 1, 
  Appearance -> {"Open"}},
 {{qlM, 8*7, "Quarantine length (in days)"}, 0, 120, 1, 
  Appearance -> {"Open"}},
 {{qcrfM, 0.25, "Quarantine contact rate fraction"}, 0, 1, 0.01, 
  Appearance -> {"Open"}},
 Delimiter,
 {{nhbrM, 2.9, "Number of hospital beds rate (per 1000 people)"}, 0, 
  100, 0.1, Appearance -> {"Open"}},
 {{nhbcrM, 0, "Number of hospital beds change rate"}, -0.5, 0.5, 
  0.001, Appearance -> {"Open"}},
 {{msprM, 200, "Medical supplies production rate"}, 0, 50000, 10, 
  Appearance -> {"Open"}},
 {{msdpM, 1.2, "Medical supplies delivery period"}, 0, 10, 0.1, 
  Appearance -> {"Open"}},
 Delimiter,
 {{displayPopulationStocks, lsPopulationKeys, 
   "Population stocks to display:"}, lsPopulationKeys, 
  ControlType -> TogglerBar},
 {{popTogetherQ, True, "Plot populations together"}, {False, True}},
 {{popDerivativesQ, False, "Plot populations derivatives"}, {False, 
   True}},
 {{popLogPlotQ, False, "LogPlot populations"}, {False, True}},
 Delimiter,
 {{displaySupplyStocks, lsSuppliesKeys, 
   "Supplies stocks to display:"}, lsSuppliesKeys, 
  ControlType -> TogglerBar},
 {{supplTogetherQ, True, "Plot supplies functions together"}, {False, 
   True}},
 {{supplDerivativesQ, False, 
   "Plot supplies functions derivatives"}, {False, True}},
 {{supplLogPlotQ, True, "LogPlot supplies functions"}, {False, 
   True}},
 Delimiter,
 {{displayMoneyStocks, lsMoneyKeys, "Money stocks to display:"}, 
  lsMoneyKeys, ControlType -> TogglerBar},
 {{moneyTogetherQ, True, "Plot money functions together"}, {False, 
   True}},
 {{moneyDerivativesQ, False, 
   "Plot money functions derivatives"}, {False, True}},
 {{moneyLogPlotQ, True, "LogPlot money functions"}, {False, True}},
 {{nPlotColumns, 1, "Number of plot columns"}, Range[5]},
 ControlPlacement -> Left, ContinuousAction -> False]

Sensitivity analysis

When making and using this kind of dynamics models it is important to see how the solutions react to changes of different parameters. For example, we should try to find answers to questions like “What ranges of which parameters bring dramatic changes into important stocks?”

Sensitivity Analysis (SA) is used to determine how sensitive is a SD model to changes of the parameters and to changes of model’s equations, [BC1]. More specifically, parameter sensitivity, which we apply below, allows us to see the changes of stocks dynamic behaviour for different sequences (and combinations) of parameter values.

Remark: This section to mirrors to a point the section with same name in [AA4], except in this notebook we are more interested in medical supplies related SA because quarantine related SA is done in [AA4].

Remark: SA shown below should be done for other stocks and rates. In order to keep this exposition short we focus on ISSP, DIP, and HP. Also, it is interesting to think in terms of “3D parameter sensitivity plots.” We also do such plots.

Evaluations by Area under the curve

For certain stocks we might be not just interested in their evolution in time but also in their cumulative values. I.e. we are interested in the so called Area Under the Curve (AUC) metric for those stocks.

There are three ways to calculate AUC for stocks of interest:

Add aggregation equations in the system of ODE’s. (Similar to the stock DIP in SEI2HR.)
- For example, in order to compute AUC for ISSP we can add to SEI2HR the equation:

(*aucISSP'[t] = ISSP[t]*)

- More details for such equation addition are given in [AA2].

Apply NIntegrate over stocks solution functions.
Apply Trapezoidal rule to stock solution function values over a certain time grid.

Below we use 1 and 3.

Variation of medical supplies delivery period

Here we calculate the solutions for a certain combination of capacities and rates:

AbsoluteTiming[
 aVarSolutions =
   Association@
    Map[
     Function[{msdpVar},
      model2 = SEI2HREconModel[t];
      model2 = 
       SetRateRules[
        model2, <|\[Kappa][MS] -> 10000, \[Kappa][HMS] -> 100, 
         mspr[HB] -> 100, msdp[HB] -> msdpVar|>];
      msdpVar -> Association[ModelNDSolve[model2, {t, 365}][[1]]]
      ],
     Union[Join[Range[0.2, 1, 0.2], Range[1, 3, 0.5]]]
     ];
 ]

(*{0.231634, Null}*)

As expected more frequent delivery results in fuller utilization of the non-occupied hospital beds:

SeedRandom[23532]
focusStock = HP;
aVals = #[focusStock][Range[0, 365]] & /@ aVarSolutions;
ListLinePlot[
 KeyValueMap[
  Callout[Tooltip[#2, #1], #1, {If[#1 < 1, RandomInteger[{120, 150}], 
      RandomInteger[{160, 260}]], Above}] &, aVals], 
 PlotLegends -> 
  SwatchLegend[Keys[aVals], 
   LegendLabel -> "Medical supplies\ndelivery period"], 
 PlotRange -> All, ImageSize -> Large, 
 PlotLabel -> focusStock[t] /. modelSEI2HREcon["Stocks"]]

Here are the corresponding AUC values:

aAUCs = TrapezoidalRule[Transpose[{Range[0, 365], #}]] & /@ aVals;
ResourceFunction["GridTableForm"][aAUCs]

BarChart[aAUCs, ChartLabels -> Keys[aAUCs], ColorFunction -> "Pastel",
  PlotLabel -> 
  Row[{focusStock[t] /. modelSEI2HREcon["Stocks"], Spacer[5], "AUC"}]]

Variation of medical supplies production rate

In order to demonstrate the effect of medical supplies production rate (mspr) it is beneficial to eliminate the hospital beds availability restriction – we assume that we have enough hospital beds for all infected severely symptomatic people.

Here we calculate the solutions for a certain combination of capacities and rates:

AbsoluteTiming[
 aVarSolutions =
   Association@
    Map[
     Function[{msprVar},
      model2 = SEI2HREconModel[t];
      model2 = 
       SetRateRules[
        model2, <|\[Kappa][MS] -> 100000, \[Kappa][HMS] -> 
          10000, \[Kappa][MSD] -> 1000, mspr[HB] -> msprVar, 
         msdp[HB] -> 1.5, mscr[ISSP] -> 0.2, mscr[TP] -> 0.001, 
         mscr[ISSP] -> 1, nhbr[TP] -> 100/1000|>];
      msprVar -> Association[ModelNDSolve[model2, {t, 365}][[1]]]
      ],
     {20, 60, 100, 200, 300, 1000, 10000}
     ];
 ]

(*{0.156794, Null}*)

Hospitalized Population

As expected we can see that with smaller production rates we get less hospitalized people:

SeedRandom[1232]
focusStock = HP;
aVals = #[focusStock][Range[0, 365]] & /@ aVarSolutions;
ListLinePlot[
 KeyValueMap[
  Callout[Tooltip[#2, #1], #1, {RandomInteger[{180, 240}], Above}] &, 
  aVals], PlotLegends -> 
  SwatchLegend[Keys[aVals], 
   LegendLabel -> "Medical supplies\nproduction rate"], 
 PlotRange -> All, ImageSize -> Large, 
 PlotLabel -> focusStock[t] /. modelSEI2HREcon["Stocks"]]

Here are the corresponding AUC values:

aAUCs = TrapezoidalRule[Transpose[{Range[0, 365], #}]] & /@ aVals;
ResourceFunction["GridTableForm"][aAUCs]

BarChart[aAUCs, ChartLabels -> Keys[aAUCs], ColorFunction -> "Pastel",
  PlotLabel -> 
  Row[{focusStock[t] /. modelSEI2HREcon["Stocks"], Spacer[5], "AUC"}]]

Medical Supplies

Here we plot the availability of MS at MS producer’s storage:

SeedRandom[821]
focusStock = MS;
aVals = #[MS][Range[0, 365]] & /@ aVarSolutions;
ListLinePlot[
 KeyValueMap[
  Callout[Tooltip[#2, #1], #1, {RandomInteger[{100, 160}], Above}] &, 
  aVals], PlotLegends -> 
  SwatchLegend[Keys[aVals], 
   LegendLabel -> "Medical supplies\nproduction rate"], 
 PlotRange -> All, ImageSize -> Large, 
 PlotLabel -> focusStock[t] /. modelSEI2HREcon["Stocks"]]

Here are the corresponding AUC values:

aAUCs = TrapezoidalRule[Transpose[{Range[0, 365], #}]] & /@ aVals;
ResourceFunction["GridTableForm"][aAUCs]

BarChart[aAUCs, 
 ChartLabels -> Map[Rotate[ToString[#], \[Pi]/6] &, Keys[aAUCs]], 
 ColorFunction -> "Pastel", 
 PlotLabel -> 
  Row[{focusStock[t] /. modelSEI2HREcon["Stocks"], Spacer[5], "AUC"}],
  ImageSize -> Medium]

Variation of delivery disruption starts

Here we compute simulation solutions for a set of delivery disruption starts using disruption length of 7 days and disruption “slowing down” factor 2:

AbsoluteTiming[
 aVarSolutions =
   Association@
    Map[
     Function[{ddsVar},
      ddsVar -> 
       Association[
        ModelNDSolve[
          SetRateRules[
           modelSEI2HREcon, <|\[Kappa][MS] -> 100000, \[Kappa][HMS] ->
              1000, mspr[HB] -> 100, ql -> 56, qsd -> 60, 
            nhbr[TP] -> 3/1000, dbp -> 1, dds -> ddsVar, ddl -> 7, 
            dsf -> 2|>], {t, 365}][[1]]]
      ],
     Append[Range[40, 120, 20], 365]
     ];
 ]

(*{0.45243, Null}*)

Note, that disruption start at day 365 means no disruption. Also, we use three hospital beds per thousand people.

Here we plot the results for HP only:

SeedRandom[009]
focusStock = HP;
aVals = #[focusStock][Range[0, 365]] & /@ aVarSolutions;
ListLinePlot[
 KeyValueMap[
  Callout[Tooltip[#2, #1], #1, {RandomInteger[{60, 140}], Bottom}] &, 
  aVals], PlotLegends -> 
  SwatchLegend[Keys[aVals], 
   LegendLabel -> "Medical supplies\ndisruption start"], 
 PlotRange -> All, ImageSize -> Large, 
 PlotLabel -> focusStock[t] /. modelSEI2HREcon["Stocks"]]

Here are the corresponding AUC values:

aAUCs = TrapezoidalRule[Transpose[{Range[0, 365], #}]] & /@ aVals;
ResourceFunction["GridTableForm"][aAUCs]

BarChart[aAUCs, 
 ChartLabels -> Map[Rotate[ToString[#], \[Pi]/6] &, Keys[aAUCs]], 
 ColorFunction -> "Pastel", 
 PlotLabel -> 
  Row[{focusStock[t] /. modelSEI2HREcon["Stocks"], Spacer[5], "AUC"}],
  ImageSize -> Medium]

Combined variability of delivery start and disruption

Here we calculate the solutions for a set of combinations of delivery periods and delivery disruption starts:

AbsoluteTiming[
 aVarSolutions =
   Association@
    Map[
     Function[{par},
      model2 = modelSEI2HREcon;
      model2 = 
       SetRateRules[
        model2, <|\[Kappa][MS] -> 100000, \[Kappa][HMS] -> 10000, 
         mspr[HB] -> 1000, dbp -> par[[1]], dds -> par[[2]], ddl -> 7,
          dsf -> 4, nhbr[TP] -> 3/1000|>];
      par -> Association[ModelNDSolve[model2, {t, 365}][[1]]]
      ],
     Flatten[Outer[List, {0.5, 1, 1.5}, {60, 100, 365}], 1]
     ];
 ]

(*{0.759922, Null}*)

As expected more frequent, less disrupted delivery brings fuller utilization of the non-occupied hospital beds:

SeedRandom[3233]
focusStock = HP;
aVals = #[focusStock][Range[0, 365]] & /@ aVarSolutions;
ListLinePlot[
 KeyValueMap[
  Callout[Tooltip[#2, ToString[#1]], 
    ToString[#1], {RandomInteger[{60, 160}], Left}] &, aVals], 
 PlotLegends -> 
  SwatchLegend[ToString /@ Keys[aVals], 
   LegendLabel -> 
    "Medical supplies\ndelivery period & disruption start"], 
 PlotRange -> All, ImageSize -> Large, 
 PlotLabel -> focusStock[t] /. modelSEI2HREcon["Stocks"]]

Here are the corresponding AUC values:

aAUCs = TrapezoidalRule[Transpose[{Range[0, 365], #}]] & /@ aVals;
ResourceFunction["GridTableForm"][aAUCs]

BarChart[aAUCs, 
 ChartLabels -> Map[Rotate[ToString[#], \[Pi]/6] &, Keys[aAUCs]], 
 ColorFunction -> "Pastel", 
 PlotLabel -> 
  Row[{focusStock[t] /. modelSEI2HREcon["Stocks"], Spacer[5], "AUC"}],
  ImageSize -> Medium]

SeedRandom[3233]
focusStock = DIP;
aVals = #[focusStock][Range[0, 365]] & /@ aVarSolutions;
ListLinePlot[
 KeyValueMap[
  Callout[Tooltip[#2, ToString[#1]], 
    ToString[#1], {RandomInteger[{60, 160}], Left}] &, aVals], 
 PlotLegends -> 
  SwatchLegend[ToString /@ Keys[aVals], 
   LegendLabel -> 
    "Medical supplies\ndelivery period & disruption start"], 
 PlotRange -> All, ImageSize -> Large, 
 PlotLabel -> focusStock[t] /. modelSEI2HREcon["Stocks"]]

ResourceFunction["GridTableForm"][Last /@ aVals]

BarChart[Last /@ aVals, 
 ChartLabels -> Map[Rotate[ToString[#], \[Pi]/6] &, Keys[aAUCs]], 
 ColorFunction -> "Pastel", 
 PlotLabel -> "Deceased Population at day 365", ImageSize -> Medium]

References

Articles

[Wk1] Wikipedia entry, “Compartmental models in epidemiology”.

[Wl2] Wikipedia entry, “Coronavirus disease 2019”.

[HH1] Herbert W. Hethcote (2000). “The Mathematics of Infectious Diseases”. SIAM Review. 42 (4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.

[BC1] Lucia Breierova, Mark Choudhari, An Introduction to Sensitivity Analysis, (1996), Massachusetts Institute of Technology.

[AA1] Anton Antonov, “Coronavirus propagation modeling considerations”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, “Basic experiments workflow for simple epidemiological models”, (2020), SystemModeling at GitHub.

[AA3] Anton Antonov, “Scaling of Epidemiology Models with Multi-site Compartments”, (2020), SystemModeling at GitHub.

[AA4] Anton Antonov, “SEI2HR model with quarantine scenarios”, (2020), SystemModeling at GitHub.

Repositories, packages

[WRI1] Wolfram Research, Inc., “Epidemic Data for Novel Coronavirus COVID-19”, WolframCloud.

[AAr1] Anton Antonov, Coronavirus propagation dynamics project, (2020), SystemModeling at GitHub.

[AAr2] Anton Antonov, “Epidemiology Compartmental Modeling Monad in R”, (2020), ECMon-R at GitHub.

[AAp1] Anton Antonov, “Epidemiology models Mathematica package”, (2020), SystemModeling at GitHub.

[AAp2] Anton Antonov, “Epidemiology models modifications Mathematica package”, (2020), SystemModeling at GitHub.

[AAp3] Anton Antonov, “Epidemiology modeling visualization functions Mathematica package”, (2020), SystemModeling at GitHub.

[AAp4] Anton Antonov, “Epidemiology modeling simulation functions Mathematica package”, (2020), SystemModeling at GitHub.

[AAp5] Anton Antonov, “System dynamics interactive interfaces functions Mathematica package”, (2020), SystemsModeling at GitHub.

Coronavirus propagation modeling tutorial presentation, OMLDS March 2020

Posted on March 30, 2020 by Anton Antonov Antonov

Last Saturday I gave a presentation at Orlando Machine Learning and Data Science (OMLDS) Meetup.

The presentation abstract heavily borrowed descriptions and plans from a similar presentation to be given for the useR! Boston Meetup in April.

This presentation was almost entirely based on frameworks, simulations, and graphics made with Mathematica. For some parts Wolfram System Modeler was used.

Here is the (main) presentation mind-map:

(Note that mind-map’s PDF has hyperlinks.)

The presentation was given online (because of COVID-19) using Zoom. The maximum number of people registered, 100. (Many were “first timers.”) Nearly 60 showed up (and stayed throughout.)

Here is a link to the video recording: https://youtu.be/odcoi9stYuY .

References

WirVsVirus 2020 hackathon participation

Posted on March 24, 2020 by Anton Antonov Antonov

Introduction

Last weekend – 2020-03-20 ÷ 2020-03-22 – I participated in the (Germany-centric) hackathon WirVsVirus. (I friend of mine who lives in Germany asked me to team up and sign up.)

Our idea proposal was accepted, listed in the dedicated overview table (see item 806). The title of our hackathon project is:

“Geo-spatial-temporal Economic Model for COVID-19 Propagation and Management in Germany”

Nearly a dozen of people enlisted to help. (We communicated through Slack.)

WebImage["https://devpost.com/software/geo-raumlich-zeitliches-\
wirtschaftsmodell-fur-covid-19"]

Multiple people helped with the discussion of ideas, directions where to find data, with actual data gathering, and related documented analysis. Of course, just discussing the proposed solutions was already a great help!

What was accomplished

Work plans

The following mind-map reflects pretty well what was planned and done:

There is also a related org-mode file with the work plan.

Data

I obtained Germany city data with Mathematica’s build-in functions and used it to heuristically derive a traveling patterns graph, [AA1].

Here is the data:

dsCityRecords = 
  ResourceFunction["ImportCSVToDataset"][
   "https://raw.githubusercontent.com/antononcube/SystemModeling/\
master/Data/dfGermanyCityRecords.csv"];
Dimensions[dsCityRecords]

(*{12538, 6}*)

Here is Geo-histogram of that data:

aCoordsToPopulations = 
  AssociationThread[
   Values /@ Normal[dsCityRecords[All, {"Lat", "Lon"}]], 
   Normal[dsCityRecords[All, "Population"]]];
grHist = GeoHistogram[aCoordsToPopulations, cellRadius, 
   ColorFunction -> (Opacity[#, Blue] &), PlotLegends -> Automatic];
If[TrueQ[renderGraphPlotsQ], grHist]

We considered a fair amount of other data. But because of the time limitations of the hackathon we had to use only the one above.

Single-site models

During the development phase I used the model SEI2R, but since we wanted to have a “geo-spatial-temporal epidemiological economics model” I productized the implementation of SEI2HR-Econ, [AAp1].

Here are the stocks, rates, and equations of SEI2HR-Econ:

Magnify[ModelGridTableForm[SEI2HREconModel[t]], 0.85]

Multi-site SEI2R (SEI2HR-Econ) over a hexagonal grid graph

I managed to follow through with a large part of the work plan for the hackathon and make multi-site scaled model that “follows the money”, [AA1]. Here is a diagram that shows the travelling patterns graph and solutions at one of the nodes:

Here is (a snapshot of) an interactive interface for studying and investigating the solution results:

For more details see the notebook [AA1]. Different parameters can be set in the “Parameters” section. Especially of interest are the quarantine related parameters: start, duration, effect on contact rates and traffic patterns.

I also put in that notebook code for exporting simulations results and programmed visualization routines in R, [AA2]. (In order other members of team to be able to explore the results.)

References

[DP1] 47_wirtschaftliche Auswirkung_Geo-spatial-temp-econ-modell, DevPost.

[WRI1] Wolfram Research, Inc., Germany city data records, (2020), SystemModeling at GitHub.

[AA1] Anton Antonov, “WirVsVirus hackathon multi-site SEI2R over a hexagonal grid graph”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, “WirVsVirus-Hackathon in R”, (2020), SystemModeling at GitHub.

[AAp1] Anton Antonov, “Epidemiology models Mathematica package”, (2020), SystemsModeling at GitHub.

Scaling of epidemiology models with multi-site compartments

Posted on March 13, 2020 by Anton Antonov Antonov

Version 1.0

Introduction

In this notebook we describe and exemplify an algorithm that allows the specification and execution geo-spatial-temporal simulations of infectious disease spreads. (Concrete implementations and examples are given.)

The assumptions of the typical compartmental epidemiological models do not apply for countries or cities that are non-uniformly populated. (For example, China, USA, Russia.) There is a need to derive epidemiological models that take into account the non-uniform distribution of populations and related traveling patterns within the area of interest.

Here is a visual aid (made with a random graph over the 30 largest cities of China):

In this notebook we show how to extend core, single-site epidemiological models into larger models for making spatial-temporal simulations. In the explanations and examples we use SEI2R, [AA2, AAp1], as a core epidemiological model, but other models can be adopted if they adhere to the model data structure of the package “EpidemiologyModels.m”, [AAp1].

From our experiments with we believe that the proposed multi-site extension algorithm gives a modeling ingredient that is hard emulate by other means within single-site models.

Definitions

Single-site: A geographical location (city, neighbourhood, campus) for which the assumptions of the classical compartmental epidemiological models hold.

Single site epidemiological model: A compartmental epidemiological model for a single site. Such model has a system of Ordinary Differential Equations (ODE’s) and site dependent initial conditions.

Multi-site area: An area comprised of multiple single sites with known traveling patterns between them. The area has a directed graph $G$ with nodes that correspond to the sites and a positive matrix $\text{\textit{tpm}}(G)$ for the traveling patterns between the sites.

Problem definition: Given (i) a single site epidemiological model $M$ , (ii) a graph $G$ connecting multiple sites, and (iii) a traveling patterns matrix $\text{\textit{tpm}}(G)$ between the nodes of $G$ derive an epidemiological model $S(M,\text{\textit{tpm}}(G))$ that simulates more adequately viral decease propagation over $G$ .

Multi-Site Epidemiological Model Extension Algorithm (MSEMEA): An algorithm that derives from a given single site epidemiological model and multi-site area an epidemiological model that can be used to simulate the geo-spatial-temporal epidemics and infectious disease spreads. (The description of MSEMEA is the main message of this notebook.)

Load packages

The epidemiological models framework used in this notebook is implemented with the packages [AAp1, AAp2, AA3]; the interactive plots functions are from the package [AAp4].

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/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelingVisualizationFunctions.m"]
Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/WL/SystemDynamicsInteractiveInterfacesFunctions.m"]

Notebook structure

The section “General algorithm description” gives rationale and conceptual steps of MSEMEA.

The next two sections of the notebook follow the procedure outline using the SEI2R model as $M$ , a simple graph with two nodes as $G$ , and both constant and time-dependent matrices for $\text{\textit{tpm}}(G)$ .

The section “Constant traveling patterns over a grid graph” presents an important test case with a grid graph that we use to test and build confidence in MSEMEA. The sub-section “Observations” is especially of interest.

The section “Time-dependent traveling patterns over a random graph” presents a nearly “real life” application of MSEMEA using a random graph and a time dependent travelling patterns matrix.

The section “Money from lost productivity” shows how to track the money losses across the sites.

The last section “Future plans” outlines envisioned (immediate) extensions work presented in this notebook.

General algorithm description

In this section we describe a modeling approach that uses different mathematical modeling approaches for (i) the multi-site travelling patterns and (ii) the single-site disease spread interactions, and then (iii) unifies them into a common model.

Splitting and scaling

The traveling between large, densely populated cities is a very different process of the usual people mingling in those cities. The usual large, dense city mingling is assumed and used in the typical compartmental epidemiological models. It seems it is a good idea to split the two processes and derive a common model.

Assume that all journeys finish within a day. We can model the people arriving (flying in) into a city as births, and people departing a city as deaths.

Let as take a simple model like SIR or SEIR and write the equation for every site we consider. This means for every site we have the same ODE’s with site-dependent initial conditions.

Consider the traveling patterns matrix $K$ , which is a contingency matrix derived from source-destination traveling records. (Or the adjacency matrix of a travelling patterns graph.) The matrix entry of $K(i,j)$ tells us how many people traveled from site $i$ to site $j$ . We systematically change the ODE’s of the sites in following way.

Assume that site $a$ had only travelers coming from site $b$ and going to site $b$ . Assume that the Total Population (TP) sizes for sites $a$ and $b$ are $N_a$ and $N_b$ respectively. Assume that only people from the Susceptible Population (SP) traveled. If the adopted single-site model is SIR, [Wk1], we take the SP equation of site $a$

$\text{SP}_a'(t)=-\frac{\beta \text{IP}_a(t) \text{SP}_a(t)}{N_a}-\text{SP}_a(t) \mu$

and change into the equation

$\text{SP}_a'(t)=-\frac{\beta \text{IP}_a(t) \text{SP}_a(t)}{N_a}-\text{SP}_a(t) \mu -\frac{K(a,b)\text{SP}_a(t)}{N_a}+\frac{K(b,a)\text{SP}_b(t)}{N_b},$

assuming that

$\frac{K(a,b)\text{SP}_a(t)}{N_a}\leq N_a ,\frac{K(b,a)\text{SP}_b(t)}{N_b}\leq N_b.$

Remark: In the package [AAp3] the transformations above are done with the more general and robust formula:

$\min \left(\frac{K(i,j)\text{SP}_i(t)}{\text{TP}_i(t)},\text{TP}_i(t)\right).$

The transformed systems of ODE’s of the sites are joined into one “big” system of ODE’s, appropriate initial conditions are set, and the “big” ODE system is solved. (The sections below show concrete examples.)

Steps of MSEMEA

MSEMEA derives a compartmental model that combines (i) a graph representation of multi-site traveling patterns with (ii) a single-site compartmental epidemiological model.

Here is a visual aid for the algorithm steps below:

Get a single-site epidemiological compartmental model data structure, .
1. The model data structure has stocks and rates dictionaries, equations, initial conditions, and prescribed rate values; see [AA2, AAp1].
Derive the site-to-site traveling patterns matrix $K$ for the sites in the graph $G$ .
For each of node of make a copy of the model and denote with .
1. In general, the models $M[i], i\in G$ have different initial conditions.
2. The models can also have different death rates, contact rates, etc.
Combine the models into the scaled model .
1. Change the equations of $M[i]$ , $i\in G$ to reflect the traveling patterns matrix $K$ .
2. Join the systems of ODE’s of $M[i]$ , $i\in G$ into one system of ODE’s.
Set appropriate or desired initial conditions for each of the populations in $S$ .
Solve the ODE’s of $S$ .
Visualize the results.

Precaution

Care should be taken when specifying the initial conditions of MSEMEA’s system of ODE’s (sites’ populations) and the traveling patterns matrix. For example, the simulations can “blow up” if the traveling patterns matrix values are too large. As it was indicated above, the package [AAp3] puts some safe-guards, but in our experiments with random graphs and random traveling patterns matrices occasionally we still get “wild” results.

Analogy with Large scale air-pollution modeling

There is a strong analogy between MSEMEA and Eulerian models of Large Scale Air-Pollution Modeling (LSAPM), [AA3, ZZ1].

The mathematical models of LSAPM have a “chemistry part” and an “advection-diffusion part.” It is hard to treat such mathematical model directly – different kinds of splittings are used. If we consider 2D LSAPM then we can say that we cover the modeling area with steer tank reactors, then with the chemistry component we simulate the species chemical reactions in those steer tanks, and with the advection-diffusion component we change species concentrations in the steer tanks (according to some wind patterns.)

Similarly, with MSEMEA we separated the travel of population compartments from the “standard” epidemiological modeling interaction between the population compartments.

Similarly to the so called “rotational test” used in LSAPM to evaluate numerical schemes, we derive and study the results of “grid graph tests” for MSEMEA.

Single site epidemiological model

Here is the SEI2R model from the package [AAp1]:

model1 = SEI2RModel[t, "InitialConditions" -> True, 
   "RateRules" -> True, 
   "TotalPopulationRepresentation" -> "AlgebraicEquation"];
ModelGridTableForm[model1]

Here we endow the SEI2R model with a (prominent) ID:

ModelGridTableForm[AddModelIdentifier[model1, 1]]

Thus we demonstrated that we can do Step 3 of MSEMEA.

Below we use ID’s that correspond to the nodes of graphs (and are integers.)

Scaling the single-site SIR model over a small complete graph

Constant travel matrices

Assume we have two sites and the following graph and matrix describe the traveling patterns between them.

Here is the graph:

gr = CompleteGraph[2, DirectedEdges -> True, GraphLayout -> "SpringElectricalEmbedding"]

And here is the traveling patterns matrix:

SeedRandom[44];
matTravel = AdjacencyMatrix[gr]*RandomInteger[{100, 1000}, {VertexCount[gr], VertexCount[gr]}];
MatrixForm[matTravel]

Note that there are much more travelers from 1 to 2 than from 2 to 1.

Here we obtain the core, single-site model (as shown in the section above):

model1 = SEI2RModel[t, "InitialConditions" -> True, 
   "RateRules" -> True, 
   "TotalPopulationRepresentation" -> "AlgebraicEquation"];

Make the multi-site compartments model with SEI2R and the two-node travel matrix using the function ToSiteCompartmentsModel of [AAp2]:

modelBig = 
  ToSiteCompartmentsModel[model1, matTravel, 
   "MigratingPopulations" -> {"Susceptible Population", "Exposed Population", "Infected Normally Symptomatic Population","Recovered Population"}];

Show the unique stocks in the multi-site model:

GetPopulationSymbols[modelBig, __ ~~ __]

(*{TP[1], SP[1], EP[1], INSP[1], ISSP[1], RP[1], MLP[1], TP[2], SP[2], 
 EP[2], INSP[2], ISSP[2], RP[2], MLP[2]}*)

From the symbolic form of the multi-model equations derive the specific equations with the adopted rate values:

ModelGridTableForm[
 KeyTake[modelBig, {"Equations"}] //. modelBig["RateRules"]]

Show the initial conditions:

RandomSample[modelBig["InitialConditions"], UpTo[12]]

(*{ISSP[2][0] == 1, TP[1][0] == 100000, EP[2][0] == 0,
  EP[1][0] == 0, SP[1][0] == 99998, RP[2][0] == 0, 
 RP[1][0] == 0, INSP[1][0] == 1, INSP[2][0] == 1, 
 TP[2][0] == 100000, SP[2][0] == 99998, 
 MLP[1][0] == 0}*)

Show the total number of equations:

Length[modelBig["Equations"]]

(*14*)

Solve the system of ODE’s of the extended model:

maxTime = 120;
AbsoluteTiming[
  aSol = Association@First@
      NDSolve[
       Join[modelBig["Equations"] //. modelBig["RateRules"], 
        modelBig["InitialConditions"]],
       GetStockSymbols[modelBig, __ ~~ "Population"],
       {t, 0, maxTime}
       ];
  ];
Length[aSol]

(*12*)

Display the solutions for each site separately:

ParametricSolutionsPlots[modelBig["Stocks"], #, None, maxTime, 
   "Together" -> True, PlotTheme -> "Detailed", 
   ImageSize -> Medium] & /@ GroupBy[Normal@aSol, #[[1, 1]] &, Association]

From the plots above we see that both sites start with total populations of $100000$ people. Because more travelers go from 1 to 2 we see that the exposed, infected, and recovered populations are larger at 2.

Time dependent travel matrices

Instead of using constant traveling patterns matrices we can use matrices with time functions as entries. It is instructive to repeat the computations above using this matrix:

SeedRandom[232]
matTravel2 = matTravel*Table[Abs[Sin[RandomReal[{0.01, 0.1}] t]], VertexCount[gr], VertexCount[gr]];
MatrixForm[matTravel2]

Here are the corresponding number of traveling people functions:

Plot[Evaluate[DeleteCases[Flatten@Normal@matTravel2, 0]], {t, 0, 120},
  PlotTheme -> "Detailed"]

Here we scale the SIR model, solve the obtained system of ODE’s, and plot the solutions:

modelBig = 
  ToSiteCompartmentsModel[model1, matTravel2, 
   "MigratingPopulations" -> {"Susceptible Population", "Exposed Population", "Infected Normally Symptomatic Population","Recovered Population"}];
aSol = Association@First@
    NDSolve[
     Join[modelBig["Equations"] //. modelBig["RateRules"], 
      modelBig["InitialConditions"]],
     GetStockSymbols[modelBig, __ ~~ ""],
     {t, 0, maxTime}
     ];
ParametricSolutionsPlots[modelBig["Stocks"], #, None, 120, 
    "Together" -> True, PlotTheme -> "Detailed", 
    ImageSize -> Medium][[1]] & /@ 
 GroupBy[Normal@KeySelect[aSol, ! MemberQ[{MLP}, Head[#]] &], #[[1, 1]] &, Association]

Note that the oscillatory nature of the temporal functions in the travelling patterns matrix are reflected in the simulation results.

Constant traveling patterns over a grid graph

In this section we do the model extension and simulation over a regular grid graph with a constant traveling patterns matrix.

Here we create a grid graph with directed edges:

{m, n} = {7, 12};
grGrid = GridGraph[{m, n}, DirectedEdges -> True, GraphLayout -> "SpringEmbedding", VertexLabels -> Automatic, ImageSize -> Large]

Note that:

There is one directed edge between any two edge-connected nodes
All horizontal edges point in one direction
All vertical edges point in one direction
The edges are directed from nodes with smaller indexes to nodes with larger indexes.

Here we make a constant traveling matrix and summarize it:

matGridTravel = AdjacencyMatrix[grGrid]*ConstantArray[1000, {VertexCount[grGrid], VertexCount[grGrid]}];
{ResourceFunction["RecordsSummary"][Flatten[matGridTravel], "All elements"][[1]], 
 ResourceFunction["RecordsSummary"][Select[Flatten[matGridTravel], # > 0 &], "Non-zero elements"][[1]], 
 MatrixPlot[matGridTravel]}

Here we scale the SEI2R model with the grid graph constant traveling matrix:

model1 = SEI2RModel[t, "InitialConditions" -> True, 
   "RateRules" -> True, 
   "TotalPopulationRepresentation" -> "AlgebraicEquation", 
   "BirthsTerm" -> True];

modelGrid = ToSiteCompartmentsModel[model1, matGridTravel, "MigratingPopulations" -> Automatic];

Change the initial conditions in the following way:

Pick initial population size per site (same for all sites)
Make a constant populations vector
At all sites except the first one put the infected populations to zero; the first site has one severely symptomatic person
Set the susceptible populations to be consistent with the total and infected populations.

maxPopulation = 10^6;
lsRPopulations = 
  ConstantArray[maxPopulation, 
   Length[GetPopulationSymbols[modelGrid, "Total Population"]]];
modelGrid =
  SetInitialConditions[
   modelGrid,
   Join[
    Join[AssociationThread[
      Through[GetPopulationSymbols[modelGrid, "Total Population"][0]],
       lsRPopulations], <|TP[1][0] -> maxPopulation|>],
    Join[Association@
      Map[#[0] -> 0 &, 
       GetPopulationSymbols[modelGrid, 
        "Infected Severely Symptomatic Population"]], <|
      ISSP[1][0] -> 1|>],
    Join[Association@
      Map[#[0] -> 0 &, 
       GetPopulationSymbols[modelGrid, 
        "Infected Normally Symptomatic Population"]], <|
      INSP[1][0] -> 0|>],
    Join[AssociationThread[
      Through[GetPopulationSymbols[modelGrid, 
         "Susceptible Population"][0]], lsRPopulations], <|
      SP[1][0] -> maxPopulation - 1|>]
    ]
   ];

Solve the system of ODE’s of the scaled model:

maxTime = 160;
AbsoluteTiming[
 aSolGrid = Association@First@
     NDSolve[
      Join[modelGrid["Equations"] //. modelGrid["RateRules"], 
       modelGrid["InitialConditions"]],
      GetStockSymbols[modelGrid, __ ~~ "Population"],
      {t, 0, maxTime}
      ];
 ]

(*{0.717229, Null}*)

Randomly sample the graph sites and display the solutions separately for each site in the sample:

Multicolumn[
 Table[
  Block[{aSol = KeySelect[aSolGrid, MatchQ[#, _Symbol[i]] &]},
   Plot[Evaluate[Map[#[t] &, Values[aSol]]], {t, 0, maxTime}, 
    PlotRange -> All, GridLines -> All, 
    PlotTheme -> "Scientific", PlotLegends -> Keys[aSol], 
    ImageSize -> 300]
   ], {i, RandomSample[Range[VertexCount[grGrid]], UpTo[9]]}],
 3]

Display solutions of the first and last site:

Multicolumn[
 Table[
  Block[{aSol = KeySelect[aSolGrid, MatchQ[#, _Symbol[i]] &]},
   Plot[Evaluate[Map[#[t] &, Values[aSol]]], {t, 0, maxTime}, 
    PlotRange -> All, GridLines -> All, 
    PlotTheme -> "Scientific", PlotLegends -> Keys[aSol], 
    ImageSize -> 300]
   ], {i, {1, VertexCount[grGrid]}}],
 3]

As expected from the graph structure, we can see in the first site plot that its total population is decreasing – nobody is traveling to the first site. Similarly, we can see in the last site plot that its total population is increasing – nobody leaves the last site.

Graph evolution visualizations

We can visualize the spatial-temporal evolution of model’s populations using sequences of graphs. The graphs in the sequences are copies of the multi-site graph each copy having its nodes colored according to the populations in the solutions steps.

Here is a sub-sequence for the total populations:

EvaluateSolutionsOverGraph[grGrid, modelGrid, "Total Population", \
aSolGrid, {1, maxTime, 12}, "NodeSizeFactor" -> 5, 
 "ColorScheme" -> "TemperatureMap", "Legended" -> True, 
 VertexLabels -> None, ImageSize -> 200]

Here is a sub-sequence for the sum of the infected populations:

EvaluateSolutionsOverGraph[grGrid, modelGrid, {"Infected Normally \
Symptomatic Population", 
  "Infected Severely Symptomatic Population"}, aSolGrid, {1, maxTime, 
  12}, "NodeSizeFactor" -> 5, 
 "ColorScheme" -> "TemperatureMap", "Legended" -> True, 
 VertexLabels -> None, ImageSize -> 200]

Here is a sub-sequence for the recovered population:

EvaluateSolutionsOverGraph[grGrid, modelGrid, "Recovered Population", \
aSolGrid, {1, maxTime, 12}, "NodeSizeFactor" -> 5, 
 "ColorScheme" -> "TemperatureMap", "Legended" -> True, 
 VertexLabels -> None, ImageSize -> 200]

Here is an animation of the sum of the infected populations:

Block[{stocks = {"Infected Normally Symptomatic Population", 
    "Infected Severely Symptomatic Population"}, 
  colorScheme = "TemperatureMap", timeStep = 4},
 Legended[
  ListAnimate[
   EvaluateSolutionsOverGraph[grGrid, modelGrid, stocks, 
    aSolGrid, {1, maxTime, timeStep}, "NodeSizeFactor" -> 6, 
    "ColorScheme" -> colorScheme, ImageSize -> 400]], 
  BarLegend[{colorScheme, 
    MinMax[EvaluateSolutionsOverGraphVertexes[grGrid, modelGrid, 
      stocks, aSolGrid, {1, maxTime, timeStep}]]}]
  ]
 ]

Curve shapes of the globally-aggregated solutions

Let us plot for each graph vertex the sum of the solutions of the two types of infected populations. Here is a sample of those graphs:

SeedRandom[1782];
ListLinePlot[#, PlotTheme -> "Detailed"] & /@ 
 KeySort[RandomSample[
   EvaluateSolutionsOverGraphVertexes[grGrid, 
    modelGrid, {"Infected Normally Symptomatic Population", 
     "Infected Severely Symptomatic Population"}, 
    aSolGrid, {1, maxTime, 4}], 12]]

We can see from the plot above that at the grid vertexes we have typical SEIR curve shapes for the corresponding infected populations.

Let us evaluate the solutions for the infected populations for over all graph vertexes and sum them. Here is the corresponding “globally-aggregated” plot:

ListLinePlot[#, PlotTheme -> "Detailed"] &@
 Total[Values[
   EvaluateSolutionsOverGraphVertexes[grGrid, 
    modelGrid, {"Infected Normally Symptomatic Population", 
     "Infected Severely Symptomatic Population"}, 
    aSolGrid, {1, maxTime, 4}]]]

We can see that the globally aggregated plot has a very different shape than the individual vertex plots. The globally aggregated plot has more symmetric look; the individual vertex plots have much steeper gradients on their left sides.

We can conjecture that a multi-site model made by MSEMEA would capture better real life situations than any single-site model. For example, by applying MSEMEA we might be more successful in our calibration attempts for the Hubei data shown (and calibrated upon) in [AA2.].

Interactive interface

With this interactive interface we see the evolution of all populations across the graph:

Manipulate[
 Block[{aSol = aSolGrid},
  DynamicModule[{gr = grGrid, imageSize = 270, factor = 3.2, 
    maxPopulation = aSolGrid[TP[VertexCount[grGrid]]][maxTime], 
    lsPops, lsRes},
   lsPops = {"Total Population", "Susceptible Population", 
     "Exposed Population", {"Infected Normally Symptomatic \
Population", "Infected Severely Symptomatic Population"}, 
     "Recovered Population"};
   If[! addInfectedPopulationsQ,
    lsPops = Flatten[lsPops]
    ];
   lsRes =
    Map[
     Labeled[
       EvaluateSolutionsOverGraph[grGrid, modelGrid, #, 
        aSolGrid, {time, time, 1}, "NodeSizeFactor" -> factor, 
        "ColorScheme" -> cf, "Legended" -> legendedQ, 
        VertexLabels -> "Name", ImageSize -> imageSize], 
       If[ListQ[#], "Infected Population", #],
       Top
       ] &, lsPops];
   Multicolumn[lsRes, 3, Dividers -> All, 
    FrameStyle -> GrayLevel[0.8]]
   ]],
 {{time, 80, "time:"}, 0, maxTime, 1, Appearance -> {"Open"}},
 {{addInfectedPopulationsQ, True, "sum infected populations:"}, {True,
    False}},
 {{cf, "TemperatureMap", "color function:"}, {"AvocadoColors", 
   "GreenBrownTerrain", "Rainbow", "TemperatureMap"}},
 {{legendedQ, False, "legended:"}, {False, True}}]

Observations

Obviously the simulations over the described grid graph, related constant traveling patterns matrix, and constant populations have the purpose to build confidence in conceptual design of MSEMEA and its implementation.

The following observations agree with our expectations for MSEMEA’s results over the “grid graph test”.

The populations plots at each site resemble the typical plots of SEI2R.
The total population at the first site linearly decreases.
The total population at the last site linearly increases.
The plots of the total populations clearly have gradually increasing gradients from the low index value nodes to the high index value nodes.
For the infected populations there is a clear wave that propagates diagonally from the low index value nodes to the high index value nodes.
1. In the direction of the general “graph flow.“
The front of the infected populations wave is much steeper (gives “higher contrast”) than the tail.
1. This should be expected from the single-site SEI2R plots.
For the recovered populations there is a clear “saturation wave” pattern that indicates that the recovered populations change from $0$ to values close to the corresponding final total populations.
The globally aggregated solutions might have fairly different shapes than the single-site solutions. Hence, we think that MSEMEA gives a modeling ingredient that is hard to replicate or emulate by other means in single-site models.

Time-dependent traveling patterns over a random graph

In this section we apply the model extension and simulation over a random graph with random time-dependent traveling patterns matrix.

SeedRandom[84];
grRandom = 
 RandomGraph[{20, 100}, DirectedEdges -> True, 
  GraphLayout -> "SpringElectricalEmbedding", 
  VertexLabels -> "Name"]

Remark: The computations in this section work with larger random graphs; we use a small graph for more legible presentation of the workflow and results. Also, the computations should clearly demonstrate the ability to do simulations with real life data.

Derive a traveling patterns matrix with entries that are random functions:

Block[{gr = grRandom},
 matRandomTravel = 
   AdjacencyMatrix[gr]*
    RandomInteger[{10, 100}, {VertexCount[gr], VertexCount[gr]}]*
    Table[Sin[RandomReal[{0.02, 0.1}] t], VertexCount[gr], 
     VertexCount[gr]];
 ]

Here is a fragment of the matrix:

Magnify[MatrixForm[
  matRandomTravel[[1 ;; 12, 
   1 ;; 12]]], 0.7]

Summarize and plot the matrix at $t=1$ :

Block[{matTravel = matRandomTravel /. t -> 1},
 {ResourceFunction["RecordsSummary"][Flatten[matTravel], 
    "All elements"][[1]], 
  ResourceFunction["RecordsSummary"][
    Select[Flatten[matTravel], # > 0 &], 
    "Non-zero elements"][[1]], 
  MatrixPlot[matTravel]}
 ]

Here we scale the SEI2R model with the random traveling matrix:

model1 = SEI2RModel[t, "InitialConditions" -> True, 
   "RateRules" -> True, 
   "TotalPopulationRepresentation" -> "AlgebraicEquation"];

modelRandom = 
  ToSiteCompartmentsModel[model1, matRandomTravel, 
   "MigratingPopulations" -> {"Susceptible Population", 
     "Exposed Population", "Infected Normally Symptomatic Population",
      "Recovered Population"}];

Change the initial conditions in the following way:

Pick maximum population size per site
Derive random populations for the sites
At all sites except the first one put the infected populations to zero; the first site has one severely symptomatic person
Set the susceptible populations to be consistent with the total and infected populations.

maxPopulation = 10^6;
lsRPopulations = 
  RandomReal[{100000, maxPopulation}, 
   Length[GetPopulationSymbols[modelRandom, "Total Population"]]];
modelRandom =
  SetInitialConditions[
   modelRandom,
   Join[
    Join[AssociationThread[
      Through[GetPopulationSymbols[modelRandom, "Total Population"][
        0]], lsRPopulations], <|TP[1][0] -> maxPopulation|>],
    Join[Association@
      Map[#[0] -> 0 &, 
       GetPopulationSymbols[modelRandom, 
        "Infected Severely Symptomatic Population"]], <|
      ISSP[1][0] -> 1|>],
    Join[Association@
      Map[#[0] -> 0 &, 
       GetPopulationSymbols[modelRandom, 
        "Infected Normally Symptomatic Population"]], <|
      INSP[1][0] -> 0|>],
    Join[AssociationThread[
      Through[GetPopulationSymbols[modelRandom, 
         "Susceptible Population"][0]], lsRPopulations], <|
      SP[1][0] -> maxPopulation - 1|>]
    ]
   ];

Here solve the obtained system of ODE’s:

maxTime = 120;
AbsoluteTiming[
 aSolRandom = Association@First@
     NDSolve[
      Join[modelRandom["Equations"] //. modelRandom["RateRules"], 
       modelRandom["InitialConditions"]],
      GetStockSymbols[modelRandom, __ ~~ "Population"],
      {t, 0, maxTime}
      ];
 ]

(*{0.208188, Null}*)

Here we plot the solutions:

ParametricSolutionsPlots[<||>, #, None, maxTime, 
    "Together" -> True, PlotTheme -> "Detailed", 
    ImageSize -> 
     Small][[1]] & /@ 
 RandomSample[
  GroupBy[Normal@
    aSolRandom, #[[1, 1]] &, 
   Association], UpTo[12]]

Graph evolution visualizations

As in the previous section we can visualize the spatial-temporal evolution of model’s populations using sequences of graphs.

Here is a globally normalized sequence:

EvaluateSolutionsOverGraph[grRandom, modelRandom, {"Recovered \
Population"}, aSolRandom, {1, maxTime, 24}, 
 "NodeSizeFactor" -> 4, "ColorScheme" -> "TemperatureMap", 
 "Normalization" -> "Global", "Legended" -> True]

Here is a locally normalized (“by vertex”) sequence:

EvaluateSolutionsOverGraph[grRandom, modelRandom, {"Recovered \
Population"}, aSolRandom, {1, maxTime, 24}, 
 "NodeSizeFactor" -> 4, "ColorScheme" -> "TemperatureMap", 
 "Normalization" -> "ByVertex", "Legended" -> True]

Money from lost productivity

The model SEI2R from [AAp1] has the stock “Money from Lost Productivity” shown as $\text{MLP}(t)$ in the equations:

ModelGridTableForm[KeyTake[SEI2RModel[t], "Equations"]]

Here are MLP plots from the two-node graph model:

ParametricSolutionsPlots[modelBig["Stocks"], #, None, 120, 
    "Together" -> True, PlotTheme -> "Detailed", 
    ImageSize -> 
     250][[1]] & /@ 
 GroupBy[Normal@
   KeySelect[aSol, MemberQ[{MLP}, Head[#]] &], #[[1,
     1]] &, Association]

Here we plot the sum of the accumulated money losses:

funcs = Values[KeySelect[aSol, MemberQ[{MLP}, Head[#]] &]];
ListLinePlot[{#, Total[Through[funcs[#]]]} & /@ Range[1, maxTime], 
 PlotTheme -> "Detailed", ImageSize -> 250]

Here is the corresponding “daily loss” (derivative):

funcs = Map[D[#[t], t] &, 
   Values[KeySelect[aSol, MemberQ[{MLP}, Head[#]] &]]];
ListLinePlot[{#, Total[funcs /. t -> #]} & /@ Range[1, maxTime], 
 PlotTheme -> "Detailed", ImageSize -> 250]

Future plans

There are multiple ways to extend the presented algorithm, MSEMEA. Here are a few most immediate ones:

Investigate and describe the conditions under which MSEMEA performs well, and under which it “blows up”
Apply MSEMEA together with single site models that have large economics parts
Do real data simulations related to the spread of COVID-19.

References

Articles, books

[Wk1] Wikipedia entry, “Compartmental models in epidemiology”.

[HH1] Herbert W. Hethcote (2000). “The Mathematics of Infectious Diseases”. SIAM Review. 42 (4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.

[AA1] Anton Antonov, “Coronavirus propagation modeling considerations”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, “Basic experiments workflow for simple epidemiological models”, (2020), SystemModeling at GitHub.

[AA3] Anton Antonov, “Air pollution modeling with gridMathematica”, (2006), Wolfram Technology Conference.

[ZZ1] Zahari Zlatev, Computer Treatment of Large Air Pollution Models. 1995. Kluwer.

Repositories, packages

[WRI1] Wolfram Research, Inc., “Epidemic Data for Novel Coronavirus COVID-19”, WolframCloud.

[AAr1] Anton Antonov, Coronavirus propagation dynamics project, (2020), SystemModeling at GitHub.

[AAp1] Anton Antonov, “Epidemiology models Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp2] Anton Antonov, “Epidemiology models modifications Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp3] Anton Antonov, “Epidemiology modeling visualization functions Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp4] Anton Antonov, “System dynamics interactive interfaces functions Mathematica package”, (2020), SystemsModeling at GitHub.

Basic experiments workflow for simple epidemiological models

Posted on March 11, 2020 by Anton Antonov Antonov

Introduction

The primary purpose of this document (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”
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.)
1. Using the parameters of the previous step.
Using an interactive interface experiment with different values of the parameters.
1. In order to form “qualitative understanding.”
Get real life data.
1. Say, for the 2019-20 coronavirus outbreak.
Attempt manual or automatic calibration of the model.
1. This step will most likely require additional data transformations and programming.
2. 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].

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

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

modelSI2R = SEI2RModel[t, "InitialConditions" -> True, "RateRules" -> True];

We can show a tabulated visualization of the model using the function ModelGridTableForm from [AAp1]:

ModelGridTableForm[modelSI2R]

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

ModelGridTableForm[KeyTake[modelSI2R, "Equations"]]

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

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

(*1*)

Here we make the births term using a birth rate that is the same as the death rate:

birthTerm = GetPopulations[modelSI2R, "Total Population"][[1]] * GetRates[modelSI2R, "Population death rate"][[1]]

(*TP[t] \[Mu][TP]*)

Here we add the births term to the equations of new model

modelSI2RNew = modelSI2R;
modelSI2RNew["Equations"] = 
  ReplacePart[
   modelSI2R["Equations"], {1, 2} -> 
    modelSI2R["Equations"][[1, 
      2]] + birthTerm];

Here we display the equations of the new model:

ModelGridTableForm[KeyTake[modelSI2RNew, "Equations"]]

Adding infected deceased population equation

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

AppendTo[modelSI2R["Equations"], 
  IDP'[t] == \[Mu][INSP]*INSP[t] + \[Mu][ISSP]*ISSP[t]];
AppendTo[modelSI2R["Stocks"], 
  IDP[t] -> "Infected Deceased Population"];
AppendTo[modelSI2R["InitialConditions"], IDP[0] == 0];

Here is how the model looks like:

ModelGridTableForm[KeyTake[modelSI2R, "Equations"]]

Parameters and actual simulation equations code

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

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

Here we set custom rates and initial conditions:

population = 58160000/400;
modelSI2R = SetRateRules[modelSI2R, <|TP[t] -> 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:

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

Simulation

Straightforward simulation for one year with using ParametricNDSolve :

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

(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

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" -> "\[CapitalDelta]", opts];
  
  lsEconPlots =
   ParametricSolutionsPlots[
    modelSI2R["Stocks"],
    KeyTake[aSol, lsEconKeys],
    {aincp, aip, spf, crisp, criap}, ndays,
    "LogPlot" -> econLogPlotQ, "Together" -> econTogetherQ, 
    "Derivatives" -> econDerivativesQ, 
    "DerivativePrefix" -> "\[CapitalDelta]", opts];
  
  lsRestPlots =
   ParametricSolutionsPlots[
    modelSI2R["Stocks"],
    KeyDrop[aSol, Join[lsPopulationKeys, lsEconKeys]],
    {aincp, aip, spf, crisp, criap}, ndays,
    "LogPlot" -> econLogPlotQ, "Together" -> econTogetherQ, 
    "Derivatives" -> econDerivativesQ, 
    "DerivativePrefix" -> "\[CapitalDelta]", opts];
  
  Multicolumn[Join[lsPopulationPlots, lsEconPlots, lsRestPlots], 
   nPlotColumns, Dividers -> All, 
   FrameStyle -> GrayLevel[0.8]]
  ],
 {{aincp, 12., "Average incubation period (days)"}, 1, 60., 1, Appearance -> {"Open"}},
 {{aip, 21., "Average infectious period (days)"}, 1, 100., 1, Appearance -> {"Open"}},
 {{spf, 0.2, "Severely symptomatic population fraction"}, 0, 1, 0.025, Appearance -> {"Open"}},
 {{crisp, 6, "Contact rate of the infected severely symptomatic population"}, 0, 30, 0.1, Appearance -> {"Open"}},
 {{criap, 3, "Contact rate of the infected normally symptomatic population"}, 0, 30, 0.1, Appearance -> {"Open"}},
 {{ndays, 90, "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]

Calibration over real data

It is important to calibrate these kind of models with real data, or at least to give a serious attempt to 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

Here is COVID-19 data taken from [WRI1] for the Chinese province Hubei:

aRealData = <|
     "RecoveredCases" -> {28, 28, 31, 
     32, 42, 45, 80, 88, 90, 141, 168, 295, 386, 522, 633, 817, 1115, 
     1439, 1795, 2222, 2639, 2686, 3459, 4774, 5623, 6639, 7862, 9128,
     10337, 11788, 11881, 15299, 15343, 16748, 18971, 20969, 23383, 
     26403, 28993, 31536, 33934, 36208, 38557, 40592, 42033, 43500, 45235}, 
    "Deaths" -> {17, 17, 24, 40, 52, 76, 125, 125, 162, 
     204, 249, 350, 414, 479, 549, 618, 699, 780, 871, 974, 1068, 
     1068, 1310, 1457, 1596, 1696, 1789, 1921, 2029, 2144, 2144, 2346,
      2346, 2495, 2563, 2615, 2641, 2682, 2727, 2761, 2803, 2835, 
     2871, 2902, 2931, 2959, 2986}, 
   "Infected" -> {399, 399, 494, 689, 964, 1302, 3349, 3341, 
     4651, 5461, 6736, 10532, 12722, 15677, 18483, 20677, 23139, 
     24881, 26965, 28532, 29659, 29612, 43437, 48175, 49030, 49847, 
     50338, 50633, 49665, 48510, 48637, 46439, 46395, 45044, 43252, 
     41603, 39572, 36829, 34617, 32610, 30366, 28174, 25904, 23972, 
     22628, 21207, 19486}|>;

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 have to be padded with a certain number of 0’s 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

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.

Clear[PadRealData];
PadRealData[aData : Association[ (_String -> _?VectorQ) ..], 
   incubationPeriod_?IntegerQ, infectiousPeriod_?IntegerQ] :=
  
  Block[{},
   Join[ConstantArray[0, incubationPeriod + infectiousPeriod], #] & /@
     aData
   ];

opts = {PlotRange -> All, PlotLegends -> None, PlotTheme -> "Detailed", PerformanceGoal -> "Speed", ImageSize -> 300};
Manipulate[
 DynamicModule[{modelSI2R = modelSI2R, lsActualEquations, aSol, 
   lsPopulationPlots, lsEconPlots, lsRestPlots},
  
  modelSI2R = SetRateRules[modelSI2R, <|TP[t] -> population|>];
  modelSI2R = 
   SetInitialConditions[
    modelSI2R, <|SP[0] -> population - 1, ISSP[0] -> 0, 
     INSP[0] -> 1|>];
  lsActualEquations = 
   Join[modelSI2R["Equations"] //. 
     KeyDrop[modelSI2R["RateRules"], lsFocusParams], 
    modelSI2R["InitialConditions"]];
  aSol =
   Association@Flatten@
     ParametricNDSolve[
      lsActualEquations, {SP, EP, INSP, RP, IDP}, {t, 0, 365}, 
      lsFocusParams];
  
  lsPopulationPlots =
   ParametricSolutionsPlots[
    modelSI2R["Stocks"],
    KeyTake[aSol, GetPopulationSymbols[modelSI2R, __ ~~ "Population"]],
    {aincp, aip, 0, criap, criap}, ndays, "Together" -> True, 
    opts];
  
  Show[lsPopulationPlots[[1]], 
   ListPlot[
    PadRealData[aRealData, Round[aincp + padOffset], 
     Round[aip + padOffset]], PlotStyle -> {Blue, Black, Red}]]
  ],
 {{population, 58160000/600, "Population"}, 58160000/1000, 58160000, 10000, Appearance -> {"Open"}},
 {{padOffset, 0, "real data padding offset"}, -100, 100, 1, Appearance -> {"Open"}},
 {{aincp, 6, "Average incubation period (days)"}, 1, 60, 1, Appearance -> {"Open"}},
 {{aip, 32, "Average infectious period (days)"}, 1, 100, 1, Appearance -> {"Open"}},
 {{criap, 0.8, "Contact rate of the infected normally symptomatic population"}, 0, 30, 0.1, Appearance -> {"Open"}},
 {{ndays, 90, "Number of days"}, 1, 365, 1, Appearance -> {"Open"}},
 ControlPlacement -> Left, ContinuousAction -> False]

Maybe good enough parameters

DynamicModule[{aincp = 6, aip = 32, criap = 0.8, ndays = 90, 
  padOffset = -8, population = 78160}, 
 DynamicModule[{modelSI2R = modelSI2R, lsActualEquations, aSol, 
   lsPopulationPlots, lsEconPlots, lsRestPlots}, 
  modelSI2R = 
   SetRateRules[modelSI2R, Association[TP[t] -> population]]; 
  modelSI2R = 
   SetInitialConditions[modelSI2R, 
    Association[SP[0] -> population - 1, ISSP[0] -> 0, 
     INSP[0] -> 1]]; 
  lsActualEquations = 
   Join[modelSI2R["Equations"] //.\[VeryThinSpace]KeyDrop[
      modelSI2R["RateRules"], lsFocusParams], 
    modelSI2R["InitialConditions"]]; 
  aSol = Association[
    Flatten[ParametricNDSolve[
      lsActualEquations, {SP, EP, INSP, RP, IDP}, {t, 0, 365}, 
      lsFocusParams]]]; 
  lsPopulationPlots = 
   ParametricSolutionsPlots[modelSI2R["Stocks"], 
    KeyTake[aSol, 
     GetPopulationSymbols[modelSI2R, __ ~~ "Population"]], {aincp, 
     aip, 0, criap, criap}, ndays, "Together" -> True, opts]; 
  Show[lsPopulationPlots[[1]], 
   ListPlot[
    PadRealData[aRealData, Round[aincp + padOffset], 
     Round[aip + padOffset]], PlotStyle -> {Blue, Black, Red}]]]]

Basic reproduction number:

Block[{aincp = 6, aip = 32, criap = 0.8, ndays = 90, padOffset = -14, population = 75000, \[Mu] = \[Mu][TP] /. modelSI2R["RateRules"]}, 
   (1/(aincp*(\[Mu] + 1/aincp)))*(criap/(1/aip))
]

(*25.5966*)

DynamicModule[{aincp = 5, aip = 26, criap = 2.3, ndays = 90, 
  padOffset = -14, population = 75000}, 
 DynamicModule[{modelSI2R = modelSI2R, lsActualEquations, aSol, 
   lsPopulationPlots, lsEconPlots, lsRestPlots}, 
  modelSI2R = 
   SetRateRules[modelSI2R, Association[TP[t] -> population]]; 
  modelSI2R = 
   SetInitialConditions[modelSI2R, 
    Association[SP[0] -> population - 1, ISSP[0] -> 0, 
     INSP[0] -> 1]]; 
  lsActualEquations = 
   Join[modelSI2R["Equations"] //.\[VeryThinSpace]KeyDrop[
      modelSI2R["RateRules"], lsFocusParams], 
    modelSI2R["InitialConditions"]]; 
  aSol = Association[
    Flatten[ParametricNDSolve[
      lsActualEquations, {SP, EP, INSP, RP, IDP}, {t, 0, 365}, 
      lsFocusParams]]]; 
  lsPopulationPlots = 
   ParametricSolutionsPlots[modelSI2R["Stocks"], 
    KeyTake[aSol, 
     GetPopulationSymbols[modelSI2R, __ ~~ "Population"]], {aincp, 
     aip, 0, criap, criap}, ndays, "Together" -> True, opts]; 
  Show[lsPopulationPlots[[1]], 
   ListPlot[
    PadRealData[aRealData, Round[aincp + padOffset], 
     Round[aip + padOffset]], PlotStyle -> {Blue, Black, Red}]]]]

Basic reproduction number:

Block[{aincp = 5, aip = 26, criap = 2.3, ndays = 90, padOffset = -14, population = 75000, \[Mu] = \[Mu][TP] /. modelSI2R["RateRules"]}, 
 criap/((aincp*(\[Mu] + 1/aincp))/aip)
]

(*59.7934*)

References

Articles

[Wk1] Wikipedia entry, “Compartmental models in epidemiology”.

[HH1] Herbert W. Hethcote (2000). “The Mathematics of Infectious Diseases”. SIAM Review. 42 (4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.

[AA1] Anton Antonov, “Coronavirus propagation modeling considerations”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, Answer of “Model calibration with phase space data”, (2019), Mathematica StackExchage.

Repositories & packages

[WRI1] Wolfram Research, Inc., “Epidemic Data for Novel Coronavirus COVID-19”, WolframCloud.

[AAr1] Anton Antonov, Coronavirus propagation dynamics project, (2020), SystemModeling at GitHub.

[AAp1] Anton Antonov, “Epidemiology models Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp2] Anton Antonov, “Epidemiology models modifications Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp3] Anton Antonov, “System dynamics interactive interfaces functions Mathematica package”, (2020), SystemsModeling at GitHub.

	EigenHorizons (@Eige… on Independent component analysis…
	Anton Antonov Antono… on Independent component analysis…
	EigenHorizons (@Eige… on Independent component analysis…
	Coronavirus propagat… on Coronavirus propagation modeli…
	Anton Antonov Antono… on SEI2HR model with quarantine…

Screenshots

References

Coronavirus

Modeling

R code/software

References

Coronavirus

Modeling

R code/software

Introduction

The models

SEI2R

SEI2HR

SEI2HR-Econ

SEI2HR-Econ’s place a development plan

Notebook structure

Load packages

SEI2HR-Econ extends SEI2HR

Equations explanations

Verbalization description of the model

Code for the equations

The infected and hospitalized populations

Total Population

Susceptible Population

Births term

Exposed Population

Infected Normally Symptomatic Population

Infected Severely Symptomatic Population

Hospitalized Population

Recovered Population

Deceased Infected Population

Hospital Beds

Hospital Medical Supplies

Medical Supplies

Medical Supplies Demand

Money for Hospital Services

Money from Lost Productivity

Quarantine scenarios

Medical supplies scenarios

Parameters and actual simulation equations code

Simulation

Different quarantine starts

Interactive interface

Sensitivity analysis

Evaluations by Area under the curve

Variation of medical supplies delivery period

Variation of medical supplies production rate

Hospitalized Population

Medical Supplies

Variation of delivery disruption starts

Combined variability of delivery start and disruption

References

Articles

Repositories, packages

References

General

Articles, blog posts

Videos

Introduction

What was accomplished

Work plans

Data

Single-site models

Multi-site SEI2R (SEI2HR-Econ) over a hexagonal grid graph

References

Introduction

Definitions

Load packages

Notebook structure

General algorithm description

Splitting and scaling

Steps of MSEMEA

Precaution

Analogy with Large scale air-pollution modeling

Single site epidemiological model

Scaling the single-site SIR model over a small complete graph

Constant travel matrices

Time dependent travel matrices

Constant traveling patterns over a grid graph

Graph evolution visualizations