Message Boards

GROUPS:

Staff Picks Biological Sciences Medical Sciences Software Development Graphics and Visualization Wolfram Language Modeling Numerical Computation

Robert Nachbar

Stochastic Epidemiology Models with Applications to the COVID-19

Robert Nachbar, Wolfram Research Inc.

Posted 1 year ago

5025 Views

11 Replies

13 Total Likes

Follow this post

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus

Stochastic simulation of SEIR model

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

(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

CompartmentalModeling`

package (download here), defining the transitions of the model, and computing the kinetic elements:

In[]:=

<<CompartmentalModeling`

In[]:=

modelName="SEIR";model=

βλ[t]

→

ℰ,ℰ

→

ℐ,ℐ

→

ℛ;forceOfInfection={λ[t_]ℐ[t]};

In[]:=

modelData=KineticCompartmentalModel[model/.forceOfInfection,t];vars=modelData["Variables"];(odes=modelData["Equations"])//Column

Out[]=

′

ℰ

[t]-ζℰ[t]+βℐ[t][t]

′

ℐ

[t]ζℰ[t]-γℐ[t]

′

ℛ

[t]γℐ[t]

′



[t]-βℐ[t][t]

We’ll parameterize the model in terms of the basic reproduction number

ℜ

, and use an incubation period (

) and duration of infectiousness (

) characteristic of COVID-19.

In[]:=

betaSoln=First@Solve

ℜ



β

,β

Out[]=

β

ℜ





In[]:=

params=<|

ℜ

5.,ββ,ζ

5.1

,γ

,

,I01|>/.betaSolnparamsN=Normal@(params//.params)

Out[]=



ℜ

5.,β

ℜ



,ζ0.196078,γ0.142857,1000,I01

Out[]=

{

ℜ

5.,β0.000714286,ζ0.196078,γ0.142857,1000,I01}

In[]:=

inits=Join[{[0]-I0,ℐ[0]I0},Thread[Through@Complement[vars,{,ℐ}][0]0]]

Out[]=

{[0]-I0+,ℐ[0]I0,ℰ[0]0,ℛ[0]0}

First, let’s see how the usual ODE solution behaves over 120 days (the code for

epidemicModelPlot

is in the initialization section at the end of the notebook).

In[]:=

tmax=120;

In[]:=

solnODE=First@NDSolve[{odes,inits}/.paramsN,vars,{t,0,tmax}];

In[]:=

epidemicModelPlot[solnODE,{0,tmax},"modelName"modelName,ImageSizeLarge]

Out[]=

	ℰ
	ℐ
	ℛ
	

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

stochasticSolve

is in the initialization section at the end of the notebook).

In[]:=

SeedRandom[1425367]solnsCME=Table[stochasticSolve[modelData["Rates"]/.paramsN,modelData["StoichiometricMatrix"],Through@modelData["Variables"][0]/.Rule@@@inits/.paramsN,modelData["Variables"],{t,0,tmax}],500];//Timing

Out[]=

{35.578,Null}

Here is a plot of the first 20 runs:

In[]:=

epidemicModelPlot[Join@@Take[solnsCME,20],{0,tmax},"modelName"modelName,ImageSizeLarge]

Out[]=

	ℰ
	ℐ
	ℛ
	

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:

In[]:=

TabView[{"median with bands"epidemicModelPlot[solnsCME,{0,tmax},"modelName""SIR",ImageSizeLarge],"median"epidemicModelPlot[solnsCME,{0,tmax},"modelName"modelName,ImageSizeLarge,"showBands"False]}]

Out[]=

median with bands

median

	ℰ
	ℐ
	ℛ
	

Bands represent 0.025 and 0.975 quantiles

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

max

we see that there are essentially two kinds of behavior:

In[]:=

Histogram[Quiet[#[tmax]/.solnsCME],{10},ChartStyleepiColor[#],PlotLabel#,PlotRange400,ImageSize472]&/@{,ℛ}//Row[#,Spacer[18]]&

Out[]=

We can partition the ensemble into an “outbreak” and “not an outbreak” subsets:

In[]:=

{outbreak,notOutbreak}=Quiet[SortBy[solnsCME,(ℛ[tmax]/.#)<50&]//SplitBy[#,(ℛ[150]/.#)<50&]&];Length/@%

Out[]=

{394,106}

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.

In[]:=

MapThread[#2epidemicModelPlot[#1,{0,tmax},"modelName"modelName,ImageSizeLarge]&,{{outbreak,notOutbreak,solnODE},{"outbreak","not outbreak","deterministic"}}]//TabView

Out[]=

outbreak

not outbreak

deterministic

	ℰ
	ℐ
	ℛ
	

Bands represent 0.025 and 0.975 quantiles

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

modelData

Association, and is shown in the table below. The key

"ComponentTransitions"

gives the data in the first column, and

"Rates"

gives that in the last column. The key

"StoichiometricMatrix"

gives the table in the middles column, whose headings are obtained from the

"Variables"

key.

In[]:=

MapThread[Join,{List/@modelData["ComponentTransitions"],RotateRight/@modelData["StoichiometricMatrix"],List/@modelData["Rates"]}]//

format table

Out[]=

	Stoichiometry
Transition		ℰ	ℐ	ℛ	Propensity
 βℐ[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:

In[]:=

state=RandomInteger[{0,1000},Length@modelData["Variables"]]time=RandomReal[{0,50}]

Out[]=

{551,554,561,940}

Out[]=

44.5866

The propensity vector



is the model rate vector evaluated at the current state and time:

In[]:=

=modelData["Rates"]/.Thread[Through[modelData["Variables"][t]]state]/.paramsN/.ttime

Out[]=

{108.039,79.1429,371.971}

The individual propensities define an empirical distribution:

In[]:=



event

=EmpiricalDistribution[(Range[Length@])]

Out[]=

DataDistribution

Type: Empirical

Data points: 3



Graphically, we can see how the propensities on the left contribute to the cumulative distribution on the right:

In[]:=

RowBarChart,

chart options

,DiscretePlotCDF[



event

,x+1],{x,0,Length[],0.01},

plot options

,Spacer[6]

Out[]=

We can use

RandomVariate

raw a random sample from the distribution to select the transition. Here is a list of 10 samples:

In[]:=

Table[RandomVariate[



event

],10]

Out[]=

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

In[]:=

oldState=state;state+=modelData["StoichiometricMatrix"]〚Echo[#,"",modelData["ComponentTransitions"]〚#〛&]&@RandomVariate[



event

]〛;TableForm[{oldState,state}//Thread,TableHeadings{modelData["Variables"],{"before","after"}}]

ℐ

→

ℛ

Out[]//TableForm=

	before	after
ℰ	551	551
ℐ	554	553
ℛ	561	562
	940	940

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:

In[]:=



Δtime

=ExponentialDistribution[Total[]]

Out[]=

ExponentialDistribution[559.154]

In[]:=

Withtmax=

Log[2]

Total[]

6,RowPlotPDF[



Δtime

,x],{x,0,tmax},

PDF options

,PlotCDF[



Δtime

,x],{x,0,tmax},ImageSize472,

CDF options

,Spacer[6]

Out[]=

Again,

RandomVariate

can be used to draw a random sample from the exponential distribution to select the time to the next event, and increment the time:

In[]:=

oldTime=time;time+=(Echo@RandomVariate[



Δtime

])

0.000653458

Out[]=

44.5873

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.

In[]:=

definePropensityFunction[rates:{__},(sMatrix_)?(MatrixQ[#1,IntegerQ]&), varsIn:{__},time_Symbol]:=Module[{mask,vars=Through@modelData["Variables"][t]},mask=Function[state,Boole[And@@NonNegative/@(state+#)]&/@sMatrix];Function[{s,t},(rates/.Thread[varss]/.timet)mask[s]]]

This code generates the propensity function for our SEIR model:

In[]:=

propensity=definePropensityFunction[modelData["Rates"]/.forceOfInfection/.paramsN,modelData["StoichiometricMatrix"],modelData["Variables"],t]

Out[]=

Function[{s$,t$},({0.196078ℰ[t],0.142857ℐ[t],0.000714286ℐ[t][t]}/.Thread[vars$265514s$]/.tt$)mask$265514[s$]]

And here is an example of its use:

In[]:=

propensity[state,time]

Out[]=

{108.039,79.,371.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

DataStructure

In[]:=

SeedRandom[1234567]stateList=Module[{time,state,idx,prop,totalProp,stateList,distEvent,distTime,result},idx=Range@Length@modelData["StoichiometricMatrix"];stateList=CreateDataStructure["DynamicArray"];state=Replace[modelData["Variables"],{-I0,ℐ->I0,_0},{1}]/.paramsN;time=0.;stateList["Append",{time,state}];While[time<120&&(totalProp=Total[prop=propensity[state,time]])>0,distEvent=EmpiricalDistribution[propidx];state+=modelData["StoichiometricMatrix"]〚RandomVariate[distEvent]〛;distTime=ExponentialDistribution[totalProp];time+=RandomVariate[distTime];stateList["Append",{time,state}];];result=Normal@stateList;stateList["DropAll"];result];//AbsoluteTiming

Out[]=

{1.12691,Null}

In[]:=

Length@stateList

Out[]=

2978

In[]:=

Last@stateList

Out[]=

{112.293,{0,0,993,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

InterpolatingFunction

objects, much like those returned by

NDSolve

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

In[]:=

{times,states}=Transpose@Normal@stateList;states=Transpose@states;

Next, we compute

InterpolatingFunction

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.

In[]:=

soln=First@{Thread[vars(Interpolation[Join@@Replace[SplitBy[Thread[{times,#}],Last],{a_,___,z_}{a},{1}],InterpolationOrder0]&/@states)]}

Out[]=

ℰInterpolatingFunction

Domain: {{0.,104.}}

Output: scalar

,ℐInterpolatingFunction

Domain: {{0.,112.}}

Output: scalar

,ℛInterpolatingFunction

Domain: {{0.,112.}}

Output: scalar

,InterpolatingFunction

Domain: {{0.,80.9}}

Output: scalar



In[]:=

epidemicModelPlot[#,{0,115},"modelName"modelName,ImageSizeLarge]&/@{soln,solnODE}//FlipView

Out[]=

	ℰ
	ℐ
	ℛ
	

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:

◼

RandomVariate

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

◼

propensity

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:

In[]:=

?stochasticSolve

Out[]=

Symbol
stochasticSolve[rates, Smatrix, initalState, vars, {t, tmin, tmax}] returns an Interpolationg function of the stochastic simulation generated with the chemical master equation method given the list of rates for the n fundemental reactions, the n*k stoichiometric matrix, the list of k state variables, and the list of time variable, start time, and end time.

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:

In[]:=

state0=Through@modelData["Variables"][0]/.Rule@@@inits/.paramsN

Out[]=

{0,1,0,999}

In[]:=

=propensity[state0,0]

Out[]=

{0.,0.142857,0.713571}

In[]:=



event

=EmpiricalDistribution[Range[Length@]]

Out[]=

DataDistribution

Type: Empirical

Data points: 3



In[]:=

RowBarChart,

chart options

,DiscretePlotCDF[



event

,x+1],{x,0,Length[],0.01},

plot options

,Spacer[6]

Out[]=

While the most likely event will be an infection, there is a good chance that it could be a recovery instead, which will eradicate the infection before it can become an outbreak. So, about one sixth (

0.142857

0.142857+0.713571

0.166806

, to be precise) of them will stop after one step, and, indeed, that is the case:

In[]:=

modelData["ComponentTransitions"]〚Table[RandomVariate[



event

],10000]〛//Counts

Out[]=



βℐ[t]

→

ℰ8330,ℐ

→

ℛ1670

In[]:=

%[ℐ

→

ℛ]

10000

//N

Out[]=

0.167

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:

In[]:=

With{prop=propensity[Through@modelData["Variables"][t],t]/.First@solnsCME},PlotAppend[prop,Total@prop]//Evaluate,{t,0,100},

options



Out[]=

	ℰ ζ → ℐ
	ℐ γ → ℛ
	 βℐ(t) → ℰ
	total

In a way, this plot give us a window into the process to design mitigation strategies. It’s quite clear that early in the process, up to about day 15, the process is almost exclusively the infection step (



βℐ(t)

→

ℰ

) and that step continues to dominate up to about day 35 That suggests that the window of opportunity for suppression is the first two weeks after the first infectious individual arrives (for this scenario), and after that time there will be an outbreak regardless of suppression measures put in place. Suppression measures include complete closure of borders, mandatory isolation at ports of entry, screening and quarantine of detected cases at ports of entry, and immediate quarantine of newly detected cases elsewhere in the population. Social distancing and lock-down orders are mitigating strategies that only lessen the effect of an outbreak.

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:

In[]:=

solnsCME⎵I0=Table[With[{paramsScenario=<|params,I0i0|>},With[{pa