Calling complete() adds differential constraints to algebraic extension of ReactionSystem

Specific Domains Modelling & Simulations
, , ,
1

Hey guys,

I’m doing some experimentation with Catalyst.jl and algebraic equations to define systems where the system volume is not constant. This can occur when the major solvent is a reactant in appreciable quantities, and it impacts all the rate laws through redefining the species concentration.

I’m having trouble converting my system to an ODEProblem() because calling the complete() function appears to add additional equatiosn to my system. Here is an example:

using Catalyst, ModelingToolkit

@independent_variables t
@parameters k1, k2, k3, rho
@species H2O(t) A(t) B(t) mH2O(t) mA(t) mB(t) V(t)


rxns = [
        Reaction(k1, [H2O], [A , B])
        Reaction(k2, [A], [B])
        Reaction(k3, [B], [H2O])
]

@named rn = ReactionSystem(rxns, combinatoric_ratelaws=false)

# set up a system of equations linking the mass of each species to its concentration
equations = [
        V ~ mH2O / rho
        H2O ~ mH2O / V
        A ~ mA / V
        B ~ mB / V
]

@named alg = ODESystem(equations, t, [mH2O mA mB], [])

rn2 = extend(alg, rn)

rn2 = complete(rn2; sys=rn2)

When I display the system ODEs, I get this:

\begin{align} \frac{\mathrm{d} \mathtt{H2O}\left( t \right)}{\mathrm{d}t} &= - \mathtt{k1} \mathtt{H2O}\left( t \right) + \mathtt{k3} B\left( t \right) \\ \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= \mathtt{k1} \mathtt{H2O}\left( t \right) - \mathtt{k2} A\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= \mathtt{k1} \mathtt{H2O}\left( t \right) + \mathtt{k2} A\left( t \right) - \mathtt{k3} B\left( t \right) \\ \frac{\mathrm{d} \mathtt{mH2O}\left( t \right)}{\mathrm{d}t} &= 0 \\ \frac{\mathrm{d} \mathtt{mA}\left( t \right)}{\mathrm{d}t} &= 0 \\ \frac{\mathrm{d} \mathtt{mB}\left( t \right)}{\mathrm{d}t} &= 0 \\ V\left( t \right) &= \frac{\mathtt{mH2O}\left( t \right)}{\mathtt{rho}} \\ \mathtt{H2O}\left( t \right) &= \frac{\mathtt{mH2O}\left( t \right)}{V\left( t \right)} \\ A\left( t \right) &= \frac{\mathtt{mA}\left( t \right)}{V\left( t \right)} \\ B\left( t \right) &= \frac{\mathtt{mB}\left( t \right)}{V\left( t \right)} \end{align}

There are now 3 additional differential equations that were not explicitly specified in the original system. I think these equations are causeing my system to be over-defined when I try to convert it to an ODEProblem:

u0 = [mA => 0.1, mB  => 0, V => 50] #mH2O is fully specified by the volume V and the density rho
ps = [k1 => 1, k2 => 3, k3 => 0.5, rho => 0.99]
oprob = ODEProblem(rn2, u0, (0.0, 5.0), ps, structural_simplify=true)

This yields:

ExtraEquationsSystemException: The system is unbalanced. There are 7 highest order derivative variables and 10 equations.
More equations than variables, here are the potential extra equation(s):
 0 ~ -H2O(t) + mH2O(t) / V(t)
 0 ~ -V(t) + mH2O(t) / rho
 0 ~ mA(t) / V(t) - A(t)
Note that the process of determining extra equations is a best-effort heuristic. The true extra equations are dependent on the model and may not be in this list.

Ideally, I want to be able to proceed with a system without these additional equations. This could happen by either

  • pruning the system after complete to remove unwanted differential equations (I don’t think this is possible because a complete system is not mutable)
  • stopping the generation of unnecessary equations from the start
  • explicitly defining these differential equations (possible, but not ideal)

The easiest option would be to prevent their generation in the first place.

Any tips on this problem?

2

I’ve solved the completeness problem -
The additional equation appear to have been added because the new variables (mH2O, mA, mB) were added using the @species macro.

Re-defining using the @variables macro eliminates the issue:

using Catalyst, ModelingToolkit

@independent_variables t
@parameters k1 k2 k3 rho
@variables mH2O(t) mA(t) mB(t) V(t)
@species H2O(t) A(t) B(t) 


rxns = [
        Reaction(k1, [H2O], [A , B])
        Reaction(k2, [A], [B])
        Reaction(k3, [B], [H2O])
]

@named rn = ReactionSystem(rxns, combinatoric_ratelaws=false)

equations = [
        V ~ mH2O / rho
        H2O ~ mH2O / V
        A ~ mA / V
        B ~ mB / V
]

@named alg = ODESystem(equations, t, [mH2O mA mB V], [rho])

rn2 = extend(alg, rn)

rn2 = complete(rn2; sys=rn2)

Now the system can be defined as fully determined:

u0 = [mA => 0.1, mB  => 0, mH2O => 49.5]
ps = [k1 => 100, k2 => 30, k3 => 0.5, rho => 0.99]
oprob = ODEProblem(rn2, u0, (0.0, 5.0), ps, structural_simplify=true)

It happens that the system defined in this example is singular, which appears during the solution step.

3

You could alternatively pass include_zero_odes = false to convert or pass the isbcspecies = true metadata when creating the species. The later is a SBML feature that tags that no equation should be generated for the species from its inclusion as a substrate or product, I.e that the species will have a separate equation given for its dynamics elsewhere.