# Microkinetic simulations yield negative concentrations using Catalyst.jl
**URL:** https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874
**Category:** Modelling & Simulations
**Tags:** question
**Created:** [November 14, 2025, 12:04pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874 "2025-11-14T12:04:07Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@lisa-schmoltzi](https://discourse.julialang.org/u/lisa-schmoltzi)
#### Post date: [November 14, 2025, 12:04pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/1 "2025-11-14T12:04:07Z")
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!
```julia-auto
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)
```
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [November 14, 2025, 1:13pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/2 "2025-11-14T13:13:28Z")
> [@lisa-schmoltzi](#):
>
> ```julia-auto
> 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](https://docs.sciml.ai/Catalyst/stable/model_creation/constraint_equations/#constraint_equations_events)
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [November 14, 2025, 1:59pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/3 "2025-11-14T13:59:29Z")
You might also try using smaller error tolerances and see if that helps.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 14, 2025, 2:00pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/4 "2025-11-14T14:00:09Z")
You can try Rodas5Pr for the residual error estimator to be more strict.
The negativity is likely not in the solution but in the interpolation of the solution.
---
### Author:  [@lisa-schmoltzi](https://discourse.julialang.org/u/lisa-schmoltzi)
#### Post date: [November 14, 2025, 2:32pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/5 "2025-11-14T14:32:14Z")
Thank you for your quick replies!
@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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 14, 2025, 2:50pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/6 "2025-11-14T14:50:43Z")
Don’t use saveat and set dense=false and see if it’s all positive
---
### Author:  [@lisa-schmoltzi](https://discourse.julialang.org/u/lisa-schmoltzi)
#### Post date: [November 14, 2025, 3:08pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/7 "2025-11-14T15:08:55Z")
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.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [November 14, 2025, 3:16pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/8 "2025-11-14T15:16:41Z")
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.
---
### Author:  [@lisa-schmoltzi](https://discourse.julialang.org/u/lisa-schmoltzi)
#### Post date: [November 14, 2025, 3:35pm UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/9 "2025-11-14T15:35:20Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 17, 2025, 4:28am UTC](https://discourse.julialang.org/t/microkinetic-simulations-yield-negative-concentrations-using-catalyst-jl/133874/10 "2025-11-17T04:28:52Z")
Okay yeah it’s just the residual accuracy.