Building a `SDEProblem` from a composed `SDESystems` does not work.

[vwiela]

I am trying to introduce a time-varying parameter following an SDE into an SDESystem obtained from Catalyst.jl ReactionSystem. However, when I want to create an SDEProblem from the composed SDESystem, I get some error about mismatching dimensions.

I created an MWE using a simple SIR-Model.

First, doing it in terms of ODEs works fine.

using DifferentialEquations
using ModelingToolkit
using Catalyst
using Plots
using StatsPlots
import ModelingToolkit: t_nounits as t, D_nounits as D

@variables begin
      β(t) = 0.01
end

@parameters begin
      γ = 0.1
end

@species begin
      S(t) = 100
      I(t) = 1
      R(t) = 0
end

rxs = [Reaction(β, [S,I], [I], [1,1],[2]),
       Reaction(γ, [I], [R])]

@named sir_simple = ReactionSystem(rxs, t)
sir_simple=complete(sir_simple)

sir_simple_ode = convert(ODESystem, sir_simple)

@parameters begin
      betak = 1.0
      betamu = 0.01
      betaeps = 0.5
end

connect_ode_sys = compose(
      ODESystem([D(sir_simple_ode.β) ~ betak * (betamu - sir_simple_ode.β)],
            t,
            name = :connect_ode_sys),
      sir_simple_ode)

connect_ode_prob = ODEProblem(complete(connect_ode_sys), [], (0.0, 10.0), [])

ode_sol = solve(connect_ode_prob, Tsit5())

plot(ode_sol,
      idxs=[sir_simple_ode.S sir_simple_ode.I sir_simple_ode.R],)

with expected output plot

However, doing the same in terms of SDESystem produces the following error.

sir_simple_sde = convert(SDESystem, sir_simple)

connect_sde_sys = compose(
      SDESystem([D(sir_simple_sde.β) ~ betak * (betamu - sir_simple_sde.β)],
            [betaeps],
            t,
            [],
            [],
            name = :connect_sde_sys),
      sir_simple_sde
)

connect_sde_prob = SDEProblem(complete(connect_sde_sys), [], (0.0, 10.0), [])

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})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/diffeqs/sdesystem.jl:646
 [11] top-level scope
    @ ~/julia_models/SIRModel.jl:67

Unfortunately, from the error messages I could not really deduce the problem.

Would appreciate any help and tips, whether this is the correct way to compose such two SDESystems to have a parameter folowing an SDE.

I am aware of the workaround in #2085, but my original system is a rather large reaction network consisting of a composition of several network_components built in Catalyst. So I do not want to write down all the equations by hand and adding diffusion terms and Brownian Motions to it by hand.

PS: Also structural_simplify does not work on the composed system and would produce the following

connect_sde = structural_simplify(connect_sde_sys)

ERROR: BoundsError: attempt to access 7-element Vector{Vector{Int64}} at index [8]
Stacktrace:
  [1] getindex
    @ ./essentials.jl:13 [inlined]
  [2] 𝑠neighbors (repeats 2 times)
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/bipartite_graph.jl:364 [inlined]
  [3] find_eq_solvables!(state::TearingState{…}, ieq::Int64, to_rm::Vector{…}, coeffs::Vector{…}; may_be_zero::Bool, allow_symbolic::Bool, allow_parameter::Bool, conservative::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/U2JaS/src/structural_transformation/utils.jl:194
  [4] find_eq_solvables!(state::TearingState{SDESystem}, ieq::Int64, to_rm::Vector{Int64}, coeffs::Vector{Int64})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/U2JaS/src/structural_transformation/utils.jl:182
  [5] linear_subsys_adjmat!(state::TearingState{SDESystem}; kwargs::@Kwargs{})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/U2JaS/src/structural_transformation/utils.jl:267
  [6] linear_subsys_adjmat!
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/structural_transformation/utils.jl:255 [inlined]
  [7] alias_eliminate_graph!(state::TearingState{SDESystem}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/alias_elimination.jl:5
  [8] alias_eliminate_graph!
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/alias_elimination.jl:4 [inlined]
  [9] alias_elimination!(state::TearingState{SDESystem}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/alias_elimination.jl:50
 [10] alias_elimination!
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/alias_elimination.jl:46 [inlined]
 [11] _structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systemstructure.jl:689
 [12] _structural_simplify!
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systemstructure.jl:674 [inlined]
 [13] #structural_simplify!#1096
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systemstructure.jl:667 [inlined]
 [14] __structural_simplify(sys::SDESystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systems.jl:83
 [15] __structural_simplify
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systems.jl:64 [inlined]
 [16] structural_simplify(sys::SDESystem, io::Nothing; simplify::Bool, split::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systems.jl:24
 [17] structural_simplify
    @ ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systems.jl:21 [inlined]
 [18] structural_simplify(sys::SDESystem)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/U2JaS/src/systems/systems.jl:21
 [19] top-level scope
    @ ~/julia_models/SIRModel.jl:69

Environment:

I used the following package versions and julia 1.10.4

[479239e8] Catalyst v14.1.0
[0c46a032] DifferentialEquations v7.13.0
[961ee093] ModelingToolkit v9.24.0
[91a5bcdd] Plots v1.40.5
[f3b207a7] StatsPlots v0.15.7
[789caeaf] StochasticDiffEq v6.66.0
[d1185830] SymbolicUtils v2.1.0
[0c5d862f] Symbolics v5.33.0

[ChrisRackauckas]

@YingboMa why are clocks in the stack trace at all?