Zero Order Extinction Reaction

Specific Domains Chemistry
,
1

Is there any way to achieve the positive_domain effect with the Catalyst.jl syntax?

In my example (decomposition of CaC2 yielding Calcium gas and solid Carbon):
CaC2 → Ca(g) + 2C(s)
I am looking for a way to “remove” Ca from the reacting system. However, I get negative mass:

  • if I override the mass kinetic rate law (using => ) to set a zero-order extinction
    dif = @reaction Dif, Ca=>0
  • or set a negative reaction rate
    dif = @reaction -Dif, 0 --> Ca

plot_1

using a positive_domain callback, the code works as expected:

plot_2

here is the full example:

using DataFrames, Plots
using Catalyst, DifferentialEquations, IfElse
using StatsBase, Optimization, OptimizationOptimJL, OptimizationOptimisers
using OptimizationBBO, OptimizationNLopt
using LsqFit

function positive_domain()
    condition(u,t,integrator) = minimum(u) < 0
    affect!(integrator) = integrator.u .= integrator.uprev
    return DiscreteCallback(condition,affect!,save_positions = (false,false)) 
end

R = 8.314 #J/mol - Universal Gas Constant
MW = [64., 40., 12.]
solids = [1, 1, 1]
@variables t 
D = Differential(t)
@parameters A E Dif
hr = 20
T_iso = 1850. 
time_step = (T_iso-1000)/hr #minus the starting time
T(t) = IfElse.ifelse(t<time_step, hr*t+1273, T_iso+273) 
k(t) = (A*1E4) * exp(-(E*1E3)/R/(T(t)))

rxn = @reaction	k(t), CaC2 --> Ca + 2C
dif = @reaction -Dif, 0 --> Ca #Ca => 0.
@named rs = ReactionSystem([rxn, dif], t)
#osys = convert(ODESystem, rs)

p = (A => 5.9E2, E => 300, Dif => 2E-2)
tspan = (0., 120.)
u0 = [:CaC2 => 1.0, :Ca => 0.0, :C => 0.0]

oprob = ODEProblem(rs, u0, tspan, p)
osol = solve(oprob, AutoTsit5(Rosenbrock23()), saveat = 1., reltol=1E-12, abstol=1E-12)#, callback=positive_domain())

plot(osol, width=2, title="Reaction System", xlabel="Time (min)", ylabel="Species (molar)")
plot!(twinx(), osol.t, T.(osol.t) .-273, legend=:bottom, label="T [C]")
1 Like
2

Yes this is something that will be very numerically sensitive and using some kind of callback is generally a good idea. The isoutofdomain solver option and PositiveDomain callback in the DiffEqCallbacks.jl library were designed for this kind of application. The one that you have there introduces error so I would recommend the two standard approaches over that.

1 Like
3

Thanks @ChrisRackauckas ! I opted for the PositiveDomain option and it seems ok.

osol = solve(oprob, Tsit5(), saveat = 1., reltol=1E-12, abstol=1E-12, callback = PositiveDomain([1.0]; save = false, abstol = 1E-9))

I have not done any Benchmark but the Optimization block I am using to fit the parameters appears to work smoother.

4

@ChrisRackauckas any idea how I can avoid my problem formulation keep pushing to the negative side?

While the below figure looks ok with all concentrations respecting the zero-barrier. Small negative perturbations (see next figure) break progressively the mass balance by pushing some compounds (gold line) to higher than the initial values.

image

image

5

Turn on the conservation law pieces?

6

are you referring to the [remove_conserved = false], like in the bellow example?

H2A_m = (1000 / 176 * 1000) / (25 / 1000) #mg/MW  ->umol/l
fct_rn = 1E-3
newparam=[0.7026, 1.35112]
rn = @reaction_network rn begin #actual reaction system
  @parameters  k_1 k_2 k_4 
  k_1*$fct_rn, cat_a --> cat_b 
  k_2 * $fct_rn * cat_b, H2A --> DHA + H2 
  (k_4*$fct_rn), DHA  --> H2A
end
vol = @reaction_network vol begin #bulk system to account for dispersion
  $newparam[2], H2 --> 0
end
full = @reaction_network begin
  $newparam[1], $(rn.H2) --> $(vol.H2)
