Hey folks,
I am trying to build a compartmental model, where some parameters are supposed to follow a SDE. But somehow can’t get the composition right.
I build the following example to illustrate the problem. First using ODEs.
using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t
using DifferentialEquations
# Define a variable infection rate parameter
@variables β(t) = 0.0
# Define Species
@species begin
S(t) = 1000
I(t) = 1
R(t) = 5
U(t) = 1
C(t) = 9999
end
@parameters begin
b = 0.00082
μi = 0.0000421
μd = 0.0000003245
d = 0.01
τ = 0.2
δ = 0.4
p = 0.4
α = 0.4
ρ = 0.9
θ = 0.8
ϕ = 100
K = 1000
bcμc = 0.00033
β_k1 = 0.4
β_k2 = 0.2
end
# Define the Reactions
#Reactions in region 1
rxs = [
Reaction(b*(S+I+R+U), nothing, [S]), # Birth: ∅ → S
Reaction(μi+μd*(S+I+R+U), [S], nothing), # Natural death: S → ∅
Reaction((1 - p) * exp(β) * (C / (C + K)), [S], [I]), # Infection: S → I
Reaction(p * α, [S], [U]), # Transition S → U
Reaction(μi+μd*(S+I+R+U), [I], nothing), # Natural death: I → ∅
Reaction(d, [I], nothing), # Disease induced death: I → ∅
Reaction(τ, [I], [R]), # Recovery: I → R
Reaction(μi+μd*(S+I+R+U), [R], nothing), # Natural death: R → ∅
Reaction((1 - p) * ρ, [R], [S]), # Transition: R → S
Reaction(ρ*p, [R], [U]), # Transition: R → U
Reaction(μi+μd*(S+I+R+U), [U], nothing), # Natural death: U → ∅
Reaction(δ, [U], [S]), # Transition: U → S
Reaction((1 - θ) * ϕ * I, nothing, [C]), # Bacteria from Infection I → C
Reaction(bcμc*C, nothing, [C]), # Birth and Death of bacteria ∅ → C
]
# Define the Reaction System
@named sir_simple = ReactionSystem(rxs, t)
cpl_sir_simple = complete(sir_simple)
@parameters begin
betak = 1.0
betamu = 0.0
betaeps = 0.5
end
connect_ode_sys = compose(
ODESystem([D(sir_simple.β) ~ betak * (betamu - sir_simple.β)],
t,
name = :connect_ode_sys),
cpl_sir_simple)
connect_ode_prob = ODEProblem(complete(connect_ode_sys), [], (0.0, 10.0), [])
ode_sol = solve(connect_ode_prob)
plot(ode_sol,
idxs = [2, 3, 4, 5]
)
This runs without any error and produces the desired output plot.
However, if I try to do vary the parameter with an SDE and have everything else as an SDESystem as well
connect_sys = compose(
SDESystem([D(sir_simple.β) ~ betak * (betamu - sir_simple.β)],
[betaeps],
t,
name = :connect_sys,
[sir_simple.β],
[betak, betamu, betaeps]),
cpl_sir_simple
)
produces an SDESystem with the correct 6 equations. But if I want to transform this into an SDEProblem
cpl_sdeprob = SDEProblem(complete(connect_sys), [], (0.0, 10.0), [])
I get the following error
ERROR: AssertionError: length(eqs) == length(eq_domain)
Stacktrace:
[1] substitute_sample_time(ci::ModelingToolkit.ClockInference{TearingState{SDESystem}}, ts::TearingState{SDESystem})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/clock_inference.jl:31
[2] infer_clocks!(ci::ModelingToolkit.ClockInference{TearingState{SDESystem}})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/clock_inference.jl:115
[3] process_DEProblem(constructor::Type, sys::SDESystem, u0map::Vector{…}, parammap::Vector{…}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, eval_module::Module, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{…}, t::Nothing, warn_initialize_determined::Bool, build_initializeprob::Bool, initialization_eqs::Vector{…}, kwargs::@Kwargs{…})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/abstractodesystem.jl:778
[4] process_DEProblem
@ ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/abstractodesystem.jl:741 [inlined]
[5] (SDEProblem{…})(sys::SDESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…}; sparsenoise::Nothing, check_length::Bool, callback::Nothing, kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/sdesystem.jl:611
[6] (SDEProblem{…})(sys::SDESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/sdesystem.jl:603
[7] (SDEProblem{true})(::SDESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/sdesystem.jl:658
[8] (SDEProblem{true})(::SDESystem, ::Vector{Any}, ::Vararg{Any})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/sdesystem.jl:657
[9] SDEProblem(::SDESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/sdesystem.jl:647
[10] SDEProblem(::SDESystem, ::Vector{Any}, ::Vararg{Any})
I am not even sure where the eq_domain is defined and what it refers to.
Maybe someone has some idea what the error exactly means and whether there is a better way to insert an SDE for a time-varying parameter into an existing model.
Thanks.