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.009561650556141641 … 1.7850304025365167 1.806542381481639]
 [0.0 0.009562161356191774 … 1.8346441116053425 1.85601506643699]
 [0.0 0.00957071911536945 … 2.063185798329754 2.0838800170790077]
 [0.0 0.009578729937247248 … 2.1052338293556505 2.1258023905718413]
 [0.0 0.009562612640450332 … 1.868826238485931 1.8900977028157497]
 [0.0 0.009562613369643841 … 1.868875768439198 1.8901470877095112]
 [0.0 0.00957678360043339 … 2.0982967315016836 2.1188860311403013]
 [0.0 0.009576477793956107 … 2.0970753346672946 2.117668286621633]
 [0.0 0.009560591757343988 … 1.6137053315391148 1.6356170666143244]
 [0.0 0.009570713786742168 … 2.0631402496431726 2.083834604506023]
 ⋮
 [0.0 0.00956807329775445 … 2.0346190821064214 2.0553985793124085]
 [0.0 0.009574584478032638 … 2.0885190501301314 2.109137573202259]
 [0.0 0.009572960545610373 … 2.079465994331637 2.100111574144705]
 [0.0 0.009565831738045398 … 1.9951458794369215 2.0160431058746666]
 [0.0 0.009570203099899652 … 2.058590136749757 2.079298076837962]
 [0.0 0.009575317747231653 … 2.0920532636726725 2.112661225705908]
 [0.0 0.009569247358926874 … 2.0489615441575815 2.069698230382221]
 [0.0 0.009564629038790731 … 1.9627294997330313 1.9837232728512775]
 [0.0 0.009571837507075593 … 2.0719670256758995 2.092635011992819]

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.009559611468855997; … ; 1.2389230392537; 1.2601008170535304;;]
 [0.0; 0.009577925591252229; … ; 2.102522373942654; 2.1230990398113856;;]

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.