support more general affects

[isaacsas]

https://github.com/augustinas1/MomentClosure.jl points out we currently can't support a geometric distribution of produced reaction products. It would be good to get a nice interface for adding this to Reactions. We could potentially do one of the following:

1. Allow substrate stoichiometries to be more general functions instead of whole numbers.
2. Add an affect kwarg that is nothing by default, but can be a user-passed function.

One issue is how this should go over into the ODE models. This make complete sense for jumps, but may need a little thought on what the right representation for ODEs is. (i.e. if using a geometric RV do we want that to propagate to a random coefficient in the ODEs, or should the ODEs be some analogous mean-field model)?

@TorkelE @ChrisRackauckas any thoughts?

[TorkelE]

If I understand right, this would allow having each reaction, instead of adding a fixed number of reactants, the number which is generated would be a stochastic variable?

Generalising stoichiometries doesn't sound like a bad thing, previously non-integer stoichiometric have been requested (which in reverse, might be easier to make for ODEs, but harder for jumps)...

Syntax wise reactants are currently Symbols, while stoichiometries are numbers. Enabling the macro to parse functions as well wouldn't be a problem. For the continuous case, would it make sense that the ODE was using the mean and somehow allow for the creation of SDEs using the distribution?

An intermediary solution would be to disable ODE/SDE generation when this kind of stuff is used (maybe even generating a jump system directly). I imagine that quite some of the base structure would be required to allow for this anyway.

(Allowing for more kwargs might also be long-term desirable)

[ChrisRackauckas]

I would think the right way to do it is to allow distributions somewhere in the macro and then have a lowering process that checks the distribution for how it lowers. Handling non-Poisson as only non-MassActionJumps would be the way to get it going until someone really needs it.

[isaacsas]

First a caveat; my comments were at the level of ReactionSystem and the Reactions it stores -- I'm more thinking about how we should encode this at the MT system level, and what our transformations to other system types should look like (or at least for ODEs/SDEs, for jumps it seems clear to me). I'm not worrying about how to handle this in the DSL yet (but it doesn't seem like it should be a big deal to allow).

@TorkelE The idea is that the number of products is not fixed. Reactants aren't really an issue since we already support custom rate laws. I thought non-integer coefficients worked at the level of ReactionSystems, but perhaps I'm misremembering.

If we're going to allow random variables for coefficients though, it seems like we should go all in and allow more general functions of state that determine the change in state after a reaction occurs. One option would just be to allow functions for each individual stoichiometric coefficient as part of a specified stoichiometry vector, but perhaps it is better for users to specify a single affect-type function that just updates the state following a reaction.

[isaacsas]

I'd like to keep ReactionSystems as generic as possible, and leave any lowering / processing to occur within them (and not at the DSL level). This keeps a symbolic and programmatic option available users that fully supports the DSL functionality.

[TorkelE]

Allowing as general things in ReactionSystems seems sensible. As long as it doesn't affect the performance of the solvers it should be fine. The case of ODE/SDEs shouldn't really be a problem (once we have decided how to translate stuff to equations, and until then we can just check if something weird has been done and give an error).

I'm less familiar with the advanced jump solvers, but again I guess one can have a check if things are "normal", and if so use what we already have, in which case it would not interfere?

[kaandocal]

While I'm not sure if this will satisfy everyone I would offer the following perspective on variable product numbers with ODE/SDE systems. If we view the ODE version of a system as a deterministic/mean-field limit we cannot directly model random variability in product numbers, and in this case the product stoichiometries should be deterministically set to their means as @TorkelE suggested. In the SDE version the randomness translates to an increase in the noise determined by the variance of the product numbers.

One way to interpret a reaction with rate r and random products is to split it up into different reactions for each possible product, each with rate r * P(product). This is equivalent for the jump formulation (you can see this by noting that the two are equivalent under the Gillespie algorithm), and from this interpretation we can recover the above behaviour for ODE and SDE systems.

One could brute-force random products by splitting up the reaction as above, but that will typically require truncation and may not be the most efficient way to do this.

[isaacsas]

Thanks for the comments. The issue I had when I took a look at this is that MT doesn’t currently work with a probability distribution from distributions.jl in it (it seemed like using such a distribution as the RV would be a nice approach). We need a way to specify random variables within MT expressions, and later detect and process them, and ideally this should work through ReactionSystems and not just the reaction_network macro. I don’t think functionality is currently in place to handle this, so it requires some work.

[isaacsas]

Maybe we can just register those distributions, I’ll have to check that.

[ChrisRackauckas]

We should just loop through Distributions.Sampleable and just register them all IMO.

[isaacsas]

We also should probably have a way to propagate RNG choice when the distribution is ultimately replaced with a call to rand.

[isaacsas]

I think we've handled the original issue here, but we could open a new issue if the current support is not general enough.