end
@named full = compose(full,[rn, vol]) #combine the two systems
osys = convert(ODESystem, full; remove_conserved = false)
p = [rn.k_1 => 11.23, rn.k_2 =>0.00346, rn.k_4 =>25.75]
time_r = (0.0, 500.0)
u0 = [rn.H2A => H2A_m, rn.DHA => 0.0, rn.H2 => 0., vol.H2 => 0.0, rn.cat_a => 10.0, rn.cat_b => 0.0]
op = ODEProblem(osys, u0, time_r, p)
sol = solve(op, Tsit5(), callback = PositiveDomain([1.0]; save = false, abstol = 1E-12))
7

@kkakosim Can you provide a MWE that produces the initial plots as well as the final plots where you get leackage?

8

I tried to keep only what demo the problem, hopefully it didn’t mess up with the physical system too much.

the problem arises when (0.0004 - k_4) gets negative values. While positive_domain() pushes all the species to positive values, it is still possible to get the negative values as in the figure below.

using ModelingToolkit, Catalyst, DifferentialEquations
using Plots
plotly()

H2A_m = 1.
rn = @reaction_network rn begin #actual reaction system
  k_2, H2A --> DHA + H2 
  (0.0004-k_4), H2A  --> DHA
end
#osys = convert(ODESystem, full; remove_conserved = true)
p = [:k_2 =>0., :k_4 =>0.00041]
time_r = (0.0, 500.0)
u0 = [:H2A => H2A_m, :DHA => 0.0, :H2 => 0.]
op = ODEProblem(rn, u0, time_r, p)
sol = solve(op, Tsit5(), callback = PositiveDomain([1.0]; save = false, abstol = 1E-12))  
plot(sol)

image
image734×488 9.77 KB

the negative values for k_4>0.00041 come into play when running optimization for fitting the parameters of the system to experimental data. I could still limit the k_4 with the [lb, ub] of the optimization loop but I wonder why the positive domain doesn’t keep the system in line.

9

Is (0.0004-k_4) supposed to ever go negative?

I looked at it, but your minimal example basically just has a negative term for the equation of DHA. I think the positive domain function is supposed to prevent steps that takes you into negative (due to numeric error or e.g. an SDE). In this case, going negative is literally the only possible steps. I.e. you have a model which says that DHA should go negative, but want to have a positive domain that prevents it. Are you sure that this is what you want?

You mention optimisation? If what happens is that you want to make an optimisation run that has a parameter k:

k, H2A  --> DHA

where the optimisation takes k to negative, then I think the optimisation is the real problem. If you are optimising these network, you should probably do it on a log scale anyway, which would naturally prevent k from going negative.

10

@Gaussia, you are right that in this particular case the issue is the negative reaction constant and there is no other way. I presumed that the “positive_domain()” would still be able, or should be able, to prevent this and thus eliminate trivial optimization solutions. Actually, without the “positive_domain()” the system would have drifted to negative values way earlier (time wise).

Can you explain the “logarithmic” implementation you suggested?

11

Basically, instead of fitting parameters for

@reaction_network begin
    p, 0 --> X
    d, X --> 0

you do

@reaction_network begin
    10^P, 0 --> X
    10^D, X --> 0

