# Compose SDESystems in ModelingToolkit

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.

Can you open an issue so this gets tracked?

I simplified the MWE and submitted an issue here.