Stochastic Network Models for Epidemiology—Application to COVID-19 Pandemic

Robert B. Nachbar, Wolfram Solutions

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Many SEIR-type epidemiology models for COVID-19 have been submitted to preprint servers since the beginning of the year, and some have been published. Nearly all these models have been based on deterministic differential equation models, which assume large, rapidly and well-mixed populations. These assumptions rarely hold for real biological populations, and other methods may be more appropriate. One frequent feature of these “published” models and observed in our own models is the need to use artificially small total population sizes to fit the reported number of cases and deaths in a given country or region. Network-based models offer an alternative to the well-mixed assumption, and stochastic methods provide an alternative to the rapid mixing assumption.

Preliminary results combining our previously developed stochastic simulation method using the Gillespie algorithm with a network architecture of subpopulations have been encouraging. In this talk we will describe our on-going efforts to fully explore the effects of network topology and various mixing scenarios on the dynamics of disease spread throughout the entire population. Comparison of the modeled outcomes with some well-characterized COVID-19 data from Europe and China will show whether or not this approach is successful, or that other modeling strategies are needed.

Spoiler: While we did implement a method to stratify a compartmental model over a network, and demonstrated deterministic and stochastic simulations, much more remains to be done to characterize the effects of network size and topology on the modeled outcomes. Therefore, this presentation is more of about methodology than application.

Outline

Background & Rationale

Epidemic Outbreaks and Pandemics

Compartmental Epidemic Models

Construction

Parametrization

Simulation

Non-Pharmaceutical Interventions (NPI)

Overview

The NPI function

Revisiting the Hubei data and model

A more realistic compartmental model

Structure

Parameterization

Simulation

Network models

Network Topologies

Examples of linear and planar networks:

Examples of nonplanar surface networks

Examples of nonplanar real-world networks

Distribute the Epidemiology Model Over a Network

A Network Model for Italy

Structure

Parameterization

Deterministic Simulation

Background & Rationale

Epidemic Outbreaks and Pandemics

\[KernelIcon];locale=(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"]);

\[KernelIcon];exampleData=fitDataWDR[locale];
Take[%,1]
Keys[%%]

Raw data:

\[KernelIcon];DateListPlot[{exampleData["CumulativeConfirmed"],exampleData["CumulativeRecovered"],exampleData["CumulativeDead"]},]

Epidemiological outcomes:

\[KernelIcon];outcomesPlot=DateListPlot[{exampleData["CumulativeConfirmed"]-exampleData["CumulativeRecovered"]-exampleData["CumulativeDead"],exampleData["NewConfirmed"],exampleData["CumulativeDead"]},]

In March of this year we attempted to fit a simple SEIRD type compartmental model to these kinds of data

We used data on the various provinces in China as the pandemic had not yet spread to many other countries

In order to fit the data we followed the example of other researchers of reducing the population to a much smaller effective population, but the results were not satisfactory

Here are example fits for the Hubei province, whose actual population is 58160000

While the curves look “OK” the parameters are really unrealistic!

(Subscript[ℜ, 0]) can be calculated from data a number of ways, and these results will not agree with those estimates

Compartmental Epidemic Models

Construction

\[KernelIcon];<<EpidemiologyModeling`

\[KernelIcon];modelSEIRD={
𝒮(Overscript[→, β λ[t]])ℰ,
ℰ(Overscript[→, ζ])ℐ,
ℐ(Overscript[→, γ])ℛ,
ℐ(Overscript[→, μ])𝒟
};

\[KernelIcon];CompartmentalModelGraph[modelSEIRD,]

\[KernelIcon];?DynamicTransmissionModel

\[KernelIcon];modelData=DynamicTransmissionModel[modelSEIRD,t,"Infectious"->ℐ,"Incidence"->{"QuarantineReduced",𝒩,𝒟}];
vars=modelData["Variables"];
(odes=modelData["Equations"])//Column
forceOfInfection=modelData["ForceOfInfection"]

Initial conditions:

\[KernelIcon];inits=#==Switch[#,
ℐ[0],I0,
𝒮[0],𝒩-I0,
_,0]&/@Through[vars[0]]

The outcomes of interest:

\[KernelIcon];prevalence=ℐ[t]
incidence=Total@Incidence[modelData,ℐ]
mortality=𝒟[t]

Compartmental Epidemic Models

Parametrization

To simulate the model we first need some parameters

The infection step 𝒮(Overscript[→, β λ[t]])ℰ is parameterized in terms of the contact rate β

The basic reproduction (Subscript[ℜ, 0]) is sometimes more useful

It can be computed using the next generation matrix method

\[KernelIcon];?NextGenerationMatrix

\[KernelIcon];ngm=Quiet@NextGenerationMatrix[odes, vars, t,{ℰ,ℐ}, forceOfInfection];
maxVal=Eigenvalues[ngm]/.𝒮[t]->𝒩//Max////;
betaSoln=First@Solve[(Subscript[ℜ, 0])==maxVal,β]

\[KernelIcon];params=<|(Subscript[ℜ, 0])->2.5,β->β,ζ->(1/(Subscript[τ, inc])),γ->(1/(Subscript[τ, rec])),μ->(1/(Subscript[τ, mort])),𝒩->QM@locale["Population"],I0->1,(Subscript[τ, inc])->6.,(Subscript[τ, rec])->7.,(Subscript[τ, mort])->12.|>/.betaSoln
paramsN=Normal@(params//.params/.𝒟[t]->0)//Echo[#,"OK?",MatchQ[Values@#,{__?Positive}]&]&

Compartmental Epidemic Models

Simulation

We need to know when the virus was first circulating

“The virus was first confirmed to have spread to Italy on 31 January 2020, when two Chinese tourists in Rome tested positive for the virus”—Wikipedia

\[KernelIcon];firstInfectionDate=(DateObject["Day: ""Fri 31 Jan 2020", DateObject[{2020,1,31},"Day","Gregorian",-4.`]]);
timeHorizon=QM@DateDifference[firstInfectionDate,Today];

Here is a simulation for the natural course of an epidemic outbreak:

This unphysical remediation to the model was a roadblock to further model fitting, and we pursued other modeling frameworks:

Stochastic simulation

Network models

Non-Pharmaceutical Interventions (NPI)

Overview

When this work was done, the pandemic was only weeks old, and we did not realize that a crucial component of the model was missing: non-pharmaceutical interventions (NPI)

