# Compose SDESystems in ModelingToolkit

**URL:** <https://discourse.julialang.org/t/compose-sdesystems-in-modelingtoolkit/117054>\
**Category:** Modelling & Simulations\
**Tags:** question, modelingtoolkit\
**Created:** [July 15, 2024, 1:38pm UTC](https://discourse.julialang.org/t/compose-sdesystems-in-modelingtoolkit/117054 "2024-07-15T13:38:08Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [July 15, 2024, 1:38pm UTC](https://discourse.julialang.org/t/compose-sdesystems-in-modelingtoolkit/117054/1 "2024-07-15T13:38:08Z")

</div>

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.

```julia
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.  
 ![Screenshot from 2024-07-15 15-33-22](https://global.discourse-cdn.com/julialang/original/3X/c/7/c70cca5ebeb0a10a9f3ea0ffceeb3413410e768b.png)

However, if I try to do vary the parameter with an SDE and have everything else as an `SDESystem` as well

```julia

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`

```julia
cpl_sdeprob = SDEProblem(complete(connect_sys), [], (0.0, 10.0), [])

```

I get the following error

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 23, 2024, 9:07pm UTC](https://discourse.julialang.org/t/compose-sdesystems-in-modelingtoolkit/117054/2 "2024-07-23T21:07:53Z")

</div>

Can you open an issue so this gets tracked?

---

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [July 24, 2024, 8:54am UTC](https://discourse.julialang.org/t/compose-sdesystems-in-modelingtoolkit/117054/3 "2024-07-24T08:54:10Z")

</div>

I simplified the MWE and submitted an issue [here](https://github.com/SciML/ModelingToolkit.jl/issues/2893).