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)?
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).
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
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…
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.
This blog post details all of the different ways that a DAE can be transformed into something solvable:
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.
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.