add infinitesimal generator or Fokker-Planck equation

[rveltz]

Hi,

It would be nice to be able to generate / solve for the law of a Markov process directly in Catalyst. In the case of Continuous time MC, it would make only sense for finite state space.

One way to this goal is to generate the infinitesimal generator. I come up with this code which is not yet working because I need to know how to instantiate the rates with parameters.

function getInfinitesimalGenerator(rn, par; smap=speciesmap(rn))
    A = zeros(Any, (numspecies(rn), numspecies(rn)))
    for (k,rx) in enumerate(reactions(rn))
        @assert length(rx.substrates) == 1
        @assert length(rx.products) == 1
        i = smap[rx.substrates[1]]
        j = smap[rx.products[1]]
        A[i, j] = rx.rate
    end
    A
end

[isaacsas]

Ideally we'd represent the full (possibly infinite) system of equations for the chemical master equation symbolically, but then specialize to a concrete system based on the initial condition or a closure condition.

[isaacsas]

Also note, rx.rate is not the full state-dependent transition rate for a multiparticle reaction like A --> B, which should be rx.rate*A(t). This is given by the associated jumpratelaw.

[isaacsas]

@rveltz once you've got a symbolic representation for the generator you should be able to create a ModelingToolkit ODESystem with the corresponding ODEs, from which you can instantiate it by calling ODEProblem on the ODESystem, see the ModelingToolkit docs. You can pass in your parameter mapping at that point.

[rveltz]

Ideally we'd represent the full (possibly infinite) system of equations for the chemical master equation symbolically, but then specialize to a concrete system based on the initial condition or a closure condition.

That seemed to difficult for me to ask for an issue :D

[ChrisRackauckas]

Yep, mentioned in https://github.com/SciML/Catalyst.jl/issues/792. It would be good to target PDESystem and try to keep it infinite, with lowering methods.

[isaacsas]

In this case it seems like what we want is a differential-difference equation system type to encode the infinite system (or can the PDE system type handle that?).

[rveltz]

Also note, rx.rate is not the full state-dependent transition rate for a multiparticle reaction like A --> B, which should be rx.rate*A(t). This is given by the associated jumpratelaw.

I was not precise enough here. What I want to do is simulate a Markov chain (MC) for which I specify the transitions using the DSL because it is convenient. The MC does not represent chemical reactions, so I want the generator to have "my" rates.

@rveltz once you've got a symbolic representation for the generator you should be able to create a ModelingToolkit ODESystem with the corresponding ODEs, from which you can instantiate it by calling ODEProblem on the ODESystem, see the ModelingToolkit docs. You can pass in your parameter mapping at that point.

That would be good but I'd be interested in having the matrix as well. Can you please give me a push on how to convert the symbolics rates to Floats given the parameters please?

[isaacsas]

Got it. I'll try to put an example together for you.

[rveltz]

Thank you!

[isaacsas]

Is this what you are looking for?

using Catalyst, ModelingToolkit, Symbolics, LinearAlgebra

rn = @reaction_network begin
    k1, S1 --> S2
    k2, S2 --> S1
    k3, S3 --> S1
end k1 k2 k3
pvals = [1.0,2.0,3.0]
pmap  = Dict(Pair.(params(rn),pvals))

osys = convert(ODESystem, rn)
eqs  = getfield.(osys.eqs,:rhs)
A = Symbolics.jacobian(eqs, osys.states)'

# method one to get a numerical generator
Avals = zeros(size(A))
for j = 1:size(A,2)
    for i = 1:size(A,1)
        Avals[i,j] = Symbolics.value(Symbolics.substitute(A[i,j],pmap))
    end
end

# method two to get a function that evaluates the numerical generator
# [1] for oop and [2] for iip...
genf = eval(build_function(A, params(rn))[1])
Avals2 = genf(pvals)

norm(Avals .- Avals2) == 0