diff --git a/ext/CatalystGraphMakieExtension.jl b/ext/CatalystGraphMakieExtension.jl index ae409b645e..a47ef7f5e5 100644 --- a/ext/CatalystGraphMakieExtension.jl +++ b/ext/CatalystGraphMakieExtension.jl @@ -1,11 +1,11 @@ module CatalystGraphMakieExtension # Fetch packages. -using Catalyst, GraphMakie -import Catalyst: lattice_plot, lattice_animation, extract_vals -import Graphs: AbstractGraph, SimpleGraph +using Catalyst, GraphMakie, Graphs, Symbolics +using Symbolics: get_variables! +import Catalyst: species_reaction_graph, incidencematgraph, lattice_plot, lattice_animation -# Creates and exports hc_steady_states function. +# Creates and exports graph plotting functions. include("CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl") - +include("CatalystGraphMakieExtension/rn_graph_plot.jl") end diff --git a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl new file mode 100644 index 0000000000..644a175b32 --- /dev/null +++ b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl @@ -0,0 +1,185 @@ +############# +# Adapted from https://github.com/MakieOrg/GraphMakie.jl/issues/52#issuecomment-1018527479 +############# + +""" + SRGraphWrap{T} + +Wrapper for the species-reaction graph containing edges for rate-dependence on species. Intended to allow plotting of multiple edges. +""" +struct SRGraphWrap{T} <: Graphs.AbstractGraph{T} + g::SimpleDiGraph{T} + rateedges::Vector{Graphs.SimpleEdge{T}} + edgeorder::Vector{Int64} +end + +# Create the SimpleDiGraph corresponding to the species and reactions +function SRGraphWrap(rn::ReactionSystem) + srg = species_reaction_graph(rn) + rateedges = Vector{Graphs.SimpleEdge{Int}}() + sm = speciesmap(rn); specs = species(rn) + + deps = Set() + for (i, rx) in enumerate(reactions(rn)) + empty!(deps) + get_variables!(deps, rx.rate, specs) + if !isempty(deps) + for spec in deps + specidx = sm[spec] + push!(rateedges, Graphs.SimpleEdge(specidx, i + length(specs))) + end + end + end + edgelist = vcat(collect(Graphs.edges(srg)), rateedges) + edgeorder = sortperm(edgelist) + SRGraphWrap(srg, rateedges, edgeorder) +end + +Base.eltype(g::SRGraphWrap) = eltype(g.g) +Graphs.edgetype(g::SRGraphWrap) = edgetype(g.g) +Graphs.has_edge(g::SRGraphWrap, s, d) = has_edge(g.g, s, d) +Graphs.has_vertex(g::SRGraphWrap, i) = has_vertex(g.g, i) +Graphs.inneighbors(g::SRGraphWrap{T}, i) where T = inneighbors(g.g, i) +Graphs.outneighbors(g::SRGraphWrap{T}, i) where T = outneighbors(g.g, i) +Graphs.ne(g::SRGraphWrap) = length(g.rateedges) + length(Graphs.edges(g.g)) +Graphs.nv(g::SRGraphWrap) = nv(g.g) +Graphs.vertices(g::SRGraphWrap) = vertices(g.g) +Graphs.is_directed(g::SRGraphWrap) = is_directed(g.g) + +function Graphs.edges(g::SRGraphWrap) + edgelist = vcat(collect(Graphs.edges(g.g)), g.rateedges)[g.edgeorder] +end + +function gen_distances(g::SRGraphWrap; inc = 0.2) + edgelist = edges(g) + distances = zeros(length(edgelist)) + for i in 2:Base.length(edgelist) + edgelist[i] == edgelist[i-1] && (distances[i] = inc) + end + distances +end + +""" + plot_network(rn::ReactionSystem; interactive=false) + +Converts a [`ReactionSystem`](@ref) into a GraphMakie plot of the species reaction graph +(or Petri net representation). Reactions correspond to small green circles, and +species to blue circles. + +Notes: +- Black arrows from species to reactions indicate reactants, and are labelled + with their input stoichiometry. +- Black arrows from reactions to species indicate products, and are labelled + with their output stoichiometry. +- Red arrows from species to reactions indicate that species is used within the + rate expression. For example, in the reaction `k*A, B --> C`, there would be a + red arrow from `A` to the reaction node. In `k*A, A+B --> C`, there would be + red and black arrows from `A` to the reaction node. +""" +# TODO: update docs for interacting with plots. The `interactive` flag sets the ability to interactively drag nodes and edges in the generated plot. Only allowed if `GLMakie` is the loaded Makie backend. +function Catalyst.plot_network(rn::ReactionSystem) + srg = SRGraphWrap(rn) + ns = length(species(rn)) + nodecolors = vcat([:skyblue3 for i in 1:ns], + [:green for i in ns+1:nv(srg)]) + ilabels = vcat(map(s -> String(tosymbol(s, escape=false)), species(rn)), + fill("", nv(srg.g) - ns)) + nodesizes = vcat([30 for i in 1:ns], + [10 for i in ns+1:nv(srg)]) + + ssm = substoichmat(rn); psm = prodstoichmat(rn) + # Get stoichiometry of reaction + edgelabels = map(Graphs.edges(srg.g)) do e + string(src(e) > ns ? + psm[dst(e), src(e)-ns] : + ssm[src(e), dst(e)-ns]) + end + edgecolors = [:black for i in 1:ne(srg)] + + num_e = ne(srg.g) + for i in 1:length(srg.edgeorder) + if srg.edgeorder[i] > num_e + edgecolors[i] = :red + insert!(edgelabels, i, "") + end + end + + graphplot(srg; + edge_color = edgecolors, + elabels = edgelabels, + elabels_rotation = 0, + ilabels = ilabels, + node_color = nodecolors, + node_size = nodesizes, + arrow_shift = :end, + arrow_size = 20, + curve_distance_usage = true, + curve_distance = gen_distances(srg) + ) +end + +""" + plot_complexes(rn::ReactionSystem; interactive=false) + + Creates a GraphMakie plot of the [`ReactionComplex`](@ref)s in `rn`. Reactions + correspond to arrows and reaction complexes to blue circles. + + Notes: + - Black arrows from complexes to complexes indicate reactions whose rate is a + parameter or a `Number`. i.e. `k, A --> B`. + - Red arrows from complexes to complexes indicate reactions whose rate + depends on species. i.e. `k*C, A --> B` for `C` a species. +""" +function Catalyst.plot_complexes(rn::ReactionSystem) + img = incidencematgraph(rn) + specs = species(rn); rxs = reactions(rn) + edgecolors = [:black for i in 1:ne(img)] + nodelabels = complexlabels(rn) + edgelabels = [repr(rx.rate) for rx in rxs] + + deps = Set() + for (i, rx) in enumerate(rxs) + empty!(deps) + get_variables!(deps, rx.rate, specs) + (!isempty(deps)) && (edgecolors[i] = :red) + end + + graphplot(img; + edge_color = edgecolors, + elabels = edgelabels, + ilabels = complexlabels(rn), + node_color = :skyblue3, + elabels_rotation = 0, + arrow_shift = :end, + curve_distance = 0.2 + ) +end + +function complexelem_tostr(e::Catalyst.ReactionComplexElement, specstrs) + if e.speciesstoich == 1 + return "$(specstrs[e.speciesid])" + else + return "$(e.speciesstoich)$(specstrs[e.speciesid])" + end +end + +# Get the strings corresponding to the reaction complexes +function complexlabels(rn::ReactionSystem) + labels = String[] + + specstrs = map(s -> String(tosymbol(s, escape=false)), species(rn)) + complexes, B = reactioncomplexes(rn) + + for complex in complexes + if isempty(complex) + push!(labels, "∅") + elseif length(complex) == 1 + push!(labels, complexelem_tostr(complex[1], specstrs)) + else + elems = map(c -> complexelem_tostr(c, specstrs), complex) + str = reduce((e1, e2) -> *(e1, " + ", e2), @view elems[2:end]; init = elems[1]) + push!(labels, str) + end + end + labels +end diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 047747da71..ad12f0fa4b 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -167,6 +167,11 @@ export hc_steady_states function make_si_ode end export make_si_ode +# GraphMakie +function plot_network end +function plot_complexes end +export plot_network, plot_complexes + ### Spatial Reaction Networks ### # Spatial reactions. diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 93bebe1888..86802fb919 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -314,6 +314,44 @@ function incidencematgraph(incidencemat::SparseMatrixCSC{Int, Int}) return graph end + +""" + species_reaction_graph(rn::ReactionSystem) + +Construct a directed simple graph where there are two types of nodes: species and reactions. +An edge from a species *s* to reaction *r* indicates that *s* is a reactant for *r*, and +an edge from a reaction *r* to a species *s* indicates that *s* is a product of *r*. By +default, the species vertices are listed first, so the first *n* indices correspond to +species nodes. + +Note: this is equivalent to the Petri net representation of a chemical reaction network. + +For example, +```julia +sir = @reaction_network SIR begin + β, S + I --> 2I + ν, I --> R +end +species_reaction_graph(sir) +""" +function species_reaction_graph(rn::ReactionSystem) + specs = species(rn) + rxs = reactions(rn) + sm = speciesmap(rn) + s = length(specs) + + edgelist = Graphs.Edge[] + for (i, rx) in enumerate(rxs) + for spec in rx.substrates + push!(edgelist, Graphs.Edge(sm[spec], s+i)) + end + for spec in rx.products + push!(edgelist, Graphs.Edge(s+i, sm[spec])) + end + end + srg = Graphs.SimpleDiGraphFromIterator(edgelist) +end + ### Linkage, Deficiency, Reversibility ### """ diff --git a/test/extensions/graphmakie.jl b/test/extensions/graphmakie.jl new file mode 100644 index 0000000000..e74b6ce48b --- /dev/null +++ b/test/extensions/graphmakie.jl @@ -0,0 +1,106 @@ +using Catalyst, GraphMakie, GLMakie, Graphs +include("../test_networks.jl") + +# Test that speciesreactiongraph is generated correctly +let + brusselator = @reaction_network begin + A, ∅ --> X + 1, 2X + Y --> 3X + B, X --> Y + 1, X --> ∅ + end + + srg = Catalyst.species_reaction_graph(brusselator) + s = length(species(brusselator)) + edgel = Graphs.Edge.([(s+1, 1), + (1, s+2), + (2, s+2), + (s+2, 1), + (s+3, 2), + (1, s+3), + (1, s+4)]) + @test all(∈(collect(Graphs.edges(srg))), edgel) + + MAPK = @reaction_network MAPK begin + (k₁, k₂),KKK + E1 <--> KKKE1 + k₃, KKKE1 --> KKK_ + E1 + (k₄, k₅), KKK_ + E2 <--> KKKE2 + k₆, KKKE2 --> KKK + E2 + (k₇, k₈), KK + KKK_ <--> KK_KKK_ + k₉, KK_KKK_ --> KKP + KKK_ + (k₁₀, k₁₁), KKP + KKK_ <--> KKPKKK_ + k₁₂, KKPKKK_ --> KKPP + KKK_ + (k₁₃, k₁₄), KKP + KKPase <--> KKPKKPase + k₁₅, KKPPKKPase --> KKP + KKPase + k₁₆,KKPKKPase --> KK + KKPase + (k₁₇, k₁₈), KKPP + KKPase <--> KKPPKKPase + (k₁₉, k₂₀), KKPP + K <--> KKPPK + k₂₁, KKPPK --> KKPP + KP + (k₂₂, k₂₃), KKPP + KP <--> KPKKPP + k₂₄, KPKKPP --> KPP + KKPP + (k₂₅, k₂₆), KP + KPase <--> KPKPase + k₂₇, KKPPKPase --> KP + KPase + k₂₈, KPKPase --> K + KPase + (k₂₉, k₃₀), KPP + KPase <--> KKPPKPase + end + srg = Catalyst.species_reaction_graph(MAPK) + @test nv(srg) == length(species(MAPK)) + length(reactions(MAPK)) + @test ne(srg) == 90 + + # Test that figures are generated properly. + f = plot_network(MAPK) + save("fig.png", f) + @test isfile("fig.png") + rm("fig.png") + f = plot_network(brusselator) + save("fig.png", f) + @test isfile("fig.png") + rm("fig.png") + + f = plot_complexes(MAPK); save("fig.png", f) + @test isfile("fig.png") + rm("fig.png") + f = plot_complexes(brusselator); save("fig.png", f) + @test isfile("fig.png") + rm("fig.png") +end + +CGME = Base.get_extension(parentmodule(ReactionSystem), :CatalystGraphMakieExtension) +# Test that rate edges are inferred correctly. We should see two for the following reaction network. +let + # Two rate edges, one to species and one to product + rn = @reaction_network begin + k, A --> B + k * C, A --> C + k * B, B --> C + end + srg = CGME.SRGraphWrap(rn) + s = length(species(rn)) + @test ne(srg) == 8 + @test Graphs.Edge(3, s+2) ∈ srg.rateedges + @test Graphs.Edge(2, s+3) ∈ srg.rateedges + # Since B is both a dep and a reactant + @test count(==(Graphs.Edge(2, s+3)), edges(srg)) == 2 + + f = plot_network(rn) + save("fig.png", f) + @test isfile("fig.png") + rm("fig.png") + f = plot_complexes(rn); save("fig.png", f) + @test isfile("fig.png") + rm("fig.png") + + # Two rate edges, both to reactants + rn = @reaction_network begin + k, A --> B + k * A, A --> C + k * B, B --> C + end + srg = CGME.SRGraphWrap(rn) + s = length(species(rn)) + @test ne(srg) == 8 + # Since A, B is both a dep and a reactant + @test count(==(Graphs.Edge(1, s+2)), edges(srg)) == 2 + @test count(==(Graphs.Edge(2, s+3)), edges(srg)) == 2 +end +