WolframAlpha.com
WolframCloud.com
All Sites & Public Resources...
Products & Services
Wolfram|One
Mathematica
Wolfram|Alpha Notebook Edition
Finance Platform
System Modeler
Wolfram Player
Wolfram Engine
WolframScript
Enterprise Private Cloud
Application Server
Enterprise Mathematica
Wolfram|Alpha Appliance
Enterprise Solutions
Corporate Consulting
Technical Consulting
Wolfram|Alpha Business Solutions
Resource System
Data Repository
Neural Net Repository
Function Repository
Wolfram|Alpha
Wolfram|Alpha Pro
Problem Generator
API
Data Drop
Products for Education
Mobile Apps
Wolfram Player
Wolfram Cloud App
Wolfram|Alpha for Mobile
Wolfram|Alpha-Powered Apps
Services
Paid Project Support
Wolfram U
Summer Programs
All Products & Services »
Technologies
Wolfram Language
Revolutionary knowledge-based programming language.
Wolfram Cloud
Central infrastructure for Wolfram's cloud products & services.
Wolfram Science
Technology-enabling science of the computational universe.
Wolfram Notebooks
The preeminent environment for any technical workflows.
Wolfram Engine
Software engine implementing the Wolfram Language.
Wolfram Natural Language Understanding System
Knowledge-based broadly deployed natural language.
Wolfram Data Framework
Semantic framework for real-world data.
Wolfram Universal Deployment System
Instant deployment across cloud, desktop, mobile, and more.
Wolfram Knowledgebase
Curated computable knowledge powering Wolfram|Alpha.
All Technologies »
Solutions
Engineering, R&D
Aerospace & Defense
Chemical Engineering
Control Systems
Electrical Engineering
Image Processing
Industrial Engineering
Mechanical Engineering
Operations Research
More...
Finance, Statistics & Business Analysis
Actuarial Sciences
Bioinformatics
Data Science
Econometrics
Financial Risk Management
Statistics
More...
Education
All Solutions for Education
Trends
Machine Learning
Multiparadigm Data Science
Internet of Things
High-Performance Computing
Hackathons
Software & Web
Software Development
Authoring & Publishing
Interface Development
Web Development
Sciences
Astronomy
Biology
Chemistry
More...
All Solutions »
Learning & Support
Learning
Wolfram Language Documentation
Fast Introduction for Programmers
Wolfram U
Videos & Screencasts
Wolfram Language Introductory Book
Webinars & Training
Summer Programs
Books
Need Help?
Support FAQ
Wolfram Community
Contact Support
Premium Support
Paid Project Support
Technical Consulting
All Learning & Support »
Company
About
Company Background
Wolfram Blog
Events
Contact Us
Work with Us
Careers at Wolfram
Internships
Other Wolfram Language Jobs
Initiatives
Wolfram Foundation
MathWorld
Computer-Based Math
A New Kind of Science
Wolfram Technology for Hackathons
Student Ambassador Program
Wolfram for Startups
Demonstrations Project
Wolfram Innovator Awards
Wolfram + Raspberry Pi
Summer Programs
More...
All Company »
Search
WOLFRAM COMMUNITY
Connect with users of Wolfram technologies to learn, solve problems and share ideas
Join
Sign In
Dashboard
Groups
People
Message Boards
Answer
(
Unmark
)
Mark as an Answer
GROUPS:
Staff Picks
Biological Sciences
Medical Sciences
Software Development
Graphics and Visualization
Wolfram Language
Modeling
Numerical Computation
13
Robert Nachbar
Stochastic Epidemiology Models with Applications to the COVID-19
Robert Nachbar, Wolfram Research Inc.
Posted
2 years ago
7448 Views
|
11 Replies
|
13 Total Likes
Follow this post
|
MODERATOR NOTE: coronavirus resources & updates:
https://wolfr.am/coronavirus
Stochastic Epidemiology Models with Applications to the COVID-19 Pandemic
As a follow-on to my earlier post about
Epidemiological Models for Influenza and COVID-19
, I’d like to show how these compartmental models can be used in a stochastic setting to study how the inherent randomness of events can affect the interpretation of modeled outcomes.
As noted before, the deterministic ODE-based models have a number of assumptions:
◼
All of the individuals of a given type are identical
◼
The population of individuals is well mixed, that is, instantaneously and uniformly mixed
◼
The number of individuals is so large that the density can be treated as a continuous variable
The last assumption rarely holds in practice for epidemiology models where the number of individuals is typically on the order of hundreds to millions. These population sizes are vastly smaller that those usually encountered in chemical kinetics, where the numbers of molecules in a beaker are on the order of
2
3
1
0
(
Avagodro's number
) or more.
So, what’s to be done? In the 1970’s, Daniel Gillespie [
1
] developed methods based on the chemical master equation to simulate chemical reactions stochastically where the last assumption, especially, does not hold. The simulation method correctly accounts for the inherent fluctuations and correlations that are ignored in the deterministic formulation. There is a nice review of the method with implementation by Higham [
2
], and a good book by Wilkinson with applications to systems biology [
3
]. Finally, we have used a hybrid deterministic-stochastic approach for viral dynamics modeling [
4
,
5
].
One of the advantages of the stochastic modeling approach is the access to distributions of individuals in the various compartments over time. Let’s look at the SEIR model.
It should be pointed out that these stochastic model simulations are not the same as agent-based models (see
Christopher
Wolfram's
recent
post
for an example). The main difference is that the latter use a uniform time step and all agents update at each step, while the former use varying length time steps and only one event occurs at each step.
Sonata—SEIR
We begin as usual by loading the
C
o
m
p
a
r
t
m
e
n
t
a
l
M
o
d
e
l
i
n
g
`
package (download
here
), defining the transitions of the model, and computing the kinetic elements:
I
n
[
]
:
=
<
<
C
o
m
p
a
r
t
m
e
n
t
a
l
M
o
d
e
l
i
n
g
`
I
n
[
]
:
=
m
o
d
e
l
N
a
m
e
=
"
S
E
I
R
"
;
m
o
d
e
l
=
β
λ
[
t
]
→
ℰ
,
ℰ
ζ
→
ℐ
,
ℐ
γ
→
ℛ
;
f
o
r
c
e
O
f
I
n
f
e
c
t
i
o
n
=
{
λ
[
t
_
]
ℐ
[
t
]
}
;
I
n
[
]
:
=
m
o
d
e
l
D
a
t
a
=
K
i
n
e
t
i
c
C
o
m
p
a
r
t
m
e
n
t
a
l
M
o
d
e
l
[
m
o
d
e
l
/
.
f
o
r
c
e
O
f
I
n
f
e
c
t
i
o
n
,
t
]
;
v
a
r
s
=
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
;
(
o
d
e
s
=
m
o
d
e
l
D
a
t
a
[
"
E
q
u
a
t
i
o
n
s
"
]
)
/
/
C
o
l
u
m
n
O
u
t
[
]
=
′
ℰ
[
t
]
-
ζ
ℰ
[
t
]
+
β
ℐ
[
t
]
[
t
]
′
ℐ
[
t
]
ζ
ℰ
[
t
]
-
γ
ℐ
[
t
]
′
ℛ
[
t
]
γ
ℐ
[
t
]
′
[
t
]
-
β
ℐ
[
t
]
[
t
]
W
e
’
l
l
p
a
r
a
m
e
t
e
r
i
z
e
t
h
e
m
o
d
e
l
i
n
t
e
r
m
s
o
f
t
h
e
b
a
s
i
c
r
e
p
r
o
d
u
c
t
i
o
n
n
u
m
b
e
r
ℜ
0
,
a
n
d
u
s
e
a
n
i
n
c
u
b
a
t
i
o
n
p
e
r
i
o
d
(
1
ζ
)
a
n
d
d
u
r
a
t
i
o
n
o
f
i
n
f
e
c
t
i
o
u
s
n
e
s
s
(
1
γ
)
c
h
a
r
a
c
t
e
r
i
s
t
i
c
o
f
C
O
V
I
D
-
1
9
.
I
n
[
]
:
=
b
e
t
a
S
o
l
n
=
F
i
r
s
t
@
S
o
l
v
e
ℜ
0
β
γ
,
β
O
u
t
[
]
=
β
γ
ℜ
0
I
n
[
]
:
=
p
a
r
a
m
s
=
<
|
ℜ
0
5
.
,
β
β
,
ζ
1
5
.
1
,
γ
1
7
.
,
3
1
0
,
I
0
1
|
>
/
.
b
e
t
a
S
o
l
n
p
a
r
a
m
s
N
=
N
o
r
m
a
l
@
(
p
a
r
a
m
s
/
/
.
p
a
r
a
m
s
)
O
u
t
[
]
=
ℜ
0
5
.
,
β
γ
ℜ
0
,
ζ
0
.
1
9
6
0
7
8
,
γ
0
.
1
4
2
8
5
7
,
1
0
0
0
,
I
0
1
O
u
t
[
]
=
{
ℜ
0
5
.
,
β
0
.
0
0
0
7
1
4
2
8
6
,
ζ
0
.
1
9
6
0
7
8
,
γ
0
.
1
4
2
8
5
7
,
1
0
0
0
,
I
0
1
}
I
n
[
]
:
=
i
n
i
t
s
=
J
o
i
n
[
{
[
0
]
-
I
0
,
ℐ
[
0
]
I
0
}
,
T
h
r
e
a
d
[
T
h
r
o
u
g
h
@
C
o
m
p
l
e
m
e
n
t
[
v
a
r
s
,
{
,
ℐ
}
]
[
0
]
0
]
]
O
u
t
[
]
=
{
[
0
]
-
I
0
+
,
ℐ
[
0
]
I
0
,
ℰ
[
0
]
0
,
ℛ
[
0
]
0
}
First, let’s see how the usual ODE solution behaves over 120 days (the code for
e
p
i
d
e
m
i
c
M
o
d
e
l
P
l
o
t
is in the initialization section at the end of the notebook).
I
n
[
]
:
=
t
m
a
x
=
1
2
0
;
I
n
[
]
:
=
s
o
l
n
O
D
E
=
F
i
r
s
t
@
N
D
S
o
l
v
e
[
{
o
d
e
s
,
i
n
i
t
s
}
/
.
p
a
r
a
m
s
N
,
v
a
r
s
,
{
t
,
0
,
t
m
a
x
}
]
;
I
n
[
]
:
=
e
p
i
d
e
m
i
c
M
o
d
e
l
P
l
o
t
[
s
o
l
n
O
D
E
,
{
0
,
t
m
a
x
}
,
"
m
o
d
e
l
N
a
m
e
"
m
o
d
e
l
N
a
m
e
,
I
m
a
g
e
S
i
z
e
L
a
r
g
e
]
O
u
t
[
]
=
ℰ
ℐ
ℛ
Nothing unusual. There is the outbreak that reaches a maximum around day 40, and nearly everyone in the population gets infected.
Now, let’s see what the stochastic simulation shows. We’ll use the same parameters and the same initial conditions, and run 500 (typically one would use more runs, but we'll keep it on the small side for this essay; we'll also use SeedRandom from time to time) simulations (the code for
s
t
o
c
h
a
s
t
i
c
S
o
l
v
e
is in the initialization section at the end of the notebook).
I
n
[
]
:
=
S
e
e
d
R
a
n
d
o
m
[
1
4
2
5
3
6
7
]
s
o
l
n
s
C
M
E
=
T
a
b
l
e
[
s
t
o
c
h
a
s
t
i
c
S
o
l
v
e
[
m
o
d
e
l
D
a
t
a
[
"
R
a
t
e
s
"
]
/
.
p
a
r
a
m
s
N
,
m
o
d
e
l
D
a
t
a
[
"
S
t
o
i
c
h
i
o
m
e
t
r
i
c
M
a
t
r
i
x
"
]
,
T
h
r
o
u
g
h
@
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
[
0
]
/
.
R
u
l
e
@
@
@
i
n
i
t
s
/
.
p
a
r
a
m
s
N
,
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
,
{
t
,
0
,
t
m
a
x
}
]
,
5
0
0
]
;
/
/
T
i
m
i
n
g
O
u
t
[
]
=
{
3
5
.
5
7
8
,
N
u
l
l
}
Here is a plot of the first 20 runs:
I
n
[
]
:
=
e
p
i
d
e
m
i
c
M
o
d
e
l
P
l
o
t
[
J
o
i
n
@
@
T
a
k
e
[
s
o
l
n
s
C
M
E
,
2
0
]
,
{
0
,
t
m
a
x
}
,
"
m
o
d
e
l
N
a
m
e
"
m
o
d
e
l
N
a
m
e
,
I
m
a
g
e
S
i
z
e
L
a
r
g
e
]
O
u
t
[
]
=
ℰ
ℐ
ℛ
We can see that there is quite a bit of variability in the individual runs, and that there is about a three week variation for the occurrence of the outbreak maximum with possibly a hint of bimodal shape. Moreover, some of the runs did not produce an outbreak. We can examine the average behavior by plotting the median and the middle 95 percentile band:
I
n
[
]
:
=
T
a
b
V
i
e
w
[
{
"
m
e
d
i
a
n
w
i
t
h
b
a
n
d
s
"
e
p
i
d
e
m
i
c
M
o
d
e
l
P
l
o
t
[
s
o
l
n
s
C
M
E
,
{
0
,
t
m
a
x
}
,
"
m
o
d
e
l
N
a
m
e
"
"
S
I
R
"
,
I
m
a
g
e
S
i
z
e
L
a
r
g
e
]
,
"
m
e
d
i
a
n
"
e
p
i
d
e
m
i
c
M
o
d
e
l
P
l
o
t
[
s
o
l
n
s
C
M
E
,
{
0
,
t
m
a
x
}
,
"
m
o
d
e
l
N
a
m
e
"
m
o
d
e
l
N
a
m
e
,
I
m
a
g
e
S
i
z
e
L
a
r
g
e
,
"
s
h
o
w
B
a
n
d
s
"
F
a
l
s
e
]
}
]
O
u
t
[
]
=
m
e
d
i
a
n
w
i
t
h
b
a
n
d
s
m
e
d
i
a
n
ℰ
ℐ
ℛ
B
a
n
d
s
r
e
p
r
e
s
e
n
t
0
.
0
2
5
a
n
d
0
.
9
7
5
q
u
a
n
t
i
l
e
s
The middle 95 percentile bands show that there is quite a bit of variability in the ensemble of solutions. Looking at the distributions of susceptible
and recovered
ℛ
at
t
t
m
a
x
we see that there are essentially two kinds of behavior:
I
n
[
]
:
=
H
i
s
t
o
g
r
a
m
[
Q
u
i
e
t
[
#
[
t
m
a
x
]
/
.
s
o
l
n
s
C
M
E
]
,
{
1
0
}
,
C
h
a
r
t
S
t
y
l
e
e
p
i
C
o
l
o
r
[
#
]
,
P
l
o
t
L
a
b
e
l
#
,
P
l
o
t
R
a
n
g
e
4
0
0
,
I
m
a
g
e
S
i
z
e
4
7
2
]
&
/
@
{
,
ℛ
}
/
/
R
o
w
[
#
,
S
p
a
c
e
r
[
1
8
]
]
&
O
u
t
[
]
=
We can partition the ensemble into an “outbreak” and “not an outbreak” subsets:
I
n
[
]
:
=
{
o
u
t
b
r
e
a
k
,
n
o
t
O
u
t
b
r
e
a
k
}
=
Q
u
i
e
t
[
S
o
r
t
B
y
[
s
o
l
n
s
C
M
E
,
(
ℛ
[
t
m
a
x
]
/
.
#
)
<
5
0
&
]
/
/
S
p
l
i
t
B
y
[
#
,
(
ℛ
[
1
5
0
]
/
.
#
)
<
5
0
&
]
&
]
;
L
e
n
g
t
h
/
@
%
O
u
t
[
]
=
{
3
9
4
,
1
0
6
}
About a fifth of the time the infection dies out before an outbreak can be established. Now, looking at the average behavior of the outbreak group, while the average behavior is very similar to that for the deterministic solution, the substantial variability remains.
I
n
[
]
:
=
M
a
p
T
h
r
e
a
d
[
#
2
e
p
i
d
e
m
i
c
M
o
d
e
l
P
l
o
t
[
#
1
,
{
0
,
t
m
a
x
}
,
"
m
o
d
e
l
N
a
m
e
"
m
o
d
e
l
N
a
m
e
,
I
m
a
g
e
S
i
z
e
L
a
r
g
e
]
&
,
{
{
o
u
t
b
r
e
a
k
,
n
o
t
O
u
t
b
r
e
a
k
,
s
o
l
n
O
D
E
}
,
{
"
o
u
t
b
r
e
a
k
"
,
"
n
o
t
o
u
t
b
r
e
a
k
"
,
"
d
e
t
e
r
m
i
n
i
s
t
i
c
"
}
}
]
/
/
T
a
b
V
i
e
w
O
u
t
[
]
=
o
u
t
b
r
e
a
k
n
o
t
o
u
t
b
r
e
a
k
d
e
t
e
r
m
i
n
i
s
t
i
c
ℰ
ℐ
ℛ
B
a
n
d
s
r
e
p
r
e
s
e
n
t
0
.
0
2
5
a
n
d
0
.
9
7
5
q
u
a
n
t
i
l
e
s
So, why are there two different regimes of behavior? Now would be a good time to examine the Gillespie Algorithm.
Adagio—Details of the Gillespie Algorithm
Elements of a stochastic simulation
The simulation is a Markov chain process, which is a random process that depends only on the current state to determine the next state. The events that occur at step of the chain are the transitions of the SEIR model, in this case. The likelihood of a particular event occurring depends on its propensity, which in turn depends on the current state, that is, the numbers of individuals in each compartment. The time between events depends on the total propensity of the system in the current state. The propensity changes from step to step, and the process ends when the maximum time of the simulation is reached, or when the total propensity is 0.
The information needed to step from the current state to the next state can be found in the
m
o
d
e
l
D
a
t
a
Association, and is shown in the table below. The key
"
C
o
m
p
o
n
e
n
t
T
r
a
n
s
i
t
i
o
n
s
"
gives the data in the first column, and
"
R
a
t
e
s
"
gives that in the last column. The key
"
S
t
o
i
c
h
i
o
m
e
t
r
i
c
M
a
t
r
i
x
"
gives the table in the middles column, whose headings are obtained from the
"
V
a
r
i
a
b
l
e
s
"
key.
I
n
[
]
:
=
M
a
p
T
h
r
e
a
d
[
J
o
i
n
,
{
L
i
s
t
/
@
m
o
d
e
l
D
a
t
a
[
"
C
o
m
p
o
n
e
n
t
T
r
a
n
s
i
t
i
o
n
s
"
]
,
R
o
t
a
t
e
R
i
g
h
t
/
@
m
o
d
e
l
D
a
t
a
[
"
S
t
o
i
c
h
i
o
m
e
t
r
i
c
M
a
t
r
i
x
"
]
,
L
i
s
t
/
@
m
o
d
e
l
D
a
t
a
[
"
R
a
t
e
s
"
]
}
]
/
/
f
o
r
m
a
t
t
a
b
l
e
O
u
t
[
]
=
S
t
o
i
c
h
i
o
m
e
t
r
y
T
r
a
n
s
i
t
i
o
n
ℰ
ℐ
ℛ
P
r
o
p
e
n
s
i
t
y
β
ℐ
[
t
]
→
ℰ
-
1
1
0
0
β
ℐ
[
t
]
[
t
]
ℰ
ζ
→
ℐ
0
-
1
1
0
ζ
ℰ
[
t
]
ℐ
γ
→
ℛ
0
0
-
1
1
γ
ℐ
[
t
]
Computing the propensity
To demonstrate the steps, we begin by assuming a current state and time:
I
n
[
]
:
=
s
t
a
t
e
=
R
a
n
d
o
m
I
n
t
e
g
e
r
[
{
0
,
1
0
0
0
}
,
L
e
n
g
t
h
@
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
]
t
i
m
e
=
R
a
n
d
o
m
R
e
a
l
[
{
0
,
5
0
}
]
O
u
t
[
]
=
{
5
5
1
,
5
5
4
,
5
6
1
,
9
4
0
}
O
u
t
[
]
=
4
4
.
5
8
6
6
The propensity vector
is the model rate vector evaluated at the current state and time:
I
n
[
]
:
=
=
m
o
d
e
l
D
a
t
a
[
"
R
a
t
e
s
"
]
/
.
T
h
r
e
a
d
[
T
h
r
o
u
g
h
[
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
[
t
]
]
s
t
a
t
e
]
/
.
p
a
r
a
m
s
N
/
.
t
t
i
m
e
O
u
t
[
]
=
{
1
0
8
.
0
3
9
,
7
9
.
1
4
2
9
,
3
7
1
.
9
7
1
}
The individual propensities define an empirical distribution:
I
n
[
]
:
=
e
v
e
n
t
=
E
m
p
i
r
i
c
a
l
D
i
s
t
r
i
b
u
t
i
o
n
[
(
R
a
n
g
e
[
L
e
n
g
t
h
@
]
)
]
O
u
t
[
]
=
D
a
t
a
D
i
s
t
r
i
b
u
t
i
o
n
T
y
p
e
:
E
m
p
i
r
i
c
a
l
D
a
t
a
p
o
i
n
t
s
:
3
Graphically, we can see how the propensities on the left contribute to the cumulative distribution on the right:
I
n
[
]
:
=
R
o
w
B
a
r
C
h
a
r
t
,
c
h
a
r
t
o
p
t
i
o
n
s
,
D
i
s
c
r
e
t
e
P
l
o
t
C
D
F
[
e
v
e
n
t
,
x
+
1
]
,
{
x
,
0
,
L
e
n
g
t
h
[
]
,
0
.
0
1
}
,
p
l
o
t
o
p
t
i
o
n
s
,
S
p
a
c
e
r
[
6
]
O
u
t
[
]
=
We can use
R
a
n
d
o
m
V
a
r
i
a
t
e
raw a random sample from the distribution to select the transition. Here is a list of 10 samples:
I
n
[
]
:
=
T
a
b
l
e
[
R
a
n
d
o
m
V
a
r
i
a
t
e
[
e
v
e
n
t
]
,
1
0
]
O
u
t
[
]
=
{
1
,
3
,
3
,
3
,
3
,
3
,
3
,
1
,
2
,
3
}
The next state is computed by incrementing the current state with stoichiometry of the selected transition:
I
n
[
]
:
=
o
l
d
S
t
a
t
e
=
s
t
a
t
e
;
s
t
a
t
e
+
=
m
o
d
e
l
D
a
t
a
[
"
S
t
o
i
c
h
i
o
m
e
t
r
i
c
M
a
t
r
i
x
"
]
〚
E
c
h
o
[
#
,
"
"
,
m
o
d
e
l
D
a
t
a
[
"
C
o
m
p
o
n
e
n
t
T
r
a
n
s
i
t
i
o
n
s
"
]
〚
#
〛
&
]
&
@
R
a
n
d
o
m
V
a
r
i
a
t
e
[
e
v
e
n
t
]
〛
;
T
a
b
l
e
F
o
r
m
[
{
o
l
d
S
t
a
t
e
,
s
t
a
t
e
}
/
/
T
h
r
e
a
d
,
T
a
b
l
e
H
e
a
d
i
n
g
s
{
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
,
{
"
b
e
f
o
r
e
"
,
"
a
f
t
e
r
"
}
}
]
»
ℐ
γ
→
ℛ
O
u
t
[
]
/
/
T
a
b
l
e
F
o
r
m
=
b
e
f
o
r
e
a
f
t
e
r
ℰ
5
5
1
5
5
1
ℐ
5
5
4
5
5
3
ℛ
5
6
1
5
6
2
9
4
0
9
4
0
The events occur randomly in time, and are therefore a
Poisson process
. It is the same process followed by radioactive decay, and follows an exponential distribution:
I
n
[
]
:
=
Δ
t
i
m
e
=
E
x
p
o
n
e
n
t
i
a
l
D
i
s
t
r
i
b
u
t
i
o
n
[
T
o
t
a
l
[
]
]
O
u
t
[
]
=
E
x
p
o
n
e
n
t
i
a
l
D
i
s
t
r
i
b
u
t
i
o
n
[
5
5
9
.
1
5
4
]
I
n
[
]
:
=
W
i
t
h
t
m
a
x
=
L
o
g
[
2
]
T
o
t
a
l
[
]
6
,
R
o
w
P
l
o
t
P
D
F
[
Δ
t
i
m
e
,
x
]
,
{
x
,
0
,
t
m
a
x
}
,
P
D
F
o
p
t
i
o
n
s
,
P
l
o
t
C
D
F
[
Δ
t
i
m
e
,
x
]
,
{
x
,
0
,
t
m
a
x
}
,
I
m
a
g
e
S
i
z
e
4
7
2
,
C
D
F
o
p
t
i
o
n
s
,
S
p
a
c
e
r
[
6
]
O
u
t
[
]
=
Again,
R
a
n
d
o
m
V
a
r
i
a
t
e
can be used to draw a random sample from the exponential distribution to select the time to the next event, and increment the time:
I
n
[
]
:
=
o
l
d
T
i
m
e
=
t
i
m
e
;
t
i
m
e
+
=
(
E
c
h
o
@
R
a
n
d
o
m
V
a
r
i
a
t
e
[
Δ
t
i
m
e
]
)
»
0
.
0
0
0
6
5
3
4
5
8
O
u
t
[
]
=
4
4
.
5
8
7
3
Coding it up
Propensity
First, we need to compute the propensity at each step, so let’s make a function that generates the model-specific propensity function. This function will take two arguments: the current state and the current time. We include the time because some models may have time-dependent rates for treatment or mitigation strategies.
For most epidemiology models, the elements of the stoichiometric matrix are 0, 1, or -1, so if one or more of the compartments for a transition are exhausted the resulting propensity for that transition is 0 and that transition cannot be selected—that’s good. However, if the stoichiometry of one of the contributing compartments is 2 or more, then a non-zero propensity will still be computed when there is only 1 individual left in that state, and if that state gets selected, the resulting new state will have a negative entry—that’s not good. So we need to define a model-specific mask function that returns a list of 1’s and 0’s according to whether or not there are sufficient individuals present in the current state to support a transition.
I
n
[
]
:
=
d
e
f
i
n
e
P
r
o
p
e
n
s
i
t
y
F
u
n
c
t
i
o
n
[
r
a
t
e
s
:
{
_
_
}
,
(
s
M
a
t
r
i
x
_
)
?
(
M
a
t
r
i
x
Q
[
#
1
,
I
n
t
e
g
e
r
Q
]
&
)
,
v
a
r
s
I
n
:
{
_
_
}
,
t
i
m
e
_
S
y
m
b
o
l
]
:
=
M
o
d
u
l
e
[
{
m
a
s
k
,
v
a
r
s
=
T
h
r
o
u
g
h
@
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
[
t
]
}
,
m
a
s
k
=
F
u
n
c
t
i
o
n
[
s
t
a
t
e
,
B
o
o
l
e
[
A
n
d
@
@
N
o
n
N
e
g
a
t
i
v
e
/
@
(
s
t
a
t
e
+
#
)
]
&
/
@
s
M
a
t
r
i
x
]
;
F
u
n
c
t
i
o
n
[
{
s
,
t
}
,
(
r
a
t
e
s
/
.
T
h
r
e
a
d
[
v
a
r
s
s
]
/
.
t
i
m
e
t
)
m
a
s
k
[
s
]
]
]
This code generates the propensity function for our SEIR model:
I
n
[
]
:
=
p
r
o
p
e
n
s
i
t
y
=
d
e
f
i
n
e
P
r
o
p
e
n
s
i
t
y
F
u
n
c
t
i
o
n
[
m
o
d
e
l
D
a
t
a
[
"
R
a
t
e
s
"
]
/
.
f
o
r
c
e
O
f
I
n
f
e
c
t
i
o
n
/
.
p
a
r
a
m
s
N
,
m
o
d
e
l
D
a
t
a
[
"
S
t
o
i
c
h
i
o
m
e
t
r
i
c
M
a
t
r
i
x
"
]
,
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
,
t
]
O
u
t
[
]
=
F
u
n
c
t
i
o
n
[
{
s
$
,
t
$
}
,
(
{
0
.
1
9
6
0
7
8
ℰ
[
t
]
,
0
.
1
4
2
8
5
7
ℐ
[
t
]
,
0
.
0
0
0
7
1
4
2
8
6
ℐ
[
t
]
[
t
]
}
/
.
T
h
r
e
a
d
[
v
a
r
s
$
2
6
5
5
1
4
s
$
]
/
.
t
t
$
)
m
a
s
k
$
2
6
5
5
1
4
[
s
$
]
]
And here is an example of its use:
I
n
[
]
:
=
p
r
o
p
e
n
s
i
t
y
[
s
t
a
t
e
,
t
i
m
e
]
O
u
t
[
]
=
{
1
0
8
.
0
3
9
,
7
9
.
,
3
7
1
.
3
}
Simulation loop
to begin the stochastic simulation, we need to instantiate the starting state and evolve it over time. We’ll capture each new state in a dynamic array
D
a
t
a
S
t
r
u
c
t
u
r
e
.
I
n
[
]
:
=
S
e
e
d
R
a
n
d
o
m
[
1
2
3
4
5
6
7
]
s
t
a
t
e
L
i
s
t
=
M
o
d
u
l
e
[
{
t
i
m
e
,
s
t
a
t
e
,
i
d
x
,
p
r
o
p
,
t
o
t
a
l
P
r
o
p
,
s
t
a
t
e
L
i
s
t
,
d
i
s
t
E
v
e
n
t
,
d
i
s
t
T
i
m
e
,
r
e
s
u
l
t
}
,
i
d
x
=
R
a
n
g
e
@
L
e
n
g
t
h
@
m
o
d
e
l
D
a
t
a
[
"
S
t
o
i
c
h
i
o
m
e
t
r
i
c
M
a
t
r
i
x
"
]
;
s
t
a
t
e
L
i
s
t
=
C
r
e
a
t
e
D
a
t
a
S
t
r
u
c
t
u
r
e
[
"
D
y
n
a
m
i
c
A
r
r
a
y
"
]
;
s
t
a
t
e
=
R
e
p
l
a
c
e
[
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
,
{
-
I
0
,
ℐ
-
>
I
0
,
_
0
}
,
{
1
}
]
/
.
p
a
r
a
m
s
N
;
t
i
m
e
=
0
.
;
s
t
a
t
e
L
i
s
t
[
"
A
p
p
e
n
d
"
,
{
t
i
m
e
,
s
t
a
t
e
}
]
;
W
h
i
l
e
[
t
i
m
e
<
1
2
0
&
&
(
t
o
t
a
l
P
r
o
p
=
T
o
t
a
l
[
p
r
o
p
=
p
r
o
p
e
n
s
i
t
y
[
s
t
a
t
e
,
t
i
m
e
]
]
)
>
0
,
d
i
s
t
E
v
e
n
t
=
E
m
p
i
r
i
c
a
l
D
i
s
t
r
i
b
u
t
i
o
n
[
p
r
o
p
i
d
x
]
;
s
t
a
t
e
+
=
m
o
d
e
l
D
a
t
a
[
"
S
t
o
i
c
h
i
o
m
e
t
r
i
c
M
a
t
r
i
x
"
]
〚
R
a
n
d
o
m
V
a
r
i
a
t
e
[
d
i
s
t
E
v
e
n
t
]
〛
;
d
i
s
t
T
i
m
e
=
E
x
p
o
n
e
n
t
i
a
l
D
i
s
t
r
i
b
u
t
i
o
n
[
t
o
t
a
l
P
r
o
p
]
;
t
i
m
e
+
=
R
a
n
d
o
m
V
a
r
i
a
t
e
[
d
i
s
t
T
i
m
e
]
;
s
t
a
t
e
L
i
s
t
[
"
A
p
p
e
n
d
"
,
{
t
i
m
e
,
s
t
a
t
e
}
]
;
]
;
r
e
s
u
l
t
=
N
o
r
m
a
l
@
s
t
a
t
e
L
i
s
t
;
s
t
a
t
e
L
i
s
t
[
"
D
r
o
p
A
l
l
"
]
;
r
e
s
u
l
t
]
;
/
/
A
b
s
o
l
u
t
e
T
i
m
i
n
g
O
u
t
[
]
=
{
1
.
1
2
6
9
1
,
N
u
l
l
}
I
n
[
]
:
=
L
e
n
g
t
h
@
s
t
a
t
e
L
i
s
t
O
u
t
[
]
=
2
9
7
8
I
n
[
]
:
=
L
a
s
t
@
s
t
a
t
e
L
i
s
t
O
u
t
[
]
=
{
1
1
2
.
2
9
3
,
{
0
,
0
,
9
9
3
,
7
}
}
At the end of the simulation, the state vector shows that there are no exposed or infectious individuals, 993 recovered, and 7 susceptible.
Result in same the form as used by NDSolve
The raw result from the iteration is a list of time and state pairs. A more convenient form would be a list of
I
n
t
e
r
p
o
l
a
t
i
n
g
F
u
n
c
t
i
o
n
objects, much like those returned by
N
D
S
o
l
v
e
. That’s not so hard to do. First, we get the times as a vector and the states as a list of vectors for each variable in the model:
I
n
[
]
:
=
{
t
i
m
e
s
,
s
t
a
t
e
s
}
=
T
r
a
n
s
p
o
s
e
@
N
o
r
m
a
l
@
s
t
a
t
e
L
i
s
t
;
s
t
a
t
e
s
=
T
r
a
n
s
p
o
s
e
@
s
t
a
t
e
s
;
Next, we compute
I
n
t
e
r
p
o
l
a
t
i
n
g
F
u
n
c
t
i
o
n
objects for each variable. Note that runs of constant state are replaced by just the first instance, and that we use an interpolation order of 0.
I
n
[
]
:
=
s
o
l
n
=
F
i
r
s
t
@
{
T
h
r
e
a
d
[
v
a
r
s
(
I
n
t
e
r
p
o
l
a
t
i
o
n
[
J
o
i
n
@
@
R
e
p
l
a
c
e
[
S
p
l
i
t
B
y
[
T
h
r
e
a
d
[
{
t
i
m
e
s
,
#
}
]
,
L
a
s
t
]
,
{
a
_
,
_
_
_
,
z
_
}
{
a
}
,
{
1
}
]
,
I
n
t
e
r
p
o
l
a
t
i
o
n
O
r
d
e
r
0
]
&
/
@
s
t
a
t
e
s
)
]
}
O
u
t
[
]
=
ℰ
I
n
t
e
r
p
o
l
a
t
i
n
g
F
u
n
c
t
i
o
n
D
o
m
a
i
n
:
{
{
0
.
,
1
0
4
.
}
}
O
u
t
p
u
t
:
s
c
a
l
a
r
,
ℐ
I
n
t
e
r
p
o
l
a
t
i
n
g
F
u
n
c
t
i
o
n
D
o
m
a
i
n
:
{
{
0
.
,
1
1
2
.
}
}
O
u
t
p
u
t
:
s
c
a
l
a
r
,
ℛ
I
n
t
e
r
p
o
l
a
t
i
n
g
F
u
n
c
t
i
o
n
D
o
m
a
i
n
:
{
{
0
.
,
1
1
2
.
}
}
O
u
t
p
u
t
:
s
c
a
l
a
r
,
I
n
t
e
r
p
o
l
a
t
i
n
g
F
u
n
c
t
i
o
n
D
o
m
a
i
n
:
{
{
0
.
,
8
0
.
9
}
}
O
u
t
p
u
t
:
s
c
a
l
a
r
I
n
[
]
:
=
e
p
i
d
e
m
i
c
M
o
d
e
l
P
l
o
t
[
#
,
{
0
,
1
1
5
}
,
"
m
o
d
e
l
N
a
m
e
"
m
o
d
e
l
N
a
m
e
,
I
m
a
g
e
S
i
z
e
L
a
r
g
e
]
&
/
@
{
s
o
l
n
,
s
o
l
n
O
D
E
}
/
/
F
l
i
p
V
i
e
w
O
u
t
[
]
=
ℰ
ℐ
ℛ
Click on the output above to flip between the stochastic and deterministic simulations.
Code optimization
The SEIR model is fairly simple and has an absorbing compartment
ℛ
(i.e., a sink), so it evaluates fairly quickly, and often terminates before the maximum time threshold is reached. If we use a population 100 times larger it will take much longer to evaluate. Running hundreds or thousands of simulations to study average behavior will take much longer. Finally, a bigger model (e.g., a gene regulatory network or a cell signalling network) will take much longer also. The simulation above took almost 1 second, so just imagine the effect of any one of the above factors on the running time!
There are two bottlenecks in the code:
◼
R
a
n
d
o
m
V
a
r
i
a
t
e
has a lot of overhead and is “slow” in this context because the distribution changes at every evaluation requiring an evaluation for every variate
◼
we can speed up the computation by going to the heart of the matter and skipping the overhead
◼
p
r
o
p
e
n
s
i
t
y
is slow because it makes use of nested functions and repeats some of the computations on every invocation
◼
we can speed up the computation by creating a function with a precomputed body
I won’t go into the details here, but the code in the initialization section at the end of the notebook has been optimized, giving a substantial speed up. It’s still top-level Wolfram Language code, however, and compilation is still and option for further improvement. Here is the usage message for the solver:
I
n
[
]
:
=
?
s
t
o
c
h
a
s
t
i
c
S
o
l
v
e
O
u
t
[
]
=
S
y
m
b
o
l
s
t
o
c
h
a
s
t
i
c
S
o
l
v
e
[
r
a
t
e
s
,
S
m
a
t
r
i
x
,
i
n
i
t
a
l
S
t
a
t
e
,
v
a
r
s
,
{
t
,
t
m
i
n
,
t
m
a
x
}
]
r
e
t
u
r
n
s
a
n
I
n
t
e
r
p
o
l
a
t
i
o
n
g
f
u
n
c
t
i
o
n
o
f
t
h
e
s
t
o
c
h
a
s
t
i
c
s
i
m
u
l
a
t
i
o
n
g
e
n
e
r
a
t
e
d
w
i
t
h
t
h
e
c
h
e
m
i
c
a
l
m
a
s
t
e
r
e
q
u
a
t
i
o
n
m
e
t
h
o
d
g
i
v
e
n
t
h
e
l
i
s
t
o
f
r
a
t
e
s
f
o
r
t
h
e
n
f
u
n
d
e
m
e
n
t
a
l
r
e
a
c
t
i
o
n
s
,
t
h
e
n
*
k
s
t
o
i
c
h
i
o
m
e
t
r
i
c
m
a
t
r
i
x
,
t
h
e
l
i
s
t
o
f
k
s
t
a
t
e
v
a
r
i
a
b
l
e
s
,
a
n
d
t
h
e
l
i
s
t
o
f
t
i
m
e
v
a
r
i
a
b
l
e
,
s
t
a
r
t
t
i
m
e
,
a
n
d
e
n
d
t
i
m
e
.
Scherzo—Flattening the curve
So, back to the SEIR model and the failure to achieve an outbreak in some runs. It all has to do with the probability of infection events happening, and therefore a sufficient number of them relative to the number of recovery events. Recall that the initial state has just one infected individual:
I
n
[
]
:
=
s
t
a
t
e
0
=
T
h
r
o
u
g
h
@
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
[
0
]
/
.
R
u
l
e
@
@
@
i
n
i
t
s
/
.
p
a
r
a
m
s
N
O
u
t
[
]
=
{
0
,
1
,
0
,
9
9
9
}
I
n
[
]
:
=
=
p
r
o
p
e
n
s
i
t
y
[
s
t
a
t
e
0
,
0
]
O
u
t
[
]
=
{
0
.
,
0
.
1
4
2
8
5
7
,
0
.
7
1
3
5
7
1
}
I
n
[
]
:
=
e
v
e
n
t
=
E
m
p
i
r
i
c
a
l
D
i
s
t
r
i
b
u
t
i
o
n
[
R
a
n
g
e
[
L
e
n
g
t
h
@
]
]
O
u
t
[
]
=
D
a
t
a
D
i
s
t
r
i
b
u
t
i
o
n
T
y
p
e
:
E
m
p
i
r
i
c
a
l
D
a
t
a
p
o
i
n
t
s
:
3
I
n
[
]
:
=
R
o
w
B
a
r
C
h
a
r
t
,
c
h
a
r
t
o
p
t
i
o
n
s
,
D
i
s
c
r
e
t
e
P
l
o
t
C
D
F
[
e
v
e
n
t
,
x
+
1
]
,
{
x
,
0
,
L
e
n
g
t
h
[
]
,
0
.
0
1
}
,
p
l
o
t
o
p
t
i
o
n
s
,
S
p
a
c
e
r
[
6
]
O
u
t
[
]
=
W
h
i
l
e
t
h
e
m
o
s
t
l
i
k
e
l
y
e
v
e
n
t
w
i
l
l
b
e
a
n
i
n
f
e
c
t
i
o
n
,
t
h
e
r
e
i
s
a
g
o
o
d
c
h
a
n
c
e
t
h
a
t
i
t
c
o
u
l
d
b
e
a
r
e
c
o
v
e
r
y
i
n
s
t
e
a
d
,
w
h
i
c
h
w
i
l
l
e
r
a
d
i
c
a
t
e
t
h
e
i
n
f
e
c
t
i
o
n
b
e
f
o
r
e
i
t
c
a
n
b
e
c
o
m
e
a
n
o
u
t
b
r
e
a
k
.
S
o
,
a
b
o
u
t
o
n
e
s
i
x
t
h
(
0
.
1
4
2
8
5
7
0
.
1
4
2
8
5
7
+
0
.
7
1
3
5
7
1
0
.
1
6
6
8
0
6
,
t
o
b
e
p
r
e
c
i
s
e
)
o
f
t
h
e
m
w
i
l
l
s
t
o
p
a
f
t
e
r
o
n
e
s
t
e
p
,
a
n
d
,
i
n
d
e
e
d
,
t
h
a
t
i
s
t
h
e
c
a
s
e
:
I
n
[
]
:
=
m
o
d
e
l
D
a
t
a
[
"
C
o
m
p
o
n
e
n
t
T
r
a
n
s
i
t
i
o
n
s
"
]
〚
T
a
b
l
e
[
R
a
n
d
o
m
V
a
r
i
a
t
e
[
e
v
e
n
t
]
,
1
0
0
0
0
]
〛
/
/
C
o
u
n
t
s
O
u
t
[
]
=
β
ℐ
[
t
]
→
ℰ
8
3
3
0
,
ℐ
γ
→
ℛ
1
6
7
0
I
n
[
]
:
=
%
[
ℐ
γ
→
ℛ
]
1
0
0
0
0
/
/
N
O
u
t
[
]
=
0
.
1
6
7
Effectiveness of early international travel restrictions
Examining how the propensity varies over the time course of the simulation offers some interesting insights. Here is a plot of the individual and total propensities from the first run of the stochastic simulation:
I
n
[
]
:
=
W
i
t
h
{
p
r
o
p
=
p
r
o
p
e
n
s
i
t
y
[
T
h
r
o
u
g
h
@
m
o
d
e
l
D
a
t
a
[
"
V
a
r
i
a
b
l
e
s
"
]
[
t
]
,
t
]
/
.
F
i
r
s
t
@
s
o
l
n
s
C
M
E
}
,
P
l
o
t
A
p
p
e
n
d
[
p
r
o
p
,
T
o
t
a
l
@
p
r
o
p
]
/
/
E
v
a
l
u
a
t
e
,
{
t
,
0
,
1
0
0
}
,
o
p
t
i
o
n
s
O
u
t
[
]
=
ℰ
ζ
→
ℐ
ℐ
γ
→
ℛ
β
ℐ
(
t
)
→
ℰ
t
o
t
a
l
I
n
a
w
a
y
,
t
h
i
s
p
l
o
t
g
i
v
e
u
s
a
w
i
n
d
o
w
i
n
t
o
t
h
e
p
r
o
c
e
s
s
t
o
d
e
s
i
g
n
m
i
t
i
g
a
t
i
o
n
s
t
r
a
t
e
g
i
e
s
.
I
t
’
s
q
u
i
t
e
c
l
e
a
r
t
h
a
t
e
a
r
l
y
i
n
t
h
e
p
r
o
c
e
s
s
,
u
p
t
o
a
b
o
u
t
d
a
y
1
5
,
t
h
e
p
r
o
c
e
s
s
i
s
a
l
m
o
s
t
e
x
c
l
u
s
i
v
e
l
y
t
h
e
i
n
f
e
c
t
i
o
n
s
t
e
p
(
β
ℐ
(
t
)
→
ℰ
)
a
n
d
t
h
a
t
s
t
e
p
c
o
n
t
i
n
u
e
s
t
o
d
o
m
i
n
a
t
e
u
p
t
o
a
b
o
u
t
d
a
y
3
5
T
h
a
t
s
u
g
g
e
s
t
s
t
h
a
t
t
h
e
w
i
n
d
o
w
o
f
o
p
p
o
r
t
u
n
i
t
y
f
o
r
s
u
p
p
r
e
s
s
i
o
n
i
s
t
h
e
f
i
r
s
t
t
w
o
w
e
e
k
s
a
f
t
e
r
t
h
e
f
i
r
s
t
i
n
f
e
c
t
i
o
u
s
i
n
d
i
v
i
d
u
a
l
a
r
r
i
v
e
s
(
f
o
r
t
h
i
s
s
c
e
n
a
r
i
o
)
,
a
n
d
a
f
t
e
r
t
h
a
t
t
i
m
e
t
h
e
r
e
w
i
l
l
b
e
a
n
o
u
t
b
r
e
a
k
r
e
g
a
r
d
l
e
s
s
o
f
s
u
p
p
r
e
s
s
i
o
n
m
e
a
s
u
r
e
s
p
u
t
i
n
p
l
a
c
e
.
S
u
p
p
r
e
s
s
i
o
n
m
e
a
s
u
r
e
s
i
n
c
l
u
d
e
c
o
m
p
l
e
t
e
c
l
o
s
u
r
e
o
f
b
o
r
d
e
r
s
,
m
a
n
d
a
t
o
r
y
i
s
o
l
a
t
i
o
n
a
t
p
o
r
t
s
o
f
e
n
t
r
y
,
s
c
r
e
e
n
i
n
g
a
n
d
q
u
a
r
a
n
t
i
n
e
o
f
d
e
t
e
c
t
e
d
c
a
s
e
s
a
t
p
o
r
t
s
o
f
e
n
t
r
y
,
a
n
d
i
m
m
e
d
i
a
t
e
q
u
a
r
a
n
t
i
n
e
o
f
n
e
w
l
y
d
e
t
e
c
t
e
d
c
a
s
e
s
e
l
s
e
w
h
e
r
e
i
n
t
h
e
p
o
p
u
l
a
t
i
o
n
.
S
o
c
i
a
l
d
i
s
t
a
n
c
i
n
g
a
n
d
l
o
c
k
-
d
o
w
n
o
r
d
e
r
s
a
r
e
m
i
t
i
g
a
t
i
n
g
s
t
r
a
t
e
g
i
e
s
t
h
a
t
o
n
l
y
l
e
s
s
e
n
t
h
e
e
f
f
e
c
t
o
f
a
n
o
u
t
b
r
e
a
k
.
In the early days of the pandemic, international travelers were introducing many more than one infected person into uninfected populations. Here is the effect of starting with 1, 2, 3, 4, or 5 initial infections:
I
n
[
]
:
=
s
o
l
n
s
C
M
E
⎵
I
0
=
T
a
b
l
e
[
W
i
t
h
[
{
p
a
r
a
m
s
S
c
e
n
a
r
i
o
=
<
|
p
a
r
a
m
s
,
I
0
i
0
|
>
}
,
W
i
t
h
[
{
p
a
r