Epidemiological Models for Influenza and COVID-19

Robert B. Nachbar
Wolfram Research, Inc.

Table of Contents

Introduction

Background and Notes

Data and Functions

Compartmental Models

Hubei Outbreak—SEIRD

Hubei Outbreak—SEIQRD

Hubei Outbreak—SEIsaRD

Summary of Results

References

Initialization

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 (Gardner 2020; Espigule-Pons 2020). There have also been some reports of epidemiological models (Götz 2020; Peng et al. 2020; Zhou et al. 2020; Jia et al. 2020; García Moreno Esteva 2020; Antonov 2020) 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

The following subsections contain useful information about the COVID-19 infection and the outbreak in Wuhan, China; it is by no means comprehensive.

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.
See this Wikipedia article for more information.

Timeline of Events

This is a timeline of major events related to the outbreak in Hubei province.
TimelinePlotLabeled[#1,Tooltip[#2,#3]]&@@@
Events
,
options

In[]:=
Out[]=

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

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

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 Initialization section 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]
→
ℐ,ℐ
γ
→
ℛ,t//Last
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[]:=
Symbol
fitDataWDR[locale, startDate] gives data for infected, recovered, and dead relative to the given date for the onset of infection for the given locale, either an administrative division or a country. For reference, cases were first reported on 31 Dec 2019.
Out[]=
Use the function to retrieve the cases in Hubei for a plot:
ListPlotfitDataWDR
Hubei, China
ADMINISTRATIVE DIVISION
,"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
[0,∞]
or
[0,1]
, respectively, during optimization.

Compartmental Models

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=
βλ[t]
→
ℐ,ℐ
γ
→
ℛ;​​GridThread@{modelSIR,{"infection","recovery"}},

In[]:=

βλ[t]
→
ℐ
infection
ℐ
γ
→
ℛ
recovery
Out[]=
Generate a graph of the SIR model:
CompartmentalModelGraphmodelSIR,

In[]:=
Out[]=
There is code in the Initialization section 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
λ[t]
is the force of infection. It is a function of all the infectious compartments in the model (in this case, just
ℐ[t]
), and thus is time dependent. The parameter
γ
is the recovery rate constant and is defined as
γ
1
τ
, where
τ
is the average duration of infection.
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]λ[t]
′
ℐ
[t]-γℐ[t]+β[t]λ[t]
′
ℛ
[t]γℐ[t]
λ[t]ℐ[t]
Out[]=
The number of infectious individuals in the population
ℐ[t]
is the prevalence of the disease. The number of individuals that become infected per unit time is called the incidence.
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},ImageSize450]​​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]+β[t]λ[t]
′
ℛ
[t]γℐ[t]
′

