Hi All,
I am playing with a basic SIR disease model and would like to make the transmission
rate (Beta) depend on the rate at which the number of disease-related deaths changes.
(e.g., when few deaths are occurring transmission rate is high because there are no social distancing precautions in place. When deaths start to occur rapidly, transmission becomes low because greater social distancing. Then when deaths become less frequent transmission rate becomes high again due to relaxing of social distancing measures).
However, I can’t figure out how to make Beta depend on the rate of deaths.
Here’s an example (that doesn’t quite do what I want):
using DifferentialEquations, Plots
#The transmission rate Beta (du[5]) is (supposed to be) a function of the rate of change of deaths (i.e., of dD/dt)
function SIR(du,u,p,t)
du[1] = dS = p[1]*u[1] - p[2]*u[1] - u[5]*u[1]*u[2]
du[2] = dI = u[5]*u[1]*u[2] - p[2]*u[2] - p[3]*u[2] - p[4]*u[2]
du[3] = dR = p[4]*u[2] - p[2]*u[3]
du[4] = dD = p[3]*u[2]
du[5] = dβ = -du[4]*p[5]
end
const β = 1e-3
u0 = [100.0, 1.0, 0.0, 0.0, β]
tspan = (0.0, 365.0)
p = [1/(30*365), 1/(85*365), 1/100, 1/40, β]
prob = ODEProblem(SIR, u0, tspan, p)
saved_values = SavedValues(Float64, Array{Float64})
cb = SavingCallback((u,t,integrator)->integrator(t,Val{1}), saved_values, saveat = 0.0:1:365.0)
sol = solve(prob, Tsit5(), saveat=1.0,callback=cb)
As you can see below, when I plot how Beta changes over time by plotting the solution
for du[5] I get a monotonic curve rather than the unimodal curve
that I expected from including du[4] in the equation:
# plot(sol,xlabel="t",ylabel="#")
plot(sol[5,:],xlabel="t",ylabel="#",label="Beta")
I expected that du[4] would give me the derivatives, similar to the callback below:
dDdt = [saved_values.saveval[i][4] for i in 2:length(sol)]
tD = saved_values.t[2:end]
plot(tD,-dDdt,xlabel="t",label="Beta")
This is all a long-winded way of asking why du[i] doesn’t give me the same thing
as integrator(t,Val{1}).
I hope what I’m asking is clear. Any help would be appreciated!