add API function to generate function for evaluating rate laws

[isaacsas]

See https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049

[isaacsas]

This would also allow users to put together stoichiometry matrix + rate law vector representations for the ODE rhs... I think such a network representation would be useful for many types of network analyses.

[yewalenikhil65]

See https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049

This is similar to once we discussed here https://github.com/SciML/Catalyst.jl/issues/269#issuecomment-745782394 I agree that it will be nice to have an API function directly. I can suggest a small example for this that I had previously done for by representing ODE RHS using stoichiometry matrix and rate law vector

using Catalyst, ModelingToolkit
rn = @reaction_network begin
           k1, A + B --> C
           k2, A --> 2*C
           k3, 2*C --> A
       end k1 k2 k3
##  generally ode RHS can be written as X' = R*rate_laws_vector
R = -substoichmat(rn)' + prodstoichmat(rn)'   # where R is transpose of net_stoich matrix

# finding RHS of ODE at some particular time instance, using solution vector obtained from solving ODEProblem say `u` at that time instance and `p` as rate constants
u = [1;2;3]; p = [0.1;0.2;0.3] same as [k1;k2;k3]

## this block can be put in suitable function
 D = [oderatelaw(rn.eqs[i]) for i in 1:nr]
 expr = build_function(D, vcat([rn.states[i] for i in 1:ns], [rn.eqs[i].rate for i in 1:nr]))
 myf = eval(expr[2]);
 out = zeros(Float32, nr)
 myf(out, vcat(u , p))   # p is array of parameters of reaction rates
  # 'out' array contains evaluated oderatelaw from the solution and parameters
##

ODE_rhs = R*out;   # X' evaluated

[yewalenikhil65]

@isaacsas Do you think its worth adding evaloderatelaw(rn::ReactionSystem) to the PR for network_api?

using Catalyst, ModelingToolkit
using DifferentialEquations

function evaloderatelaw(rn::ReactionSystem)
    D = oderatelaw.(rn.eqs)
    expr = build_function(D, [rn.states; rn.ps])
    myf = eval(expr[2]);
    return myf
end

rn = @reaction_network begin
    k1, 2*x1 --> x2
    k2, x1 --> x3
    k3, x3 --> x4
    k4, x2 + x4 --> x5
end k1 k2 k3 k4
ns = numspecies(rn);   nr = numreactions(rn);
R = -substoichmat(rn)' + prodstoichmat(rn)';

k = Float32[0.1, 0.2, 0.13, 0.3];
tspan = (0.f0, 20.f0);      datasize = 41;
tsteps = range(tspan[1], tspan[2], length = datasize)
u0 = Float32[1; 0.5 ;0.; 0.; 0.]
oprob = ODEProblem(rn, u0, tspan, k)
sol = solve(oprob, Tsit5(), saveat = tsteps)

## ideal derivative or ODE-RHS
Dsol = similar(Array(sol))
du = copy(u0);
for j in 1:datasize
    oprob.f(du, sol.u[j], k, tsteps[j])
    Dsol[:,j] =  du
end


## checking using evaloderatelaw()
D_sol = similar(Dsol);
for j in 1:datasize
    myf = evaloderatelaw(rn)
    out = zeros(Float32,nr);
    myf(out, [sol.u[j] ; k] )
    D_sol[:,j] = R*out;
end

# check the entries of Dsol and D_sol.. they should be almost same

[isaacsas]

Yes, we'll add something like that soon. Maybe even next week. I think we actually want to have a system type that encapsulates the rate law vector combined with the net stoichiometry matrix, with all the assorted transformations.

[isaacsas]

Something like that is needed for many kinds of reaction network analysis.

[yewalenikhil65]

Yes, we'll add something like that soon. Maybe even next week. I think we actually want to have a system type that encapsulates the rate law vector combined with the net stoichiometry matrix, with all the assorted transformations.

I would feel glad to contribute for this, if you could explain the features required using some example ?

[isaacsas]

I was thinking we should add a new system type that ReactionSystems can be converted into, and which just stores a vector of the rate law expressions and the stoichiometry matrix. This is a relatively simple representation, but one that some users might prefer. It would be easy to convert to it from normal ReactionSystems, but hard to go the other direction in any optimal way. It would also be straightforward to add conversions to ODEs/SDEs/Jumps for it, though that is a lot of redundancy we would be adding...

[yewalenikhil65]

I was thinking we should add a new system type that ReactionSystems can be converted into, and which just stores a vector of the rate law expressions and the stoichiometry matrix. This is a relatively simple representation, but one that some users might prefer. It would be easy to convert to it from normal ReactionSystems, but hard to go the other direction in any optimal way. It would also be straightforward to add conversions to ODEs/SDEs/Jumps for it, though that is a lot of redundancy we would be adding...

Somthing like RateLawMapping ?

And a way to

1. [ratelawvector, netstoich] = convert(rn::ReactionSystem, RateLawMapping)
2. osys = convert(ODESystem, RateLawMapping) # thats simply netstoich*ratelawvector and similar for SDEs, ,,.. care should be taken so that we use oderatelaw and jumpratelaw

I think this should be fairly easy to do.. Yeah, i understand why u say it would be difficult to go back to ReactionSystem..

[yewalenikhil65]

