In this document we give usage examples for the functions of the package, “SystemDynamicsModelGraph.m”, [AAp1]. The package provides functions for making dependency graphs for the stocks in System Dynamics (SD) models. The primary motivation for creating the functions in this package is to have the ability to introspect, proofread, and verify the (typical) ODE models made in SD.
A more detailed explanation is:
For a given SD system S of Ordinary Differential Equations (ODEs) we make Mathematica graph objects that represent the interaction of variable dependent functions in S.
Those graph objects give alternative (and hopefully convenient) way of visualizing the model of S.
Load packages
The following commands load the packages [AAp1, AAp2, AAp3]:
When the second argument given to ModelDependencyGraph is Automatic the stocks in the equations are heuristically found with the function ModelHeuristicStocks:
Rapid specification of Machine Learning (ML) workflows using natural language commands.
The easiest things to automate with ML are ML workflows.
This presentation demonstrates that with natural language interfaces to ML algorithms.
Motivation
Assume that:
We want to create conversation agents that help Data Science (DS) and ML practitioners to quickly create first, initial versions of different DS and ML workflows for different programming languages and related packages.
We expect that the initial versions of programming code are tweaked further. (In order to produce desired outcomes in the application area of interest.)
The workflows considered
In this presentation we focus on these three ML areas:
Quantile Regression (QR)
Latent Semantic Analysis (LSA)
Recommendations
First example data
Assume we have a data frame with temperature data. Here is a summary:
summary(dfTemperatureData)
## Time Temperature
## Min. :3.629e+09 Min. : 4.72
## 1st Qu.:3.661e+09 1st Qu.:19.67
## Median :3.692e+09 Median :23.33
## Mean :3.692e+09 Mean :22.28
## 3rd Qu.:3.724e+09 3rd Qu.:26.11
## Max. :3.755e+09 Max. :31.17
ggplot2::ggplot(dfTemperatureData) + ggplot2::geom_point( ggplot2::aes( x = Time, y = Temperature ) )
Quantile Regression workflow: first example
Here is a Quantile Regression (QR) workflow specification:
qrmon2 <-
eval( expr =to_QRMon_R_command(
"create from dfTemperatureData; compute quantile regression with 12 knots and probabilities 0.25, 0.5, and 0.75; show date list plot with date origin 1900-01-01;", parse =TRUE) )
How it is done?
For a given ML domain (like QR or LSA) we create two types of Domain Specific Languages (DSL’s):
a software monad (i.e. programming language pipeline package) and
a DSL that is a subset of a spoken language.
These two DSL’s are combined: the latter translates natural language commands into the former.
By executing those translations we interpret commands of spoken DSL’s into ML computational results.
Note, that we assume that there is a separate system that converts speech into text.
Development cycle
Here is a clarification diagram:
MonadicMakingOfMLConversationalAgents
Grammars and parsers
For each type of workflow is developed a specialized DSL translation Raku module.
Each Raku module:
Has grammars for parsing a sequence of natural commands of a certain DSL
Translates the parsing results into corresponding software monad code
Different programming languages and packages can be the targets of the DSL translation.
(At this point are implemented DSL-translators to Python, R, and Wolfram Language.)
## Warning in QRMonSetData(res, data): The argument data is expected to be a data frame with columns: { Regressor, Value }.
## Warning in QRMonSetData(res, data): Proceeding by renaming the first columm "Time" as "Regressor" and renaming the second columm "Temperature" as "Value".
Latent Semantic Analysis workflow: translation
Here is how a sequence of natural commands that specifies a LSA workflow is translated into code for a LSA software monad:
lsaExpr <-
to_LSAMon_R_command(
"create from textHamlet; make document term matrix with automatic stop words and without stemming; apply lsi functions global weight function idf, local term weight function none, normalizer function cosine; extract 12 topics using method SVD, max steps 120, and min number of documents per term 2; show thesaurus table for ghost and grave;", parse =FALSE )
lsaExpr
lsamon2 <- eval(expr =parse( text =paste(lsaExpr)))
## Warning in NonNegativeMatrixFactorization::NearestWords(lsaObj$H, word, : More that one column name corresponds to the search word; the first match is used.
Consider the making of a recommender over the Titanic data:
dfTitanic[sample(1:nrow(dfTitanic), 12), ]
## id passengerClass passengerAge passengerSex passengerSurvival
## 1225 id.1225 3rd 20 male died
## 443 id.443 2nd 20 male died
## 761 id.761 3rd 30 male survived
## 835 id.835 3rd 30 male died
## 706 id.706 3rd -1 male died
## 339 id.339 2nd 30 male died
## 515 id.515 2nd 0 male survived
## 10 id.10 1st 70 male died
## 579 id.579 2nd 30 male died
## 1248 id.1248 3rd -1 female survived
## 673 id.673 3rd -1 male died
## 1205 id.1205 3rd 20 male died
## id passengerClass passengerAge passengerSex passengerSurvival
## id.1 : 1 1st:323 Min. :-1.00 female:466 died :809
## id.10 : 1 2nd:277 1st Qu.:10.00 male :843 survived:500
## id.100 : 1 3rd:709 Median :20.00
## id.1000: 1 Mean :23.55
## id.1001: 1 3rd Qu.:40.00
## id.1002: 1 Max. :80.00
## (Other):1303
Recommender workflow
Here are recommender workflow specification and evaluation results:
smrmon2 <-
eval( expr =to_SMRMon_R_command(
"create from dfTitanic; apply the LSI functions inverse document frequency, term frequency, and cosine; compute the top 6 recommendations for the profile female=1, 30=1; extend recommendations with dfTitanic; show pipeline value", parse =TRUE ) )
The approach taken in the design and implementation of the natural language commands interpreters can handle misspellings:
smrmon2 <-
eval( expr =to_SMRMon_R_command(
"create from dfTitanic; aply the LSI functions inverse document frequency, term frequency, and cosine; compute the top 6 recomendations for the profle female=1, 30=1; extend recommendations with dfTitanic; show pipeline value" ) )
## [1] "Possible misspelling of 'apply' as 'aply'."
## [2] "Possible misspelling of 'recommendations' as 'recomendations'."
## [3] "Possible misspelling of 'profile' as 'profle'."
## Score Index id passengerClass passengerAge passengerSex passengerSurvival
## 1 2 1 id.1 1st 30 female survived
## 2 2 33 id.1027 3rd 30 female survived
## 3 2 68 id.1059 3rd 30 female died
## 4 2 72 id.1062 3rd 30 female survived
## 5 2 99 id.1087 3rd 30 female died
## 6 2 108 id.1095 3rd 30 female survived
Not just ML workflows
Obviously this approach can be used for any type of computational workflows.
Here is an example of an Epidemiology Modeling workflow:
ecmCommands <-
'create with the model susceptible exposed infected two hospitalized recovered; assign 100000 to the susceptible population; set infected normally symptomatic population to be 0; set infected severely symptomatic population to be 1; assign 0.56 to contact rate of infected normally symptomatic population; assign 0.58 to contact rate of infected severely symptomatic population; assign 0.1 to contact rate of the hospitalized population; simulate for 240 days; plot populations results;'
(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.)
About This Data The CSV file and charts on this site show a relative volume of directions requests per country/region or city compared to a baseline volume on January 13th, 2020. We define our day as midnight-to-midnight, Pacific time. Cities represent usage in greater metropolitan areas and are stably defined during this period. In many countries/regions and cities, relative volume has increased since January 13th, consistent with normal, seasonal usage of Apple Maps. Day of week effects are important to normalize as you use this data. Data that is sent from users’ devices to the Maps service is associated with random, rotating identifiers so Apple doesn’t have a profile of your movements and searches. Apple Maps has no demographic information about our users, so we can’t make any statements about the representativeness of our usage against the overall population.
Observations
The observations listed in this subsection are also placed under the relevant statistics in the following sections and indicated with “Observation”.
The directions request volumes reference date for normalization is 2020-01-13 : all the values in that column are 100.
From the community clusters of the nearest neighbor graphs (derived from the time series of the normalized driving directions requests volume) we see that countries and cities are clustered in expected ways. For example, in the community graph plot corresponding to “{city, driving}” the cities Oslo, Copenhagen, Helsinki, Stockholm, and Zurich are placed in the same cluster. In the graphs corresponding to “{city, transit}” and “{city, walking}” the Japanese cities Tokyo, Osaka, Nagoya, and Fukuoka are clustered together.
In the time series plots the Sundays are indicated with orange dashed lines. We can see that from Monday to Thursday people are more familiar with their trips than say on Fridays and Saturdays. We can also see that on Sundays people (on average) are more familiar with their trips or simply travel less.
Apple mobile data was provided in this WWW page: https://www.apple.com/covid19/mobility , [APPL1]. (The data has to be download from that web page – there is an “agreement to terms”, etc.)
We can visualize the data using heat-map plots. Here we use the package “HeatmapPlot.m”, [AAp2].
Remark: Using the contingency matrices prepared for the heat-map plots we can do further analysis, like, finding correlations or nearest neighbors. (See below.)
Observation: From the community clusters of the nearest neighbor graphs (derived from the time series of the normalized driving directions requests volume) we see that countries and cities are clustered in expected ways. For example in the community graph plot corresponding to “{city, driving}” the cities Oslo, Copenhagen, Helsinki, Stockholm, and Zurich are placed in the same cluster. In the graphs corresponding to “{city, transit}” and “{city, walking}” the Japanese cities Tokyo, Osaka, Nagoya, and Fukuoka are clustered together.
Time series analysis
Time series
In this section for each date we sum all cases over the region-transportation pairs, make a time series, and plot them.
Remark: In the plots the Sundays are indicated with orange dashed lines.
Observation: In the time series plots the Sundays are indicated with orange dashed lines. We can see that from Monday to Thursday people are more familiar with their trips than say on Fridays and Saturdays. We can also see that on Sundays people (on average) are more familiar with their trips or simply travel less.
“Forecast”
He we do “forecast” for code-workflow demonstration purposes – the forecasts should not be taken seriously.
(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.)
“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].
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:
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
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:
0z7878rbvh3ut
Note that with Max we specified that the total population cannot be less than .
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.
1gercwr5szizt
Because we hospitalize the severely infected people only instead of the term
09xmdu8xq7e0p
we have the terms
0dkc6td55qv7p
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.
083x1htqp0kkx
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.
01cfzpu1jv3o4
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.
0zhs2d6lt6gz8
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:
1c35mhvi92ffe
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.
04y15eqwki71f
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.
1ikrv4a51h9xo
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.)
1is45mamhsayd
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, .
The MS producer has limited capacity for delivery, .
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:
1uxby6ipa49dw
Here is the corresponding HMS equation:
0ri4ohdddjxle
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, .
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:
1dy38wd13b0h6
Here is the MS formula with the MS delivery term replaced with “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, :
(*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:
0xc7vujwb2npw
Medical Supplies Demand
The stock Medical Supplies Demand (MSD) simply accumulates the MS demand derived from population stocks and their corresponding mscr:
0ut8smszrobx2
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.
0evzoim1tjwsj
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.
1fhewr6izy46n
Quarantine scenarios
In order to model quarantine scenarios we use piecewise constant functions for the contact rates and .
Remark: Other functions can be used, like, functions derived through some statistical fitting.
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].)
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:
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:
[HH1] Herbert W. Hethcote (2000). “The Mathematics of Infectious Diseases”. SIAM Review. 42 (4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.
“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:
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
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:
1cg8j8gglqwi1
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 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
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:
0peb290e1lx2d
Note that with Max we specified that the total population cannot be less than .
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.
0h1zwgod2xs4g
Because we hospitalize the severely infected people only instead of the term
09xmdu8xq7e0p
we have the terms
0dkc6td55qv7p
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.
0mrj0gt13ppsn
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.
0hy96yb7pg654
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.
19hml5ravi0a8
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.
0ippnyp2hf90k
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.
042iy3oy8o454
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.
1glekhgj1v0e7
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.)
1is45mamhsayd
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.
0evzoim1tjwsj
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.
1133o7bm975sr
Quarantine scenarios
In order to model quarantine scenarios we use piecewise constant functions for the contact rates and .
Remark: Other functions can be used, like, functions derived through some statistical fitting.
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:
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.
Remark: The AUC measure for a stock is indicated with the prefix “∫”. For example AUC for ISSP is marked with “∫ISSP”.
Ranges
Below we use the following sets of quarantine starts and quarantine durations.
[HH1] Herbert W. Hethcote (2000). “The Mathematics of Infectious Diseases”. SIAM Review. 42 (4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.
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.
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.)
Yesterday in one of the forums I frequent it was announced that New York Times has published COVID-19 data on GitHub. I decided to make a Mathematica notebook that gives data links and related code for data ingestions. (And rudimentary data analysis.)
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.)
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:
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:
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.)
Mathematica for prediction algorithms 
