Analysis of Interventions On The SEIRHD Epidemic Model

First, let's implement the classic SEIRHD epidemic model with ModelingToolkit:

using EasyModelAnalysis
@variables t
Dₜ = Differential(t)
@variables S(t)=0.9 E(t)=0.05 I(t)=0.01 R(t)=0.2 H(t)=0.1 D(t)=0.01
@variables T(t)=0.0 η(t)=0.0 cumulative_I(t)=0.0
@parameters β₁=0.6 β₂=0.143 β₃=0.055 α=0.003 γ₁=0.007 γ₂=0.011 δ=0.1 μ=0.14
eqs = [T ~ S + E + I + R + H + D
    η ~ (β₁ * E + β₂ * I + β₃ * H) / T
    Dₜ(S) ~ -η * S
    Dₜ(E) ~ η * S - α * E
    Dₜ(I) ~ α * E - (γ₁ + δ) * I
    Dₜ(cumulative_I) ~ I
    Dₜ(R) ~ γ₁ * I + γ₂ * H
    Dₜ(H) ~ δ * I - (μ + γ₂) * H
    Dₜ(D) ~ μ * H];
@named seirhd = ODESystem(eqs)
seirhd = structural_simplify(seirhd)
prob = ODEProblem(seirhd, [], (0, 110.0))
sol = solve(prob)
plot(sol)

Let's solve a few problems:

Provide a forecast of cumulative Covid-19 cases and deaths over the 6-week period from May 1 – June 15, 2020 under no interventions, including 90% prediction intervals in your forecasts. Compare the accuracy of the forecasts with true data over the six-week timespan.

get_uncertainty_forecast(prob, [cumulative_I], 0:100, [β₁ => Uniform(0.0, 1.0)], 6 * 7)

42-element Vector{Matrix{Float64}}:
 [0.0 0.009564252345864444 … 1.9499407904437587 1.970972575929712]
 [0.0 0.009569629416527661 … 2.053001086536485 2.0737257134297336]
 [0.0 0.009566415276867532 … 2.007475144599528 2.0283356336885263]
 [0.0 0.009566967961926801 … 2.0176520560361904 2.0384822131333964]
 [0.0 0.009572663061224564 … 2.0775911483562473 2.0982423338670846]
 [0.0 0.009558925179451944 … 0.7333687466788062 0.7470811083961493]
 [0.0 0.009578918846540824 … 2.1058389404809104 2.126405692201991]
 [0.0 0.009558561994882857 … 0.4641156018031798 0.4707705502425733]
 [0.0 0.009573075690805018 … 2.080171872262306 2.1008153420054883]
 [0.0 0.009569824263714772 … 2.054960357433404 2.0756791347032206]
 ⋮
 [0.0 0.009570838931048983 … 2.0641999828234487 2.0848911710506033]
 [0.0 0.009564415537710211 … 1.9556708614569496 1.9766856200733856]
 [0.0 0.009561044995116511 … 1.7027743069567096 1.7245021976341823]
 [0.0 0.009568490894947256 … 2.040090712469655 2.060853887126654]
 [0.0 0.009564615021334799 … 1.9622804232602151 1.9832755318070965]
 [0.0 0.00956942478285921 … 2.0508712042699977 2.071602189740902]
 [0.0 0.009568300171756744 … 2.037647180145349 2.0584176442421978]
 [0.0 0.0095674354006629 … 2.025317636522279 2.0461248729309127]
 [0.0 0.009560102340903591 … 1.4712120589589712 1.4932141743088414]

plot_uncertainty_forecast(prob, [cumulative_I], 0:100, [β₁ => Uniform(0.0, 1.0)], 6 * 7)

get_uncertainty_forecast_quantiles(prob, [cumulative_I], 0:100, [β₁ => Uniform(0.0, 1.0)],
    6 * 7)

2-element Vector{Matrix{Float64}}:
 [0.0; 0.009558986381174731; … ; 0.7843100691758393; 0.7991724108711555;;]
 [0.0; 0.009576794394062718; … ; 2.098335534353523; 2.118924717788651;;]

plot_uncertainty_forecast_quantiles(prob, [cumulative_I], 0:100, [β₁ => Uniform(0.0, 1.0)],
    6 * 7)

Based on the forecasts, do we need additional interventions to keep cumulative Covid deaths under 6000 total? Provide a probability that the cumulative number of Covid deaths will stay under 6000 for the next 6 weeks without any additional interventions.