I was thinking we should add a new system type that ReactionSystems can be converted into, and which just stores a vector of the rate law expressions and the stoichiometry matrix. This is a relatively simple representation, but one that some users might prefer. It would be easy to convert to it from normal ReactionSystems, but hard to go the other direction in any optimal way. It would also be straightforward to add conversions to ODEs/SDEs/Jumps for it, though that is a lot of redundancy we would be adding...

hey @isaacsas , we can discuss this later also after we improve https://github.com/SciML/Catalyst.jl/pull/359 Just noting here it down for future reference

Expressing a new system(which can be converted to ReactionSystem) as one which just stores stoichiometric matrices and rate law expressions will be fairly easy from complex_stoich_mat and complex_incidence_mat from https://github.com/SciML/Catalyst.jl/pull/359. if user specifies Z and B matrices and also rate law expressions, we can write a function similar to loadrxnetworks from https://github.com/SciML/ReactionNetworkImporters.jl/blob/master/src/parsing_routines_matrixnetworks.jl I will do it quick

[yewalenikhil65]

I was thinking we should add a new system type that ReactionSystems can be converted into, and which just stores a vector of the rate law expressions and the stoichiometry matrix. This is a relatively simple representation, but one that some users might prefer. It would be easy to convert to it from normal ReactionSystems, but hard to go the other direction in any optimal way. It would also be straightforward to add conversions to ODEs/SDEs/Jumps for it, though that is a lot of redundancy we would be adding...

hey @isaacsas , we can discuss this later also after we improve #359 Just noting here it down for future reference

Expressing a new system(which can be converted to ReactionSystem) as one which just stores stoichiometric matrices and rate law expressions will be fairly easy from complex_stoich_mat and complex_incidence_mat from #359. if user specifies Z and B matrices and also rate law expressions, we can write a function similar to loadrxnetworks from https://github.com/SciML/ReactionNetworkImporters.jl/blob/master/src/parsing_routines_matrixnetworks.jl I will do it quick

hi, Sir @isaacsas , example for this.. let me know yours thoughts whenever you are free

rns = @reaction_network begin
    k₁ , ∅ --> X₁
    ( k₂/(1 + X₁*X₂ + X₃*X₄ ), k₃/(1 + X₁*X₂ + X₃*X₄ )), 2X₁ + X₂ ↔ 3X₃ + X₄
    k₄, X₄ --> ∅
end k₁ k₂ k₃ k₄

Z = [0 1 2 0 0;
      0 0 1 0 0;
      0 0 0 3 0;
      0 0 0 1 1]    # Can also be obtained from  complex_stoich_mat(rns), once https://github.com/SciML/Catalyst.jl/pull/359 is merged
B = [-1 0 0 1;
      1 0 0 0;
      0 -1 1 0;
      0 1 -1 0;
      0 0 0 -1]     # Can also be obtained from complex_incidence_mat(rns) , once https://github.com/SciML/Catalyst.jl/pull/359

r = reactions(rns)
rateExprs = [r[i].rate for i in 1:numreactions(rns)]

# reconstructing the same ReactionSystem from Z,B and rateExprs, params(rns)
@parameters t
myspecies = [] # empty list of species
isempty(myspecies) && (myspecies = [funcsym(:X,i)(t) for i = 1:numspecies(rns)])



function LoadRxNetwork(rateexprs::AbstractVector, Z::AbstractMatrix,
                        B::AbstractMatrix, myspecies::AbstractVector,
                        pars::AbstractVector)

    numspecs, numcomp = size(Z)
    @assert all(>=(0),Z)  # All entries in Z should be >= 0
    @assert numcomp == size(B, 1)
    @assert all(∈([-1,0,1]),B)   # All entries in  B should be -1, 0 Or 1
    numrxs = size(B, 2)

    rn = Vector{Reaction{Any,Int64}}(undef,numrxs)     # create the network
    sub_indices = argmin(B, dims=1)     # cartesian indices of substrate complexes
    prod_indices = argmax(B, dims=1)    # cartesian indicies of products complexes

    for i ∈ 1:numrxs
        subStoichInd = getindex.(findall(!iszero, Z[:,sub_indices[i][1]]))
        prodStoichInd = getindex.(findall(!iszero, Z[:,prod_indices[i][1]]))

        if subStoichInd == Int64[] &&  prodStoichInd != Int64[]
            rn[i] = Reaction(rateexprs[i], nothing, myspecies[prodStoichInd],
                        nothing, Z[prodStoichInd,prod_indices[i][1]])

        elseif subStoichInd != Int64[] &&  prodStoichInd == Int64[]
            rn[i] = Reaction(rateexprs[i], myspecies[subStoichInd], nothing,
                        Z[subStoichInd,sub_indices[i][1]], nothing)
        else
            rn[i] = Reaction(rateexprs[i], myspecies[subStoichInd],
                                myspecies[prodStoichInd],
                        Z[subStoichInd,sub_indices[i][1]], Z[prodStoichInd,prod_indices[i][1]])
        end
    end

    rs = ReactionSystem(rn, t, myspecies, pars)
end
prn = LoadRxNetwork( rateExprs,Z,B,myspecies,params(rns) )
using Test
@test all(prn == rns)

[isaacsas]

Adding a new importer that takes the reaction complex representation and builds a ReactionSystem would be useful and could be added to ReactionNetworkImporters.

[yewalenikhil65]

Adding a new importer that takes the reaction complex representation and builds a ReactionSystem would be useful and could be added to ReactionNetworkImporters.

Done at https://github.com/SciML/ReactionNetworkImporters.jl/pull/53