Epidemiological Models for Influenza and COVID-19
Epidemiological Models for Influenza and COVID-19
Robert B. Nachbar
Wolfram Research, Inc.
Table of Contents
Table of Contents
Introduction
Introduction
The COVID-19 outbreak, initially in China and now throughout the world, has captured the interest of a large number of organizations and individuals alike. Some effort has been spent to model (or at least visualize) the geographic spread (; ). There have also been some reports of epidemiological models (; ; ; ; ; ) that have been developed in an effort to estimate parameters that can be used to project the severity of the outbreak, its duration and the mortality rate.
This report aims to put some of these modeling efforts into perspective so that conclusions and predictions can be better understood. We will be using compartmental models that allow one to describe the flow of individuals from one health state to another. We will attempt to employ just the right kinds of compartments and connections that are supported by the available data.
We explore several epidemiological models for the COVID-19 outbreak and use data from various provinces to estimate model parameters. We primarily want to see how similar or dissimilar the outbreaks are to each other, point out the shortcomings and make some suggestions for improvements.
Background and Notes
Background and Notes
The following subsections contain useful information about the COVID-19 infection and the outbreak in Wuhan, China; it is by no means comprehensive.
Summary
Summary
Symptom reports for COVID-19 first occurred in Wuhan, capital and largest city in China’s Hubei province, in late December. COVID-19 is caused by the SARS-CoV-2 virus, which appears to have jumped to humans from an animal at the Huanan Seafood Market (though the exact animal of origin is unclear).
The incubation period for COVID-19 is usually around five days, with reports varying from one to 14 days (and in one case, 27 days). The median time from development of symptoms to death between six and 41 days, with a median of 14 days. Available data from South Korea suggests approximately a 14-day recovery period as well.
Timeline of Events
Timeline of Events
This is a timeline of major events related to the outbreak in Hubei province.
TimelinePlotLabeled[#1,Tooltip[#2,#3]]&@@@,
In[]:=
Out[]=
Data and Functions
Data and Functions
This section contains the computational resources required for replicating the results in this paper, including data sources and custom functions for cleaning and processing the data.
Epidemic Data
Epidemic Data
Time series for reported cases, deaths and recoveries from each province are available from the Wolfram Data Repository. This can be imported using ResourceData.
Import and view the epidemic data:
covid19Data=ResourceData["Epidemic Data for Novel Coronavirus COVID-19"]
In[]:=
Out[]=
Functions
Functions
To simplify this workflow, I have made use of a custom modeling package and written several functions. The following are brief descriptions of these functions and their uses. (Full code for these functions can be found in the of this notebook.)
The CompartmentalModeling package is used for computations directly dealing with compartmental models. Given a list of transitions for a model, KineticComparmentalModel computes the equations for a model, and CompartmentalModelGraph generates a graph of the model. For clarity, we use a form of schematic diagram known as a Petri net.
Graph a basic SIR model:
CompartmentalModelGraphℐ,ℐℛ,
βλ[t]
→
γ
→
In[]:=
Out[]=
Generate the equations for the SIR model:
KineticCompartmentalModelℐ,ℐℛ,t//Last
βλ[t]
→
γ
→
In[]:=
{[t]-γℐ[t]+β[t]λ[t],[t]γℐ[t],[t]-β[t]λ[t]}
′
ℐ
′
ℛ
′
Out[]=
The fitDataWDR function returns reported infections, deaths and recoveries for a given locale and start date. It includes a "makeCorrectionsForHubei" option for removing the jump in confirmed cases in Hubei that occurred when the testing method changed.
Get usage info for the function:
?fitDataWDR
In[]:=
Out[]=
Use the function to retrieve the cases in Hubei for a plot:
ListPlotfitDataWDR,"1 Jan 2020",
In[]:=
Out[]=
A sumSquaredError function is defined to give the sum of squared errors for a model, given a set of parameters. I have also included additional helper functions toLog, fromLog, toLogit and fromLogit that are useful for passively constraining parameters to the domain or , respectively, during optimization.
[0,∞]
[0,1]
Compartmental Models
Compartmental Models
Basic Compartmental Models: SIR, SEIR
Basic Compartmental Models: SIR, SEIR
We will be using compartmental models, which have had numerous applications in biology, ecology, chemistry and medicine. For our purposes, each compartment represents a group of individuals in the same health state, for example susceptible or infectious. The connections between compartments indicate the direction and rate of movement from one health state to another.
One of the simplest compartmental models for epidemiology has three compartments: susceptible, infectious, and recovered or SIR. Susceptible individuals come in contact with infectious individuals and become infected. After some period of time, infectious individuals recover, are not longer infectious and have permanent immunity. The process can be described as a set of transitions.
Starting with a list of transitions and their labels, we can use CompartmentalModelGraph to visualize this model.
Define and display transitions for the SIR model:
modelSIR=ℐ,ℐℛ;GridThread@{modelSIR,{"infection","recovery"}},
βλ[t]
→
γ
→
In[]:=
βλ[t] → | infection |
ℐ γ → | recovery |
Out[]=
Generate a graph of the SIR model:
CompartmentalModelGraphmodelSIR,
In[]:=
Out[]=
There is code in the at the end of the notebook for generating these graphs.
Each compartment becomes a time-dependent variable in the model. The coefficient is the transmission rate constant, and is the force of infection. It is a function of all the infectious compartments in the model (in this case, just ), and thus is time dependent. The parameter is the recovery rate constant and is defined as , where is the average duration of infection.
β
λ[t]
ℐ[t]
γ
γ
1
τ
τ
This transmission model describes the rate of change of each of the compartments, which can be modeled mathematically as a system of ordinary differential equations (ODEs).
KineticCompartmentalModel[modelSIR,t]//Last//Append[#〚{3,1,2}〛,λ[t]ℐ[t]]&//Column[#,Left,BaseStyle{18},Spacings{Automatic,{{Automatic},1,Automatic}}]&closeCellGroup[]
In[]:=
′ |
′ ℐ |
′ ℛ |
λ[t]ℐ[t] |
Out[]=
The number of infectious individuals in the population is the prevalence of the disease. The number of individuals that become infected per unit time is called the incidence.
ℐ[t]
Many epidemiologic models include population demographics, that is, birth and natural death. For comparison, this is the SIR model with demographics.
Λ ↦ | birth |
βλ[t] → | infection |
ℐ γ → | recovery |
μ ⇥ | death |
ℐ μ ⇥ | death |
ℛ μ ⇥ | death |
Out[]=
CompartmentalModelGraph[modelSIRwithDemographics,GraphLayout{"LayeredDigraphEmbedding","RootVertex"Λ,"Orientation"Left},PerformanceGoal"Quality",BaseStyle{18},ImageSize450]closeCellGroup[]
In[]:=
Out[]=
odesSIRwithDemographics=KineticCompartmentalModel[modelSIR,t]//Last;Column[Append[odesSIRwithDemographics,λ[t]ℐ[t]],Left,BaseStyle{18},Spacings{Automatic,{{Automatic},1,Automatic}}]closeCellGroup[]
In[]:=
′ ℐ |
′ ℛ |
′ |
λ[t]ℐ[t] |
Out[]=
The shape for in the schematic diagram is different because its units are , whereas all the other rates are .
Λ
persons
time
1
time
From the base SIR model, we add an exposed compartment to form an SEIR, which is the base model used in this exploration.
βλ[t] → | infection |
ℰ ζ → | incubation |
ℐ γ → | recovery |
Out[]=
CompartmentalModelGraph[modelSEIR,GraphLayout{"LayeredDigraphEmbedding","RootVertex","Orientation"Left},PerformanceGoal"Quality",BaseStyle{18},ImageSizeLarge]closeCellGroup[]
In[]:=
Out[]=
Additional Compartments and Assumptions
Additional Compartments and Assumptions
More complex models will be considered throughout this report, adding new compartments to the SEIR model to account for dead, quarantined and infectious (symptomatic/asymptomatic) populations. We’ll apply each of these models to the Hubei outbreak, hoping to find one with a reasonable structure that recapitulates the reported data.
Each alternative will be tested using full and effective population sizes. Full population sizes will come from Wolfram|Alpha; values of , , , and will be tested for effective population sizes (represented as ).
50×
3
10
100×
3
10
500×
3
10
1×
6
10
50×
6
10
This paper generally assumes:
No birth or natural mortality due to the short time span of the coverage of the model. Death due to disease will be included.
◼
No migration into or out of the population.
◼
An initial count of 10 for infected individuals (). (Du et al. [] estimated the number to be 7.78 in Wuhan.)
I
0
◼
First infections occurring on December 20, 2019, or December 30, 2019 (with the first confirmed reports). We don’t know with certainty when the first infection occurred. Du et al. () assumed December 1, 2019.
◼
No variation due to age or gender. Even though the severity of the infection and risk of death have been reported to be age and gender dependent, the available data for modeling does not have sufficient detail to support such stratification.
◼
Standard incidence, which means that the force of infection is normalized by the size of the population (the alternative is mass action incidence, which does not normalize).
◼
All current cases are reported. Underreporting of cases is a serious concern, primarily due to the infection being asymptomatic in many individuals, especially children (). This paper makes no attempt to extrapolate and only uses numbers directly from the reports.
◼
This work was carried out over a number of days, with the available Wolfram Data Repository data keeping pace. However, to assure a consistent analysis and interpretation of results, we’ll work with data for the period of January 21, 2020, to February 26, 2020.
◼
More recent data can be used to test the models.
◼
SEIRD Model
SEIRD Model
We start by introducing a compartment to account for those who have already died from the disease. An effective population size with standard incidence will be used. This assumes that individuals in the population cannot come in contact with all other individuals in the population (in contrast to models for influenza and SARS []). may be adjusted for individual computations; initially we’ll use 100,000 individuals.
Adding the dead compartment brings the model to five compartments and four transitions.
Define the transitions for the SEIRD model:
modelSEIRD=ℰ,ℰℐ,ℐℛ,ℐ;
βλ[t]
→
ζ
→
γ
→
μ
→
In[]:=
Label the transitions and show them in a grid:
GridThread@{modelSEIRD,{"infection","incubation","recovery","death"}},
In[]:=
βλ[t] → | infection |
ℰ ζ → | incubation |
ℐ γ → | recovery |
ℐ μ → | death |
Out[]=
Display the graph for the model:
CompartmentalModelGraphmodelSEIRD,
In[]:=
Out[]=
Next we define and solve the ODEs for the model.
Generate the ODEs for the model:
{varsSEIRD,odesSEIRD}=Values@KineticCompartmentalModel[modelSEIRD,t]〚{1,-1}〛;
In[]:=
Display the ODEs in a column:
Column[odesSEIRD]
In[]:=
′ |
′ ℰ |
′ ℐ |
′ ℛ |
′ |
Out[]=
Define for the SEIRD model:
λ(t)
forceOfInfectionSEIRD=λ[t_];
ℐ[t]
-[t]
In[]:=
Solve for each variable:
{susceptibleSEIRD,exposedSEIRD,infectedSEIRD,recoveredSEIRD,deadSEIRD}={,ℰ,ℐ,ℛ,}/.ParametricNDSolve[Join[odesSEIRD/.forceOfInfectionSEIRD,{[0]-I0,ℰ[0]0,ℐ[0]I0,ℛ[0]0,[0]0}],varsSEIRD,{t,0,500},{,I0,β,ζ,γ,μ}];
In[]:=
SEIQRD Model
SEIQRD Model
An SEIQRD model uses an additional quarantined compartment. Quarantined individuals do not contribute to the force of infection. The duration of circulation of infected individuals is constant even though we know it was initially infinitely long and later became very short as isolation and quarantine measures were instituted; however, no detailed information about when changes occurred is available. The duration of infection and the mortality rate for circulating and quarantined infectious individuals is the same.
[t]
The additional quarantine transition is between infection and recovery for a total of six compartments and seven transitions.
Define the transitions for the SEIQRD model:
modelSEIQRD=ℰ,ℰℐ,ℐ,ℐℛ,ℛ,ℐ,;
βλ[t]
→
ζ
→
δ
→
γ
→
γ
→
μ
→
μ
→
In[]:=
Label the transitions and show them in a grid:
GridThread@{modelSEIQRD,{"infection","incubation","quarantine","recovery","recovery","death","death"}},
In[]:=
βλ[t] → | infection |
ℰ ζ → | incubation |
ℐ δ → | quarantine |
ℐ γ → | recovery |
γ → | recovery |
ℐ μ → | death |
μ → | death |
Out[]=
Display the graph for the model:
CompartmentalModelGraphmodelSEIQRD,
In[]:=
Out[]=
Now we can define and solve the ODEs for the model.
Generate the variables and ODEs for the SEIQRD model:
{varsSEIQRD,odesSEIQRD}=Values@KineticCompartmentalModel[modelSEIQRD,t]〚{1,-1}〛;
In[]:=
Display the ODEs in a column:
Column[odesSEIQRD]
In[]:=
′ |
′ ℰ |
′ ℐ |
′ |
′ ℛ |
′ |
Out[]=
Define for the model:
λ(t)
forceOfInfectionSEIQRD=λ[t_];
ℐ[t]
-[t]-[t]
In[]:=
Solve for each variable:
{susceptibleSEIQRD,exposedSEIQRD,infectedSEIQRD,quarantinedSEIQRD,recoveredSEIQRD,deadSEIQRD}={,ℰ,ℐ,,ℛ,}/.ParametricNDSolve[Join[odesSEIQRD/.forceOfInfectionSEIQRD,{[0]-I0,ℰ[0]0,ℐ[0]I0,[0]0,ℛ[0]0,[0]0}],varsSEIQRD,{t,0,500},{,I0,β,ζ,γ,δ,μ}];
In[]:=
SEIsaRD Model
SEIsaRD Model
Finally, we add symptomatic and asymptomatic infectious compartments (ignoring quarantine) for an SEIsaRD model. The duration of infection for symptomatic and asymptomatic infectious individuals is the same. We have two different recovered compartments for bookkeeping. Death due to disease occurs only for symptomatic infectious individuals.
This model has a total of seven compartments and seven transitions.
Define the transitions for the SEIsaRD model:
modelSEIsaRD=ℰ,ℰ,ℰ,,,;
βλ[t]
→
ζρ
→
ℐ
"s"
ζ(1-ρ)
→
ℐ
"a"
ℐ
"s"
γ
→
ℛ
"s"
ℐ
"a"
γ
→
ℛ
"a"
ℐ
"s"
μ
→
In[]:=
Label the transitions and show them in a grid:
GridThread@{modelSEIsaRD,{"infection","incubation, symptomatic","incubation, asymptomatic","recovery, symptomatic","recovery, asymptomatic","death, symptomatic"}},
In[]:=
βλ[t] → | infection |
ℰ ζρ → ℐ s | incubation, symptomatic |
ℰ ζ(1-ρ) → ℐ a | incubation, asymptomatic |
ℐ s γ → ℛ s | recovery, symptomatic |
ℐ a γ → ℛ a | recovery, asymptomatic |
ℐ s μ → | death, symptomatic |
Out[]=
Display the graph for the model:
CompartmentalModelGraphmodelSEIsaRD,
In[]:=
Out[]=
Now we can define and solve the ODEs for the model.
Generate the variables and ODEs for the SEIsaRD model:
{varsSEIsaRD,odesSEIsaRD}=Values@KineticCompartmentalModel[modelSEIsaRD,t]〚{1,-1}〛;
In[]:=
Display the ODEs in a column:
Column[odesSEIsaRD]
In[]:=
′ ℐ s |
′ ℰ |
′ |
′ ℐ a ℐ a |
′ ℐ s ℐ s ℐ s |
′ ℛ a ℐ a |
′ ℛ s ℐ s |
Out[]=
Define for the model:
λ(t)
forceOfInfectionSEIsaRD=λ[t_][t]+[t];
ℐ
"s"
ℐ
"a"
-[t]
In[]:=
Solve for each variable:
{susceptibleSEIsaRD,exposedSEIsaRD,infectedSymptomaticSEIsaRD,infectedAsymptomaticSEIsaRD,recoveredSymptomaticSEIsaRD,recoveredAsymptomaticSEIsaRD,deadSEIsaRD}={,ℰ,,,,,}/.ParametricNDSolve[Join[odesSEIsaRD/.forceOfInfectionSEIsaRD,{[0]-I0,ℰ[0]0,[0]ρI0,[0](1-ρ)I0,[0]0,[0]0,[0]0}],varsSEIsaRD,{t,0,500},{,I0,β,ζ,ρ,γ,μ}];
ℐ
"s"
ℐ
"a"
ℛ
"s"
ℛ
"a"
ℐ
"s"
ℐ
"a"
ℛ
"s"
ℛ
"a"
In[]:=
Hubei Outbreak—SEIRD
Hubei Outbreak—SEIRD
In this section, we apply the SEIRD model to the Hubei outbreak with two different start dates and two effective population sizes, then compare the results from each fit.
Functions
Functions
Here, several functions are defined to compute and visualize the fit for the SEIRD model in Hubei. Note that fitErrorSEIRD expects the parameters to be log transformed and uses fromLog to back transform before using them in the functions for the model variables.
Function for computing the sum of squared errors between the data and the SEIRD model:
ClearAll[fitErrorSEIRD]fitErrorSEIRD[fitData_, pop_, I0_, (logBeta_)?NumberQ, (logZeta_)?NumberQ, (logGamma_)?NumberQ, (logMu_)?NumberQ] := Module[params = MapThread[#1 #2[#3]&, {β, ζ, γ, μ}, {fromLog, fromLog, fromLog, fromLog}, {logBeta, logZeta, logGamma, logMu}], sumSquaredError[ fitData, infectedSEIRD[pop, I0, β, ζ, γ, μ], recoveredSEIRD[pop, I0, β, ζ, γ, μ], deadSEIRD[pop, I0, β, ζ, γ, μ], params ] ]
In[]:=
Function for plotting the data and the SEIRD model:
ClearAll[plotFitSEIRD]plotFitSEIRD[fitData_, {pop_, I0_, beta_, zeta_, gamma_, mu_}, opts:OptionsPattern[{Plot, ListPlot}]] := Module[{tmax = 100, ymax = 40000, scale = 50000/pop}, Show[ Plot[ scale*susceptibleSEIRD[pop, I0, beta, zeta, gamma, mu][t], exposedSEIRD[pop, I0, beta, zeta, gamma, mu][t], infectedSEIRD[pop, I0, beta, zeta, gamma, mu][t], recoveredSEIRD[pop, I0, beta, zeta, gamma, mu][t], deadSEIRD[pop, I0, beta, zeta, gamma, mu][t] // Evaluate, {t, 0, tmax}, PlotRangePadding -> {Automatic, {Scaled[0.03], Scaled[0.03]}}, PlotLegends -> {N[scale]*, ℰ, ℐ, ℛ, }, PlotLabel -> StringForm[" = ``, SSE = ``", pop, fitErrorSEIRD[fitData, pop, I0, Log[beta], Log[zeta], Log[gamma], Log[mu]]], FrameLabel -> {"time", "# people"}, PlotStyle -> ColorData["Crayola"]["Green"], ColorData["Crayola"]["Sunglow"], ColorData["Crayola"]["Scarlet"], ColorData["Crayola"]["NavyBlue"], ColorData["Crayola"]["OuterSpace"] , Evaluate@FilterRules[{opts}, Options[Plot]], PlotRange -> {{0, tmax}, {0, ymax}}, GridLines -> {None, None} ], ListPlot[fitData, PlotStyle -> ColorData["Crayola"]["Scarlet"], ColorData["Crayola"]["NavyBlue"], ColorData["Crayola"]["OuterSpace"] , Evaluate@FilterRules[{opts}, Options[ListPlot]] ] ] ]
In[]:=
Outbreak Start—December 20, 2019
Outbreak Start—December 20, 2019
Fitting Data
Fitting Data
We use the fitDataWDR function described above to fit a model to the data for Hubei, with an assumed initial infection date of December 20, 2019.
Compute a fit for the data:
fitData1=fitData=fitDataWDR,"20 Dec 2019",;
In[]:=
Plot the fit:
ListPlotfitData,
In[]:=
Out[]=
Optimization—100000
Optimization—
100000
We assume the effective population size is 100,000.
fitErrorSEIRD[fitData1,100000,10,Sequence@@toLog@{0.32,0.60,0.072,0.0028}]
In[]:=
1.48125×
10
10
Out[]=
randomSampleWithError=With[{delta=Log[10.]{-1,+1},init=Log@{0.32,0.60,0.072,0.0028}},{Thread[{β,ζ,γ,μ}Exp@{##}],fitErrorSEIRD[fitData,100000,10,##]}&@@@Prepend[Transpose[RandomReal[#+delta,50000]&/@init],init]]//SortBy[#,Last]&;//Timing
In[]:=
{524.2,Null}
Out[]=
Take[SortBy[Last]@randomSampleWithError,15]//tfn
In[]:=
1 | 2 | |
1 | {β0.767774,ζ0.0631035,γ0.0317259,μ0.00390758} | 1.48908× 9 10 |
2 | {β0.770991,ζ0.0621624,γ0.036652,μ0.000819925} | 1.735 |