Support for (rate) parameters that are distributions

[TorkelE]

I have been looking a bit at https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/parametric_stoichiometry/#Gene-expression-with-randomly-produced-amounts-of-protein and I think it would be good to

• [ ] Enable an easier syntax to generate the model.
• [ ] Enable an easier syntax to simulate the model.
• [ ] Enable any parameters to be set to a distribution.

Easier syntax to construct models with randomly distributed parameters.

It would be neat if one could do something like:

using Catalyst, Distributions
p = Uniform(0.8,1.2)
@reaction_network begin
    $p, 0 --> X
    d, X --> 0 
end

or even

using Catalyst, Distributions
p = Uniform
@reaction_network begin
    $p(0.8,1.2), 0 --> X
    d, X --> 0 
end

(possibly also having parameters in the distribution, e.g. p(p_min,p_max))

Maybe we cannot go all the way here, but further along would be cool. I'm not sure if this is something that would make sense, but if possible, we could try to register some common distributions.

Enable an easier syntax to simulate the model

I tried to make the example work with the simplified syntax:

using Catalyst, Distributions, JumpProcesses, Plots
burstyrn = @reaction_network burstyrn begin
    k₊, G₋ --> G₊
    k₋*P^2, G₊ --> G₋
    kₚ, G₊ --> G₊ + $m*P
    γₚ, P --> ∅
end
u0 = [:G₊ => 1, :G₋ => 0, :P => 1]
tspan = (0., 6.0) 
p = [:k₊ => 0.05, :k₋ => 0.001, :kₚ => 2288.5714285714284, :γₚ => 1, :b => 70]
dprob = DiscreteProblem(burstyrn, u0, tspan, u0)
jprob = JumpProblem(burstyrn, dprob, Direct())
sol = solve(jprob, SSAStepper())
plot(sol, idxs=:P, legend = false, xlabel = "time", ylabel = "P(t)")

but this fail. It seems like something we could make work.

General parameters as distributions.

Ultimately, it seems it would be cool if one could just declare a normal model, and then at the end designate a parameter according to some distribution:

using Catalyst, Distributions, OrdinaryDiffEq, Plots
rn = @reaction_network begin
    (p,d), 0 <--> X
end
u0 = [:X => 0.0]
tspan = (0.0, 10.0)
ps = [:p => Uniform(0.8,1.2), :d => 0.1]
oprob = ODEProblem(rn, u0, tspan, ps)
sol = solve(oprob)
plot(sol)

Then each time oprob is simulated, a new p is drawn from its distribution. This would be an incredibly good way to model extrinsic noise.

[TorkelE]

(since this is mostly a ModelingToolkit issue, I have opened an issue there as well: https://github.com/SciML/ModelingToolkit.jl/issues/2293)

[isaacsas]

One issue that has to be considered is if the random parameter is only supposed to be sampled once at the start of a simulation, and then it is fixed, vs. if it is randomly sampled each time it is evaluated. Bursty transcription needs to be resampled during the course of a simulation to model the different number of transcripts each time transcription occurs.

[TorkelE]

I agree that there are two types of distributions here. For the stoichiometric coefficients, these should remain distributions. Parameters in rate, they should be drawn at the beginning.

[isaacsas]

I don't see why a user might not want to have such terms in rates, or coupled equations either potentially. And having a random stoichiometric coefficient fixed at the start of a simulation could also make sense (for example, to model a cell with a random number of plasmid copies).