up

Author: vyudu

src/network_analysis.jl

@@ -195,8 +195,8 @@ end
 """
     laplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false)
 
-    Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian, and Φ is the [`massactionvector`](@ref). A_k is an n-by-n matrix, where n is the number of complexes, where A_{ij} = k_{ij} if a reaction exists between the two complexes and 0 otherwise. 
-    Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via parammap. 
+    Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian, and Φ is the [`massactionvector`](@ref). A_k is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. 
+    Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
 """
 function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) 
     D = incidencemat(rn; sparse); K = fluxmat(rn, pmap; sparse)
@@ -212,19 +212,13 @@ end
 function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) 
     rates = reactionrates(rn)
 
-    if !isempty(pmap)
-        length(pmap) != length(parameters(rn)) && error("Incorrect number of parameters were specified.")
-        pmap = symmap_to_varmap(rn, pmap)
-        pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
-        rates = [substitute(rate, pmap) for rate in rates]
-    end
-
+    !isempty(pmap) && (rates = substitutevals(rn, pmap, parameters(rn), rates))
     rcmap = reactioncomplexmap(rn); nc = length(rcmap); nr = length(rates)
     mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)
-    K = if sparse
-        fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates)
+    if sparse
+        return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates)
     else
-        fluxmat(Matrix{mtype}, rcmap, rates)
+        return fluxmat(Matrix{mtype}, rcmap, rates)
     end
 end
 
@@ -265,47 +259,53 @@ function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false)
     fluxmat(rn, pdict; sparse)
 end
 
+# Helper to substitute values into a (vector of) symbolic expressions. The syms are the symbols to substitute and the symexprs are the expressions to substitute into.
+function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs) 
+    length(map) != length(syms) && error("Incorrect number of parameter-value pairs were specified.")
+    map = symmap_to_varmap(rn, map)
+    map = Dict(ModelingToolkit.value(k) => v for (k, v) in map)
+    vals = [substitute(expr, map) for expr in symexprs]
+end
+
 """
-    massactionvector(rn::ReactionSystem, scmap = Dict())
+    massactionvector(rn::ReactionSystem, scmap = Dict(); combinatoric_ratelaws = true)
 
     Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be A*B, and for the complex 2X + Y, the corresponding entry would be X^2*Y. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the [`networklaplacian`](@ref). This utility returns Φ(x).
     Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap.
+    If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws).
 """
-function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) 
+function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = true)
     r = numreactions(rn); rxs = reactions(rn)
     sm = speciesmap(rn); specs = species(rn)
 
     if !all(r -> ismassaction(r, rn), rxs)
         error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.")
     end
 
-    if !isempty(scmap)
-        length(scmap) != length(specs) && error("Incorrect number of species concentrations were specified.")
-        scmap = symmap_to_varmap(rn, scmap)
-        scmap = Dict(ModelingToolkit.value(k) => v for (k, v) in scmap)
-        specs = [substitute(s, scmap) for s in specs]
-    end
+    !isempty(scmap) && (specs = substitutevals(rn, scmap, specs, specs))
 
     vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Num : eltype(specs)
     Φ = Vector{vtype}()
     rcmap = reactioncomplexmap(rn)
     for comp in keys(reactioncomplexmap(rn))
         subs = map(ce -> getfield(ce, :speciesid), comp)
         stoich = map(ce -> getfield(ce, :speciesstoich), comp)
-        push!(Φ, prod(vtype[specs[s]^α for (s, α) in zip(subs, stoich)]))
+        maprod = prod(vtype[specs[s]^α for (s, α) in zip(subs, stoich)])
+        combinatoric_ratelaws && (maprod /= prod(map(factorial, stoich)))
+        push!(Φ, maprod)
     end
 
     Φ 
 end
 
-function massactionvector(rn::ReactionSystem, scmap::Tuple) 
+function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = true) 
     sdict = Dict(scmap)
-    massactionvector(rn, sdict)
+    massactionvector(rn, sdict; combinatoric_ratelaws)
 end
 
-function massactionvector(rn::ReactionSystem, scmap::Vector) 
+function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = true) 
     sdict = Dict(scmap)
