Making GraphComplexNetwork work

[yewalenikhil65]

Catalyst.jl/src/graphs.jl

Line 286 in 8aa0a0a

function Graph(rn::ReactionSystem)

• A little observation while I am trying this, according to

  Catalyst.jl/src/graphs.jl

  Line 289 in 8aa0a0a

  statenodes = [Node(string(s.f.name), Attributes(:shape => "circle", :color => "#6C9AC3")) for s in specs]

  
  my code

    rxs = reactions(rn);
    specs = species(rn);
    complexes,B = reactioncomplexes(rn);
    fun = rcel -> specs[rcel.speciesid]*rcel.speciesstoich
    compfun(rc) = rc == Catalyst.ReactionComplex{Int64}[] ? 0 : sum(fun, rc)  
    compnodes = [Node(string(compfun(rc)),Attributes(:shape => "circle",
                      :color => "#6C9AC3")) for rc in complexes]

Our complexes created do not have .f.name fields as in elements of specs which seems to be required in edgify function

Catalyst.jl/src/graphs.jl

Line 246 in 8aa0a0a

val = String(p[1].f.name)

,So i have modified edgify function as follows

function myedgify(δ, i, reverse::Bool)
    attr = Attributes()
    return map(δ) do p
        val = string(p[1])   # modified here,. instead of String(p[1].f.name)
      weight = "$(p[2])"
      attr = Attributes(:label => weight, :labelfontsize => "6")
      return Edge(reverse ? ["rx_$i", "$val"] :
                            ["$val", "rx_$i"], attr)
    end
end

• Do we need transnodes ? I think not.

    transnodes = [Node(string("rx_$i"), Attributes(:shape => "point", :color => "#E28F41", :width => ".1")) for (i, r) in enumerate(rxs)]

• Rough way to create edges for connecting complexes as per reactions

  comps = [compfun(rc) for rc in complexes]
   edges = map(enumerate(rxs)) do (i, r)
      subcomp = comps[argmin(B[:,i])]  # returns complex in substrate participating in i'th reaction
      prodcomp = comps[argmax(B[:,i])] # returns complex in product participating in i'th reaction

#  we don't need the weights of edges really, but for now i have assigned -1 as weight from substrate complex to transnode and  1 as weight from transnode to product complex
      vcat(myedgify(zip([subcomp],[-1]), i, false),    
           myedgify(zip([prodcomp], [1]), i, true))  
    end

We don't need the weights of edges really, but for now i have assigned -1 as weight from substrate complex to transnode and 1 as weight from transnode to product complex (as entries in B matrix)

This doesn't work yet due to some issue . Following is full function with all changes explained above.

function myedgify(δ, i, reverse::Bool)
    attr = Attributes()
    return map(δ) do p
        val = string(p[1])
      weight = "$(p[2])"
      attr = Attributes(:label => weight, :labelfontsize => "6")
      return Edge(reverse ? ["rx_$i", "$val"] :
                            ["$val", "rx_$i"], attr)
    end
end

function GraphComplexNetwork(rn::ReactionSystem)
    rxs   = reactions(rn);
    specs = species(rn);
    complexes,B = reactioncomplexes(rn);
    fun = rcel -> specs[rcel.speciesid]*rcel.speciesstoich
    compfun(rc) = rc == Catalyst.ReactionComplex{Int64}[] ? 0 : sum(fun, rc)
    compnodes = [Node(string(compfun(rc)),Attributes(:shape => "circle",
                      :color => "#6C9AC3")) for rc in complexes]

    transnodes = [Node(string("rx_$i"), Attributes(:shape => "point", :color => "#E28F41", :width => ".1")) for (i, r) in enumerate(rxs)]    # may not need this, but kept fr 

    stmts = vcat(compnodes,transnodes)
    comps = [compfun(rc) for rc in complexes]
    edges = map(enumerate(rxs)) do (i, r)
      subcomp = comps[argmin(B[:,i])]
      prodcomp = comps[argmax(B[:,i])]
      vcat(myedgify(zip([subcomp],[-1]), i, false),
           myedgify(zip([prodcomp], [1]), i, true))
    end
  #  es = myedgifyrates(rxs, Xₜₑₘₚ)   # temporarily not implementing edgifyrates
  #  (!isempty(es)) && push!(edges, es)

    stmts2 = Vector{Statement}()
    append!(stmts2, stmts)
    append!(stmts2, collect(flatten(edges)))
    g = Digraph("G", stmts2; graph_attrs=graph_attrs, node_attrs=node_attrs,
                                                  edge_attrs=edge_attrs)
    return g
end

Let's take an example Example

rn = @reaction_network begin
    k₁, A --> B
    k₂, B + B --> C + B
    k₃, B + C --> A + C
end k₁ k₂ k₃

using this example for running line-by-line statements in GraphComplexNetwork results in stmts2 as

