From 92d6437ef94a3087616140e538ea5cea4c118f35 Mon Sep 17 00:00:00 2001 From: Nikhil Yewale Date: Tue, 20 Apr 2021 02:49:17 +0530 Subject: [PATCH] Update networkapi.jl updating networkapi.jl with new API functions based on intermediate complexes in reactions --- src/networkapi.jl | 146 ++++++++++++++++++++++++++++++++-------------- 1 file changed, 103 insertions(+), 43 deletions(-) diff --git a/src/networkapi.jl b/src/networkapi.jl index 460791f280..807f93702a 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -184,61 +184,121 @@ end """ - complexstoichmat(rn; smap=speciesmap(rn)) -Returns, -Z --> complex composition matrix/complex stoichiometeric matrix -B --> Incidence matrix -Δ --> Outgoing matrix -kd --> Array of rate of reactions -complexes_mat --> unique intermediate complexes that define the system, - -We can define a laplacian matrix from these returned arrays from complexstoichmat as, -KD = diagm((map(Num, kd))) -L = B*KD*Δ', and rewrite the ODESystem as x' = -Z*L*exp.(Z'*log.(states(rn))) - -and/or also check that, netstoichmat(rn) == Z*B +Function reaction_complexes(rn) +Returns complexes of reaction network as +`complexes` and, `complexes_mat` , both are useful representations """ -function complexstoichmat(rn; smap=speciesmap(rn)) +function reaction_complexes(rn) + complexes_mat = []; complexes = []; r = reactions(rn); - nr = numreactions(rn); ns = numspecies(rn); - - complexes_mat = []; - for i ∈ 1:numreactions(rn) + for i in 1:numreactions(rn) push!(complexes_mat, [r[i].substrates r[i].substoich]) push!(complexes_mat, [r[i].products r[i].prodstoich]) - end - complexes_mat = unique(complexes_mat) - nc = length(complexes_mat) # Number of unique complexes - - # complex composition matrix/complex stoichiometeric matrix - Z = zeros(Int, ns, nc); - for i ∈ 1:nc - for j in 1:length(complexes_mat[i][:,1]) - Z[smap[complexes_mat[i][:,1][j]] ,i] = complexes_mat[i][:,2][j] + if r[i].substrates == Any[] + push!(complexes, 0) + push!(complexes, sum(r[i].products.*r[i].prodstoich)) + elseif r[i].products == Any[] + push!(complexes, sum(r[i].substrates.*r[i].substoich)) + push!(complexes, 0) + else + push!(complexes, sum(r[i].substrates.*r[i].substoich)) + push!(complexes, sum(r[i].products.*r[i].prodstoich)) end end + return [complexes_mat[unique(i -> complexes[i] , eachindex(complexes))] unique(complexes)] +end - B = zeros(Int, nc, nr); # Incidence matrix +""" +Returns rates of reactions +""" +function reaction_rates(rn) kd = []; - for i ∈ 1:nr - for j ∈ 1:nc - if isequal([r[i].substrates r[i].substoich], complexes_mat[j]) - B[j, i] = -1; # found the reactant as complexes - continue; - end - if isequal([r[i].products r[i].prodstoich], complexes_mat[j]) - B[j, i] = 1; # found the product as complexes - continue; - end - end + r =reactions(rn) + for i in 1:numreactions(rn) push!(kd, r[i].rate) end + return kd +end + +""" +complex_stoich_matrix(rn; cmp_mat = reaction_complexes(rn)[:,1],smap = speciesmap(rn)) +Returns a complex stooichiometric matrix..all entries are non-negative + +""" +function complex_stoich_matrix(rn; cmp_mat = reaction_complexes(rn)[:,1], + smap = speciesmap(rn)) + Z = zeros(Int, numspecies(rn), length(cmp_mat)); + for i in 1:length(cmp_mat) + for j in 1:length(cmp_mat[i][:,1]) + Z[smap[cmp_mat[i][:,1][j]] ,i] = cmp_mat[i][:,2][j] + end + end + return Z +end + +""" +complex_incidence_matrix(rn; cmp_mat = reaction_complexes(rn)[:,2]) +Returns a complex incidence matrix.. defined as, +Bᵢⱼ = -1, if the ith complex is the substrate of the jth reaction, + = 1, if the ith complex is the product of the jth reaction, + = 0, otherwise + +""" +function complex_incidence_matrix(rn; cmp_mat = reaction_complexes(rn)[:,2]) + r = reactions(rn); + B = zeros(Int, length(cmp_mat), numreactions(rn)); + for i in 1:numreactions(rn) + if r[i].substrates == Any[] && r[i].products != Any[] + sub_index = findall(x -> isequal(x ,0) ,cmp_mat)[1] + + prod_index = findall(x -> isequal(x, sum(r[i].products.*r[i].prodstoich)) + ,cmp_mat)[1] + + elseif r[i].substrates != Any[] && r[i].products == Any[] + sub_index = findall(x -> isequal(x, sum(r[i].substrates.*r[i].substoich)) + ,cmp_mat)[1] + prod_index = findall(x -> isequal(x ,0) ,cmp_mat)[1] + + elseif (r[i].substrates != Any[] && r[i].products != Any[]) + sub_index = findall(x -> isequal(x, sum(r[i].substrates.*r[i].substoich)) + ,cmp_mat)[1] + + prod_index = findall(x -> isequal(x, sum(r[i].products.*r[i].prodstoich)) + ,cmp_mat)[1] + end + B[sub_index, i] = -1; + B[prod_index, i] = 1; + end + return B +end - ## Outgoing matrix - Δ = copy(B) - Δ[Δ.==1] .= 0 - return Z, B , Δ ,kd ,complexes_mat +""" +complex_outgoing_matrix(rn; cmp_mat = reaction_complexes(rn)[:,2]) +Returns a complex_outgoing_matri defined as, +Δᵢⱼ = 0 , if Bᵢⱼ = 1 + = Bᵢⱼ , otherwise +""" +function complex_outgoing_matrix(rn, cmp_mat = reaction_complexes(rn)[:,2]) + r = reactions(rn); + Δ = zeros(Int, length(cmp_mat), numreactions(rn)); + for i in 1:numreactions(rn) + if r[i].substrates == Any[] + sub_index = findall(x -> isequal(x ,0) ,cmp_mat)[1] + else + sub_index = findall(x -> isequal(x, sum(r[i].substrates.*r[i].substoich)) + ,cmp_mat)[1] + end + Δ[sub_index, i] = -1; + end + return Δ end +""" +We can define a laplacian matrix from these returned arrays from complexstoichmat as, +KD = diagm((map(Num, kd))) +L = B*KD*Δ', and rewrite the ODESystem as x' = -Z*L*exp.(Z'*log.(states(rn))) +and/or also check that, netstoichmat(rn) == Z*B +""" + ######################## reaction network operators #######################