In mid January, China shut down all travel within and across its borders

People were required to stay at home unless absolutely necessary

Special hospitals and quarantine centers were built or established to isolate infectious persons

We knew that these measures were taking place, but modeling them in detail seemed to be, and actually is, a formidable task

We also didn’t realize how profoundly and immediately these measures can change the course of an epidemic outbreak; we were focused on understanding and modeling the natural history of the disease SARS-CoV-2

Modeling their combined effects is feasible, as they all reduce the contact rate β, and it therefore becomes a function of time β(t)

Because β∝(Subscript[ℜ, 0]), we can think of modulating (Subscript[ℜ, 0]) for which there is a clear target:

When (Subscript[ℜ, 0])>1 the outbreak will persist

When (Subscript[ℜ, 0])<1 the outbreak will dissipate

Here are some examples for NPI functions:

Because the mitigation measures taken in different countries, provinces, and states proceeded in stages and with varying degrees of compliance, we’ll use a combination of the latter two

Non-Pharmaceutical Interventions (NPI)

Revisiting the Hubei data and model

NYT, 6 August 2020:

JAN. 23
Wuhan, a city of more than 11 million, was cut off by the Chinese authorities.

The Chinese authorities closed off Wuhan by canceling planes and trains leaving the city, and suspending buses, subways and ferries within it. At this point, at least 17 people had died and more than 570 others had been infected, including in Taiwan, Japan, Thailand, South Korea and the United States.

https://www.nytimes.com/article/coronavirus-timeline.html

\[KernelIcon];tau=QM@DateDifference["20 Dec 2019","23 Jan 2020"]

