# Using Catalyst ODEproblem and algebraic equations

**URL:** <https://discourse.julialang.org/t/using-catalyst-odeproblem-and-algebraic-equations/128375>\
**Category:** Chemistry\
**Tags:** question, package\
**Created:** [April 24, 2025, 2:26pm UTC](https://discourse.julialang.org/t/using-catalyst-odeproblem-and-algebraic-equations/128375 "2025-04-24T14:26:50Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![AnthonySoive](https://avatars.discourse-cdn.com/v4/letter/a/a88e57/32.png) [@AnthonySoive](https://discourse.julialang.org/u/AnthonySoive)\
**Post date:** [April 24, 2025, 2:26pm UTC](https://discourse.julialang.org/t/using-catalyst-odeproblem-and-algebraic-equations/128375/1 "2025-04-24T14:26:50Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [April 24, 2025, 11:32pm UTC](https://discourse.julialang.org/t/using-catalyst-odeproblem-and-algebraic-equations/128375/2 "2025-04-24T23:32:03Z")

</div>

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.)

```julia
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⁻])

```

---

<div class="post-metadata">

**Author:** ![AnthonySoive](https://avatars.discourse-cdn.com/v4/letter/a/a88e57/32.png) [@AnthonySoive](https://discourse.julialang.org/u/AnthonySoive)\
**Post date:** [April 25, 2025, 12:23pm UTC](https://discourse.julialang.org/t/using-catalyst-odeproblem-and-algebraic-equations/128375/3 "2025-04-25T12:23:11Z")

</div>

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.

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

---

<div class="post-metadata">

**Author:** ![AnthonySoive](https://avatars.discourse-cdn.com/v4/letter/a/a88e57/32.png) [@AnthonySoive](https://discourse.julialang.org/u/AnthonySoive)\
**Post date:** [April 25, 2025, 12:37pm UTC](https://discourse.julialang.org/t/using-catalyst-odeproblem-and-algebraic-equations/128375/4 "2025-04-25T12:37:55Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [April 25, 2025, 12:50pm UTC](https://discourse.julialang.org/t/using-catalyst-odeproblem-and-algebraic-equations/128375/5 "2025-04-25T12:50:33Z")

</div>

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