Algebraic equations/observables with HybridProblems

[kaandocal]

Dear Catalyst team,

I'm trying to cook up a simple model of the E. coli tetracycline response and am running into a small issue trying to define observables/pseudo-species. Take the following system:

using Catalyst

Tc_input(t) = 100 * (t >= 30)

rn = @reaction_network begin 
    @species R(t) 
    @variables Tc(t)
    @parameters D_Tc α K_Tc h_R K_R Tc_ext(::Real)::Real
    
    @equations begin
        D(Tc) ~ D_Tc * (Tc_ext(t) - Tc) 
        R_f ~ R / (1 + Tc/K_Tc)
    end

    h_R / (1 + R_f/K_R), 0 --> R
end 


u0 = [ :R => 0, :Tc => 0 ]
p = [  :D_Tc => 0.015,  :α => 1,  :K_Tc => 0.1,  :h_R => 0.4,  :K_R => 2, :Tc_ext => Tc_input ]
tspan = (0, 180.)

HereTc is a variable, R is a species and R_f is a variable that is a function of R and Tc. This compiles to an ODEProblem just fine:

oprob = ODEProblem(rn, u0, tspan, p; mtkcompile=true)

But trying to convert this to a HybridProblem fails because of the algebraic equation:

hprob = HybridProblem(rn, u0, tspan, p; mtkcompile=true)

ERROR: Initial condition underdefined. Some are missing from the variable map.
Please provide a default (`u0`), initialization equation, or guess
for the following variables:

R_f(t)

