SEI2HR model with quarantine scenarios

Version 0.8
Anton Antonov
​MathematicaForPrediction at WordPress​
​SystemModeling at GitHub​
March, April 2020

Introduction

The epidemiology compartmental model, [Wk1], presented in this notebook -- SEI2HR -- deals with the left-most and middle rectangles in this diagram:
In[]:=
ImageResize[Import["https://github.com/antononcube/SystemModeling/raw/master/Projects/Coronavirus-propagation-dynamics/Diagrams/Coronavirus-propagation-simple-dynamics.jpeg"],900]
Out[]=
“SEI2HR” stands for “Susceptible, Exposed, Infected two, Hospitalized, Recovered” (populations.)
In this notebook we also deal with quarantine scenarios.
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 are introduced and explained in the notebook [AA2]. SEI2R differs from the classical SEIR model , [Wk1, HH1], with the following elements:
1. Two separate infected populations: one is "severely symptomatic", the other is "normally symptomatic"
2. The monetary equivalent of lost productivity due to infected or died people is tracked.

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
    • ◼
  • Two types of infected populations: normally symptomatic and severely symptomatic
    • ◼
  • Hospitalized population
    • ◼
  • 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 are tracked
    • ◼
  • The money for hospital services are tracked

    • SEI2HR’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:
    Out[]=

    Notebook structure

    The rest of notebook has the following sequence of sections:
    ◼
  • Package load section
    • ◼
  • SEI2HR structure in comparison of SEI2R
    • ◼
  • Explanations of the equations of SEI2HR
    • ◼
  • Quarantine scenario modeling preparation
    • ◼
  • Parameters and initial conditions setup
    • ◼
  • Populations, hospital beds, quarantine scenarios
    • ◼
  • Parametric simulation solution
    • ◼
  • Interactive interface
    • ◼
  • Sensitivity analysis

    • Load packages

    The epidemiological models framework used in this notebook is implemented with the packages [AAp1, AAp2, AA3]; many of the plot functions are from the package [AAp4].
    In[]:=
    Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModels.m"];Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelModifications.m"];Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/Projects/Coronavirus-propagation-dynamics/WL/EpidemiologyModelingVisualizationFunctions.m"];Import["https://raw.githubusercontent.com/antononcube/SystemModeling/master/WL/SystemDynamicsInteractiveInterfacesFunctions.m"];

    SEI2HR extends SEI2R

    The model SEI2HR is an extension of the model SEI2R, [AA2].
    Here is SEI2R:
    In[]:=
    reprTP="AlgebraicEquation";​​lsModelOpts={"Tooltips"True,TooltipStyle{BackgroundYellow,CellFrameColorGray,FontSize20}};​​modelSEI2R=SEI2RModel[t,"InitialConditions"True,"RateRules"True,"TotalPopulationRepresentation"reprTP];​​ModelGridTableForm[modelSEI2R,lsModelOpts]
    Out[]=
    Stocks
    #
    Symbol
    Description
    1
    TP[t]
    Total Population
    2
    SP[t]
    Susceptible Population
    3
    EP[t]
    Exposed Population
    4
    INSP[t]
    Infected Normally Symptomatic Population
    5
    ISSP[t]
    Infected Severely Symptomatic Population
    6
    RP[t]
    Recovered Population
    7
    MLP[t]
    Money of Lost Productivity
    ,Rates
    #
    Symbol
    Description
    1
    μ[TP]
    Population death rate
    2
    μ[INSP]
    Infected Normally Symptomatic Population death rate
    3
    μ[ISSP]
    Infected Severely Symptomatic Population death rate
    4
    sspf[SP]
    Severely Symptomatic Population Fraction
    5
    β[INSP]
    Contact rate for the normally symptomatic population
    6
    β[ISSP]
    Contact rate for the severely symptomatic population
    7
    aip
    Average infectious period
    8
    aincp
    Average incubation period
    9
    lpcr[ISSP,INSP]
    Lost productivity cost rate (per person per day)
    ,Equations
    #
    Equation
    1
    ′
    SP
    [t]-
    INSP[t]
    SP[t]
    β[INSP]
    TP[t]
    -
    ISSP[t]
    SP[t]
    β[ISSP]
    TP[t]
    -
    SP[t]
    μ[TP]
    2
    ′
    EP
    [t]
    INSP[t]
    SP[t]
    β[INSP]
    TP[t]
    +
    ISSP[t]
    SP[t]
    β[ISSP]
    TP[t]
    -
    EP[t]
    1
    aincp
    +
    μ[TP]
    3
    ′
    INSP
    [t]-
    INSP[t]
    aip
    +
    EP[t]
    (1-
    sspf[SP]
    )
    aincp
    -
    INSP[t]
    μ[INSP]
    4
    ′
    ISSP
    [t]-
    ISSP[t]
    aip
    +
    EP[t]
    sspf[SP]
    aincp
    -
    ISSP[t]
    μ[ISSP]
    5
    ′
    RP
    [t]
    INSP[t]
    +
    ISSP[t]
    aip
    -
    RP[t]
    μ[TP]
    6
    ′
    MLP
    [t]
    lpcr[ISSP,INSP]
    (-
    RP[t]
    -
    SP[t]
    +
    TP[t]
    )
    7
    TP[t]
    Max[0,
    EP[t]
    +
    INSP[t]
    +
    ISSP[t]
    +
    RP[t]
    +
    SP[t]
    ]
    ,RateRules
    #
    Symbol
    Value
    1
    μ[TP]
    1
    45625
    2
    μ[ISSP]
    0.035
    aip
    3
    μ[INSP]
    0.01
    aip
    4
    β[ISSP]
    0.15
    5
    β[INSP]
    0.15
    6
    aip
    26
    7
    aincp
    6
    8
    sspf[SP]
    0.2
    9
    lpcr[ISSP,INSP]
    1
    ,InitialConditions
    #
    Equation
    1
    SP[0]99998
    2
    EP[0]0
    3
    ISSP[0]1
    4
    INSP[0]1
    5
    RP[0]0
    6
    MLP[0]0
    7
    TP[0]100000
    
    Here is SEI2HR:
    Here are the “differences” between the two models:

    Equations explanations

    In this section we provide rationale of the model equations of SEI2HR.
    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.
    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. If there are enough hospital beds the severely symptomatic infected are hospitalized. 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 deaths from infection are tracked (accumulated.) 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:
    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 dictionary
    Here are the mortality rates from the SEI2HR dictionary
    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. The stock Total Population has the following (obvious) algebraic equation:
    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:
    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.
    • ◼
  • If the new ISSP people are more than the Available Hospital Beds (AHB) we take as many as AHB.
    • ◼
  • Hospitalized people have the same average infection period (aip).
    • ◼
  • Hospitalized (severely symptomatic) people have their own mortality rate.
    • 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.

    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.)

    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

    Remark: Other functions can be used, like, functions derived through some statistical fitting.
    Here is an example plot :
    To model quarantine with a piecewise constant function we use the following parameters:

    Parameters and actual simulation equations code

    Here are the parameters we want to experiment with (or do calibration with):
    Here we set custom rates and initial conditions:
    Here is the system of ODE’s we use to do parametrized simulations:

    Simulation

    Solutions

    Straightforward simulation for one year using ParametricNDSolve :

    Example evaluation

    Here are the parameters of a stock solution:
    Here we replace the parameters with concrete rate values (kept in the model object):
    Here is an example evaluation of a solution using the parameter values above:

    Interactive interface

    Using the interface in this section we can interactively see the effects of changing the focus parameters.

    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 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: The sensitivity analysis shown below should be done for other stocks and rates. In order to keep this exposition short we focus on ISSP, DIP, and HP.
    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:
    1
    .
    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:
    • ◼
  • More details for such equation addition are given in [AA2].
    • 2
    .
    Apply NIntegrate over stocks solution functions.
    3
    .
    Apply Trapezoidal rule to stock solution function values over a certain time grid.
    Below we use 1 and 3.

    Ranges

    Below we use the following sets of quarantine starts and quarantine durations.
    Note that putting the quarantine start to be at day 365 means “no quarantine.”

    Number of infected people

    Quarantine starts sensitivity

    Note that the plots and tabulated differences with “no quarantine” indicate that there is a very narrow range to choose an effective quarantine start.

    Quarantine duration sensitivity

    Number of deceased people

    Quarantine starts sensitivity

    Number of hospitalized people

    Quarantine starts sensitivity

    Infected Severely Symptomatic Population stock integral with respect to quarantine start and length

    In this section the 3D plot of AUC of ISSP is calculated using Trapezoidal rule.

    Deceased Infected Population stock with respect to quarantine start and length

    We can see from SEI2HR’s equations that DIP is already an AUC type of value. We can just plot the DIP values at the time horizon (one year.)

    Hospitalized Population stock integral with respect to quarantine start and length

    In this section the 3D plot of AUC of HP is calculated using Trapezoidal rule.

    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.

    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.

    Project management files

    [AAo1] Anton Antonov, WirVsVirus-Hackathon-work-plan.org, (2020), SystemsModeling at GitHub.
    [AAo2] Anton Antonov, WirVsVirus-hackathon-Geo-spatial-temporal-model-mind-map, (2020), SystemsModeling at GitHub.