Miscomputation when using a stoichiometric coefficient != 1

[odurif1]

version=Catalyst v13.5.1

The result is incorrect when computing a simple chemical reaction 2A -> B

It seems to be an issue when using a stoichiometric coefficient != 1

Example

rn = @reaction_network begin
    k1, 2A --> B
end

k = [
    :k1 => 2.1e-13
]

u0 = [
    :A => 2.8e10,
    :B => 0
]

tspan = (0.0,10.0)
op    = ODEProblem(rn, u0, tspan, k)
sol   = solve(op)

Bellow a plot comparing this numerical result with the analytical solution: $\ f(t)=\dfrac{1}{(2 k_1 t+1/[A]_0)}$

[isaacsas]

Have you read https://docs.sciml.ai/Catalyst/stable/introduction_to_catalyst/introduction_to_catalyst/#Reaction-rate-laws-used-in-simulations

Catalyst assumes microscopic rates by default, so your rate might differ by a factor of two from what you are expecting. You can turn this off as described in that section.

[odurif1]

Thank you, it is exactly a factor two, and my issue is solved using _combinatoricratelaws=false

ODEProblem(rn, u0, tspan, k; combinatoric_ratelaws=false)

I did not knew this subtility between microscopic rates and transition rate . Out of curiosity, why is _combinatoricratelaws=false is set to true by default ?

[isaacsas]

We decided to use the physical interpretation that is common for stochastic jump process models, our most microscopic supported model, by default. The rate has a clear physical meaning at that level, i.e. the probability per time one random set of substrates can react. We then used this interpretation for all mathematical model types to ensure they are comparable.

You can also disable this as a keyword argument to ReactionSystem, unfortunately we haven’t added support for disabling it in the DSL yet. @TorkelE we should probably add a DSL option to set this keyword argument to ensure the DSL fully supports all ReactionSystem features.

[TorkelE]

Yes, adding a DSL option for this is quick work, I will do it once the other changes are done.

@odurif1 Think of it like this. Imagine that you have three hypothetical molecules x, y1, and y2. You then have two reactions: k, x + x --> X2 k, y1 + y2 --> Y2 Here, the binding sites where x binds x and y1 binds y2 are identical, so these reactions have identical rates.

Now imagine that you put N x molecules in one container. In another container, you put N y1 molecules and N y2 molecules. In this case, Y2 should be produced at a higher rate than X2 (since there is a lot more of its substrates). This is only the case if you actually adjust the rates for the combinatorics of the reactions, which is why we do it.

Practically, this adjustment is often omitted. However, some would do it, so in either case some people will miss it. In the end we decided to use the approach which is technically more correct (while also trying to raise awareness that this adjustment is a thing, and technically should be carried out).

[odurif1]

Thanks a lot for your help. I am making a mental picture but personally, I found this very confusing, certainly because I am not used to this approach.

In chemical kinetics the system $2 A \rightarrow A_2$ is written as:

$\frac{d[A]}{dt} = - 2 k [A]^2$ following the reactant or $\frac{d[A_2]}{dt} = k [A]^2$ following the product

It is written this way because the physical interpretation is here straightforward, 2 $A$ are consumed when one $A _2$ is produced. However, the rate coefficient $k$ is the same whether you consider the reactant or the product because it is the same reaction.

In the microscopic approach employed here, the rate coefficient has to be twice higher when considering the reactant compared to the product. Thus we would have to consider not one rate coefficient for this reaction but two differents, which makes things, to my point of view, confusing.

But let me thank you for this wonderful and easy-to-use Julia package!

[isaacsas]

If you think about the rate law you are using, it can't make sense when populations are sufficiently low (for example, when there is only 1 molecule of species A). This is because it is an approximation, typically thought of as arrising when populations are sufficiently large or in related thermodynamic limits (though there are other ways to derive it).

If one considers the more fundamental representation of keeping track of the exact number of A's, call this $A(t)$, one then needs a rate law that is zero when there is only one A molecule. At the microscopic level, one assumes there is an inherent probability per time, $k$, that a pair of species A molecules randomly distributed within the domain may react. Assuming the A molecules are independent and well-mixed, the rate law, or propensity as it is called in the stochastic modeling literature, is then

k \frac{A (A-1)}{2}

i.e. the rate for a pair to react multiplied by the number of distinct pairs. This represents the total probability per time the reaction will occur within the system.

If one converts to concentration units, i.e. units of number per volume for simplicity, this becomes $\tfrac{k V}{2} [A] ([A] - \tfrac{1}{V})$. If we imagine a thermodynamic limit where $[A]$ is held fixed but the number of A molecules and the volume are increased, this becomes well-approximated by $\hat{k} [A]^2$ where $\hat{k} = k V/2$. This would correspond to the rate you are using (up to whether you want units of Molars or number concentrations).