Scenario 2: Limiting Hospitalizations

Generate the Model and Dataset

using EasyModelAnalysis, Optimization, OptimizationMOI, NLopt, Plots, Random
Random.seed!(12345)

@variables t
Dₜ = Differential(t)
@variables S(t)=0.97 E(t)=0.02 I(t)=0.01 R(t)=0.0 H(t)=0.0 D(t)=0.0
@variables T(t)=10000.0 η(t)=0.0 cumulative_I(t)=0.0
@parameters β₁=0.06 β₂=0.015 β₃=0.005 α=0.003 γ₁=0.007 γ₂=0.001 δ=0.2 μ=0.04
eqs = [T ~ S + E + I + R + H + D
       η ~ (β₁ * E + β₂ * I + β₃ * H)
       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.0, 60.0), saveat = 1.0)
sol = solve(prob)
u60 = sol[:, end]
plot(sol)

Model Analysis

Parameterize model either using data from the previous two months (October 28th – December 28th, 2021), or with relevant parameter values from the literature.

data = [I => sol[I], R => sol[R], H => sol[H], D => sol[D]]
prior_mean = [0.06, 0.015, 0.005, 0.003, 0.007, 0.001, 0.2, 0.04]
prior_sd = [0.006, 0.0015, 0.0005, 0.0003, 0.0007, 0.0001, 0.02, 0.004]
p = [β₁, β₂, β₃, α, γ₁, γ₂, δ, μ]
p_priors = Pair.(p,
    Truncated.(Normal.(prior_mean, prior_sd), prior_mean - 3 * prior_sd,
        prior_mean + 3 * prior_sd))
tsave = collect(0.0:1.0:60.0)
fit = bayesian_datafit(prob, p_priors, tsave, data, noise_prior = InverseGamma(10, 0.1))

8-element Vector{Pair{Num, Vector{Float64}}}:
 β₁ => [0.058980512410633734, 0.060709019773107685, 0.059936197872994146, 0.059005901263172914, 0.05919372853055342, 0.05909521570248349, 0.0600656399323654, 0.061005687613304535, 0.06198160899956821, 0.05748887405234092  …  0.06367307575277482, 0.06322698620938744, 0.057443643851258816, 0.058045922952013954, 0.0606594089671135, 0.05835016816713921, 0.058509552737268794, 0.05999524068735448, 0.059017618705764614, 0.05806689660981153]
 β₂ => [0.015142496318686138, 0.016865564810250793, 0.013854987359653242, 0.01430583825709135, 0.015385787693635748, 0.013901361142549013, 0.016800204450601916, 0.016674289008912317, 0.013428243107202444, 0.01607609588907482  …  0.013301892335897555, 0.017097189159200082, 0.011037654690079834, 0.012815702061808508, 0.018141933428469717, 0.018406153643079577, 0.012482058009669124, 0.01436744124663044, 0.013739729232611174, 0.015095973542836712]
 β₃ => [0.005303731520807747, 0.004970368572924244, 0.004376004609818733, 0.004417492318113291, 0.0047252223272287765, 0.004233510942598392, 0.005408319769722632, 0.00531510784067321, 0.0054194841916741015, 0.00499678244289874  …  0.005631571326975174, 0.004452337411085278, 0.004588077404222818, 0.005556966670840406, 0.004564071652396834, 0.004781929822737369, 0.005680925328497973, 0.0042979467868497385, 0.004557499431516839, 0.005509743395437111]
  α => [0.0030689659611853226, 0.0029293307708221005, 0.002994481099420599, 0.0030860701664391496, 0.0030746308089188512, 0.003128055690288789, 0.0029681310068029403, 0.002900196606288852, 0.002864856168387103, 0.0031975881996634337  …  0.002728646386193508, 0.0027851937759614874, 0.003263866410586809, 0.0032573127145358173, 0.002905742011287106, 0.003160800896734668, 0.003142754444368257, 0.0029862036955391755, 0.003098571819494564, 0.003175939858786616]
 γ₁ => [0.00647001077085299, 0.00696102075495229, 0.00715240446375618, 0.006771922064891163, 0.007273125394372118, 0.0067300086276892955, 0.006871021402260889, 0.006399855279191031, 0.007681490565133257, 0.0064310088739329965  …  0.006536628213925486, 0.00776707598199531, 0.007448510587106107, 0.007588225397568336, 0.006421258767237342, 0.007384119163878126, 0.006010732418939788, 0.006431905483780899, 0.006029766415477084, 0.007201998572066841]
 γ₂ => [0.0010655984904385895, 0.000986566488568434, 0.0009871370269153175, 0.0009525925570688357, 0.0011605499135209777, 0.0010608785679746642, 0.0010471177564019972, 0.0010950327116036273, 0.0009239121770743016, 0.0010287301747286743  …  0.0011133056036631297, 0.000857903598171144, 0.0010099978221655586, 0.0009673184851765511, 0.0010814247934741918, 0.0009315134854186095, 0.0010278439705771578, 0.0010064639554860485, 0.0010071233226898561, 0.0010708304571543548]
  δ => [0.20340341231722536, 0.19775025187797526, 0.20525413679893378, 0.2054372508123316, 0.19829186386815506, 0.1996907754419417, 0.20248001513307454, 0.20192547972551142, 0.20023281482046992, 0.19834496104042654  …  0.2052230196832906, 0.2032833670031723, 0.1978781658246187, 0.19681803625680655, 0.20109068752662343, 0.19675556127901336, 0.19499377604844698, 0.19863408957546694, 0.20475543103564264, 0.19349270881223396]
  μ => [0.039770153626576536, 0.04031129302443778, 0.03938864239950999, 0.040017502339673675, 0.03993530873618308, 0.04045195014214249, 0.039952437370086524, 0.039961958776015306, 0.040115106841463785, 0.03920586986880458  …  0.03990442350630812, 0.040548271048057374, 0.04055600487213483, 0.040107408388911095, 0.03969660396539464, 0.04022614865784575, 0.03992044617493397, 0.040170497563029484, 0.040294704470686504, 0.039247157481951817]