[t]-β[t]λ[t]
λ[t]ℐ[t]
Out[]=
The shape for
Λ
in the schematic diagram is different because its units are
persons
time
, whereas all the other rates are
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},ImageSizeLarge]​​closeCellGroup[]
In[]:=
Out[]=

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
50×
3
10
,
100×
3
10
,
500×
3
10
,
1×
6
10
and
50×
6
10
will be tested for effective population sizes (represented as

).
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 (
    I
    0
    ). (Du et al. [2020] estimated the number to be 7.78 in Wuhan.)
    • ◼
  • 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. (2020) 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 (Cai et al. 2020). 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

    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 [Martcheva 2015]).  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:
    GridThread@{modelSEIRD,{"infection","incubation","recovery","death"}},
    
    In[]:=
    
    βλ[t]
    →
    ℰ
    infection
    ℰ
    ζ
    →
    ℐ
    incubation
    ℐ
    γ
    →
    ℛ
    recovery
    ℐ
    μ
    →
    
    death
    Out[]=
    Display the graph for the model:
    CompartmentalModelGraphmodelSEIRD,
    
    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[]:=
    ′
    
    [t]μℐ[t]
    ′
    ℰ
    [t]-ζℰ[t]+β[t]λ[t]
    ′
    ℐ
    [t]ζℰ[t]-γℐ[t]-μℐ[t]
    ′
    ℛ
    [t]γℐ[t]
    ′
    
    [t]-β[t]λ[t]
    Out[]=
    Define
    λ(t)
    for the SEIRD model:
    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

    An SEIQRD model uses an additional quarantined compartment. Quarantined individuals
    [t]
    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.
    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:
    GridThread@{modelSEIQRD,{"infection","incubation","quarantine","recovery","recovery","death","death"}},
    
    In[]:=
    
    βλ[t]
    →
    ℰ
    infection
    ℰ
    ζ
    →
    ℐ
    incubation
    ℐ
    δ
    →
    
    quarantine
    ℐ
    γ
    →
    ℛ
    recovery
    
    γ
    →
    ℛ
    recovery
    ℐ
    μ
    →
    
    death
    
    μ
    →
    
    death
    Out[]=
    Display the graph for the model:
    CompartmentalModelGraphmodelSEIQRD,
    
    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[]:=
    ′
    
    [t]μℐ[t]+μ[t]
    ′
    ℰ
    [t]-ζℰ[t]+β[t]λ[t]
    ′
    ℐ
    [t]ζℰ[t]-γℐ[t]-δℐ[t]-μℐ[t]
    ′
    
    [t]δℐ[t]-γ[t]-μ[t]
    ′
    ℛ
    [t]γℐ[t]+γ[t]
    ′
    
    [t]-β[t]λ[t]
    Out[]=
    Define
    λ(t)
    for the model:
    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

    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:
    GridThread@{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:
    CompartmentalModelGraphmodelSEIsaRD,
    
    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[]:=
    ′
    
    [t]μ
    ℐ
    s
    [t]
    ′
    ℰ
    [t]-ζ(1-ρ)ℰ[t]-ζρℰ[t]+β[t]λ[t]
    ′
    
    [t]-β[t]λ[t]
    ′
    ℐ
    a
    [t]ζ(1-ρ)ℰ[t]-γ
    ℐ
    a
    [t]
    ′
    ℐ
    s
    [t]ζρℰ[t]-γ
    ℐ
    s
    [t]-μ
    ℐ
    s
    [t]
    ′
    ℛ
    a
    [t]γ
    ℐ
    a
    [t]
    ′
    ℛ
    s
    [t]γ
    ℐ
    s
    [t]
    Out[]=
    Define
    λ(t)
    for the model:
    forceOfInfectionSEIsaRD=λ[t_]
    ℐ
    "s"
    [t]+
    ℐ
    "a"
    [t]
    -[t]
    ;
    In[]:=
    Solve for each variable:
    {susceptibleSEIsaRD,exposedSEIsaRD,infectedSymptomaticSEIsaRD,infectedAsymptomaticSEIsaRD,recoveredSymptomaticSEIsaRD,recoveredAsymptomaticSEIsaRD,deadSEIsaRD}={,ℰ,
    ℐ
    "s"
    ,
    ℐ
    "a"
    ,
    ℛ
    "s"
    ,
    ℛ
    "a"
    ,}/.ParametricNDSolve[Join[odesSEIsaRD/.forceOfInfectionSEIsaRD,{[0]-I0,ℰ[0]0,
    ℐ
    "s"
    [0]ρI0,
    ℐ
    "a"
    [0](1-ρ)I0,
    ℛ
    "s"
    [0]0,
    ℛ
    "a"
    [0]0,[0]0}],varsSEIsaRD,{t,0,500},{,I0,β,ζ,ρ,γ,μ}];
    In[]:=

    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

    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

    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
    Hubei, China
    ADMINISTRATIVE DIVISION
    ,"20 Dec 2019",
    ;
    In[]:=
    Plot the fit:
    ListPlotfitData,
    
    In[]:=
    Out[]=

    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,