Using Catalyst ODEproblem and algebraic equations

Specific Domains Chemistry
,
1

Dear all,
I am trying to add algebraic equations to a Catalyst reaction system. I read old post on this subject where It was not possible to do with ODEproblem function. Is it still the case? For example, I’ve tried this small example :

using ModelingToolkit: t_nounits as t, D_nounits as D

@parameters k
@species CaOH₂(t) Ca²⁺(t) H⁺(t) ℍ(t)
diss_A = @reaction k, CaOH₂ + 2H⁺ => Ca²⁺ + 2ℍ
@named sys = ReactionSystem([diss_A], t)

sys_complete = complete(sys)

odesys = convert(ODESystem, sys_complete)

@variables OH⁻(t)
eqs = [
    equations(odesys)...,
    OH⁻ ~ 1e-14 / H⁺  # Équation algébrique
]

daesys = ODESystem(eqs, t; name=:extended_system)

#simpsys = structural_simplify(daesys)
simpsys = complete(daesys)

u0 = [CaOH₂ => 1.0, Ca²⁺ => 0.0, H⁺ => 1e-7, ℍ => 1.0, OH⁻ => 1e-7]  # Car B(0) = C(0) + 1
tspan = (0.0, 10.0)
params = [k => 1.0]

prob = ODAEProblem(simpsys, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, vars=[CaOH₂, Ca²⁺, H⁺, ℍ, OH⁻], labels=["CaOH₂" "Ca²⁺" "H⁺" "ℍ" "OH⁻"])```
However, I get the error message:
The LHS cannot contain nondifferentiated variables

If I use structural_simplify, algebraic equations are removed.
Is there any other way to do that?

Thanks
2

Good news, in Catalyst 15 you can just write your full model in the DSL. (And I assume you wanted mass action kinetics so that the rate law is k * \text{CaOH₂} * (\text{H⁺})^2 and not k, so I’ve changed the arrow in your model and turned off combinatoric ratelaws, which if left on use k * \text{CaOH₂} * (\text{H⁺})^2/2.)

using Catalyst, OrdinaryDiffEqTsit5, Plots

# note I use --> to get the mass action rate law
rn = @reaction_network begin
          @combinatoric_ratelaws false
           k, CaOH₂ + 2H⁺ --> Ca²⁺ + 2ℍ
           @equations begin
               OH⁻ ~ 1e-14 / H⁺
           end
       end
osys = structural_simplify(convert(ODESystem, rn))

# note, you shouldn't give an initial value for the algebraic variable, OH^-, 
# as that will make the system overdetermined -- ModelingToolkit will solve 
# for its initial value from the algebraic equation
# I also changed some of these values to get a non-trivial plot
u0 = [:CaOH₂ => 1.0, :Ca²⁺ => 0.0, :H⁺ => 1.0, :ℍ => 1.0]
tspan = (0.0, 10.0)
params = [:k => 1.0]
oprob = ODEProblem(osys, u0, tspan, params)
sol = solve(oprob, Tsit5())
plot(sol, idxs = [:CaOH₂, :Ca²⁺, :H⁺, :ℍ, :OH⁻])
3 Likes
3

Thank you so much. This is effectively a great news to write all the system in the DSL.
Yes, of course, in my example k is not constant. It depends on several parameters and species.

k_CaOH₂ * A_CaOH₂ * abs(1 - (a_Ca²⁺ / (a_H⁺)^2 / K_CaOH₂)) * (CaOH₂ >= 0) * (H⁺ >= 0),  CaOH₂ + 2H⁺ => Ca²⁺ + 2ℍ

where k_CaOH₂ is the rate constant, A the specific reactive surface (which depends on CaOH₂), a the activity and K the equilibrium contant. That was the reason why I noticed a “simple” arrow.

4

Thank you also for your responsiveness and the improvements you made to the code.

5

Glad that helped, and that the arrow use was intentional. Sometimes people use it not realizing it disables mass action type kinetics.