Question 1

Forecast Covid cases and hospitalizations over the next 3 months under no interventions.

prob = remake(prob; u0 = u60,
    p = Pair.(getfield.(fit, :first), mean.(getfield.(fit, :second))))
forecast_threemonths = solve(prob, tspan = (0.0, 90.0), saveat = 1.0)
plot(forecast_threemonths)

Question 2

Based on the forecast, do we need interventions to keep total Covid hospitalizations under a threshold of 3000 on any given day? If there is uncertainty in the model parameters, express the answer probabilistically, i.e., what is the likelihood or probability that the number of Covid hospitalizations will stay under this threshold for the next 3 months without interventions?

need_intervention = maximum(forecast_threemonths[H]) > 0.05

true

post_mean = mean.(getfield.(fit, :second))
post_sd = sqrt.(var.(getfield.(fit, :second)))
trunc_min = post_mean .- 3 * post_sd
trunc_max = post_mean .+ 3 * post_sd
post_trunc = Truncated.(Normal.(post_mean, post_sd), trunc_min, trunc_max)
posterior = Pair.(getfield.(fit, :first), post_trunc)
prob_violating_threshold(prob, posterior, [H > 0.05])

0.03123388506212185

Question 3

Assume a consistent policy of social distancing/masking will be implemented, resulting in a 50% decrease from baseline transmission. Assume that we want to minimize the time that the policy is in place, and once it has been put in place and then ended, it can't be re-implemented. Looking forward from “today’s” date of Dec. 28, 2021, what are the optimal start and end dates for this policy, to keep projections below the hospitalization threshold over the entire 3-month period? How many fewer hospitalizations and cases does this policy result in?

function f(ts, p = nothing)
    tstart = ts[1]
    tstop = ts[2]
    tstop - tstart
end

