Combinatoric rate law scaling

[paulflang]

In the docstring for Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem) I read:

- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate
law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
`combinatoric_ratelaws=false` then the ratelaw is `k*S^2`, i.e. the scaling factor is
ignored.

I am thinking that most users who set combinatoric_ratelaws=true, would want k*S*(S-1)/2!. Users who set combinatoric_ratelaws=false probably want to stick to k*S^2 as is.

[isaacsas]

I don't think there is any single convention that is always followed, but many ODE models replace S-1 with S in my experience (especially if they are working in units of concentration). We tried to go with this choice for consistency with this definition. The scaling is even more confusing as some authors would use 2*k*S^2 and some would use k*S^2 for the total rate one loses species S (i.e. the stoichiometric coefficient times the rate law). We could definitely add an option to use S*(S-1) though if that is something that is needed for your applications.

[paulflang]

Hi @isaacsas ! Thanks for getting back. To my knowledge there are clear conventions. There are microscopic/stochastic and macroscopic/deterministic rate laws. So for mass action reactions

• Microscopic rate laws are a constant times the possible combinations of particles colliding. This is (n choose k).
• Macroscopic rate laws are a constant time the product of reactant concentrations.

Since the division by 2 stems from n choose 2 = n!/((n-k)!*k!) there is no point in dividing by 2! if using n^2 in the numerator. Therefore, I believe any user who sets combinatoric_ratelaws=true wants to use k * n choose k and any user who sets combinatoric_ratelaws=false wants to use k * n^2.

I just checked some literature to confirm that and found Eqn 20, 21 here, and Table 1 here.

To avoid confusions with stoichiometries, reaction velocities in deterministic and stochastic models described in units of reactions per time. To determine how the abundance of the reactants and products change, this number is multiplied by their stoichiometry.

[isaacsas]

Isn't losing the "-1" consistent with the rigorous mean-field limit of the stochastic model? If we assume the stochastic rate law is k*A*(A-1)/2 which we agree on, the large-population mean-field limit assumes k = khat / V where V is a system size parameter (usually volume times Avogadro's number). Convert to units of Molar concentration, C_A(t) = A(t)/V. In that case the mean-field (large-population) limit gives d C_A / dt = -khat (C_A)^2 = -2 * (khat * (C_A)^2/2) (the -1 becomes a -1/V which approaches zero as V -> infinity). But please do let me know if that is not the rigorous limit; perhaps I am misremembering the result in the homodimerization case or perhaps it is not covered by Kurtz's work. Formally, reintroducing V to convert back to units of "number of" we get the large-population mean-field approximation for the number of A as dA/dt = - 2 * (khat/V * A^2/2) = -2 * (k*A^2/2), i.e. twice the current ODE rate law. So our current default gives consistency between the stochastic model and the rigorous large-population mean-field ODE model.

With all that said, I can certainly see that one might prefer to think of the ODE model as a moment approximation that replaces averages of products with products of averages in the stochastic model, and hence keep the "-1" as you suggest. Similarly, we might want to offer an option to use the jump rate laws within the CLE too (though I don't know if this might make it more prone to going negative than it already is).

It might also make sense for us to consider changing the ODE defaults to combinatoric_ratelaws=false. At least one other user has brought this up in the past. (Though we've had the current convention for a long time, so this would be a big and potentially subtle breaking change for anyone that has been using Catalyst to model ODEs with the current scaling.)

Maybe the best solution would be to offer all three and require users to knowingly choose which they want (or perhaps we allow passing a function to set mass-action rate laws during convert, and offer all three as premade options macroscopic(), meanfield() and microscopic()). This would at least break user code in the sense that it would just not work till they actively choose a rate law function to use.

@TorkelE @ChrisRackauckas any thoughts?

[isaacsas]

Also, this has been a long weekend and I'm a bit out of it, so if any of that is confusing feel-free to let me know!

[isaacsas]

It also looks like we use combinatoric_ratelaws in an analogous manner for the jump rate laws, allowing users to indicate they don't want the constant factorial terms (perhaps because they have already included them in their rates).

[paulflang]

@isaacsas yes, you are completly right. sorry for my confusion and thanks for the discussion. it does make sense to loose the -1 here. And I also think it does make sense to have combinatoric_ratelaws default to true considering that ReactionSystems use microscopic parameters. Regarding the parameter V needed to convert from microscopic to to macroscopic systems, I commented on March 8 on https://github.com/SciML/CellMLToolkit.jl/issues/15

[isaacsas]

Let me move this over to Catalyst, and keep it open there. I think generalizing our scalings to handle the three options I mentioned above is worth thinking about longer-term at least.