Time-varying parameter specification in DiscreteProblem?

Specific Domains Modelling & Simulations
1

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)
2

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

3

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
4

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