function g(res, ts, p = nothing)
    tstart = ts[1]
    tstop = ts[2]
    start_intervention = (t == tstart) => [β₁ ~ β₁ / 2, β₂ ~ β₂ / 2, β₃ ~ β₃ / 2]
    stop_intervention = (t == tstop) => [β₁ ~ β₁ * 2, β₂ ~ β₂ * 2, β₃ ~ β₃ * 2]
    @named opttime_sys = ODESystem(eqs, t;
        discrete_events = [
            start_intervention,
            stop_intervention
        ])
    opttime_sys = structural_simplify(opttime_sys)
    prob = ODEProblem(opttime_sys, [], [0.0, 90.0])
    prob = remake(prob; u0 = u60)
    sol = solve(prob, saveat = 0.0:1.0:90.0, tstops = [tstart, tstop])
    hospitalizations = sol(0.0:1.0:90.0, idxs = H)
    if SciMLBase.successful_retcode(sol.retcode)
        res .= vcat(hospitalizations, tstop - tstart)
    else
        res .= Inf
    end
end

optf = OptimizationFunction(f, Optimization.AutoFiniteDiff(), cons = g)
optprob = Optimization.OptimizationProblem(optf, [0.0, 90.0], lb = [0.0, 0.0],
    ub = [90.0, 90.0],
    lcons = vcat(fill(-Inf, 91), 0.0),
    ucons = vcat(fill(0.05, 91), Inf))
min_intervention_timespan = solve(optprob,
    OptimizationMOI.MOI.OptimizerWithAttributes(NLopt.Optimizer,
        "algorithm" => :GN_ORIG_DIRECT,
        "maxtime" => 60.0))
min_intervention_timespan.u

2-element Vector{Float64}:
 15.000000000000002
 35.0

res = zeros(92)
g(res, min_intervention_timespan.u)
maximum(res[1:91])

0.04995936120428489

Question 4

Assume there is a protocol to kick in mitigation policies when hospitalizations rise above 80% of the hospitalization threshold (i.e. 80% of 3000). When hospitalizations fall back below 80% of the threshold, these policies expire.

When do we expect these policies to first kick in?

What is the minimum impact on transmission rate these mitigation policies need to have the first time they kick in, to (1) ensure that we don't reach the hospitalization threshold at any time during the 3-month period, and (2) ensure that the policies only need to be implemented once, and potentially expired later, but never reimplemented? Express this in terms of change in baseline transmission levels (e.g. 10% decrease, 50% decrease, etc.).

function f(reduction_rate, p = nothing)
    reduction_rate = reduction_rate[1]
end
function g(res, reduction_rate, p = nothing)
    reduction_rate = reduction_rate[1]
    root_eqs = [H ~ 0.05 * 0.8]
    affect = [
        β₁ ~ β₁ * (1 - reduction_rate),
        β₂ ~ β₂ * (1 - reduction_rate),
        β₃ ~ β₃ * (1 - reduction_rate)
    ]
    @named mask_system = ODESystem(eqs, t; continuous_events = root_eqs => affect)
    mask_system = structural_simplify(mask_system)
    prob = ODEProblem(mask_system, [], (0.0, 90.0))
    prob = remake(prob; u0 = u60)
    sol = solve(prob, saveat = 0.0:1.0:90.0)
    hospitalizations = sol(0.0:1.0:90.0, idxs = H)
    if SciMLBase.successful_retcode(sol.retcode)
        res .= hospitalizations
    else
        res .= Inf
    end
end
optf = OptimizationFunction(f, Optimization.AutoFiniteDiff(), cons = g)
optprob = OptimizationProblem(optf, [0.0], lb = [0.0], ub = [1.0], lcons = fill(-Inf, 91),
    ucons = fill(0.05, 91))
min_intervention_strength = solve(optprob,
    OptimizationMOI.MOI.OptimizerWithAttributes(NLopt.Optimizer,
        "algorithm" => :GN_ORIG_DIRECT,
        "maxtime" => 60.0))
min_intervention_strength.u

1-element Vector{Float64}:
 0.9074074074074074

res = zeros(91)
g(res, min_intervention_strength.u)
maximum(res[1:91])

0.049689990075674755

Question 5

Now assume that instead of NPIs, the Board wants to focus all their resources on an aggressive vaccination campaign to increase the fraction of the total population that is vaccinated. What is the minimum intervention with vaccinations required in order for this intervention to have the same impact on cases and hospitalizations, as your optimal answer from question 3? Depending on the model you use, this may be represented as an increase in total vaccinated population, or increase in daily vaccination rate (% of eligible people vaccinated each day), or some other representation.

This requires an additional model structure, not very interesting to showcase the SciML stack.