diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md index 4a84198ee9..fb56fc0d3b 100644 --- a/docs/src/steady_state_functionality/nonlinear_solve.md +++ b/docs/src/steady_state_functionality/nonlinear_solve.md @@ -75,13 +75,13 @@ nl_prob = NonlinearProblem(two_state_model, u_guess, p; remove_conserved = true) nothing # hide ``` here it is important that the quantities used in `u_guess` correspond to the conserved quantities we wish to use. E.g. here the conserved quantity $X1 + X2 = 3.0 + 1.0 = 4$ holds for the initial condition, and will hence also hold in the computed steady state as well. We can now find the steady states using `solve` like before: -```@example steady_state_solving_claws + We note that the output only provides a single value. The reason is that the actual system solved only contains a single equation (the other being eliminated with the conserved quantity). To find the values of $X1$ and $X2$ we can [directly query the solution object for these species' values, using the species themselves as inputs](@ref simulation_structure_interfacing_solutions): -```@example steady_state_solving_claws + ## [Finding steady states through ODE simulations](@id steady_state_solving_simulation) The `NonlinearProblem`s generated by Catalyst corresponds to ODEs. A common method of solving these is to simulate the ODE from an initial condition until a steady state is reached. Here we do so for the dimerisation system considered in the previous section. First, we declare our model, initial condition, and parameter values. diff --git a/ext/CatalystGraphMakieExtension.jl b/ext/CatalystGraphMakieExtension.jl index a47ef7f5e5..f6f6665813 100644 --- a/ext/CatalystGraphMakieExtension.jl +++ b/ext/CatalystGraphMakieExtension.jl @@ -1,7 +1,7 @@ module CatalystGraphMakieExtension # Fetch packages. -using Catalyst, GraphMakie, Graphs, Symbolics +using Catalyst, GraphMakie, Graphs, Symbolics, SparseArrays using Symbolics: get_variables! import Catalyst: species_reaction_graph, incidencematgraph, lattice_plot, lattice_animation diff --git a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl index 64e3afe50b..4420124907 100644 --- a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl +++ b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl @@ -3,18 +3,22 @@ ############# """ - SRGraphWrap{T} + MultiGraphWrap{T} -Wrapper for the species-reaction graph containing edges for rate-dependence on species. Intended to allow plotting of multiple edges. +Wrapper intended to allow plotting of multiple edges. This is needed in the following cases: +- For the species-reaction graph, multiple edges can exist when a reaction depends on some species for its rate, and if that species is produced by the reaction. +- For the complex graph, multiple edges can exist between a pair of nodes if there are multiple reactions between the same complexes. This might include a reversible pair of reactions - we might have three total edges if one reaction is reversible, and we have a separate reaction going from one complex to the other. + +`gen_distances` sets the distances between the edges that allows multiple to be visible on the plot at the same time. """ -struct SRGraphWrap{T} <: Graphs.AbstractGraph{T} +struct MultiGraphWrap{T} <: Graphs.AbstractGraph{T} g::SimpleDiGraph{T} - rateedges::Vector{Graphs.SimpleEdge{T}} + multiedges::Vector{Graphs.SimpleEdge{T}} edgeorder::Vector{Int64} end -# Create the SimpleDiGraph corresponding to the species and reactions -function SRGraphWrap(rn::ReactionSystem) +# Create the SimpleDiGraph corresponding to the species and reactions, the species-reaction graph +function MultiGraphWrap(rn::ReactionSystem) srg = species_reaction_graph(rn) rateedges = Vector{Graphs.SimpleEdge{Int}}() sm = speciesmap(rn); specs = species(rn) @@ -32,29 +36,60 @@ function SRGraphWrap(rn::ReactionSystem) end edgelist = vcat(collect(Graphs.edges(srg)), rateedges) edgeorder = sortperm(edgelist) - SRGraphWrap(srg, rateedges, edgeorder) + MultiGraphWrap(srg, rateedges, edgeorder) +end + +# Automatically set edge order if not supplied +function MultiGraphWrap(g::SimpleDiGraph{T}, multiedges::Vector{Graphs.SimpleEdge{T}}) where T + edgelist = vcat(collect(Graphs.edges(g)), multiedges) + edgeorder = sortperm(edgelist) + MultiGraphWrap(g, multiedges, 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] +Base.eltype(g::MultiGraphWrap) = eltype(g.g) +Graphs.edgetype(g::MultiGraphWrap) = edgetype(g.g) +Graphs.has_edge(g::MultiGraphWrap, s, d) = has_edge(g.g, s, d) +Graphs.has_vertex(g::MultiGraphWrap, i) = has_vertex(g.g, i) +Graphs.inneighbors(g::MultiGraphWrap{T}, i) where T = inneighbors(g.g, i) +Graphs.outneighbors(g::MultiGraphWrap{T}, i) where T = outneighbors(g.g, i) +Graphs.ne(g::MultiGraphWrap) = length(g.multiedges) + length(Graphs.edges(g.g)) +Graphs.nv(g::MultiGraphWrap) = nv(g.g) +Graphs.vertices(g::MultiGraphWrap) = vertices(g.g) +Graphs.is_directed(::Type{<:MultiGraphWrap}) = true +Graphs.is_directed(g::MultiGraphWrap) = is_directed(g.g) + +function Graphs.edges(g::MultiGraphWrap) + edgelist = vcat(collect(Graphs.edges(g.g)), g.multiedges)[g.edgeorder] end -function gen_distances(g::SRGraphWrap; inc = 0.2) +function gen_distances(g::MultiGraphWrap; inc = 0.4) edgelist = edges(g) distances = zeros(length(edgelist)) - for i in 2:Base.length(edgelist) - edgelist[i] == edgelist[i-1] && (distances[i] = inc) + edgedict = Dict(edgelist[1] => [1]) + for (i, e) in enumerate(@view edgelist[2:end]) + if edgelist[i] != edgelist[i+1] + edgedict[e] = [i+1] + else + push!(edgedict[e], i+1) + end + end + + for (edge, inds) in edgedict + if haskey(edgedict, Edge(dst(edge), src(edge))) + distances[inds[1]] != 0. && continue + inds_ = edgedict[Edge(dst(edge), src(edge))] + + len = length(inds) + length(inds_) + sp = -inc/2*(len-1) + ep = sp + inc*(len-1) + dists = collect(sp:inc:ep) + distances[inds] = dists[1:length(inds)] + distances[inds_] = -dists[length(inds)+1:end] + else + sp = -inc/2*(length(inds)-1) + ep = sp + inc*(length(inds)-1) + distances[inds] = collect(sp:inc:ep) + end end distances end @@ -77,7 +112,7 @@ Notes: red and black arrows from `A` to the reaction node. """ function Catalyst.plot_network(rn::ReactionSystem) - srg = SRGraphWrap(rn) + srg = MultiGraphWrap(rn) ns = length(species(rn)) nodecolors = vcat([:skyblue3 for i in 1:ns], [:green for i in ns+1:nv(srg)]) @@ -104,16 +139,16 @@ function Catalyst.plot_network(rn::ReactionSystem) 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) + 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 @@ -130,27 +165,46 @@ end depends on species. i.e. `k*C, A --> B` for `C` a species. """ function Catalyst.plot_complexes(rn::ReactionSystem) - img = incidencematgraph(rn) + img = incidencematgraph(rn); D = incidencemat(rn; sparse=true) specs = species(rn); rxs = reactions(rn) - edgecolors = [:black for i in 1:ne(img)] - nodelabels = complexlabels(rn) + edgecolors = [:black for i in 1:length(rxs)] edgelabels = [repr(rx.rate) for rx in rxs] deps = Set() + edgelist = Vector{Graphs.SimpleEdge{Int}}() + rows = rowvals(D) + vals = nonzeros(D) + + # Construct the edge order for reactions. for (i, rx) in enumerate(rxs) + inds = nzrange(D, i) + val = vals[inds] + row = rows[inds] + (sub, prod) = val[1] == -1 ? (row[1], row[2]) : (row[2], row[1]) + push!(edgelist, Graphs.SimpleEdge(sub, prod)) + 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 + # Resolve differences between reaction order and edge order in the incidencematgraph. + rxorder = sortperm(edgelist); edgelist = edgelist[rxorder] + multiedges = Vector{Graphs.SimpleEdge{Int}}() + for i in 2:length(edgelist) + isequal(edgelist[i], edgelist[i-1]) && push!(multiedges, edgelist[i]) + end + img_ = MultiGraphWrap(img, multiedges) + + graphplot(img_; + edge_color = edgecolors[rxorder], + elabels = edgelabels[rxorder], + ilabels = complexlabels(rn), + node_color = :skyblue3, + elabels_rotation = 0, + arrow_shift = :end, + curve_distance_usage = true, + curve_distance = gen_distances(img_) ) end diff --git a/test/extensions/graphmakie.jl b/test/extensions/graphmakie.jl index ecba9d3f7e..a7fe1d1487 100644 --- a/test/extensions/graphmakie.jl +++ b/test/extensions/graphmakie.jl @@ -1,4 +1,4 @@ -using Catalyst, GraphMakie, CairoMakie, Graphs +using Catalyst, GraphMakie, CairoMakie, Graphs, SparseArrays include("../test_networks.jl") # Test that speciesreactiongraph is generated correctly @@ -74,11 +74,11 @@ let k * C, A --> C k * B, B --> C end - srg = CGME.SRGraphWrap(rn) + srg = CGME.MultiGraphWrap(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 + @test Graphs.Edge(3, s+2) ∈ srg.multiedges + @test Graphs.Edge(2, s+3) ∈ srg.multiedges # Since B is both a dep and a reactant @test count(==(Graphs.Edge(2, s+3)), edges(srg)) == 2 @@ -96,10 +96,88 @@ let k * A, A --> C k * B, B --> C end - srg = CGME.SRGraphWrap(rn) + srg = CGME.MultiGraphWrap(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 + +function test_edgeorder(rn) + # The initial edgelabels in `plot_complexes` is given by the order of reactions in reactions(rn). + D = incidencemat(rn; sparse=true); img = incidencematgraph(rn) + rxs = reactions(rn) + edgelist = Vector{Graphs.SimpleEdge{Int}}() + rows = rowvals(D) + vals = nonzeros(D) + + for (i, rx) in enumerate(rxs) + inds = nzrange(D, i) + val = vals[inds] + row = rows[inds] + (sub, prod) = val[1] == -1 ? (row[1], row[2]) : (row[2], row[1]) + push!(edgelist, Graphs.SimpleEdge(sub, prod)) + end + + rxorder = sortperm(edgelist); sortededgelist = edgelist[rxorder] + multiedges = Vector{Graphs.SimpleEdge{Int}}() + for i in 2:length(sortededgelist) + isequal(sortededgelist[i], sortededgelist[i-1]) && push!(multiedges, sortededgelist[i]) + end + img_ = CGME.MultiGraphWrap(img, multiedges) + + # Label iteration order is given by edgelist[rxorder. Actual edge drawing iteration order is given by edges(g) + @test edgelist[rxorder] == Graphs.edges(img_) + return rxorder +end + +# Test edge order for complexes. +let + # Multiple edges + rn = @reaction_network begin + k1, A --> B + (k2, k3), C <--> D + k4, A --> B + end + rxorder = test_edgeorder(rn) + edgelabels = [repr(rx.rate) for rx in reactions(rn)] + # Test internal order of labels is preserved + @test edgelabels[rxorder][1] == "k1" + @test edgelabels[rxorder][2] == "k4" + + # Multiple edges with species dependencies + rn = @reaction_network begin + k1, A --> B + (k2, k3), C <--> D + k4, A --> B + hillr(D, α, K, n), C --> D + k5*B, A --> B + end + rxorder = test_edgeorder(rn) + edgelabels = [repr(rx.rate) for rx in reactions(rn)] + labels = ["k1", "k4", "k5*B(t)", "k2", "Catalyst.hillr(D(t), α, K, n)", "k3"] + @test edgelabels[rxorder] == labels + + rs = @reaction_network begin + ka, Depot --> Central + (k12, k21), Central <--> Peripheral + ke, Central --> 0 + end + test_edgeorder(rs) + + rn = @reaction_network begin + (k1, k2), A <--> B + k3, C --> B + (α, β), (A, B) --> C + k4, B --> A + (k5, k6), B <--> A + k7, B --> C + (k8, k9), C <--> A + (k10, k11), (A, C) --> B + (k12, k13), (C, B) --> A + end + rxorder = test_edgeorder(rn) + edgelabels = [repr(rx.rate) for rx in reactions(rn)] + @test edgelabels[rxorder][1:3] == ["k1", "k6", "k10"] +end