Scaling of Epidemiology Models with Multi-site Compartments

Accendo Data LLC

Introduction

In this notebook we describe and exemplify an algorithm that allows the specification and execution of 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 nonuniformly populated (for example, China, USA and Russia). There is a need to derive epidemiological models that take into account the nonuniform 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, 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 (ODEs) 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

tpm(G)

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

tpm(G)

between the nodes of

, derive an epidemiological model

S(M,tpm(G))

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

In[]:=

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/SystemDynamicsInteractiveInterfacesFunctions.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

, a simple graph with two nodes as

and both constant and time-dependent matrices for

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 subsection “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 traveling 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 traveling 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 than 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) to 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 ODEs 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

K(i,j)

tells us how many people traveled from site

to site

. We systematically change the ODEs 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

′

(t)-

(t)

(t)μ,

(

)

and change into the equation

′

(t)-

(t)

(t)μ-

K(a,b)

(t)

K(b,a)

(t)

(

)

assuming that

K(a,b)

(t)

≤

K(b,a)

(t)

≤

(

)

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

min

K(i,j)

(t)

(t).

(

)

The transformed systems of ODEs of the sites are joined into one “big” system of ODEs, 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:

CompoundExpression[

]

Out[]=

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

, make a copy of the model

and denote with

M[i]

In general, the models

M[i],i∈G

have different initial conditions.

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

Combine the models

M[i],i∈G

into the scaled model

Change the equations of

M[i]

i∈G

to reflect the traveling patterns matrix

Join the systems of ODEs of

M[i]

i∈G

into one system of ODEs.

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

Solve the ODEs of

Visualize the results.

Precaution

Care should be taken when specifying the initial conditions of MSEMEA’s system of ODEs (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 in some safeguards, 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]:

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

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

Below we use IDs 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:

And here is the traveling patterns matrix:

Note that there are many 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):

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

Show the unique stocks in the multi-site model:

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

Show the initial conditions:

Show the total number of equations:

Solve the system of ODEs of the extended model:

Display the solutions for each site separately:

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

Here are the corresponding number of traveling people functions:

Here we scale the SIR model, solve the obtained system of ODEs and plot the solutions:

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:

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:

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

Change the initial conditions in the following ways:

◼

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.

Solve the system of ODEs of the scaled model:

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

Display solutions of the first and last site:

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 the 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 subsequence for the total populations:

Here is a subsequence for the sum of the infected populations:

Here is a subsequence for the recovered population:

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

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:

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 a 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:

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.

The wave propagates 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.

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:

Here is a fragment of the matrix:

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

Change the initial conditions in the following ways:

◼

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.

Here solve the obtained system of ODEs:

Here we plot the solutions:

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:

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

Money from Lost Productivity

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

Here we plot the sum of the accumulated money losses:

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

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

Do real-data simulations related to the spread of COVID-19.

References

Articles and Books

[Wk1] Wikipedia, "Compartmental Models in Epidemiology".

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

[AA1] Antonov, Anton. (2020). "Coronavirus Propagation Modeling Considerations." SystemModeling at GitHub.

[AA2] Antonov, Anton. (2020). "Basic Experiments Workflow for Simple Epidemiological Models." SystemModeling at GitHub.

[AA3] Antonov, Anton. (2006). "Air Pollution Modeling with gridMathematica." Wolfram Technology Conference.

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

Repositories and Packages

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

[AAr1] Antonov, Anton. (2020). “Coronavirus Propagation Dynamics Project.” SystemModeling at GitHub.

[AAp1] Antonov, Anton. (2020). "Epidemiology Models Mathematica Package." SystemModeling at GitHub.

[AAp2] Antonov, Anton. (2020). "Epidemiology Models Modifications Mathematica Package." SystemsModeling at GitHub.

[AAp3] Antonov, Anton. (2020). "Epidemiology Modeling Visualization Functions Mathematica Package." SystemModeling at GitHub.

[AAp4] Antonov, Anton. (2020). "System Dynamics Interactive Interfaces Functions Mathematica Package." SystemModeling at GitHub.