Stacktrace:
  [1] varmap_to_vars(varmap::ModelingToolkitBase.AtomicArrayDict{…}, vars::Vector{…}; tofloat::Bool, use_union::Bool, container_type::Type, buffer_eltype::Type, toterm::Function, check::Bool, allow_symbolic::Bool, is_initializeprob::Bool, substitution_limit::Int64, missing_values::ModelingToolkitBase.MissingGuessValue.var"typeof(MissingGuessValue)")
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/systems/problem_utils.jl:390
  [2] varmap_to_vars
    @ ~/.julia/packages/ModelingToolkitBase/QKi0f/src/systems/problem_utils.jl:334 [inlined]
  [3] __process_SciMLProblem(constructor::Type, sys::System, op::Dict{…}; floatT::Type, u0Type::Type, u0_eltype::Type, build_initializeprob::Bool, implicit_dae::Bool, t::Int64, guesses::Dict{…}, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, eval_expression::Bool, eval_module::Module, fully_determined::Nothing, check_initialization_units::Bool, tofloat::Bool, u0_constructor::typeof(identity), p_constructor::typeof(identity), check_length::Bool, symbolic_u0::Bool, warn_cyclic_dependency::Bool, circular_dependency_max_cycle_length::Int64, circular_dependency_max_cycles::Int64, initsys_mtkcompile_kwargs::@NamedTuple{}, substitution_limit::Int64, use_scc::Bool, time_dependent_init::Bool, algebraic_only::Bool, missing_guess_value::ModelingToolkitBase.MissingGuessValue.var"typeof(MissingGuessValue)", allow_incomplete::Bool, is_initializeprob::Bool, is_steadystateprob::Bool, return_operating_point::Bool, compiler_options::CompilerOptions, init_compiler_options::CompilerOptions, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/systems/problem_utils.jl:1888
  [4] process_SciMLProblem(constructor::Any, sys::System, op::Any; u0_eltype::Nothing, u0_constructor::Function, p_constructor::Function, symbolic_u0::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/systems/problem_utils.jl:1782
  [5] (ODEProblem{…})(sys::System, op::Any, tspan::Tuple{…}; callback::Any, check_length::Bool, eval_expression::Bool, expression::Type, eval_module::Module, check_compatibility::Bool, _skip_events::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/odeproblem.jl:116
  [6] ODEProblem
    @ ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/odeproblem.jl:105 [inlined]
  [7] (JumpProblem{…})(sys::System, op::Dict{…}, tspan::Tuple{…}; check_compatibility::Bool, eval_expression::Bool, eval_module::Module, checkbounds::Bool, cse::Bool, aggregator::JumpProcesses.NullAggregator, callback::Nothing, rng::Nothing, save_positions::Tuple{…}, kwargs::@Kwargs{})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/jumpproblem.jl:27
  [8] JumpProblem
    @ ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/jumpproblem.jl:1 [inlined]
  [9] JumpProblem
    @ ./none:-1 [inlined]
 [10] JumpProblem
    @ ./none:-1 [inlined]
 [11] HybridProblem(rs::ReactionSystem{…}, u0::Vector{…}, tspan::Tuple{…}, p::Vector{…}; name::Symbol, physical_scales::Nothing, default_scale::Catalyst.PhysicalScale.T, combinatoric_ratelaws::Bool, expand_catalyst_funs::Bool, save_positions::Tuple{…}, checks::Bool, mtkcompile::Bool, use_jump_ratelaws::Bool, kwargs::@Kwargs{})
    @ Catalyst ~/.julia/packages/Catalyst/MJhlE/src/reactionsystem_conversions.jl:1458
 [12] top-level scope
    @ ~/code/src/model.jl:108
Some type information was truncated. Use `show(err)` to see complete types.

Sounds like I have to provide an initial guess for R_f. But then I get an error because of the algebraic constraint:

ERROR: The LHS cannot contain nondifferentiated variables. Please run `mtkcompile` or use the DAE form.\nGot R_f(t) ~ R(t) / (1 + Tc(t) / K_Tc)
Stacktrace:
  [1] error(s::LazyString)
    @ Base ./error.jl:44
  [2] check_operator_variables(eqs::Vector{Equation}, ::Type{Differential})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/utils.jl:595
  [3] generate_rhs(sys::System; implicit_dae::Bool, scalar::Bool, expression::Type, wrap_gfw::Type, eval_expression::Bool, eval_module::Module, override_discrete::Bool, cachesyms::Nothing, compiler_options::CompilerOptions, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/systems/codegen.jl:85
  [4] (ODEFunction{…})(sys::System; u0::Any, p::Any, tgrad::Bool, jac::Bool, t::Int64, eval_expression::Bool, eval_module::Module, sparse::Bool, steady_state::Bool, checkbounds::Bool, sparsity::Bool, analytic::Any, simplify::Bool, cse::Bool, initialization_data::Any, expression::Type, check_compatibility::Bool, nlstep::Bool, nlstep_compile::Bool, nlstep_scc::Bool, optimize::Nothing, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/odeproblem.jl:32
  [5] __process_SciMLProblem(constructor::Type, sys::System, op::Dict{…}; floatT::Type, u0Type::Type, u0_eltype::Type, build_initializeprob::Bool, implicit_dae::Bool, t::Int64, guesses::Dict{…}, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, eval_expression::Bool, eval_module::Module, fully_determined::Nothing, check_initialization_units::Bool, tofloat::Bool, u0_constructor::typeof(identity), p_constructor::typeof(identity), check_length::Bool, symbolic_u0::Bool, warn_cyclic_dependency::Bool, circular_dependency_max_cycle_length::Int64, circular_dependency_max_cycles::Int64, initsys_mtkcompile_kwargs::@NamedTuple{}, substitution_limit::Int64, use_scc::Bool, time_dependent_init::Bool, algebraic_only::Bool, missing_guess_value::ModelingToolkitBase.MissingGuessValue.var"typeof(MissingGuessValue)", allow_incomplete::Bool, is_initializeprob::Bool, is_steadystateprob::Bool, return_operating_point::Bool, compiler_options::CompilerOptions, init_compiler_options::CompilerOptions, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/systems/problem_utils.jl:1944
  [6] process_SciMLProblem(constructor::Any, sys::System, op::Any; u0_eltype::Nothing, u0_constructor::Function, p_constructor::Function, symbolic_u0::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/systems/problem_utils.jl:1782
  [7] (ODEProblem{…})(sys::System, op::Any, tspan::Tuple{…}; callback::Any, check_length::Bool, eval_expression::Bool, expression::Type, eval_module::Module, check_compatibility::Bool, _skip_events::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/odeproblem.jl:116
  [8] ODEProblem
    @ ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/odeproblem.jl:105 [inlined]
  [9] (JumpProblem{…})(sys::System, op::Dict{…}, tspan::Tuple{…}; check_compatibility::Bool, eval_expression::Bool, eval_module::Module, checkbounds::Bool, cse::Bool, aggregator::JumpProcesses.NullAggregator, callback::Nothing, rng::Nothing, save_positions::Tuple{…}, kwargs::@Kwargs{})
    @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/jumpproblem.jl:27
 [10] JumpProblem
    @ ~/.julia/packages/ModelingToolkitBase/QKi0f/src/problems/jumpproblem.jl:1 [inlined]
 [11] JumpProblem
    @ ./none:-1 [inlined]
 [12] JumpProblem
    @ ./none:-1 [inlined]
 [13] HybridProblem(rs::ReactionSystem{…}, u0::Vector{…}, tspan::Tuple{…}, p::Vector{…}; name::Symbol, physical_scales::Nothing, default_scale::Catalyst.PhysicalScale.T, combinatoric_ratelaws::Bool, expand_catalyst_funs::Bool, save_positions::Tuple{…}, checks::Bool, mtkcompile::Bool, use_jump_ratelaws::Bool, kwargs::@Kwargs{})
    @ Catalyst ~/.julia/packages/Catalyst/MJhlE/src/reactionsystem_conversions.jl:1458
 [14] top-level scope
    @ ~/code/eric/src/model.jl:108
Some type information was truncated. Use `show(err)` to see complete types.

I could change R_f to be an observable, but that it limiting and I run into other issues. What is the best way to do this with HybridProblems?

[ChrisRackauckas]

Yes that value needs a guess.

[isaacsas]

Even ignoring that, I don’t think there are working solvers for DAEs + jumps? So wouldn’t this require MTK to run simplification on models with jumps to try to eliminate all algebraic equations (which it doesn’t do)?

[ChrisRackauckas]

The DAE solvers should work fine.

[isaacsas]

I don't think there are any correctness tests for DAEs + jumps. So one should probably carefully test that if they run through the results are correct.

[kaandocal]

Can MTK simplify this system simply by eliding the variable R_f? I can replace R_f by its value in the model, but it would be nice to have observables be a bit more powerful, e.g. being able to nest them!

[isaacsas]

MTK could in theory extend simplification to include jumps, but currently I don't believe they are supported at all. Someone would have to dig into the simplification code in MTKBase to add support, and I don't think it would be trivial (since you'd need to handle updating jump expressions, and possibly to convert jump types, i.e. from mass action to constant rate, if they depend on an elided variable).

[ChrisRackauckas]

Why would it have to extend? It would just carry forward the state information and observed, which it should be doing already?

[isaacsas]

If you move a state to be an observable and it is also used in a mass action jump then MTK would need to convert the mass action jump to a constant rate jump that can properly access values of observables. Mass action jumps can only access state variables within u.

[isaacsas]

But yes, if DAEs + jump solvers work then I guess it should go through if the system isn't modified (Catalyst should already construct the jump with the appropriate type, so if the system is not modified/reduced in MTK I think it should be ok for such a solver).