-    massactionvector(rn, sdict)
+    massactionvector(rn, sdict; combinatoric_ratelaws)
 end
 
 @doc raw"""

test/network_analysis/network_properties.jl

@@ -350,11 +350,44 @@ let
 
     Φ = Catalyst.massactionvector(MAPK)
     specs = species(MAPK)
-    @test isequal(Φ[1], MAPK.KKK * MAPK.E1)
-    @test isequal(Φ[6], MAPK.KKK * MAPK.E2)
-    @test isequal(Φ[8], MAPK.KK_KKK_)
+    truevec = [MAPK.KKK * MAPK.E1,
+               MAPK.KKKE1,
+               MAPK.KKK_ * MAPK.E1,
+               MAPK.KKK_ * MAPK.E2,
+               MAPK.KKKE2,
+               MAPK.KKK * MAPK.E2,
+               MAPK.KK * MAPK.KKK_,
+               MAPK.KK_KKK_,
+               MAPK.KKP * MAPK.KKK_,
+               MAPK.KKPKKK_,
+               MAPK.KKPP * MAPK.KKK_,
+               MAPK.KKP * MAPK.KKPase,
+               MAPK.KKPKKPase,
+               MAPK.KKPPKKPase,
+               MAPK.KK * MAPK.KKPase,
+               MAPK.KKPP * MAPK.KKPase,
+               MAPK.KKPP * MAPK.K,
+               MAPK.KKPPK,
+               MAPK.KKPP * MAPK.KP, 
+               MAPK.KPKKPP,
+               MAPK.KPP * MAPK.KKPP,
+               MAPK.KP * MAPK.KPase,
+               MAPK.KPKPase, 
+               MAPK.KKPPKPase, 
+               MAPK.K * MAPK.KPase,
+               MAPK.KPP * MAPK.KPase,
+              ]
+    @test isequal(Φ, truevec)
 
     K = fluxmat(MAPK)
+    # Construct matrix from incidence matrix
+    mat = zeros(Num, 30, 26); D = incidencemat(MAPK)
+    rates = reactionrates(MAPK)
+    for (i, col) in enumerate(eachcol(D))
+        sub = findfirst(==(-1), col)
+        mat[i, sub] = rates[i]
+    end
+    @test isequal(K, mat)
     @test isequal(K[1, 1], MAPK.k₁)
     @test all(==(0), K[1, 2:end])
     @test isequal(K[2, 2], MAPK.k₂) 
@@ -400,23 +433,75 @@ let
     @test all(iszero, simplify(numeqs - Y*A_k*Φ))
 end
 
-# Test handling for weird complexes. 
+# Test handling for weird complexes and combinatoric rate laws. 
 let
     rn = @reaction_network begin
         k1, 2X + Y + 3Z --> ∅
+        (k2, k3), 2Y + 2Z <--> 3X
+    end
+
+    Φ = Catalyst.massactionvector(rn); specs = species(rn)
+    crvec = [rn.X^2/2 * rn.Y * rn.Z^3/6,
+             1.,
+             rn.Y^2/2 * rn.Z^2/2,
+             rn.X^3/6]
+    @test isequal(Φ, crvec)
+    ncrvec = [rn.X^2 * rn.Y * rn.Z^3,
+              1.,
+              rn.Y^2 * rn.Z^2,
+              rn.X^3]
+    Φ_2 = Catalyst.massactionvector(rn; combinatoric_ratelaws = false)
+    @test isequal(Φ_2, ncrvec)
+
+    # Test that the ODEs generated are the same.
+    eqs = Catalyst.assemble_oderhs(rn, specs)
+    S = netstoichmat(rn); Y = complexstoichmat(rn); K = fluxmat(rn); A_k = laplacianmat(rn)
+    @test all(iszero, simplify(eqs - S*K*Φ))
+    @test all(iszero, simplify(eqs - Y*A_k*Φ))
+
+    eq_ncr = Catalyst.assemble_oderhs(rn, specs; combinatoric_ratelaws = false)
+    @test all(iszero, simplify(eq_ncr - S*K*Φ_2))
+    @test all(iszero, simplify(eq_ncr - Y*A_k*Φ_2))
+
+    # Test that the ODEs with rate constants are the same.
+    k = rand(rng, numparams(rn))
+    ratevec = collect(zip(parameters(rn), k))
+    ratemap = Dict(ratevec)
+    K = fluxmat(rn, ratemap); A_k = laplacianmat(rn, ratemap)
+
+    numeqs = similar(eqs)
+    for i in 1:length(eqs)
+        numeqs[i] = substitute(eqs[i], ratemap)
+    end
+    @test_broken all(iszero, simplify(numeqs - S*K*Φ))
+    @test_broken all(iszero, simplify(numeqs - Y*A_k*Φ))
+
+    numeqs_ncr = similar(eq_ncr)
+    for i in 1:length(eq_ncr)
+        numeqs_ncr[i] = substitute(eq_ncr[i], ratemap)
     end
-    Φ = Catalyst.massactionvector(rn)
-    @test isequal(Φ[1], rn.X^2 * rn.Y * rn.Z^3)
-    @test isequal(Φ[2], 1)
+    @test all(iszero, simplify(numeqs_ncr - S*K*Φ_2))
+    @test all(iszero, simplify(numeqs_ncr - Y*A_k*Φ_2))
 
+    # Test that handling of species concentrations is correct.
     u0vec = [:X => 3., :Y => .5, :Z => 2.]
     u0map = Dict(u0vec)
     u0tup = Tuple(u0vec)
 
     Φ = Catalyst.massactionvector(rn, u0vec)
-    @test isequal(Φ[1], 36.)
-    Φ = Catalyst.massactionvector(rn, u0tup)
+    @test isequal(Φ[1], 3.)
+    Φ = Catalyst.massactionvector(rn, u0tup; combinatoric_ratelaws = false)
     @test isequal(Φ[1], 36.)
     Φ = Catalyst.massactionvector(rn, u0map)
-    @test isequal(Φ[1], 36.)
+    @test isequal(Φ[1], 3.)
+
+    # Test full simplification.
+    u0map = symmap_to_varmap(rn, u0map)
+    numeqs = [substitute(eq, u0map) for eq in numeqs]
+    @test isapprox(numeqs, S*K*Φ)
+    @test isapprox(numeqs, Y*A_k*Φ)
+
+    numeqs_ncr = [substitute(eq, u0map) for eq in numeqs_ncr]
+    @test isapprox(numeqs_ncr - S*K*Φ_2)
+    @test isapprox(numeqs_ncr - Y*A_k*Φ_2)
 end