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

Specific Domains Modelling & Simulations
, ,
1

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

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

2

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

3

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

4

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