\[KernelIcon];With[{npi=({{\[Piecewise], ({{1, 0≤t<34}, {(E^(-0.08` (-34+t))), 34≤t<∞}, {0, True}})}})},
DateListPlot[Table[npi,{t,0,120}],"20 Dec 2019",]
]

Including this simple NPI function in the model made a huge improvement:

But, there is still room for more improvement

(Subscript[ℜ, 0]) is still too large

the infectious period is too long

the mortality rate is too low

A more realistic compartmental model

Before we move on to distributing the model over a network, we should construct a more realistic model

Structure

The CDC website gives a fairly detailed description of the natural history, or timeline,  of COVID-19

https://www.cdc.gov/coronavirus/2019-ncov/hcp/planning-scenarios.html

The the data given for these events allow them to be placed on a timeline for the available age groups

A plausible compartmental can be constructed:

Susceptible: 𝒮

Exposed: ℰ

Asymptomatic, mild symptomatic, and severe symptomatic infectious compartments: (Subscript[ℐ, a]), (Subscript[ℐ, m]), (Subscript[ℐ, s])

Hospitalized and intensive care unit (effectively quarantined) compartments: (Subscript[𝒬, m]), (Subscript[𝒬, s])

Asymptomatic, mild symptomatic, and severe symptomatic recovered: (Subscript[ℛ, a]), (Subscript[ℛ, m]), (Subscript[ℛ, s])

Dead: 𝒟

\[KernelIcon];modelName="COVID-19";
model={
𝒮(Overscript[→, β λ[t]])ℰ,
ℰ(Overscript[→, α ζ])(Subscript[ℐ, a]),
ℰ(Overscript[→, (1-α)(1-σ)ζ])(Subscript[ℐ, m]),
ℰ(Overscript[→, (1-α)σ ζ])(Subscript[ℐ, s]),
(Subscript[ℐ, a])(Overscript[→, (Subscript[γ, a])])(Subscript[ℛ, a]),
(Subscript[ℐ, m])(Overscript[→, η])(Subscript[𝒬, m]),(Subscript[𝒬, m])(Overscript[→, (Subscript[γ, m])])(Subscript[ℛ, m]),
(Subscript[ℐ, s])(Overscript[→, (1-δ)η])(Subscript[𝒬, s]),(Subscript[𝒬, s])(Overscript[→, (Subscript[γ, s])])(Subscript[ℛ, s]),
(Subscript[ℐ, s])(Overscript[→, δ (Subscript[μ, i])])𝒟,
(Subscript[𝒬, s])(Overscript[→, δ (Subscript[μ, q])])𝒟
};

\[KernelIcon];modelGraph=CompartmentalModelGraph[model,]//

\[KernelIcon];modelData=DynamicTransmissionModel[model,t,"Infectious"->{(Subscript[ℐ, a]),(Subscript[ℐ, m]),(Subscript[ℐ, s])},"Incidence"->{"QuarantineReduced",𝒩,{(Subscript[𝒬, m]),(Subscript[𝒬, s]),𝒟}}];
vars=modelData["Variables"];
(odes=modelData["Equations"])//Column
forceOfInfection=modelData["ForceOfInfection"]

A more realistic compartmental model

Parameterization

\[KernelIcon];ngm=Quiet@NextGenerationMatrix[odes, vars, t,{ℰ,(Subscript[ℐ, a]),(Subscript[ℐ, m]),(Subscript[ℐ, s])}, forceOfInfection];
maxVal=Eigenvalues[ngm]/.𝒮[t]->𝒩//Max////;
betaSoln=First@Solve[(Subscript[ℜ, 0])==maxVal,β]

The parameter values are those recommended by the CDC

(Subscript[ℜ, 0]) and (Subscript[τ, mort]) were determined by optimization using the prevalence and mortality data

\[KernelIcon];params=<|(Subscript[ℜ, 0])->6.2,β->β,ζ->(1/(Subscript[τ, inc])),η->(1/(Subscript[τ, hosp])),(Subscript[γ, a])->(1/(Subscript[τ, rec])),(Subscript[γ, m])->(1/(Subscript[τ, relHosp])),(Subscript[γ, s])->(1/(Subscript[τ, relICU])),(Subscript[μ, i])->(1/(Subscript[τ, mort])),(Subscript[μ, q])->(1/((Subscript[τ, mort])-(Subscript[τ, hosp]))),𝒩->QM@locale["Population"](*5 (10^7)*),α->0.40,σ->0.39,δ->0.12,I0->1,(Subscript[τ, inc])->6.,(Subscript[τ, hosp])->5.6,(Subscript[τ, rec])->7.,(Subscript[τ, relHosp])->3.9,(Subscript[τ, relICU])->12.,(Subscript[τ, mort])->20.0|>/.betaSoln
paramsN=Normal@(params//.params/.𝒟[t]->0)/.{β[t]->β}//Echo[#,"OK?",MatchQ[Values@#,{__?Positive}]&]&

A more realistic compartmental model

Simulation

\[KernelIcon];inits=Join[{𝒮[0]==𝒩-I0,(Subscript[ℐ, m])[0]==I0},Thread[Through@Complement[vars,{𝒮,(Subscript[ℐ, m])}][0]==0]]

We already know the date of the first infection

\[KernelIcon];firstInfectionDate=(DateObject["Day: ""Fri 31 Jan 2020", DateObject[{2020,1,31},"Day","Gregorian",-4.`]]);
timeHorizon=QM@DateDifference[firstInfectionDate,Today];

To define the NPI function, we also need to know when mitigation policies were put in place, and when they were relaxed

“On 8 March 2020, Prime Minister Giuseppe Conte expanded the quarantine to all of Lombardy and 14 other northern provinces, and on the following day to all of Italy, placing more than 60 million people in lockdown.”—Wikipedia

“In May, many restrictions were gradually eased, and on 3 June freedom of movement across regions and other European countries was restored.”—Wikipedia

The parameters for the NPI function were determined by optimization using the prevalence and mortality data

\[KernelIcon];solnODE=First@NDSolve[{odes,inits}/.forceOfInfection/.betaSoln/.{(Subscript[ℜ, 0])->npi (Subscript[ℜ, 0])}/.paramsN,vars,{t,0,timeHorizon}];

These are the outcomes we want to compare against the data

\[KernelIcon];prevalence=(Subscript[ℐ, m])[t]+(Subscript[ℐ, s])[t]+(Subscript[𝒬, m])[t]+(Subscript[𝒬, s])[t];
incidence=Total@Values@Incidence[modelData,{(Subscript[ℐ, m]),(Subscript[ℐ, s])}]//Simplify;
deaths=𝒟[t];

\[KernelIcon];FlipView[{modelPlot=DateListPlot[Transpose@Table[{prevalence,incidence,deaths},{t,0,timeHorizon}]/.paramsN/.solnODE,firstInfectionDate,],Show[outcomesPlot,modelPlot,PlotRange->{0,130000}]}]
closeCellGroup[]

Not bad, however, some of the parameters are outside or near the edges of the ranges suggested by the CDC:

({{"Parameter", "Value", "CDC Median", "CDC IQR"}, {(Subscript[ℜ, 0]), 6.2`, 2.5`, "\!\(\*FormBox[TagBox[FormBox[\"2.`\", TraditionalForm], TraditionalForm, Rule[Editable, True]], TraditionalForm]\) – \!\(\*FormBox[TagBox[FormBox[\"4.`\", TraditionalForm], TraditionalForm, Rule[Editable, True]], TraditionalForm]\)"}, {(Subscript[τ, mort]), 20.`, 15.`, "\!\(\*FormBox[TagBox[FormBox[\"9.`\", TraditionalForm], TraditionalForm, Rule[Editable, True]], TraditionalForm]\) – \!\(\*FormBox[TagBox[FormBox[\"24.`\", TraditionalForm], TraditionalForm, Rule[Editable, True]], TraditionalForm]\)"}})

Network models

The network concept for our models has these features:

Each vertex (node) represents a homogeneous community population that is well mixed

An independent epidemiological model operates at each vertex

The directed edges between vertices represent migration of individuals

Thus, the number of individuals in each vertex can change over time

The flow is bidirectional (a directed edge for each direction) and follows mass action kinetics

The network topology (number of vertices and their connections) is static over time

In principle, the flow rate constants along a pair of edges could go to 0, effectively disconnecting the two vertices

The aggregate behavior of the interconnected, independent models represents the aggregate epidemiology of a whole country, province, or state.

Mathematically, the model has been stratified by community

Network Topologies

There is a whole zoo of topologies from which to choose, from highly structured and symmetric to random

Random networks with desirable properties are typically used to model real world phenomena

Regular networks offer the opportunity to study the effects of particular network properties and how they scale with network size

Examples of linear and planar networks:

({{({{}, {(Row[(Subscript[N, v])==6, (Subscript[N, e])==5, ])}}), ({{}, {(Row[(Subscript[N, v])==6, (Subscript[N, e])==6, ])}}), ({{}, {(Row[(Subscript[N, v])==38, (Subscript[N, e])==49, ])}})}, {({{}, {(Row[(Subscript[N, v])==20, (Subscript[N, e])==31, ])}}), ({{}, {(Row[(Subscript[N, v])==20, (Subscript[N, e])==43, ])}}), ({{}, {(Row[(Subscript[N, v])==40, (Subscript[N, e])==102, ])}})}})

Examples of nonplanar surface networks

({{({{}, {(Row[(Subscript[N, v])==54, (Subscript[N, e])==108, ])}}), ({{}, {(Row[(Subscript[N, v])==54, (Subscript[N, e])==162, ])}}), ({{}, {(Row[(Subscript[N, v])==100, (Subscript[N, e])==294, ])}})}})

Examples of nonplanar real-world networks

({{"Barabasi Albert Graph Distribution", "Watts Strogatz Graph Distribution"}, {, }, {(Row[(Subscript[N, v])==60, (Subscript[N, e])==174, ]), (Row[(Subscript[N, v])==60, (Subscript[N, e])==120, ])}})

The Barabasi-Albert graph distribution is good for generating scale-free networks, such as the internet, citation networks, and some social networks (Wikipedia)

The Watts-Strogatz graph distribution is good for networks with small-world properties, such as short average path lengths and high clustering (Wikipedia)

Distribute the Epidemiology Model Over a Network

By way of example, we can use the simple three-compartment SIR model and a simple 4-node linear network

Here are the basic steps for constructing a network model

First, construct the model for a single node

\[KernelIcon];modelSIR={𝒮(Overscript[→, β λ[t]])ℐ,ℐ(Overscript[→, γ])ℛ};
\[KernelIcon];modelDataSIR=DynamicTransmissionModel[modelSIR,t,"Infectious"->{ℐ},"Incidence"->{"Standard",𝒩}];
modelDataSIR["Equations"]//Column

\[KernelIcon];ngm=Quiet@NextGenerationMatrix[modelDataSIR["Equations"], modelDataSIR["Variables"], t,{ℐ}, modelDataSIR["ForceOfInfection"]];
maxVal=Eigenvalues[ngm]/.ℛ[t]->0//Max;
betaSolnSIR=First@Solve[Echo[(Subscript[ℜ, 0])==maxVal],β]

\[KernelIcon];paramsSIR=<|(Subscript[ℜ, 0])->2.5,β->β,γ->(1/7.),𝒩->(10^3),I0->1|>/.betaSolnSIR
paramsSIRN=paramsSIR//.paramsSIR//Normal

Then construct the network as a Graph

\[KernelIcon];nNodes=4;
networkGraph=linearNetwork[nNodes,ImageSize->Large,BaseStyle->{28}]

Distribute the model over the network and update the list of parameters

\[KernelIcon];?distributeModel

\[KernelIcon];networkModelData=distributeModel[modelDataSIR,networkGraph,𝓂,
"StratifiedParameters"->{β,𝒩},
"AncillaryExpressions"->betaSolnSIR]/.𝓂[__]:>𝓂;
Keys@%

\[KernelIcon];networkVars=networkModelData["Variables"];
networkForceOfInfection=networkModelData["ForceOfInfection"];
networkOdes=networkModelData["Equations"];
networkBetaSoln=networkModelData["AncillaryExpressions"];

\[KernelIcon];networkParams=Join[
\[LeftAssociation](Subscript[ℜ, 0])->2.5,γ->1/7.,I0->1,𝓂->(1/90.)\[RightAssociation],
Association@Table[β[v]->β[v],{v,Range@nNodes}],
Association@Table[𝒩[v]->1000,{v,Range@nNodes}]/N@nNodes
]/.networkBetaSoln
networkParamsN=networkParams//.networkParams//Normal

Finally, define the initial conditions and simulate

\[KernelIcon];networkInits=With[{nonZero=Join[{𝒮[1][0]==𝒩[1]-I0,ℐ[1][0]==I0},𝒮[#][0]==𝒩[#]&/@Range[2,nNodes]]},Join[nonZero,Thread[Complement[Through@networkVars[0],First/@nonZero]==0]]]

\[KernelIcon];networkSolnODE=First@NDSolve[{networkOdes,networkInits}/.networkBetaSoln/.networkForceOfInfection/.networkParamsN,networkVars,{t,0,timeHorizon}];

The following figure compares the single-node model to the linear network model from their structures through their simulations

The same rate constants are used in each model

The same total population is used in each model

Clearly, distributing the model over a network can yield qualitatively different results

Moreover, stochastic simulation demonstrates that an infection might never take hold to produce an outbreak

A Network Model for Italy

Structure

Recall the model structure:

\[KernelIcon];modelGraph

We can use a randomly generated network, although the choice of the generation method is subjective

The Watts-Strogatz method may be appropriate because the nodes will have a small number of connections  and will better reflect a population under restrictions

\[KernelIcon];nNodes=100;

\[KernelIcon];SeedRandom[944397];
Row[{networkWattsStrogatz=RandomGraph[WattsStrogatzGraphDistribution[nNodes,0.3]],degreePlot[networkWattsStrogatz,]},Spacer[36]]

The Barabasi-Albert method will yield a few very highly connected nodes, which can act as super spreaders

\[KernelIcon];SeedRandom[944397];
Row[{networkBarabasiAlbert=RandomGraph[BarabasiAlbertGraphDistribution[nNodes,3]],degreePlot[networkBarabasiAlbert,]},Spacer[36]]

Note that the parameters 𝒷 and 𝒩 need to be stratified, as well as the equation for β in terms of (Subscript[ℜ, 0])

For simplicity, we assume that the migration rate constant is the same along all network edges

networkModelData=distributeModel[modelData,networkWattsStrogatz,𝓂,
"StratifiedParameters"->{β,𝒩},
"AncillaryExpressions"->betaSoln]/.𝓂[__]:>𝓂;
Length@networkModelData["Variables"]

\[KernelIcon];networkVars=networkModelData["Variables"];
networkForceOfInfection=networkModelData["ForceOfInfection"];
networkOdes=networkModelData["Equations"];
networkBetaSoln=networkModelData["AncillaryExpressions"];

A Network Model for Italy

Parameterization

\[KernelIcon];networkParams=Join[
<|(Subscript[ℜ, 0])->6.2,ζ->(1/(Subscript[τ, inc])),η->(1/(Subscript[τ, hosp])),(Subscript[γ, a])->(1/(Subscript[τ, rec])),(Subscript[γ, m])->(1/(Subscript[τ, relHosp])),(Subscript[γ, s])->(1/(Subscript[τ, relICU])),(Subscript[μ, i])->(1/(Subscript[τ, mort])),(Subscript[μ, q])->(1/((Subscript[τ, mort])-(Subscript[τ, hosp]))),α->0.40,σ->0.39,δ->0.12,I0->1,(Subscript[τ, inc])->6.,(Subscript[τ, hosp])->5.6,(Subscript[τ, rec])->7.,(Subscript[τ, relHosp])->3.9,(Subscript[τ, relICU])->12.,(Subscript[τ, mort])->20.0,𝓂->(1/90.)|>,
Association@Table[β[v]->β[v],{v,Range@nNodes}],
Association@Table[𝒩[v]->QM@locale["Population"],{v,Range@nNodes}]/N@nNodes
]/.networkBetaSoln;
networkParamsN=Normal@(networkParams//.networkParams/.𝒟[_][t]->0)//Echo[#,"OK?",MatchQ[Values@#,{__?Positive}]&]&;

A Network Model for Italy

Deterministic Simulation

\[KernelIcon];networkInits=With[{nonZero=Join[{𝒮[1][0]==𝒩[1]-I0,(Subscript[ℐ, m])[1][0]==I0},𝒮[#][0]==𝒩[#]&/@Range[2,nNodes]]},Join[nonZero,Thread[Complement[Through@networkVars[0],First/@nonZero]==0]]];

\[KernelIcon];networkSolnODE=First@NDSolve[{networkOdes,networkInits}/.networkBetaSoln/.networkForceOfInfection/.{(Subscript[ℜ, 0])->npi (Subscript[ℜ, 0])}/.networkParamsN,networkVars,{t,0,timeHorizon}];

Here is a plot of all the compartments, except susceptible 𝒮

epidemicModelPlot[FilterRules[networkSolnODE,Except[𝒮[_]]],{0,timeHorizon},"modelName"->modelName,PlotRange->{0,30000},"scale"->1,PlotLegends->None,]

Aggregate the compartments by health state

\[KernelIcon];aggregatedNetworkSolnODE=aggregateNetworkResult[networkSolnODE];

\[KernelIcon];epidemicModelPlot[FilterRules[aggregatedNetworkSolnODE,Except[𝒮]],{0,timeHorizon},"modelName"->modelName,PlotRange->{0,600000},"scale"->1,]

Compare the outcomes:

aggregated network model vs. data

aggregated network model vs. non-network model

\[KernelIcon];With[{networkModelPlot=DateListPlot[Transpose@Table[{prevalence,incidence,deaths},{t,0,timeHorizon}]/.networkParamsN/.aggregatedNetworkSolnODE,firstInfectionDate,]},
FlipView[{
Show[networkModelPlot,outcomesPlot,PlotRange->{0,All}],Show[networkModelPlot,modelPlot/.c_RGBColor:>Sequence[c,Dashed],PlotRange->{0,All}]
}]//Labeled[#,Style["Watts-Strogatz Network",36,Bold,Gray],Top]&
]
closeCellGroup[]

Summary

Artificially small population sizes lead to unrealistic optimized parameter values

Non-pharmaceutical intervention is a crucial part of COVID-19 epidemiology models

A model has been proposed with

asymptomatic, mild, and severe infectious compartments

hospitalized (quarantined) compartments

disease induced death compartment

non-pharmaceutical intervention

The proposed model successfully recapitulates the observed prevalence and death

The model has been extended by stratification over a community network

The detailed network model behavior strongly depends on the network structure

The aggregated network model behavior tends to be similar to the non-network model behavior, but more work needs to be done

Resources

Packages

CloudObject["https://www.wolframcloud.com/obj/rnachbar/Published/CompartmentalModeling.wl"]

CloudObject["https://www.wolframcloud.com/obj/rnachbar/Published/EpidemiologyModeling.wl"]

Wolfram Community Posts

Epidemiological Models for Influenza and COVID-19

Stochastic Epidemiology Models with Applications to the COVID-19

Web resources

COVID-19 Data & Resources

Resources for COVID-19 and the SARS-CoV-2 Virus

Thank You!

Initialization

SetDirectory@NotebookDirectory[];

SetOptions[ListPlot,{Axes->False,Frame->True}];
SetOptions[DateListPlot,{Axes->False,Frame->True,Joined->False}];
$PlotTheme = {"Scientific", "VibrantColor"};

color[i_]:=ColorData[{"Indexed","VibrantColor"}][i]

Unprotect[Missing]; 
UpValues[Missing] = {}; 
Missing /: (m_Missing) + ___ := m
Missing /: (m_Missing)*___ := m
Missing /: (m_Missing)^_ := m
Protect[Missing];

QM = QuantityMagnitude;

closeCellGroup[] := (SelectionMove[EvaluationCell[], Next,Cell]; 
FrontEndExecute[FrontEndToken["SelectionCloseUnselectedCells"]])

<<EpidemiologyModeling`
<<NPI`

fitDataWDR

Note: some countries started reporting reporting data by administrative division in May,

(* time is measured in Days *)
toModelTime[t0_DateObject, date_DateObject] := 
	QM[DateDifference[t0, date]]
toDataDate[t0_DateObject, time:(_Integer|_Real)] := 
	DatePlus[t0, time]

ResourceUpdate["Epidemic Data for Novel Coronavirus COVID-19"];

ClearAll[fitDataWDR]
fitDataWDR::usage = "fitDataWDR[locale] gives an Associaiton of the reported daily cumulative \
confirmed, recovered, and dead cases and computed daily new confirmed, recovered, and dead cases \
for the given locale, either an administrative division or a country."; 

fitDataWDR::locnf = "Locale `` not found."; 

$covid19Data = ResourceData["Epidemic Data for Novel Coronavirus COVID-19"]; 

Options[fitDataWDR] = {
	"dateRange" -> All
	}; 

fitDataWDR[locale:Entity["AdministrativeDivision", _] | Entity["Country", _], opts:OptionsPattern[]] := 
	Module[{cumConfirmed, cumRecovered, cumDead, confirmed, recovered, dead, dateRange}, 
		Switch[locale, 
			Entity["AdministrativeDivision", _], 
				If[!MemberQ[$covid19Data[All, "AdministrativeDivision"], locale], 
					Message[fitDataWDR::locnf, locale]; 
					Return[$Failed]
					], 
			Entity["Country", _], 
				If[!MemberQ[$covid19Data[All, "Country"], locale], 
					Message[fitDataWDR::locnf, locale]; 
					Return[$Failed]
					]
			]; 
		
		dateRange = OptionValue["dateRange"]; 
		dateRange = Switch[dateRange, 
			All, 
				$covid19Data[[1, "ConfirmedCases"]] /@ {"FirstDate", "LastDate"}, 
			{_, _}, 
				dateRange = (Quiet[DateObject[#1, "Instant"]] & ) /@ 
					(dateRange //. {{All, b_} :> {$covid19Data[[1,"ConfirmedCases"]]["FirstDate"], b}, 
						{a_, All} :> {a, $covid19Data[[1,"ConfirmedCases"]]["LastDate"]}}); 
				If[MatchQ[dateRange, {DateObject[{__}, _String, _String, _Real], DateObject[{__}, _String, _String, _Real]}], 
					dateRange, 
					Message[fitDataWDR::optvg, "dateRange", OptionValue["dateRange"], "All or a list of two dates"]; 
					Return[$Failed]
					], 
			_, 
				Message[fitDataWDR::optvg, "dateRange", OptionValue["dateRange"], "All or a list of two dates"]; 
				Return[$Failed]
			](* // Echo[#, "dateRange"]&*); 
		
		Switch[locale, 
			Entity["AdministrativeDivision", _], 
				cumConfirmed = $covid19Data[SelectFirst[#AdministrativeDivision == locale & ], "ConfirmedCases"]; 
				cumRecovered = $covid19Data[SelectFirst[#AdministrativeDivision == locale & ], "RecoveredCases"]; 
				cumDead = $covid19Data[SelectFirst[#AdministrativeDivision == locale & ], "Deaths"], 
			Entity["Country", _], 
				cumConfirmed = Plus @@ Normal[$covid19Data[Select[#Country == locale & ], 
					"ConfirmedCases"] /. Missing["NotAvailable"] -> 0]; 
				cumRecovered = Plus @@ Normal[$covid19Data[Select[#Country == locale & ], 
					"RecoveredCases"] /. Missing["NotAvailable"] -> 0]; 
				cumDead = Plus @@ Normal[$covid19Data[Select[#Country == locale & ], 
					"Deaths"] /. Missing["NotAvailable"] -> 0]
			]; 
		
		cumConfirmed = TimeSeriesWindow[cumConfirmed, dateRange]; 
		cumRecovered = TimeSeriesWindow[cumRecovered, dateRange]; 
		cumDead = TimeSeriesWindow[cumDead, dateRange]; 
		
		(*Normal /@ *)AssociationThread[
			{"CumulativeConfirmed", "CumulativeRecovered", "CumulativeDead", 
				"NewConfirmed", "NewRecovered", "NewDead"}, 
			Join[{cumConfirmed, cumRecovered, cumDead}, 
				Differences /@ {cumConfirmed, cumRecovered, cumDead}]
			]
		]

Figure for epidemic outbreak

DynamicModule[{susc,exp,inf,rec,dead,funcs,pop,startDate,tmax},
{susc,exp,inf,rec,dead}={𝒮,ℰ,ℐ,ℛ,𝒟}/.ParametricNDSolve[{odes/.forceOfInfection/.betaSoln,inits}/.FilterRules[paramsN,Except[𝒩]],vars,{t,0,timeHorizon},{𝒩}];
funcs={prevalence,incidence,mortality}/.paramsN/.{𝒮->susc[𝒩],ℰ->exp[𝒩],ℐ->inf[𝒩],ℛ->rec[𝒩],𝒟->dead[𝒩]};
pop=QM@locale["Population"];
startDate=firstInfectionDate;
tmax=timeHorizon;

TabView[{
"Full"->Column[
{DateListPlot[Thread@Table[funcs/.𝒩->pop,{t,0,tmax}],startDate,Joined->True,PlotMarkers->None,PlotLabel->locale["Name"],FrameLabel->{"date","# individuals"},PlotLegends->Placed[{"Prevalence","Incidence","Mortality"},After],PlotRange->All,PlotStyle->{epiColor[ℐ],epiColor[ℰ],epiColor[𝒟]},ImageSize->12 72,BaseStyle->{28},AspectRatio->1/(2 GoldenRatio)],
Show[outcomesPlot,ImageSize->12 72,BaseStyle->{28},AspectRatio->1/(2 GoldenRatio)]
}
],
"Various"->Manipulate[
Show[
DateListPlot[Thread@Table[funcs/.𝒩->10^log𝒩,{t,0,tmax}],startDate,Joined->True,PlotMarkers->None,PlotLabel->locale["Name"],FrameLabel->{"date","# individuals"},PlotLegends->Placed[{"Prevalence","Incidence","Mortality"},After],PlotRange->All,GridLines->{None,Automatic},PlotStyle->{epiColor[ℐ],epiColor[ℰ],epiColor[𝒟]},ImageSize->12 72,BaseStyle->{28}],
First@outcomesPlot
],

{{log𝒩,Log10@pop,Dynamic@Style[StringForm["𝒩 = ``",NumberForm[Round[10^log𝒩],DigitBlock->3]],28]},Log10@10000,Log10@pop,(Log10@pop-Log10@10000)/30}
]
}]
]

Figure for different NPIs

With[{
options=Sequence[ImageSize->Large,Axes->True,AxesOrigin->{0,0},AxesLabel->{"t"},AxesStyle->Thick,
Frame->False,
Ticks->{None,{{0,0},{0.6,1},{1,"(Subscript[ℜ, 0])"}}},TicksStyle->36,
PlotRange->{0,1.25},PlotRangePadding->Scaled[.05],
PlotStyle->Thick],
r0={0.6,Directive[Blue]},
low={0.25,Directive[Gray, Dashed]},
low2={0.5,Directive[Gray, Dashed]}
},
Row[{
Plot[0.25+0.75(1-UnitStep[t-2]),{t,0,10},GridLines->{None,{r0,low}},options],
Plot[0.25+0.75(1-UnitStep[t-2])+0.75 1/3(UnitStep[t-6]),{t,0,10},GridLines->{None,{r0,low,low2}},options],
Plot[0.25+0.75(E^(-0.4 t)),{t,0,10},GridLines->{None,{r0,low}},options]
},Spacer[36]]
]

Parameter grid

Grid[{
{"Parameter","Value","CDC Median","CDC IQR"},
{(Subscript[ℜ, 0]),6.2,2.5,StringForm["`` – ``",TraditionalForm@2.0,TraditionalForm@4.0]},
{(Subscript[τ, mort]),20.0,15.0,StringForm["`` – ``",TraditionalForm@9.0,TraditionalForm@24.0]}
},
Dividers->{{True,True,{LightGray},True},{True,True,{LightGray},True}},
BaseStyle->{ShowStringCharacters->False}
]//TraditionalForm

epidemicModelPlot

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

epiColor[𝒮|(Subscript[𝒮, _])|𝒮[_]|(Subscript[𝒮, _])[_]] = ColorData["Crayola"]["Green"];
epiColor[ℐ|(Subscript[ℐ, _])|ℐ[_]|(Subscript[ℐ, _])[_]] = ColorData["Crayola"]["Scarlet"];
epiColor[ℛ|(Subscript[ℛ, _])|ℛ[_]|(Subscript[ℛ, _])[_]] = ColorData["Crayola"]["Blue"];
epiColor[ℰ|(Subscript[ℰ, _])|ℰ[_]|(Subscript[ℰ, _])[_]] = ColorData["Crayola"]["Sunglow"];
epiColor[𝒬|(Subscript[𝒬, _])|𝒬[_]|(Subscript[𝒬, _])[_]] = ColorData["Crayola"]["Fuchsia"];
epiColor[𝒜|(Subscript[𝒜, _])|𝒜[_]|(Subscript[𝒜, _])[_]] = ColorData["Crayola"]["Salmon"];
epiColor[𝒟|(Subscript[𝒟, _])|𝒟[_]|(Subscript[𝒟, _])[_]] = ColorData["Crayola"]["Outerspace"];
epiColor[_] = ColorData["Crayola"]["Gray"];

ClearAll[epidemicModelPlot];

epidemicModelPlot::usage = "epidemicModelPlot[solution, {tmin, tmax}] plots each variable in solution, \
which has the form {var1 -> fun1, var2 -> fun2, ...} where funi is an InterpolatingFunction. 
epidemicModelPlot[{soln1, soln2, ...}, {tmin, tmax}] plots the median and 0.025 to 0.975 \
quantile band for each variable."; 

Options[epidemicModelPlot] = {
	"modelName" -> None, 
	"plotType" -> "Linear", 
	"scale" -> Automatic, 
	"showBands" -> True
	}; 

(* TODO: Tooltips *)

epidemicModelPlot[soln:{(_ -> (_InterpolatingFunction | _Function))..}, {tmin_, tmax_}, 
	opts:OptionsPattern[{epidemicModelPlot, Plot}]] := 
	Module[{vars, t, colors, plotType, plotFunc, scale, modelName, ymin, ymax, 
			gridLines, plotOptions, stochasticQ, singleVariableQ, plotLabel}, 
		modelName = OptionValue["modelName"]; 
		plotType = OptionValue["plotType"]; 
		plotFunc = Switch[plotType, 
			"Linear", Plot, 
			"Log", LogPlot, 
			_, 
				Message[epidemicModelPlot::optvg, "plotType", plotType, "\"Linear\" or \"Log\""];
				Return[$Failed]
			]; 
		
		vars = Keys @ soln; 
		singleVariableQ = SameQ @@ vars; 
		stochasticQ = MatchQ[InterpolatingFunctionInterpolationOrder /@ Values @ soln, 
			{{0}..}];
		
		colors = If[singleVariableQ,
			Automatic, 
		  (* else *)
			epiColor /@ vars
			]; 
		
		ymax = If[singleVariableQ, 
			InterpolatingFunctionValuesOnGrid /@ Values @ soln // Max, 
		  (* else *)
			(Through@vars[0] /. soln // Total) / Max[Counts[vars]]
			]; 
		ymin = Switch[plotFunc, LogPlot, 10^(-Log10[ymax]/10), _, 0]; 
		
		scale = OptionValue["scale"]; 
		scale = Switch[scale, 
			Automatic, 
				Switch[plotType, 
					"Linear", 10^(3 Floor[Log10[ymax] / 3]), 
					"Log", 1
					], 
			_Integer?Positive, 
				scale, 
			_, 
				Message[epidemicModelPlot::optvg, "scale", scale, "Automatic or a positive integer"];
				Return[$Failed]
			]; 

		gridLines = Switch[plotType, 
			"Linear", None, 
			"Log", {None, {1}}
			]; 
		
		plotLabel = If[stochasticQ, "Stochastic", "Deterministic"]; 
		plotLabel =  plotLabel <> 
			If[MatchQ[modelName, None], 
				"", 
				ToString@StringForm[" ``", modelName]
				];
		plotLabel =  plotLabel <> 
			If[singleVariableQ, 
				ToString@StringForm[", `` (n = ``)", First@vars, Length@vars], 
				""
				];
		
		plotOptions = Sequence @@ FilterRules[{opts}, Options[plotFunc]]; 
		
		Quiet @ plotFunc[(*If[singleVariableQ, 
				Max[ymin/2, (First@vars)[t]] /. List /@ soln, 
			  (* else *)
				Through@vars[t] /. soln
				]*)MapThread[Tooltip, {Max[ymin /scale / 2, # / scale]& /@ Through[Values[soln][t]], Keys[soln]}] // Evaluate,
			{t, tmin, tmax},
			PlotStyle -> colors,
			Frame -> True, 
			FrameLabel -> {"time (d)", If[scale == 1, 
					"cases", 
					StringForm["# individuals (``)", Switch[scale, 10^3, "thousands", 10^6, "millions", _, scale]]
					]},
			PlotLabel -> plotLabel, 
			PlotRangePadding -> {Scaled[0.02], Scaled[0.03]}, 
			GridLines -> gridLines, 
			Evaluate @ plotOptions, 
			PlotLegends -> If[singleVariableQ, None, Union @ vars],
			PlotRange -> {{tmin, tmax}, {ymin, ymax} / scale}
			]
		]

epidemicModelPlot[solns:{{(_ -> (_InterpolatingFunction | _Function))..}..}, {tmin_, tmax_}, 
	opts:OptionsPattern[{epidemicModelPlot, Plot}]] := 
	Module[{vars, t, colors, plotType, plotFunc, scale, smodelName, ymin, ymax, 
			gridLines, plotOptions, stochasticQ, 
			median, bands, n, filling}, 
		modelName = OptionValue["modelName"]; 
		plotType = OptionValue["plotType"]; 
		plotFunc = Switch[plotType, 
			"Linear", Plot, 
			"Log", LogPlot, 
			_, 
				Message[epidemicModelPlot::optvg, "plotType", plotType, "\"Linear\" or \"Log\""];
				Return[$Failed]
			]; 
		filling = OptionValue["showBands"]; 
		If[MatchQ[filling, Except[True | False]], 
			Message[epidemicModelPlot::opttf, "showBands", filling]; 
			Return[$Failed]
			]; 
		
		vars = Keys @ First @ solns; 
		n = Length@vars; 
		stochasticQ = MatchQ[InterpolatingFunctionInterpolationOrder /@ Values @ First @ solns, 
			{{0}..}];
		median = Median /@ Thread@(Through@vars[t] /. solns) // Quiet;
		
		colors = epiColor /@ vars; 
		If[filling, 
			bands = Outer[Quantile, Thread@(Through@vars[t] /. solns), 
				{0.025, 0.975}, 1] // Quiet(* // (Echo[#, "", Dimensions]&)*);
			filling = MapThread[First@#1 -> {{Last@#1}, #2} &, 
				{Partition[n + Range@(2 n), 2], Opacity[0.5, #] & /@ colors}];
			colors = Join[colors, ConstantArray[None, 2 n]];, 
		  (* else *)
			bands = {};
			filling = {};
			]; 
		
		ymax = Through@vars[0] /. solns // Map[Total, #]& // Max; 
		ymin = Switch[plotFunc, LogPlot, 10^(-Log10[ymax]/10), _, 0]; 
		
		scale = OptionValue["scale"]; 
		scale = Switch[scale, 
			Automatic, 
				Switch[plotType, 
					"Linear", 10^(3 Floor[Log10[ymax] / 3]), 
					"Log", 1
					], 
			_Integer?Positive, 
				scale, 
			_, 
				Message[epidemicModelPlot::optvg, "scale", scale, "Automatic or a positive integer"];
				Return[$Failed]
			]; 

		gridLines = Switch[plotType, 
			"Linear", None, 
			"Log", {None, {1}}
			]; 
		
		plotOptions = Sequence @@ FilterRules[{opts}, Options[plotFunc]]; 
		
		Quiet @ plotFunc[Flatten[Flatten@Join[median, bands]] / scale // Evaluate,
			{t, tmin, tmax},
			PlotStyle -> colors, 
			Filling -> filling, 
			PlotLegends -> vars,
			Frame -> True, 
			FrameLabel -> {"time (d)", If[scale == 1, 
					"cases", 
					StringForm["# individuals (``)", Switch[scale, 10^3, "thousands", 10^6, "millions", _, scale]]
					]},
			PlotLabel -> StringForm["Median `` ``(n = ``)", 
				If[stochasticQ, "Stochastic", "Deterministic"], 
				If[MatchQ[modelName, None], "", modelName<>" "], 
				Length@solns], 
			PlotRangePadding -> {Scaled[0.02], Scaled[0.03]}, 
			GridLines -> gridLines, 
			Evaluate @ plotOptions, 
			PlotRange -> {{tmin, tmax}, {ymin, ymax} / scale}
			] // 
			If[MatchQ[filling, {}], 
				#, 
				  (* else *)
				Labeled[#, 
					Style["Bands represent 0.025 and 0.975 quantiles", Gray, 10, 
						FontFamily -> "Sanserif"], 
					Bottom]
				]&
		]

linearNetwork

ClearAll[linearNetwork]
linearNetwork[n_Integer?Positive,opts:OptionsPattern[Graph]]:=
Module[{v=Range[n],e},
e=UndirectedEdge@@@Partition[v,2,1];
Graph[v,e,opts,VertexLabels->Automatic]
]

distributeModel

ClearAll[distributeModel]

distributeModel::usage = "distributeModel[modelData, graph, migrationRate] returns a new model data Association \
for the input model data distributed over the input graph. The symbol migrationRate is used to construct \
rate constants for the transitions between like compartments in at different nodes. Undirected edges in graph \
are converted to opposing directed edges."; 

Options[distributeModel] = {
	"StationaryVariables" -> None, 
	"StratifiedParameters" -> None, 
	"AncillaryExpressions" -> None
	}; 

distributeModel[modelData_, (network_Graph)/;UndirectedGraphQ[network] && EdgeCount[network]>0, 
		migrationRate_, opts:OptionsPattern[]] := 
	Module[{v, e}, 
		v = VertexList[network]; 
		e = Join[DirectedEdge @@@ EdgeList[network], 
			DirectedEdge @@@ Reverse /@ EdgeList[network]]; 
		distributeModel[modelData, Graph[v, e], migrationRate, opts]
		]

distributeModel[modelData_Association, (network_Graph)?DirectedGraphQ, 
		(migrationRate_)?( !ValueQ[#1] && !NumericQ[#1] & ), opts:OptionsPattern[]] := 
	Module[{statVars, stratParams, ancilExpr, model, vars, t, force, forceVars, vertices, edges, 
			varRules, paramRules, allRules, migration, networkModel}, 
		statVars = Flatten[{OptionValue["StationaryVariables"]} /. None -> {}]; 
		stratParams = Flatten[{OptionValue["StratifiedParameters"]} /. None -> {}]; 
		ancilExpr = Flatten[{OptionValue["AncillaryExpressions"]} /. None -> {}]; 
		
		model = modelData["ComponentTransitions"]; 
		vars = modelData["Variables"]; 
		t = modelData["Time"]; 
		force = modelData["ForceOfInfection"]; 
		forceVars = Keys[force] /. (v_)[modelData["Time"]] :> v; 
		
		vertices = VertexList[network]; 
		edges = EdgeList[network]; 
		
		varRules = MapThread[Join, {
			Table[Thread[vars -> Through[vars[v]]], {v, vertices}], 
			Table[Thread[forceVars -> Through[forceVars[v]]], {v, vertices}]}]; 
		paramRules = MapThread[Join, {
			Table[Thread[stratParams -> Through[stratParams[v]]], {v, vertices}]}]; 
		allRules = MapThread[Join, {varRules, paramRules}]; 
		
		migration = Outer[Transition[#1[First[#2]], migrationRate @@ #2, #1[Last[#2]]] & , 
			Complement[vars, statVars], edges]; 
		
		networkModel = Join[Join @@ (model /. allRules), Join @@ migration]; 
		
		Association[
			Quiet[DynamicTransmissionModel[networkModel, t]], 
			Association[
				"ForceOfInfection" -> Join @@ (force /. allRules), 
				"AncillaryExpressions" -> Join @@ (ancilExpr /. allRules)
				]
			](* // Timing // Echo[#, "DynamicTransmissionModel timing", First]& // Last*)
		]

distributeModel[modelData,linearNetwork[1],𝓂]

aggregateNetworkResult

ClearAll[aggregateNetworkResult]

aggregateNetworkResult[networkSolns:{{(_[_Integer] -> _InterpolatingFunction)..}..}] := 
	aggregateNetworkResult /@ networkSolns

aggregateNetworkResult[networkSoln:{(_[_Integer] -> _InterpolatingFunction)..}] := 
	Normal[GroupBy[Through[Keys[networkSoln][t$]], Head /* Head, 
		With[{total = Total[#1]}, Function[t, total]] & ]] /. Dispatch @ networkSoln

aggregateNetworkResult[networkResults:{__Association}] := 
	aggregateNetworkResult /@ networkResults

aggregateNetworkResult[networkResult:_Association] /; MatchQ[Keys@networkResult, {(_[_Integer])..}]:= 
	GroupBy[Keys[networkResult], Head, Total[#1] & ] /. Dispatch @ Normal @ networkResult

degreePlot

ClearAll[degreePlot]
degreePlot[g_Graph,opts:OptionsPattern[BarChart]]:=
Module[{δ=VertexDegree[g],degreeCounts,total},
degreeCounts=Counts@δ;
degreeCounts=KeySort@Join[AssociationThread[(Range@@MinMax@Keys@degreeCounts)->0],degreeCounts];
total=Total@degreeCounts;

BarChart[Values@degreeCounts/total,
Axes->False,
Frame->{{True,False},{True,False}},
FrameLabel->{"vertex degree","fraction (%)"},ChartLabels->Keys@degreeCounts,
LabelingFunction->(Placed[total #/.{0->""},Above]&),
opts]
]