# Catalyst.jl - what is the appropriate way to model equilibrium (i.e. pKa?)

**URL:** <https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149>\
**Category:** Modelling & Simulations\
**Tags:** question, package, catalyst\
**Created:** [June 23, 2025, 11:44pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149 "2025-06-23T23:44:54Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [June 23, 2025, 11:44pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/1 "2025-06-23T23:44:54Z")

</div>

Hello again,

I’m working to add some additional equilibrium conditions to a reaction network, particularly pKa.

I’m inclined to think that [algebraic constraint equations](https://docs.sciml.ai/Catalyst/stable/model_creation/constraint_equations/) are the appropriate choice here assuming that the equilibrium constraint is enforced at every timestep for a particular interaction (i.e. the equilibrium reaction is much faster than other system dynamics).

Am I on the right track, or is there a better way to implement this?

I’m thinking something like:

```julia
eq_rn = @reaction_network begin

    @species H_cat(t), OH_an(t), H2O(t)
    # other reactions in the network with some time dependent kinetics
    7.9, HA --> H_cat + A_an
    1.2, B_an + H_cat --> HB
    (3.1, 5.7), 2H_cat + C_dian <--> H2C
    
    @equations begin
        # strictly enforce the pKa of water - now OH_an is an observable dependent variable
        14 ~ (H_cat(t) * OH_an(t)) / H2O(t)
    end
end

```

Perhaps this should be modeled instead with a normal bidirectional reaction (\<–\>) using a diffusion limit estimation for the forward rate constant (k1) and a multiple of the pKa for the reverse rate constant (k2 = k1 / pKa)?

---

<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:** [June 25, 2025, 1:22pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/2 "2025-06-25T13:22:00Z")

</div>

Your model here looks underspecified in that the number of equations you will have is 4 (3 ODEs and one algebraic equation), but you have way more species than that.

Catalyst should work fine for adding algebraic constraints as equations – you can then use structural\_simplify in MTK to handle reduction and getting a representation that is good to use with DAE solvers – but you do need to make sure that you are specifying a consistent, well-posed set of equations collectively (i.e. each species/variable in the final model can be associated with one independent ODE or algebraic equation).

---

<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:** [June 25, 2025, 1:59pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/3 "2025-06-25T13:59:24Z")

</div>

If you want to enforce some type of additional constraint beyond a well-posed set of model equations, you could consider using DiffEqCallbacks that project to the desired manifolds: [Manifold Projection · DiffEqCallbacks.jl](https://docs.sciml.ai/DiffEqCallbacks/stable/projection/)

---

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [June 25, 2025, 7:25pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/4 "2025-06-25T19:25:22Z")

</div>

Thank you for that feedback, thats very helpful.

In connection with this - is there a preferred / natural way to add these algebraic constraints to an existing reaction network system?

I am thinking of behavior similar to the ReactionNetworkImporters [MatrixNetwork](https://docs.sciml.ai/ReactionNetworkImporters/stable/#Loading-a-matrix-representation) tools.

It would be ideal to be able to add algebraic constraints via some matrix representation or other similar approach. Then reaction networks could be constructed from external files (like a .csv) that are read in as matrices and passed to create a `ReactionSystem`.

If this behavior is not allowed, we would be restricted to specifying algebraic constraints ‘manually’ in the Julia code rather than being able to generate these constraints programmatically. I don’t know enough about the underlying mechanics of the @reactionsystem macro to have a good feel for what would be allowed or the natural way to do this…

---

<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:** [June 27, 2025, 2:38pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/5 "2025-06-27T14:38:12Z")

</div>

You could always use the symbolic interface and generate any symbolic algebraic constraints you want programmatically after reading your CSV.

Generally using the DSL is not a good workflow for generating models from other input sources (it is basically a wrapper over the symbolic interface itself). ReactionNetworkImporters, for example, just uses the symbolic interface under the hood to map matrices to ReactionSystems.

---

<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:** [June 27, 2025, 3:11pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/6 "2025-06-27T15:11:18Z")

</div>

This blog post details all of the different ways that a DAE can be transformed into something solvable:

> **[Machine learning with hard constraints: Neural Differential-Algebraic...](https://www.stochasticlifestyle.com/machine-learning-with-hard-constraints-neural-differential-algebraic-equations-daes-as-a-general-formalism/)**
>
> We recently released a new manuscript Semi-Explicit Neural DAEs: Learning Long-Horizon Dynamical Systems with Algebraic Constraints where we showed a way to develop neural networks where any arbitrary constraint function can be directly imposed...

ModelingToolkit, and thus by extension Catalyst, doesn’t handle them all yet, but it’s getting there. Probably in two months we’ll have the stuff to handle overdetermined equation sets and generate manifold callbacks etc. It’s pretty active work.

Until then, using the manifold projection callback yourself is probably the quickest way that tends to work pretty well.

> [@spackman](#):
>
> In connection with this - is there a preferred / natural way to add these algebraic constraints to an existing reaction network system?

Note that the difficult part about this is not adding the algebraic constraints but finding the equations to delete. ModelingToolkit’s DAE algorithms currently cannot handle an overdetermined set, so if you add an algebraic condition you need to remove one of the differential equations. It should make sense though: if `x^2 + y^2 ~ 1` then have an `x'` and a `y'` equation is overdetermined, since solving for `x` determines `y` via the constraint and so you shouldn’t have both. But the current algorithm will just complain there’s too many equations and you need to go back and either remove the `D(x) ~ ...` equation or `D(y) ~ ...` equation. Removing the wrong one could introduce a singularity, etc. etc. so it is generally hard but we should make a heuristic for it, but we just don’t have this part yet automated.

---

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [July 21, 2025, 4:06pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/7 "2025-07-21T16:06:57Z")

</div>

Ah, I see…  
I will be very interested to see this work when it is complete, it sounds like it coule be exactly what I need!

---

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [July 21, 2025, 4:16pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/8 "2025-07-21T16:16:27Z")

</div>

> [@isaacsas](#):
>
> Generally using the DSL is not a good workflow for generating models from other input sources (it is basically a wrapper over the symbolic interface itself).

Thanks for this insight- can you provide an clarity on what is meant by the symbolic interface and the acronym ‘DSL’ (perhaps ‘domain specific language’)?

I am concered that an approach using the symbolic interface would require ‘manually’ writing out all of the differential equations in julia code, rather than being able to abstract them away to a .csv file like with ReactionNetworkImporters.

My ideal solution would allow for the complete differential algebraic equations (both network and algebraic constraints) to be defined by an external file, and any ‘manual’ edits to the julia code for a specific system would only be in specifying the system timesteps / solution times, or additional event criteria. In this way, any code I write could be easily extended to different systems by simply swapping out the external file to update the network topology or kinetic parameters.

---

<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:** [July 24, 2025, 11:28am UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-equilibrium-i-e-pka/130149/9 "2025-07-24T11:28:16Z")

</div>

DSL = domain specific language, which in Catalyst is using the `@reaction_network` macro and its custom notations to specify a model. Symbolic interface means manually building each `Reaction`, variable, species, and `Equation`, and then manually passing them into a `ReactionSystem`. The latter is usually an easier approach when trying to build a system programmatically via code (vs. specifying a full model directly in a file).