I use Catalyst.jl to model kinetic data of a chemical reaction network.
Somehow, for some specific sets of initial concentrations, the resulting time-concentration course results in negative concentrations. Can anyone explain why this happens or where the error in my code is? At t = 9270 s, the reactant sub2 is added.
Attached is the plot where the unexpected behavior occurs.
When every parameter of the simulation is kept the same, except for the initial concentration of :prod, which changed from 0.00059 (behavior as expected) to 0.00058 (weird behavior). Other parameters can be changed too to “switch” between the two types of results.
Thank you!
using Catalyst, DifferentialEquations, Optim
import DifferentialEquations as DE
using Plots
rn = @reaction_network begin
(k1, k_1), sub1 + cat <--> catI
(k2, k_2), catI + sub2 <--> cat + prod
end
# initial conditions
conds = Dict(:cat => 0.00096, :prod => 0.00058, :sub2 => 0.0, :sub1 => 0.015, :catI => 0.0)
timespan = (0.0, 15000.0)
# parameters for reaction network that fit the experimental data (not given in this example) best
pmap = Dict(:k_1 => 8.824955678175811e-9, :k_2 => 8396.690567743815, :k1 => 0.07834400385602105, :k2 => 73536.72395735912)
# time at which :sub2 is added in the experiment
addtime = 9270
# with a concentration of
sub2conc = 0.00042288949
# of species :sub2 which is index of the rn:
allspecies = [Symbol(replace(string(ispec), "(t)" => "")) for ispec in species(rn)]
idxsub2 = findfirst(==(:sub2), allspecies)
# simulate reaction network
affect!(integrator) = integrator.u[idxsub2] += sub2conc
cb = DE.PresetTimeCallback(addtime, affect!)
oprob = ODEProblem(rn, conds, timespan, pmap)
sol = DE.solve(oprob, Rodas5P(), saveat=10, callback = cb)
plot(sol, lw = 4)
This may be a stiffness issue with your system. Solving your example with KenCarp47 instead gives non-negative solutions. @ChrisRackauckas might be better able to suggest what solver you should use.
As an aside, I’d suggest encoding your event as a symbolic event. ModelingToolkit no longer guarantees that it will map the ordering in species to the order in u as far as I am aware, so you might not actually be updating the variable you think you are updating (and even if it works now it could break in the future). See Coupled ODEs, Algebraic Equations, and Events · Catalyst.jl
@isaacsas Using KenCarp47 solves the problem, the simulation is now as expected.
Thanks for the hint on how to code the event properly, I’ll change that in my code.
@ChrisRackauckas Did I get it correct that the interpolation is necessary because I entered the keyword “saveat”? I did this to calculate a loss function between two simulations of different reaction networks. Is there any way to circumvent that or should I stick to an ODE solver like Rodas5Pr or KenCarp47 instead?
If I don’t use saveat, how can I calculate the loss function between two different simulations? I did this using the same grid for both and I don’t find a solution without doing that.
Chris is just suggesting turning saveat off and dense = false as a test to debug what is going on with Rodas5p, not that it is what you should ultimately do. Then you can look at the values in sol.u to see if there are any that are negative.
Ah sorry, I got that wrong.
Using Rodas5Pr works in all cases (using saveat or using dense=false).
With Rodas5P using saveat or dense=false, both show the same plot with imaginary frequencies and have negative values.
Just to prevent confusion, I used either saveat or dense=false.