julia> stmts2
14-element Vector{Statement}:
 Node("A(t)", OrderedDict{Symbol, AttributeValue}(:shape => "circle", :color => "#6C9AC3"))
 Node("B(t)", OrderedDict{Symbol, AttributeValue}(:shape => "circle", :color => "#6C9AC3"))
 Node("2B(t)", OrderedDict{Symbol, AttributeValue}(:shape => "circle", :color => "#6C9AC3"))
 Node("B(t) + C(t)", OrderedDict{Symbol, AttributeValue}(:shape => "circle", :color => "#6C9AC3"))
 Node("A(t) + C(t)", OrderedDict{Symbol, AttributeValue}(:shape => "circle", :color => "#6C9AC3"))
 Node("rx_1", OrderedDict{Symbol, AttributeValue}(:shape => "point", :color => "#E28F41", :width => ".1"))
 Node("rx_2", OrderedDict{Symbol, AttributeValue}(:shape => "point", :color => "#E28F41", :width => ".1"))
 Node("rx_3", OrderedDict{Symbol, AttributeValue}(:shape => "point", :color => "#E28F41", :width => ".1"))
 Edge(NodeID[NodeID("A(t)", "", ""), NodeID("rx_1", "", "")], OrderedDict{Symbol, AttributeValue}(:label => "-1", :labelfontsize => "6"))
 Edge(NodeID[NodeID("rx_1", "", ""), NodeID("B(t)", "", "")], OrderedDict{Symbol, AttributeValue}(:label => "1", :labelfontsize => "6"))
 Edge(NodeID[NodeID("2B(t)", "", ""), NodeID("rx_2", "", "")], OrderedDict{Symbol, AttributeValue}(:label => "-1", :labelfontsize => "6"))
 Edge(NodeID[NodeID("rx_2", "", ""), NodeID("B(t) + C(t)", "", "")], OrderedDict{Symbol, AttributeValue}(:label => "1", :labelfontsize => "6"))
 Edge(NodeID[NodeID("B(t) + C(t)", "", ""), NodeID("rx_3", "", "")], OrderedDict{Symbol, AttributeValue}(:label => "-1", :labelfontsize => "6"))
 Edge(NodeID[NodeID("rx_3", "", ""), NodeID("A(t) + C(t)", "", "")], OrderedDict{Symbol, AttributeValue}(:label => "1", :labelfontsize => "6"))

julia> GraphComplexNetwork(rn)
Error: <stdin>: syntax error in line 5 near '('
Error showing value of type Graph:
ERROR: failed process: Process(`dot -Tsvg`, ProcessExited(1)) [1]

Stacktrace:
  [1] pipeline_error
    @ ./process.jl:525 [inlined]
  [2] open(f::var"#43#46"{Graph}, cmds::Cmd, args::IOContext{IOBuffer}; kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:write,), Tuple{Bool}}})
    @ Base ./process.jl:400
  [3] run_graphviz(io::IOContext{IOBuffer}, graph::Graph; prog::Nothing, format::String)
    @ Main ~/nikhil/s rao reduction of chem reaction/model_reduction_julia/myGraphs.jl:121
  [4] show(io::IOContext{IOBuffer}, #unused#::MIME{Symbol("image/svg+xml")}, graph::Graph)
    @ Main ~/nikhil/s rao reduction of chem reaction/model_reduction_julia/myGraphs.jl:133
  [5] show(io::IOContext{IOBuffer}, m::String, x::Graph)
    @ Base.Multimedia ./multimedia.jl:111
  [6] displayinplotpane(x::Graph)
    @ Atom ~/.julia/packages/Atom/4PVOT/src/display/showdisplay.jl:67
  [7] display(d::Atom.JunoDisplay, x::Graph)
    @ Atom ~/.julia/packages/Atom/4PVOT/src/display/showdisplay.jl:118
  [8] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:328
  [9] #invokelatest#2
    @ ./essentials.jl:708 [inlined]
 [10] invokelatest
    @ ./essentials.jl:706 [inlined]
 [11] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:247
 [12] (::REPL.var"#40#41"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:231
 [13] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:462
 [14] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:229
 [15] (::REPL.var"#do_respond#61"{Bool, Bool, Atom.var"#246#247"{Module}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:798
 [16] #invokelatest#2
    @ ./essentials.jl:708 [inlined]
 [17] invokelatest
    @ ./essentials.jl:706 [inlined]
 [18] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/LineEdit.jl:2441
 [19] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:1126
 [20] (::REPL.var"#44#49"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ./task.jl:411

[yewalenikhil65]

It seems, Node name as a string should not contain characters like ( or + etc. Need to think of alternative representation of complexes ,
like A(t)+B(t) complex to be represented as A₊B or something similar

[yewalenikhil65]

Answered and rectified at SciML#380