This notebook is an excerpt from the set of full notebooks which can be found on the Wolfram Cloud, starting with: EpidemiologicalModelsForInfluenzaAndCOVID-19--part_1.nb

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 [JHU][JEP]. There have also been some reports of epidemiological models [TG][PYZ][ZCW][JDL][EGE][AA] 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 has to main goals. 1) It 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. As the noted statistician G. E. P. Box admonished us, “All models are wrong, but some models are useful.” It is also important know the assumptions on which the models are predicated. 2) It demonstrates the breadth of the Wolfram Language which makes these analyses relatively straight forward, from the retrieval of data from the web, to modeling and data fitting, to exposition and presentation, all in a single interactive document.
It is instructive to first examine two influenza outbreaks that occurred in 1978 in two different boarding schools. The two populations were well defined in terms of size and the assumption of rapid and uniform mixing. These two examples demonstrate that the reported retrospective data can be analyzed satisfactorily by mechanistically different models, and that limitations of the models sometimes preclude the explanation of all the observations.
Next 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 the outbreaks are to each other, point out the short comings, and make some suggestions for improvements.

Compartmental models

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 indication 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

βλ[t]
→
ℐ
infection
ℐ
γ
→
ℛ
recovery
Schematically, the model looks like this
For clarity, we will use a slightly different form of schematic diagram known as a Petri net, which looks like this
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).
′

[t]-β[t]λ[t]
′
ℐ
[t]-γℐ[t]+β[t]λ[t]
′
ℛ
[t]γℐ[t]
λ[t]ℐ[t]
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
′

[t]Λ-μ[t]-β[t]λ[t]
′
ℐ
[t]-γℐ[t]-μℐ[t]+β[t]λ[t]
′
ℛ
[t]γℐ[t]-μℛ[t]
λ[t]ℐ[t]
The shape for
Λ
in the schematic diagram is different because its units are
persons
time
, whereas all the other rates are
1
time
.
Because the duration of the influenza and COVID-19 outbreaks that are discussed here are so short (weeks or months), we can make the assumption that births and natural deaths can be ignored.
These compartmental models and their ODE representation were first developed by Kermack and McKendrick [KM] almost 100 years ago. One of the assumptions is that the population of individuals is well mixed, that is, every individual has equal likelihood to come into contact with every other individual. This assumption rarely holds in practice, but useful results can be obtained nonetheless. Other assumptions to keep in mind are that all individuals in each class are identical, and that recovery confers life-long perfect immunity.

Basic reproduction number

The rate constants i these models can be used to estimate various epidemiological parameters, such as average recover time
1
γ
. Another parameter that is often used to characterize epidemics is the basic reproduction number, or R0. It represents the number of individuals a single infected individual will infect in a completely susceptible population during her/his lifetime. If the value is less than 1, the outbreak will die out. If it is greater than 1, then the epidemic will persist. A very good description of this parameter, its uses, and what various values greater than 1 mean was recently published [EY].

Delay differential equations

While delay differential equations (DDE) are often used in population models for ecology, and many other areas as well, they can lead to unphysical artefacts. For example, shown below is the SIR model solution (solid curves) and the solution for the corresponding DDE model (dashed curves). While the total population size (
[t][t]+ℐ[t]+ℛ[t]
) is constant, the infected compartment in the DDE model (dashed red) has several excursions below 0 and the recovered compartment (dashed blue) has corresponding excursions above . Adding a scaling factor to those terms in the DDEs will only minimize and not eliminate the issue. The oscillations are also of concern.
While DDE models have been used for epidemiological models [AV], they must be used with care.

Parameter optimization

All of the parameters in these model need to remain in the range
(0,∞)
or
(0,1)
, so constraints should be used. Wolfram Language has many efficient functions for optimization with constraints, and for many systems they work well. From many years of experience with these kinds of models (epidemiologic, pharmacokinetic, and viral dynamics), a simpler passive constraint method (transformation) with unconstrained optimization works better.
Values for parameters that should fall between 0 and ∞ are log-transformed and the parameters in the ODEs are back-transformed with
Exp
. Similarly, for parameters whose values should be between 0 and 1 we use the logit transformation and its inverse. To facilitate coding, the functions
toLog
,
fromLog
,
toLogit
, and
fromLogit
are defined in the initialization section.

toLog
,
fromLog

toLog[x_] := Log[x]
fromLog[y_] := Exp[y]
fromLog[toLog[x]]x
Out[]=
True

toLogit
,
fromLogit