Now, -inf < P, D < + inf (so you don’t have to worry about bounds. Also, these parameters typically exist more on a log scale, typically improving performance a lot. Packages like PEtab.jl (recommended) fixes this automatically for you.

I think positive_domain is mainly for cases like SDEs where you might jump out of the domain, but @ChrisRackauckas would need to clarify this.

1 Like
12

No it’s for cases with ODEs where you’re solving at a higher tolerance but want to enforce the solution stays positive, given that the true solution is actually positive. That last part is very crucial though, and IIRC it even has a big warning in the docs that it only works if the true solution has to be positive, and thus it doesn’t apply in this case.

1 Like
13

Can you be more specific on this? the reaction rate may get negative values but the physical problem requires conservation and a positive true solution.

14

If the derivative is zero when a value is zero then it will not have a positive solution. The tooling to retain positivity will fail in this case.

1 Like
15

I never thought of these problems in such a “computational” way. So is there any practice to enforce a positive solution while allowing for zero values? both possible for a reaction-system.

16

So, a well-defined reaction system should never be able to go negative. If you have a model where this is the case, it is most likely ill-defined. I.e. in your minimal example you have a negative rate, which simply is not a thing for these models (unless you are using them in a domain where this is the case?).

Exactly why do your model permit negative values? My impression is that it is because of your optimisation procedure scanning these values? In this case the problem seems to be that your optimiser checks incorrect models, and the solution is to prevents this (i.e. putting bounds on the optimisation, or use log-scale). This should also speed up optimisation time since you avoid checking models which you know at the start are not true.

If the cause is something else though, please elaborate and we can see if there is another solution.

17

You are right:

in this case, as well as many similar, often, I am scanning in an exploratory mode and may have reactions which are obsolete - ideally the parameter fitting should grant them trivial e.g., -->zero rate. Also, a negative rate is not necessarily a “infeasible” because in this exploratory mode may imply the “reverse” reaction. Another challenge is that most of the times the objective function compares with one of the compounds or another indirect quantity (like total mass) which might smoothly lay in the positive domain and without noise.

I will try to create a MWE on the same, for the sake of curiosity, this is what happened to me:

  • in this case. I was fitting the parameters on H2, with good agreement, until I realized that it was converging to infeasible system/solution (i.e., negative concentrations for other compounds).
  • in another case, I was fitting three reactions over a very wide range of temperatures because none of the activation energies is known, e.g.,
    k1(T), A + B → C (production of C)
    k2(T), C + A → D (side reaction of C)
    k3(T), C → E (degradation of C)
    • all or some of the parameters of the three Arrhenius-type k(T) equations (2 per equation) are fitted against the remaining mass of the mixture (i.e, mass of A+C)
    • with the specific optimization bounds being very wide, especially the 3rd reaction (degradation) pushes the compound C to negative concentrations.
    • despite that the objective function yielded an excellent agreement because it is indeed unaware of the rest of the issues.

Sidenote, I was not aware of PEtab.jl - looks interesting and mature enough. I am trying to move there.

Finally, understanding that all this is in volunteering mode, if you are interested in seeing more MWE please let me know and I will work on one providing the experimental data as well, because without them a MWE may be misleading :smiley:

18

An extensive introduction in using PEtab + Catalyst can be found here: https://docs.sciml.ai/Catalyst/dev/inverse_problems/petab_ode_param_fitting/
(in addition to the package’s own documentation). I highly recommend it for these kinds of optimisation.

Not an expert, but as long as the constant A is positive in the Arrhenius equation, isn’t it always positive?

Finally, for a reaction

k, X --> Y

the reaction

-k, X --> Y

does not represent the reverse reaction. I.e. the first one yields the ODE:

dX/dt = -k*X
dY/dt = +k*X

while the second one is:

dX/dt = k*X
dY/dt = -k*X

(which totally can go negative, and to infinity, simultaneously. Two things that should not happen)
The reverse reactition

k, Y --> X

instead yield the ODE

dX/dt = +k*Y
dY/dt = -k*Y

A zero value for k should still be possible and yield no problem, but you can include this in your bounds (and probably you’d check if k < small number, just round it to 0. I still really do not see why you’d want to allow reactions with negative rates.

19

Yes Correct.

that’s why I had reverse in quotes, it doesn’t represent the reverse reaction as we per the definition but it does have the reverse effect of generating X instead of Y, though not with a first-order dependency to Y but to X, like an autocatalytic reaction.

Based on all the discussion, I updated the system in my previous posts:

  k_1*$fct_rn, cat_a --> cat_b 
  k_2 * $fct_rn * cat_b, H2A --> DHA + H2 
  (k_4*$fct_rn), DHA  --> H2A

to this scheme:

$k_1 * $fct_rn * cat_a, cat_a --> cat_b 
(10^k_2) * cat_b, H2A --> DHA + H2 
($k_deg * 10^k_4), H2A --> DHA

which as you mentioned, cannot go negative and does converge fast enough. In addition, when the last step is trivial/negligible it will yield a very high negative k_4

20

Iäd still protest that -k, X --> Y represent the “reverse” reaction (using quotation), as this still suggests that it represents a reaction in some forms (from what I can tell, it really does not). I.e. you cannot actually write down a reaction that that has the same effect on the system as -k, X --> Y does (hint, the ODE suggests that Y is being consumed, however, that means that it has to exist on the reaction lhs, which means it would affect the rate, which it do not).

Your updated scheme look good. I don’t know what all your parameters mean, but you probably want to log scale all you want to estimate (unless you already have done it. If you plan to spend lots of time doing this, I’d definitely have a look at using PEtab instead though.