# Time-varying parameter specification in DiscreteProblem?

I have a simple DiscreteProblem, based on the SIR epidemic model. I’d like to pass a time-varying parameter to the model, to reduce \beta by a factor v(t). As this is a discrete-time map, I’d like to define v as a vector and just retrieve value at a discrete-time index, as DataInterpolations.jl seems like overkill. Any suggestions on how best to implement this?

using OrdinaryDiffEq
using Plots

@inline function rate_to_proportion(r::Float64,t::Float64)
1-exp(-r*t)
end

function sir_map!(unp1,un,p,t)
(S,I,C) = un
(β,γ,δt) = p
infection = rate_to_proportion(β*I,δt)*S
recovery = rate_to_proportion(γ,δt)*I
@inbounds begin
unp1[1] = S-infection
unp1[2] = I+infection-recovery
unp1[3] = C+infection
end
nothing
end

tspan = (0.0, 40.0)
u0 = [0.99,0.01,0.0] # S,I,C
p = [0.5,0.25,0.1] # β,γ,δt
prob = DiscreteProblem(sir_map!,u0,tspan,p)
sol = solve(prob,FunctionMap(),dt=p[3])
plot(sol)


You can @symbolic_register a function that takes t and returns the correct value

Here are some ideas:

1. I would simply pass the vector v as a parameter, e.g.,
v = exp.(.- (0:0.1:40) ./ 10)  # Note: Covers same times as your simulation
p = (0.5,0.25,1//10,v)  # β₀,γ,δt,v

2. In the function, access v at the right place, i.e.,
• either converting time to index
• or running the simulation in discrete steps in the first place

Here is code illustrating both options:

@inline function rate_to_proportion(r::Real,t::Real)  # Float not general enough for Rational time
1-exp(-r*t)
end

function sir_map_v1!(unp1,un,p,t)
(S,I,C) = un
(β₀,γ,δt,v) = p
β = β₀ * v[Int(t // δt) + 1]  # Note: Rational to prevent rounding issues
... # cut some lines here
end

function sir_map_v2!(unp1,un,p,t)
(S,I,C) = un
(β₀,γ,δt,v) = p
β = β₀ * v[t]  # Note: t is actually steps here
... # cut some lines here
end

# Run option 1
tspan = (0//1, 40//1)
p = (0.5,0.25,1//10,v)  # Rational δt
prob = DiscreteProblem(sir_map_v1!,u0,tspan,p)
sol = solve(prob,FunctionMap(),dt=p[3])

# Run option 2
tspan = (1, 401)  # Discrete steps matching indices of v
p = (0.5,0.25,0.1,v)
prob = DiscreteProblem(sir_map_v2!,u0,tspan,p)
sol = solve(prob,FunctionMap(),dt=1)  # discrete stepping, i.e, dt isa Int


Thanks! v2 is cleaner for my purposes. I didn’t know that specifying tspan and dt would force DiscreteProblem to take integer steps.