Lattice reaction network

[TorkelE]

Creates a simple interface that distributes a reaction system across a connected graph, and allows reactions between connected nodes.

Currently added MetaGraphs to dependencies. The plotting interface also required Plots and GraphRecipes.

Definitely things in the code that could be made smoother. Also plan to bring stuff more in line with the rest of MTK, e.g. make LatticeReactionSystem <: MT.AbstractTimeDependentSystem, and sort out SpatialReaction in some nice way. But I figured it would be good to have something up first.

Example of usage:

using Catalyst
using OrdinaryDiffEq
using Graphs, MetaGraphs
using Plots

# Make the base network (a brusselator).
rs = @reaction_network base_network begin
    A, ∅ → X
    1, 2X + Y → 3X
    B, X → Y
    1, X → ∅
end A B;

# Declare a set of spatial reactions.
# Rate followed by substrates (compartment1,comaprtment2) and products (compartment1,comaprtment2).
# Stoichiomtery currently given using species=>stoichiometry, but can use somethign better.
srs = [SpatialReaction((@parameters D)[1],([:X=>1],[]),([],[:X=>1]))]

# Creates a graph on which to simualte network (now just a grid)
lattice = MetaGraph(Graphs.SimpleGraphs.grid([20,20]))
# Meta graoh allows setting properties in each node. here :X denotes initial concentration of :X.
# If no :X is set here, whatever X value is used in ODEProblem input will be used.
foreach(v -> (set_prop!(lattice,v,:X,10*rand());set_prop!(lattice,v,:Y,10*rand())), vertices(lattice))

# Bundles everything into a LatticeReactionSystem.
lrs = LatticeReactionSystem(rs,srs,lattice)

# Creates an ODEProblem
u0 =[] #Could set X and Y value, but since all nodes have one, not needed,
tspan = (0.0,100.0)
p = [:D => 0.1, :A=>1.0, :B=>4.0]
oprob = ODEProblem(lrs,u0,tspan,p)

# Simulates the problem.
@time sol = solve(oprob,QNDF(),saveat=0.1)

# The "plot_spatial_sol" function allows us to make an animation.
@time anim = plot_spatial_sol(sol,lattice,2;samp_freq=10,nodesize=0.2)
@time gif(anim, "anim.gif", fps = 2)

[isaacsas]

I was traveling today, but will take a look tomorrow. Looks cool!

[TorkelE]

Stochastic systems work as well (haven't tried SDE yet though):

dprob = DiscreteProblem(lrs,u0,tspan,p)
@time jprob = JumpProblem(lrs,dprob,RSSACR(),save_positions=(false,false))
@time sol = solve(jprob,SSAStepper(),saveat=0.1)

@time anim = plot_spatial_sol(sol,lattice,2;samp_freq=10,nodesize=0.2)
@time gif(anim, "anim.gif", fps = 10)

[isaacsas]

So one thing we need to be consider with having new system types is that it means we'll have to reimplement any compositional modeling type functionality (subsystems, flatten, extend, compose, all the various querying functions like (states, species, reactions, etc)). They'd have to be generalized for spatial systems anyways, but it is something for us to keep in mind.

[TorkelE]

Yes. I should be in Boston by the end of the month. Hopefully we could meet up some time to discuss through what's the best way of setting everything up.

[TorkelE]

I wanted to have rates dependant on both edge and contabased parameters. The former should be quite easy (will add), but it was not obvious what would be the best way to write a general rate with a mix of parameters from two adjacent cells.

[isaacsas]

Yeah, specifying rate expressions is tricky. One really wants a stencil type language to specify them, which I think is available in Symbolics but not really been advertised. In JumpProcesses we just allowed a few factored forms, like D[s] * L[i,j] where s is a species index and L[i,j] could be a discrete Laplacian, or more generally L[s,i,j].

[isaacsas]

(above i and j denote destination and source vertices)

[TorkelE]

Update:

• MetaGraph Dependency removed.
• Initial conditions/parameters for individual compartments/connections are now designated by passing vectors in u0/p.
• Parameter can now be specified for individual connections. Parameters specified in the initial reaction system are expanded across all vertices, and parameters for the spatial reactions are expanded across all edges.
• Now takes a DiGraph as input, permitting different connections in both directions between cells (this also simplifies the internal generation of the network). Graphs that are inputted are converted to digraphs (additional support for this, e.g. when using graph having the same parameter in both directions) needs to be implemented.
• The expanded reaction network is now never saved, but only generated when an O/S/Discrete/JumpProblem is created.
• Can now designate xs and ys vectors with the position of all compartments in the plot.

New workflow:

# Fetch packages
using Catalyst
using OrdinaryDiffEq
using Graphs
using Plots

# Declare the base network (a brusselator), the lattice, and the spatial reactions..
rs = @reaction_network base_network begin
    A, ∅ → X
    1, 2X + Y → 3X
    B, X → Y
    1, X → ∅
end A B;
srs = [SpatialReaction((@parameters D)[1],([:X=>1],[]),([],[:X=>1]))]
lattice = Graphs.SimpleGraphs.grid([20,20])

# Bundle into a LatticeReactionSystem.
lrs = LatticeReactionSystem(rs,srs,lattice)

# Set initial conditions and parameter values.
u0 = [:X => 10*rand(length(vertices(lattice))), :Y => 10*rand(length(vertices(lattice)))]
p = [:A => 1.0, :B => 4.0, :D => 0.2]

# Simulate system.
@time oprob = ODEProblem(lrs,u0,(0.0,100.0),p)
@time sol = solve(oprob,QNDF();saveat=0.1)

# Plot animation.
@time anim = plot_spatial_sol(sol,lattice,2;samp_freq=4,nodesize=0.2)
@time gif(anim, "anim.gif", fps=15)

Currently, the 4 runtimes are

5.148575 seconds (22.17 M allocations: 1.436 GiB, 5.01% gc time)
3.662192 seconds (2.31 k allocations: 24.032 MiB)
68.967249 seconds (1.17 G allocations: 29.839 GiB, 8.03% gc time, 0.19% compilation time)
2.301520 seconds (507 allocations: 34.781 KiB)

with the most intense part being the actual plotting. I think that maybe Makie might be better, since it allows updating the graph, and not replotting each frame in full. I will look into this.

[TorkelE]

Given that we have https://github.com/SciML/Catalyst.jl/pull/571, I'm closing this.