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

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 with nodes that correspond to the sites and a positive matrix for the traveling patterns between the sites.

Problem definition: Given (i) a single site epidemiological model , (ii) a graph connecting multiple sites, and (iii) a traveling patterns matrix between the nodes of derive an epidemiological model that simulates more adequately viral decease propagation over .

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

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 , a simple graph with two nodes as , and both constant and time-dependent matrices for .

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 , which is a contingency matrix derived from source-destination traveling records. (Or the adjacency matrix of a travelling patterns graph.) The matrix entry of tells us how many people traveled from site to site . We systematically change the ODE’s of the sites in following way.

Assume that site had only travelers coming from site and going to site . Assume that the Total Population (TP) sizes for sites and are and 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

and change into the equation

assuming that

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

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

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 for the sites in the graph .

For each of node of make a copy of the model and denote with .

In general, the models have different initial conditions.

The models can also have different death rates, contact rates, etc.

Combine the models into the scaled model .

Change the equations of , to reflect the traveling patterns matrix .

Join the systems of ODE’s of , into one system of ODE’s.

Set appropriate or desired initial conditions for each of the populations in .

Solve the ODE’s of .

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.

From the plots above we see that both sites start with total populations of 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:

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.

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:

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:

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.

In the direction of the general “graph flow.“

The front of the infected populations wave is much steeper (gives “higher contrast”) than the tail.

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

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:

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

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

Using the parameters of the previous step.

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

In order to form “qualitative understanding.”

Get real life data.

Say, for the 2019-20 coronavirus outbreak.

Attempt manual or automatic calibration of the model.

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

Only manual calibration is shown in this notebook.

Load packages of the framework

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

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.

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

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:

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.

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