toLogit[x_] := Log[x/(1 - x)]
fromLogit[y_] := Exp[y]/(1 + Exp[y])
fromLogit[toLogit[x]]x//Simplify[#,Assumptionsx∈Reals]&
Out[]=
True
Influenza
There are two textbook examples of influenza outbreaks in boarding schools. They are retrospective reports because the studies were unplanned, however, from a modeling perspective they present several advantages:
◼
  • The populations were closed (no migration in or out)
    • ◼
  • There were no births or deaths, so models without demographics can be used
    • ◼
  • The populations were well mixed as the students were attending classes together
    • ◼
  • There is a daily record of the number of students with influenza

    • The data

    Outbreak 1

    Dataset

    ◼
  • From a textbook problem, Martcheva. pp. 126-139 [MM][A]
    • …, in January-February 1978, an epidemic of influenza occurred in a boarding school in the north of England. The boarding school housed a total of 763 boys, who were at risk during the epidemic. On January 22, three boys were sick.
    512 boys (67%) spent between three and seven days away from class, …
    influenzaData1=
    influenza data; Brit. Med. J. (1978)
    ;​​Keys@%
    Out[]=
    {NumberBoys,NumberConfinedToBed,NumberConvelescent,NumberNotInClass}
    DateListPlot[KeyTake[influenzaData1,{"NumberConfinedToBed","NumberConvelescent"}],FrameLabel{"date","# boys"},ImageSizeMedium,PlotMarkers{
    ,
    ,
    ,
    ,
    }]//Rasterize
    Out[]=

    Fitting data

    firstDate1=First@First@influenzaData1["NumberConfinedToBed"]​​lastDate1=First@Last@influenzaData1["NumberConfinedToBed"]
    Out[]=
    Day: Sun 22 Jan 1978
    Out[]=
    Day: Sat 4 Feb 1978
    fitData1={toModelTime[firstDate1,#1],#2}&@@@influenzaData1["NumberConfinedToBed"]//Echo[#,"",Multicolumn]&;
    »
    {0,0}
    {4,223}
    {8,191}
    {12,12}
    {1,6}
    {5,294}
    {9,124}
    {13,4}
    {2,26}
    {6,258}
    {10,68}
    ​
    {3,74}
    {7,237}
    {11,27}
    ​
    populationSize1=influenzaData1["NumberBoys"]
    Out[]=
    763

    Outbreak 2

    Dataset

    ◼
  • From another textbook problem, Martcheva, pp. 144-145 [MM][HJR]
    • The West Country English Boarding School housed 578 boys. An epidemic of influenza began on 15 January 1978.
    One hundred and sixty-six boys (29 per cent) were treated in the sick bays by bed rest and aspirin. Certainly this is not the total number with influenza, as the older boys had their own rooms in which some remained, treating themselves and avoiding detection and supervision as a result of the dislocated curriculum.
    influenzaData2=
    influenza data; J. Royal College Gen. Practitioners (1980)
    ;​​Keys@%
    Out[]=
    {NumberBoys,NumberWithInfluenza,NumberInSickBays}
    DateListPlot[influenzaData2["NumberWithInfluenza"],FrameLabel{"date","number of boys with influenza"},AspectRatio66/156,ImageSizeLarge,PlotMarkers{
    ,
    ,
    ,
    ,
    }]//Rasterize
    Out[]=

    Fitting data

    firstDate2=First@First@influenzaData2["NumberWithInfluenza"]//DatePlus[#,-1]&​​lastDate2=First@Last@influenzaData2["NumberWithInfluenza"]
    Out[]=
    Day: Sun 15 Jan 1978
    Out[]=
    Day: Wed 15 Feb 1978
    fitData2={toModelTime[firstDate2,#1],#2}&@@@influenzaData2["NumberWithInfluenza"]//Echo[#,"",Multicolumn]&;
    »
    {1,2}
    {6,15}
    {12,55}
    {18,30}
    {23,11}
    {28,4}
    {2,5}
    {8,31}
    {13,58}
    {19,23}
    {24,12}
    {29,2}
    {3,10}
    {9,42}
    {15,52}
    {20,19}
    {25,9}
    {30,1}
    {4,12}
    {10,45}
    {16,42}
    {21,13}
    {26,7}
    {31,1}
    {5,14}
    {11,53}
    {17,40}
    {22,14}
    {27,5}
    ​
    populationSize2=influenzaData2["NumberBoys"]
    Out[]=
    578

    SEIR Model

    ◼
  • SEIR: SIR model with exposed compartment
    • ◼
  • No birth or natural mortality due to short time span of coverage of model
    • ◼
  • Mass action incidence
    • ◼
  • These are the transitions:
    • 
    βλ[t]
    →
    ℰ
    infection
    ℰ
    ζ
    →
    ℐ
    incubation
    ℐ
    γ
    →
    ℛ
    recovery
    ◼
  • These are the ODEs
    • forceOfInfectionSEIR={λ[t_]ℐ[t]};​​{varsSEIR,odesSEIR}=Values@KineticCompartmentalModel
    βλ[t]
    →
    ℰ,ℰ
    ζ
    →
    ℐ,ℐ
    γ
    →
    ℛ,t{1,-1};​​Column[odesSEIR]
    Out[]=
    ′
    ℰ
    [t]-ζℰ[t]+β[t]λ[t]
    ′
    ℐ
    [t]ζℰ[t]-γℐ[t]
    ′
    ℛ
    [t]γℐ[t]
    ′
    
    [t]-β[t]λ[t]
    In[]:=
    {susceptibleSEIR,infectedSEIR,exposedSEIR,recoveredSEIR}={,ℐ,ℰ,ℛ}/.ParametricNDSolve[Join[odesSEIR/.forceOfInfectionSEIR,{[0]-I0,ℰ[0]0,ℐ[0]I0,ℛ[0]0}],Head/@varsSEIR,{t,0,100},{,I0,β,ζ,γ}];

    First dataset

    Manual fit

    Out[]=
    ​
    β = 0.0123027 (infection)
    ζ = 0.635387 (incubation)
    γ = 0.456576 (recovery)
    t
    max
    y
    max
    
    ℰ
    ℐ
    ℛ

    Statistical fit

    fit1=fit=NonlinearModelFit[fitData1,infectedSEIR[populationSize1,1,β,ζ,γ][t],{{β,0.012},{ζ,0.61},{γ,0.43}},t]
    Out[]=
    FittedModel
    InterpolatingFunction
    Domain: {{0.,100.}}
    Output: scalar
    [t]
    
    fit["EstimatedVariance"]
    Out[]=
    291.969
    fit["ParameterTable"]
    Out[]=
    Estimate
    Standard Error
    t-Statistic
    P-Value
    β
    0.0104591
    0.00133223
    7.85086
    7.80853×
    -6
    10
    ζ
    0.677183
    0.0882803
    7.67082
    9.71785×
    -6
    10
    γ
    0.474842
    0.0184577
    25.7259
    3.53682×
    -11
    10
    fitWithDataPlot[fit,{firstDate1,lastDate1}]
    Out[]=
    residualsPlot[fit]
    Out[]=

    Sensitivity analysis

    modelSensitivityPlot[infectedSEIR[populationSize1,1,β,ζ,γ],fit,{{0,14},{-25,400}},{0.004,0.25,0.20}]
    Out[]=
    β
    ζ
    γ
    observed
    calculated
    negative sensitivity
    positive sensitivity

    Total incidence & epidemic size

    incidence=NIntegrate[ζexposedSEIR[populationSize1,1,β,ζ,γ][t]/.fit["BestFitParameters"],{t,0,100}]
    Out[]=
    762.
    incidence
    populationSize1
    Out[]=
    0.998689
    ◼
  • Nearly every boy is predicted to become infected, whereas only 512 boys (67%) were reported to have been confined to bed
    • ◼
  • Graphically, we can plot the epidemic size
    -[t]
    to show this
    • epidemicSizePlot[susceptibleSEIR[populationSize1,1,β,ζ,γ],fit,populationSize1,influenzaData1["NumberNotInClass"],{{0,15},{0,800}}]
    Out[]=

    Second dataset

    [see full notebook for elided content]

    SEIQR Model

    Since the data report the “number of boys confined to bed”, they are actually quarantined or isolated and not freely circulating among the student body
    ◼
  • SEIQR: SEIR model with quarantined compartment
    • ◼
  • No birth or natural mortality due to short time span of coverage of model
    • ◼
  • Mass action incidence
    • ◼
  • These are the transitions:
    • Out[]=
    
    βλ[t]
    →
    ℰ
    infection
    ℰ
    ζ
    →
    ℐ
    incubation
    ℐ
    δ
    →
    
    quarantine
    ℐ
    γ
    →
    ℛ
    recovery
    
    γ
    →
    ℛ
    recovery
    ◼
  • These are the ODEs
    • forceOfInfectionSEIQR={λ[t_]ℐ[t]};​​{varsSEIQR,odesSEIQR}=Values@KineticCompartmentalModel
    βλ[t]
    →
    ℰ,ℰ
    ζ
    →
    ℐ,ℐ
    δ
    →
    ,ℐ
    γ
    →
    ℛ,
    γ
    →
    ℛ,t{1,-1};​​Column[odesSEIQR]
    Out[]=
    ′
    ℰ
    [t]-ζℰ[t]+β[t]λ[t]
    ′
    ℐ
    [t]ζℰ[t]-γℐ[t]-δℐ[t]
    ′
    
    [t]δℐ[t]-γ[t]
    ′
    ℛ
    [t]γℐ[t]+γ[t]
    ′
    
    [t]-β[t]λ[t]
    In[]:=
    {susceptibleSEIQR,exposedSEIQR,infectedSEIQR,quarantinedSEIQR,recoveredSEIQR}={,ℰ,ℐ,,ℛ}/.ParametricNDSolve[Join[odesSEIQR/.forceOfInfectionSEIQR,{[0]-I0,ℰ[0]0,ℐ[0]I0,[0]0,ℛ[0]0}],Head/@varsSEIQR,{t,0,100},{,I0,β,ζ,γ,δ}];

    First dataset

    [see full notebook for elided content]

    Second dataset

    Manual fit

    Out[]=
    ​
    β = 0.009 (infection)
    ζ = 5.76 (incubation)
    γ = 0.21 (recovery)
    δ = 4.2 (quarantine)
    t
    max
    y
    max
    scale
    ℰ
    ℐ
    
    scaleℛ

    Statistical fit

    fit4=fit=NonlinearModelFit[fitData2,quarantinedSEIQR[populationSize2,1,β,ζ,γ,δ][t],{{β,0.0090},{ζ,3.90},{γ,0.21},{δ,4.2}},t]
    Out[]=
    FittedModel
    InterpolatingFunction
    Domain: {{0.,100.}}
    Output: scalar
    [t]
    
    fit["EstimatedVariance"]
    Out[]=
    5.31158
    fit["ParameterTable"]
    Out[]=
    Estimate
    Standard Error
    t-Statistic
    P-Value
    β
    0.0124396
    0.0310227
    0.400983
    0.691839
    ζ
    1.96899
    2.5298
    0.778318
    0.443686
    γ
    0.317253
    0.098527
    3.21996
    0.00353842
    δ
    5.30268
    13.6043
    0.389781
    0.7
    fitWithDataPlot[fit,{firstDate2,lastDate2}]
    Out[]=
    residualsPlot[fit]
    Out[]=

    Sensitivity analysis

    modelSensitivityPlot[quarantinedSEIQR[populationSize2,1,β,ζ,γ,δ],fit,{{0,32},{-5,80}},{0.0005,0.6,0.08,0.25}]
    Out[]=
    β
    ζ
    γ
    δ
    observed
    calculated
    negative sensitivity
    positive sensitivity

    Total incidence & epidemic size

    incidence=NIntegrate[ζexposedSEIQR[populationSize2,1,β,ζ,γ,δ][t]/.fit["BestFitParameters"],{t,0,100}]
    Out[]=
    234.078
    incidence
    populationSize2
    Out[]=
    0.40498
    epidemicSizePlot[susceptibleSEIQR[populationSize2,1,β,ζ,γ,δ],fit,populationSize2,influenzaData2["NumberInSickBays"],{{0,40},{0,600}}]
    Out[]=
    ◼
  • Again, overestimating the epidemic size, but much improved over the model without quarantine

    • Summary

    Parameter values and statistics

    In[]:=
    parameters=Keys@#["BestFitParameters"]&/@{fit1,fit2,fit3,fit4};
    Out[]=
    Parameter Tables
    SEIR
    SEIQR
    dataset 1
    Estimate
    Standard Error
    t-Statistic
    P-Value
    β
    0.0104591
    0.00133223
    7.85086
    7.80853×
    -6
    10
    ζ
    0.677183
    0.0882803
    7.67082
    9.71785×
    -6
    10
    γ
    0.474842
    0.0184577
    25.7259
    3.53682×
    -11
    10
    Estimate
    Standard Error
    t-Statistic
    P-Value
    β
    0.0270749
    0.0320636
    0.844411
    0.418183
    ζ
    1.02765
    0.880404
    1.16725
    0.270189
    γ
    0.421429
    0.0214521
    19.6451
    2.55636×
    -9
    10
    δ
    5.61746
    14.3029
    0.392749
    0.702746
    dataset 2
    Estimate
    Standard Error
    t-Statistic
    P-Value
    β
    0.0102136
    0.000826261
    12.3612
    2.16619×
    -12
    10
    ζ
    0.190368
    0.0182367
    10.4387
    8.60582×
    -11
    10
    γ
    0.844692
    0.0274663
    30.7538
    5.64834×
    -22
    10
    Estimate
    Standard Error
    t-Statistic
    P-Value
    β
    0.0124396
    0.0310227
    0.400983
    0.691839
    ζ
    1.96899
    2.5298
    0.778318
    0.443686
    γ
    0.317253
    0.098527
    3.21996
    0.00353842
    δ
    5.30268
    13.6043
    0.389781
    0.7
    Out[]=
    Parameter Correlation Matrices
    SEIR
    SEIQR
    dataset 1
    β
    ζ
    γ
    β
    1.
    -0.962581
    -0.330638
    ζ
    -0.962581
    1.
    0.313405
    γ
    -0.330638
    0.313405
    1.
    β
    ζ
    γ
    δ
    β
    1.
    0.702698
    0.0289383
    0.962664
    ζ
    0.702698
    1.
    -0.470858
    0.868825
    γ
    0.0289383
    -0.470858
    1.
    -0.157721
    δ
    0.962664
    0.868825
    -0.157721
    1.
    dataset 2
    β
    ζ
    γ
    β
    1.
    -0.957492
    -0.446813
    ζ
    -0.957492
    1.
    0.615198
    γ
    -0.446813
    0.615198
    1.
    β
    ζ
    γ
    δ
    β
    1.
    -0.998185
    0.97948
    0.999961
    ζ
    -0.998185
    1.
    -0.989368
    -0.997621
    γ
    0.97948
    -0.989368
    1.
    0.977696
    δ
    0.999961
    -0.997621
    0.977696
    1.

    Goodness of fit

    Out[]=
    Model
    2
    R
    (higher is better)
    SEIR
    SEIQR
    dataset 1
    0.99004
    0.989476
    dataset 2
    0.983136
    0.994303
    Out[]=
    Model AIC (lower is better)
    SEIR
    SEIQR
    dataset 1
    124.203
    127.309
    dataset 2
    166.057
    136.725
    Out[]=
    Model likelihood (higher is better)
    SEIR
    SEIQR
    dataset 1
    1.
    0.2116
    dataset 2
    4.27227×
    -7
    10
    1.
    ◼
  • For the first influenza dataset, the SEIR model is somewhat better than the SEIQR model
    • ◼
  • For the second influenza dataset, the SEIQR model is much better than the SEIR model

    • Epidemiological parameters

    Out[]=
    Average incubation period
    1
    ζ
    (d)
    SEIR
    SEIQR
    dataset 1
    1.47671
    0.973095
    dataset 2
    5.25297
    0.507874
    Out[]=
    Average duration of infection
    1
    γ
    (d)
    SEIR
    SEIQR
    dataset 1
    2.10597
    2.37288
    dataset 2
    1.18386
    3.15206
    Out[]=
    Average duration of circulation
    1
    δ
    (h)
    SEIQR
    dataset 1
    4.27239
    dataset 2
    4.52601
    ◼
  • The basic reproduction number
    R
    0
    was computed from formulas in Brauer [FB]:
    • ◼
  • SEIR model
    • Out[]=
    Basic reproduction number (​
    R
    0
    ​)
    SEIR
    dataset 1
    16.8063
    dataset 2
    6.98887
    ◼
  • SEIQR model
    • Out[]=
    Basic reproduction number ​
    R
    γ
    ​
    SEIQR
    dataset 1
    8.72622
    dataset 2
    4.27398
    ◼
  • Biggerstaff et al. [BCR] report much smaller values for the basic reproduction number for influenza:
    • Out[]=
    Outbreak
    n
    Median
    R
    0
    1918 pandemic
    51
    1.8
    1957 pandemic
    7
    1.65
    1968 pandemic
    7
    1.8
    2009 pandemic
    78
    1.46
    seasonal influenza
    47
    1.28

    General observations & conclusions

    ◼
  • Retrospective data without full details of collection can be interpreted in different ways
    • ◼
  • Recapitulating heterogenous data (e.g., prevalence and epidemic size) can be difficult
    • ◼
  • Different models may do better for different datasets for the same disease

    • Suggestions for follow-up and improvements

    ◼
  • Recent reading [MM][FT] suggests that standard incidence should be used for SEIQR models
    • ◼
  • Investigate the effect on model fitting
    • ◼
  • investigate the effect on the basic reproduction number
    • ◼
  • Why is
    R
    0
    computed to be so much higher than reported?
    • ◼
  • Try Jacobian method [MM]
    • ◼
  • Try next generation method [MM]
    • COVID-19

    Epidemic data from Wolfram Data Repository (WDR)

    ◼
  • The Johns Hopkins University data are available from the Wolfram Data Repository, and are described at https://www.wolframcloud.com/obj/resourcesystem/published/DataRepository/resources/Epidemic-Data-for-Novel-Coronavirus-COVID-19
    • ResourceObject["Epidemic Data for Novel Coronavirus COVID-19"]
    Out[]=
    ResourceObject
    Name: Epidemic Data for Novel Coronavirus COVID-19
    »
    Type: DataResource
    Description:
    Estimated cases of novel coronavirus (COVID-19, formerly known as 2019-nCoV) inf...
    
    ◼
  • The function
    fitDataWDR
    (defined in the Initialization section at the end of the notebook) provides convenient access to the data for modeling
    • ◼
  • Note: a correction to the confirmed cases data for Hubei is applied [BK]
    • ?fitDataWDR
    Out[]=
    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.
    Options@fitDataWDR
    Out[]=
    {makeCorrectionForHubeiFalse,dateRangeAll}
    ListPlotfitDataWDR
    Anhui, China
    ADMINISTRATIVE DIVISION
    ,"1 Jan 2020",PlotLabel"Anhui",FrameLabel{"time (d)","# individuals"},PlotLegends{"confirmed"-("recovered"+"dead"),"recovered","dead"}//Rasterize
    Out[]=
    ListPlotfitDataWDR
    China
    COUNTRY
    ,"1 Jan 2020",PlotLabel"Mainland China",FrameLabel{"time (d)","# individuals"},PlotLegends{"confirmed"-("recovered"+"dead"),"recovered","dead"}//Rasterize
    Out[]=
    ◼
  • Get the population size for the top three provinces
    • WolframAlpha["population Hubei China","Result"]
    Out[]=
    58160000
    people
    WolframAlpha["population Guangdong China","Result"]
    Out[]=
    104303132
    people
    WolframAlpha["population Henan China","Result"]
    Out[]=
    94023567
    people

    Models

    ◼
  • We’ll explore a number of models, hoping to find one with a reasonable structure that recapitulates the reported data
    • ◼
  • SEIRD: SEIR model with dead compartment
    • ◼
  • SEIQRD: SEIQR model with dead compartment
    • ◼
  • SEIsaRD: SEIR model with symptomatic and asymptomatic infectious compartments
    • ◼
  • We’ll explore all the the models for one location, then apply the “best” model to the other locations
    • ◼
  • We’ll examine alternative models with full and effective population sizes
    • ◼
  • Full population sizes will come from Wolfram|Alpha
    • ◼
  • Effective population size
    
    to be examined:
    50×
    3
    10
    ,
    100×
    3
    10
    ,
    500×
    3
    10
    ,
    1×
    6
    10
    ,
    50×
    6
    10
    • ◼
  • We’ll examine alternative models with mass action or standard incidence
    • ◼
  • General assumptions:
    • ◼
  • No birth or natural mortality due to short time span of coverage of model
    • ◼
  • Death due to disease will be included
    • ◼
  • There is no migration into or out of the population
    • ◼
  • We don’t know with certainty when the first infection occurred
    • ◼
  • Du et al. assumed 1 Dec 2019 [DWC]
    • ◼
  • We will examine datasets assuming the first infections occurred on 20 Dec 2019, and 30 Dec 2019 (when the first reported infections)
    • ◼
  • The initial number of infected individuals
    I0
    is
    10
    ; Du et al. estimated the number to be 7.78 in Wuhan [DWC]
    • ◼
  • 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 stratification by age and gender
    • ◼
  • Under reporting of cases is a serious concern, primarily due to the infection being asymptomatic in many individuals, especially children [CXL]; we will assume that all cases are reported
    • ◼
  • This work was carried out over a number of days, with the available WDR data keeping pace. However to assure a consistent analysis and interpretation of results we’ll work with data for the period 21 Jan 2020 to 26 Feb 2020.
    • ◼
  • More recent data can be used to test the models.

    • Hubei outbreak—SEIRD

    SEIRD Model with standard incidence

    ◼
  • SEIRD: SEIR model with dead compartment
    • ◼
  • 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; this is in contrast to models for influenza and SARS [MM]
    • ◼
  •  may be adjusted for individual locations; initially we'll use
    100000
    individuals
    • ◼
  • These are the transitions:
    • Out[]=
    
    βλ[t]
    →
    ℰ
    infection
    ℰ
    ζ
    →
    ℐ
    incubation
    ℐ
    γ
    →
    ℛ
    recovery
    ℐ
    μ
    →
    
    death
    ◼
  • These are the ODEs
    • forceOfInfectionSEIRD=λ[t_]
    ℐ[t]
    -[t]
    ;​​{varsSEIRD,odesSEIRD}=Values@KineticCompartmentalModel
    βλ[t]
    →
    ℰ,ℰ
    ζ
    →
    ℐ,ℐ
    γ
    →
    ℛ,ℐ
    μ
    →
    ,t{1,-1};​​Column[odesSEIRD]
    Out[]=
    ′
    
    [t]μℐ[t]
    ′
    ℰ
    [t]-ζℰ[t]+β[t]λ[t]
    ′
    ℐ
    [t]ζℰ[t]-γℐ[t]-μℐ[t]
    ′
    ℛ
    [t]γℐ[t]
    ′
    
    [t]-β[t]λ[t]
    In[]:=
    {susceptibleSEIRD,exposedSEIRD,infectedSEIRD,recoveredSEIRD,deadSEIRD}={,ℰ,ℐ,ℛ,}/.ParametricNDSolve[Join[odesSEIRD/.forceOfInfectionSEIRD,{[0]-I0,ℰ[0]0,ℐ[0]I0,ℛ[0]0,[0]0}],Head/@varsSEIRD,{t,0,500},{,I0,β,ζ,γ,μ}];

    outbreak start 20 Dec 2019

    Fitting data

    t0=DateObject["20 Dec 2019"]
    Out[]=
    Day: Fri 20 Dec 2019
    ◼
  • We use the option
    "makeCorrectionForHubei"
    remove the jump in confirmed cases in Hubei that occurred when the testing method changed, which at the time seemed like a reasonable data preprocessing step. Unfortunately, for the data after about 15 March 2020 this practice leads to an artefact that gives negative values.
    • fitData1=fitData=fitDataWDR
    Hubei, China
    ADMINISTRATIVE DIVISION
    ,t0,"makeCorrectionForHubei"True,"dateRange"{All,"26 Feb 2020"};​​Dimensions/@%
    Out[]=
    {{36,2},{36,2},{36,2}}
    ListPlot[fitData,PlotLabel"Hubei",FrameLabel{"time (d)","# individuals"},PlotLegends{"confirmed"-("recovered"+"dead"),"recovered","dead"}]//Rasterize
    Out[]=

    Manual fit

    Out[]=
    ​
    effective population size
    50000
    100000
    500000
    1000000
    5000000
    β = 0.328095 (infection)
    ζ = 0.6 (incubation)
    γ = 0.072 (recovery)
    μ = 0.0028 (death)
    t
    max
    y
    max
    0.5
    ℰ
    ℐ
    ℛ
    

    Optimization—
    100000

    fitErrorSEIRD[fitData,100000,10,logBeta,logZeta,logGamma,logMu]/.Thread[{logBeta,logZeta,logGamma,logMu}Log@{0.32,0.60,0.072,0.0028}]
    Out[]=
    1.48125×
    10
    10
    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
    Out[]=
    {535.716,Null}
    Take[SortBy[Last]@randomSampleWithError,15]//tfn
    Out[]//TableForm=
    1
    2
    1
    {β0.783142,ζ0.0617432,γ0.0230569,μ0.00744279}
    1.60102×
    9
    10
    2
    {β0.80396,ζ0.0608227,γ0.0325909,μ0.00106566}
    1.69269×
    9
    10
    3
    {β0.68954,ζ0.0715426,γ0.0328715,μ0.00365816}
    1.74281×
    9
    10
    4
    {β0.781349,ζ0.0605596,γ0.0261109,μ0.0173146}
    2.01264×
    9
    10
    5
    {β0.662634,ζ0.0703356,γ0.0341459,μ0.00139895}
    2.0374×
    9
    10
    6
    {β0.754485,ζ0.0674952,γ0.0360232,μ0.000922667}
    2.03869×
    9
    10
    7
    {β0.65945,ζ0.0699785,γ0.0337624,μ0.00174329}
    2.08808×
    9
    10
    8
    {β0.824407,ζ0.0690598,γ0.0351951,μ0.0103829}
    2.14372×
    9
    10
    9
    {β0.734833,ζ0.0683832,γ0.0203583,μ0.00999414}
    2.14916×
    9
    10
    10
    {β0.778671,ζ0.0628293,γ0.0281187,μ0.000346511}
    2.15915×
    9
    10
    11
    {β0.791376,ζ0.0618965,γ0.0257106,μ0.00152217}
    2.18647×
    9
    10
    12
    {β0.733267,ζ0.0722545,γ0.0462848,μ0.0089711}
    2.2431×
    9
    10
    13
    {β0.690317,ζ0.0657344,γ0.0366964,μ0.00152161}
    2.24614×
    9
    10
    14
    {β0.524454,ζ0.108293,γ0.0281095,μ0.0143855}
    2.29754×
    9
    10
    15
    {β0.794401,ζ0.0703482,γ0.044285,μ0.0175903}
    2.3318×
    9
    10
    Cases[randomSampleWithError,{_,e_}/;Log10[e]<10.3e]//Histogram
    Out[]=
    Cases[randomSampleWithError,{p_,e_}/;Log10[e]<9.75Append[Take[Values@p,2],Log10[e]]]//ListPointPlot3D[MapThread[Tooltip,{#,Range@Length@#}]]&
    Out[]=
    (opts=NMinimize[fitErrorSEIRD[fitData,100000,10,logBeta,logZeta,logGamma,logMu],MapThread[{#1,#2-0.1,#2+0.1}&,{{logBeta,logZeta,logGamma,logMu},Log@#}],Method{"NelderMead","PostProcess"False},MaxIterations500]&/@Values/@randomSampleWithError〚{1,48,38,22,27},1〛)//Timing
    Out[]=
    31.7781,6.71022×
    8
    10
    ,{logBeta0.371418,logZeta-3.41883,logGamma-3.78253,logMu-5.27678},2.08648×
    9
    10
    ,{logBeta-1.60938,logZeta6.60015,logGamma-3.33887,logMu-4.45224},2.0863×
    9
    10
    ,{logBeta-1.61179,logZeta22.7617,logGamma-3.33755,logMu-4.49022},2.08587×
    9
    10
    ,{logBeta-1.61005,logZeta10.3016,logGamma-3.33896,logMu-4.45934},2.08585×
    9
    10
    ,{logBeta-1.61006,logZeta37.4288,logGamma-3.33893,logMu-4.45944}
    In[]:=
    opts={{6.710218035671642`*^8,{logBeta0.37141796551238854`,logZeta-3.4188303507765982`,logGamma-3.782525253430217`,logMu-5.2767786841232915`}},{2.0864811216891184`*^9,{logBeta-1.6093779748474608`,logZeta6.600152877185301`,logGamma-3.338867980122241`,logMu-4.452243377967799`}},{2.0862993932778804`*^9,{logBeta-1.6117883830332087`,logZeta22.761748916252788`,logGamma-3.3375455082279317`,logMu-4.490216567952472`}},{2.085865955403772`*^9,{logBeta-1.6100471264040541`,logZeta10.301614036432737`,logGamma-3.3389626645594963`,logMu-4.4593361850039095`}},{2.0858508310300074`*^9,{logBeta-1.6100566406571377`,logZeta37.42881230452504`,logGamma-3.338934608497214`,logMu-4.459435844887798`}}};
    Out[]=
    (fitParams=Thread[{β,ζ,γ,μ}#]&/@(Last/*Values/*Exp)/@opts)//Column
    Out[]=
    {β1.44979,ζ0.0327507,γ0.0227651,μ0.00510886}
    {β0.200012,ζ735.208,γ0.0354771,μ0.0116524}
    β0.19953,ζ7.67895×
    9
    10
    ,γ0.035524,μ0.0112182
    {β0.199878,ζ29780.6,γ0.0354737,μ0.01157}
    β0.199876,ζ1.7994×
    16
    10
    ,γ0.0354747,μ0.0115689
    ◼
  • The table below shows the estimates for the epidemiological parameters
    • ◼
  • the average duration of infection takes into account both recovery and death
    • Out[]=
    average incubation
    period (d)
    average duration of
    infection (d)
    average recovery
    time (d)
    average life span
    while infected (d)
    30.5337
    35.8757
    43.9268
    195.738
    0.00136016
    21.2181
    28.1872
    85.8193
    1.30226×
    -10
    10
    21.3939
    28.1499
    89.1407
    0.0000335789
    21.2568
    28.1899
    86.4301
    5.55742×
    -17
    10
    21.2569
    28.1891
    86.4387
    ◼
  • The last four fits have unrealistically short incubation periods (
    1
    ζ
    )

    • Optimization—
    50000

    fitErrorSEIRD[fitData,50000,10,logBeta,logZeta,logGamma,logMu]/.Thread[{logBeta,logZeta,logGamma,logMu}Log@{0.25,0.81,0.025,0.0028}]
    Out[]=
    7.55454×
    8
    10
    randomSampleWithError=With[{delta=Log[10.]{-1,+1},init=Log@{0.32,0.60,0.072,0.0028}},{Thread[{β,ζ,γ,μ}Exp@{##}],fitErrorSEIRD[fitData,50000,10,##]}&@@@Prepend[Transpose[RandomReal[#+delta,50000]&/@init],init]​​]//SortBy[#,Last]&;//Timing
    Out[]=
    {544.733,Null}
    Take[randomSampleWithError,15]//tfn
    Out[]//TableForm=
    1
    2
    1
    {β0.298171,ζ0.345907,γ0.0200451,μ0.00148306}
    4.23445×
    8
    10
    2
    {β0.270818,ζ0.457373,γ0.0184726,μ0.000699187}
    4.66637×
    8
    10
    3
    {β0.288063,ζ0.378503,γ0.0206175,μ0.00339612}
    4.71975×
    8
    10
    4
    {β0.278524,ζ0.446642,γ0.018432,μ0.00171766}
    4.74215×
    8
    10
    5
    {β0.30724,ζ0.351347,γ0.0225475,μ0.00164229}
    5.02901×
    8
    10
    6
    {β0.258096,ζ0.49521,γ0.0188345,μ0.000323483}
    5.07881×
    8
    10
    7
    {β0.327702,ζ0.257677,γ0.0178705,μ0.000443775}
    5.07981×
    8
    10
    8
    {β0.29057,ζ0.384673,γ0.0185249,μ0.00638532}
    5.19525×
    8
    10
    9
    {β0.448074,ζ0.157858,γ0.0188528,μ0.00434991}
    5.41922×
    8
    10
    10
    {β0.2721,ζ0.428386,γ0.0167909,μ0.000365572}
    5.42847×
    8
    10
    11
    {β0.212864,ζ1.95119,γ0.0185973,μ0.00103472}
    5.50506×
    8
    10
    12
    {β0.222443,ζ1.33832,γ0.016305,μ0.00463509}
    5.55485×
    8
    10
    13
    {β0.466086,ζ0.1453,γ0.0190575,μ0.000498626}
    5.56683×
    8
    10
    14
    {β0.296613,ζ0.337273,γ0.0215859,μ0.00302788}
    5.61315×
    8
    10
    15
    {β0.33093,ζ0.247425,γ0.0165961,μ0.000523316}
    5.67203×
    8
    10
    Cases[randomSampleWithError,{_,e_}/;Log10[e]<10.3e]//Histogram
    Out[]=
    Cases[randomSampleWithError,{p_,e_}/;Log10[e]<9.75Append[Take[Values@p,2],Log10[e]]]//ListPointPlot3D[MapThread[Tooltip,{#,Range@Length@#}]]&
    Out[]=
    (opts=NMinimize[fitErrorSEIRD[fitData,50000,10,logBeta,logZeta,logGamma,logMu],MapThread[{#1,#2-0.1,#2+0.1}&,{{logBeta,logZeta,logGamma,logMu},Log@#}],Method{"NelderMead","PostProcess"False},MaxIterations500]&/@Values/@randomSampleWithError〚{1,88,12,17,19,51},1〛)//Timing
    Out[]=
    26.8058,4.12079×
    8
    10
    ,{logBeta-1.19747,logZeta-1.06882,logGamma-3.93667,logMu-5.92011},4.12079×
    8
    10
    ,{logBeta-1.19747,logZeta-1.0688,logGamma-3.93668,logMu-5.92001},4.12079×
    8
    10
    ,{logBeta-1.1974,logZeta-1.06899,logGamma-3.93667,logMu-5.92009},4.41355×
    8
    10
    ,{logBeta-1.20461,logZeta-1.07757,logGamma-3.89677,logMu-7.49488},4.12084×
    8
    10
    ,{logBeta-1.19732,logZeta-1.06887,logGamma-3.93702,logMu-5.9105},4.12079×
    8
    10
    ,{logBeta-1.19747,logZeta-1.06882,logGamma-3.93668,logMu-5.92009}
    In[]:=
    opts={{4.1207928762936115`*^8,{logBeta-1.1974666762663535`,logZeta-1.068822223739908`,logGamma-3.936672407434169`,logMu-5.920107052365711`}},{4.120792884212842`*^8,{logBeta-1.1974729141480858`,logZeta-1.0687976562149881`,logGamma-3.9366814230483707`,logMu-5.920011783270811`}},{4.1207929010304075`*^8,{logBeta-1.1974033926543777`,logZeta-1.0689916272672186`,logGamma-3.9366739559771093`,logMu-5.920091291282083`}},{4.4135482987447315`*^8,{logBeta-1.2046113419789002`,logZeta-1.0775655387259848`,logGamma-3.8967679093486955`,logMu-7.494876420961075`}},{4.120836777009539`*^8,{logBeta-1.1973209653024746`,logZeta-1.068868272540231`,logGamma-3.937021469031884`,logMu-5.910498943464927`}},{4.120792876819955`*^8,{logBeta-1.1974665599902918`,logZeta-1.0688222836071968`,logGamma-3.9366789629760923`,logMu-5.920091970674995`}}};
    Out[]=
    (fitParams=Thread[{β,ζ,γ,μ}#]&/@(Last/*Values/*Exp)/@opts)//Column
    Out[]=
    {β0.301958,ζ0.343413,γ0.019513,μ0.00268491}
    {β0.301956,ζ0.343421,γ0.0195129,μ0.00268517}
    {β0.301977,ζ0.343355,γ0.019513,μ0.00268496}
    {β0.299808,ζ0.340423,γ0.0203074,μ0.000555925}
    {β0.302002,ζ0.343397,γ0.0195062,μ0.00271083}
    {β0.301958,ζ0.343413,γ0.0195129,μ0.00268495}
    ◼
  • The table below shows the estimates for the epidemiological parameters
    • ◼
  • the average duration of infection takes into account both recovery and death
    • Out[]=
    average incubation
    period (d)
    average duration of
    infection (d)
    average recovery
    time (d)
    average life span
    while infected (d)
    2.91195
    45.0492
    51.2478
    372.452
    2.91188
    45.049
    51.2482
    372.416
    2.91244
    45.0492
    51.2479
    372.446
    2.93752
    47.9309
    49.243
    1798.8
    2.91208
    45.0105
    51.2657
    368.89
    2.91195
    45.0494
    51.2481
    372.446
    ◼
  • The first three and last two fits are essentially the same
    • ◼
  • The fourth fit under-predicts the deaths

    • outbreak start 30 Dec 2019

    Fitting data

    t0=DateObject["30 Dec 2019"]
    Out[]=
    Day: Mon 30 Dec 2019
    fitData2=fitData=fitDataWDR
    Hubei, China
    ADMINISTRATIVE DIVISION
    ,t0,"makeCorrectionForHubei"True,"dateRange"{All,"26 Feb 2020"};​​Dimensions/@%
    Out[]=
    {{36,2},{36,2},{36,2}}
    ListPlot[fitData,PlotLabel"Hubei",FrameLabel{"time (d)","# individuals"},PlotLegends{"confirmed"-("recovered"+"dead"),"recovered","dead"}]//Rasterize
    Out[]=
    [see full notebook for elided content]

    Comparisons

    Out[]//TableForm=
    startDate
    populationSize
    β
    ζ
    γ
    μ
    SSE
    20 Dec 2019
    50000
    0.301958
    0.343413
    0.019513
    0.00268491
    4.12079×
    8
    10
    20 Dec 2019
    100000
    1.44979
    0.0327507
    0.0227651
    0.00510886
    6.71022×
    8
    10
    30 Dec 2019
    50000
    0.380397
    0.412503
    0.0215062
    0.00353125
    2.01761×
    8
    10
    30 Dec 2019
    100000
    2.32889
    0.0323497
    0.0237666
    0.00537906
    5.24727×
    8
    10
    Out[]=
    ◼
  • increasing the population size makes ℐ peak broader
    • ◼
  • shifting t0 to a later date (closer to the data) does not have a noticeable effect
    • ◼
  • in all cases, the susceptible compartment goes to 0 and the epidemic size is 100% of the population

    • Hubei outbreak—SEIQRD

    SEIQRD Model with standard incidence

    ◼
  • SEIQRD: SEIQR model with dead 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
    • ◼
  • These are the transitions:
    • Out[]=
    
    βλ[t]
    →
    ℰ
    infection
    ℰ
    ζ
    →
    ℐ
    incubation
    ℐ
    δ
    →
    
    quarantine
    ℐ
    γ
    →
    ℛ
    recovery
    
    γ
    →
    ℛ
    recovery
    ℐ
    μ
    →
    
    death
    
    μ
    →
    
    death
    ◼
  • These are the ODEs
    • forceOfInfectionSEIQRD=λ[t_]
    ℐ[t]
    -[t]-[t]
    ;​​{varsSEIQRD,odesSEIQRD}=Values@KineticCompartmentalModel
    βλ[t]
    →
    ℰ,ℰ
    ζ
    →
    ℐ,ℐ
    δ
    →
    ,ℐ
    γ
    →
    ℛ,
    γ
    →
    ℛ,ℐ
    μ
    →
    ,
    μ
    →
    ,t{1,-1};​​Column[odesSEIQRD]
    Out[]=
    ′
    
    [t]μℐ[t]+μ[t]
    ′
    ℰ
    [t]-ζℰ[t]+β[t]λ[t]
    ′
    ℐ
    [t]ζℰ[t]-γℐ[t]-δℐ[t]-μℐ[t]
    ′
    
    [t]δℐ[t]-γ[t]-μ[t]
    ′
    ℛ
    [t]γℐ[t]+γ[t]
    ′
    
    [t]-β[t]λ[t]
    In[]:=
    {susceptibleSEIQRD,exposedSEIQRD,infectedSEIQRD,quarantinedSEIQRD,recoveredSEIQRD,deadSEIQRD}={,ℰ,ℐ,,ℛ,}/.ParametricNDSolve[Join[odesSEIQRD/.forceOfInfectionSEIQRD,{[0]-I0,ℰ[0]0,ℐ[0]I0,[0]0,ℛ[0]0,[0]0}],Head/@varsSEIQRD,{t,0,500},{,I0,β,ζ,γ,δ,μ}];

    outbreak start 20 Dec 2019, 30 Dec 2019

    Optimization—
    500000
    ,
    100000
    ,
    50000

    [see full notebook for elided content]

    Comparisons

    Out[]//TableForm=
    startDate
    populationSize
    β
    ζ
    γ
    δ
    μ
    SSE
    20 Dec 2019
    50000
    0.813
    0.376419
    0.0156619
    0.372657
    0.00224936
    4.93066×
    8
    10
    20 Dec 2019
    100000
    3.15091
    1.68687
    0.0206171
    2.78794
    0.00429265
    7.18726×
    8
    10
    20 Dec 2019
    500000
    16.4034
    6.02909
    0.0195967
    16.0677
    0.00383133
    1.40148×
    9
    10
    30 Dec 2019
    50000
    2.0433
    0.345355
    0.0191734
    1.06171
    0.0024311
    2.72508×
    8
    10
    30 Dec 2019
    100000
    26.3555
    1.66974
    0.0218454
    24.0305
    0.00438372
    4.19801×
    8
    10
    30 Dec 2019
    500000
    7.90897
    45.8682
    0.0208074
    7.765
    0.00404869
    9.93811×
    8
    10
    Out[]=
    ◼
  • increasing the population size makes  and ℛ rise more slowly, and decreases the epidemic size
    • ◼
  • shifting t0 to a later date (closer to the data) appears to make the  peak narrower, and reverses the population size trends for β and δ

    • Hubei outbreak—SEIsaRD

    SEIsaRD Model

    ◼
  • SEIsaRD: SEIR model with symptomatic and asymptomatic infectious compartments
    • ◼
  • Standard incidence
    • ◼
  • The duration of infection for symptomatic and asymptomatic infectious individuals is the same
    • ◼
  • We have two different recovered compartments for book keeping
    • ◼
  • death due to disease occurs only for symptomatic infectious individuals
    • ◼
  • These are the transitions:
    • Out[]=
    
    βλ[t]
    →
    ℰ
    infection
    ℰ
    ζ
    ρ
    s
    →
    ℐ
    s
    incubation, symptomatic
    ℰ
    ζ(1-
    ρ
    s
    )
    →
    ℐ
    a
    incubation, asymptomatic
    ℐ
    s
    γ
    →
    ℛ
    s
    recovery, symptomatic
    ℐ
    a
    γ
    →
    ℛ
    a
    recovery, asymptomatic
    ℐ
    s
    μ
    →
    
    death, symptomatic
    ◼
  • These are the ODEs
    • forceOfInfectionSEIsaRD=λ[t_]
    ℐ
    "s"
    [t]+
    ℐ
    "a"
    [t]
    -[t]
    ;​​{varsSEIsaRD,odesSEIsaRD}=Values@KineticCompartmentalModel
    βλ[t]
    →
    ℰ,ℰ
    ζρ
    →
    ℐ
    "s"
    ,ℰ
    ζ(1-ρ)
    →
    ℐ
    "a"
    ,
    ℐ
    "s"
    γ
    →
    ℛ
    "s"
    ,
    ℐ
    "a"
    γ
    →
    ℛ
    "a"
    ,
    ℐ
    "s"
    μ
    →
    ,t{1,-1};​​Column[odesSEIsaRD]
    Out[]=
    ′
    
    [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]
    In[]:=
    {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}],Head/@varsSEIsaRD,{t,0,500},{,I0,β,ζ,ρ,γ,μ}];

    outbreak start 20 Dec 2019

    Optimization—
    100000
    ,
    50000

    [see full notebook for elided content]

    Comparisons

    Out[]//TableForm=
    startDate
    populationSize
    β
    ζ
    ρ
    γ
    μ
    SSE
    20 Dec 2019
    50000
    0.201151
    2.68809
    0.990566
    0.0128232
    0.00454096
    3.69765×
    8
    10
    20 Dec 2019
    100000
    0.238916
    1.19482
    0.509441
    0.0175598
    0.00424525
    2.23757×
    8
    10
    Out[]=
    ◼
  • Increasing the population size increases the recovery rate
    • ◼
  • in all cases, the susceptible compartment goes to 0 and the epidemic size is 100% of the population

    • Summary

    ◼
  • All three models (SEIRD, SEIQRD, and SEIsaRD) exhibit long exponential tails for the ℐ or  peak
    • ◼
  • The recovered compartment never recovers as fast as the data
    • ◼
  • The deaths compartment is usually fit well
    • ◼
  • Most of the models gave an epidemic size of 100 %
    • ◼
  • More work is needed

    • Suggestions for improvements & new studies

    ◼
  • Coburn et al. [CWB] describe an isolation compartment in equilibrium with the susceptible compartment
    • ◼
  • This compartment effectively reduces the susceptible compartment size and should allow mass action incidence to be used
    • ◼
  • Jia et al. [JDL] use a similar compartment for their COVID-19 model, as well as mass action incidence
    • ◼
  • They only show confirmed cases against the diagnosed compartment
    D
    in the manuscript; plotting all the data against
    D
    ,
    R
    , and the sum of the two death sinks shows that the model is inadequate
    • ◼
  • New model under consideration
    • ◼
  • The three circled compartments represent the reported data:  ≡ confirmed cases,
    ℛ
    d
    ≡ recovered,  ≡ deaths
    • ◼
  • Are the data for Hubei too idiosyncratic for modeling? Are the data from another province or country (e.g., South Korea) more amenable to modeling?
    • References
    [JHU] “Mapping 2019-nCoV”, https://systems.jhu.edu/research/public-health/ncov/
    [JEP] J. Espigule-Pons “Mapping Novel Coronavirus COVID-19 Outbreak”, https://community.wolfram.com/groups/-/m/t/1868945
    [TG] T. Götz, “First attempts to model the dynamics of the coronavirus outbreak 2020”, https://arxiv.org/pdf/2002.03821.pdf
    [PYZ] L. Peng, W. Yang, D. Zhang, C. Zhuge, L. Hong “Epidemic analysis of COVID-19 in China by dynamical modeling”, https://www.medrxiv.org/content/10.1101/2020.02.16.20023465v1
    [ZCW] Y. Zhou, Z. Chen, X. Wu, Z. Tian, L. Cheng, L. Ye “The Outbreak Evaluation of COVID-19 in Wuhan District of China”, https://arxiv.org/pdf/2002.09640.pdf
    [JDL] J. Jia, J. Ding, S. Liu, G. Liao, J. Li, B. Duan, G. Wang, R. Zhang “Modeling the Control of COVID-19: Impact of
    Policy Interventions and Meteorological Factors”, https://arxiv.org/pdf/2003.02985.pdf
    [EGE] E. G. M E. “An SEIR like model that fits the coronavirus infection data”, https://community.wolfram.com/groups/-/m/t/1888335
    [AA] A. Antonov “Basic experiments workflow for simple epidemiological models”, https://community.wolfram.com/groups/-/m/t/1895686
    [KM] W.O. Kermack, A. G. McKendrick “Contributions to the Mathematical Theory of Epidemics—I” reprinted in Bull. Math. Biol., 53, 33-55 (1991). https://link.springer.com/article/10.1007/BF02464423
    [EY] E. Yong “The Deceptively Simple Number Sparking Coronavirus Fears”, The Atlantic, 28 Jan 2020, https://www.theatlantic.com/science/archive/2020/01/how-fast-and-far-will-new-coronavirus-spread/605632/
    [AV] J. Arino, P. van den Driessche “Time delays in Epidemic Models; Modeling and Numerical Considerations” in “Delay Differential Equations and Applications”, O. Arino (ed.) Springer, 2006.
    [FB] F. Brauer “Reproduction numbers and final size relations”, https://www.fields.utoronto.ca/programs/scientific/10-11/drugresistance/emergence/fred1.pdf
    [BCR] M. Biggerstaff, S. Cauchemez, C. Reed, M. Gambhir, L. Finelli “Estimates of the reproduction number for seasonal, pandemic, and zoonotic influenza: a systematic review of the literature” BMC Infectious Diseases, 14, 480 (2014), http://www.biomedcentral.com/1471-2334/14/480
    [MM] M. Martcheva “An introduction to mathematical epidemiology” Springer, 2015.
    [A] Anonymous, Anonymous, Brit. Med. J., 1978, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1603269/pdf/brmedj00115-0064.pdf
    [HJR] H. J. Rose “The use of amantadine and influenza vaccine in a type A influenza epidemic in a boarding school”, Journal of Royal College of General Practitioners, 30, 619-621 (1980). PubMedCentral
    [FT] Z. Feng, H. R. Thieme “Recurrent Outbreaks of Childhood Diseases Revisited: The Impact of Isolation”, Math. Biosciences, 128, 93-130 (1995). https://doi.org/10.1016/0025-5564(94)00069-C
    [BK] S. Boseley, L. Kuo “Huge rise in coronavirus cases casts doubt over scale of epidemic”, The Guardian, 13 Feb 2020, https://www.theguardian.com/world/2020/feb/13/huge-rise-coronavirus-cases-raises-doubts-scale-epidemic-china
    [DWC] Z. Du, L. Wang, S. Cauchemex, X. Xu, X. Wang, B. J. Cowling, L. A. Meyers “Risk for Transportation of 2019 Novel Coronavirus (COVID-19) from Wuhan to Cities in China”, https://doi.org/10.1101/2020.01.28.20019299
    [CXL] J. Cai, J. Xu, D. Lin, Z. Yang, L. Xu, Z, Qu, Y. Zhang, H. Zhang, R. Jia, P. Liu, X. Wang, Y. Ge, A. Xia, H. Tian, H. Chang, C. Wang, J. Li, J. Wang, M. Zheng “A Case Series of children with 2019 novel coronavirus infection: clinical and epidemiological features”, Clinical Infectious Diseases, https://doi.org/10.1093/cid/ciaa198
    [CWB] B. J. Coburn, B. G. Wagner, S. Blower “Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1)”, BMC Medicine, 7, (2009), http://www.biomedcentral.com/1741-7015/7/30
    Initialization
    [see full notebook for elided content]

    CITE THIS NOTEBOOK

    Epidemiological models for Influenza and COVID-19​
    by Robert Nachbar​
    Wolfram Community, STAFF PICKS, March, 2020
    ​https://community.wolfram.com/groups/-/m/t/1896178