A simple interface for general spatial reaction networks

[kaandocal]

Motivation

Many modellers work with spatial reaction networks, e.g. ones defined on a graph or on a spatial grid, which can be seen as a special type of graph. These can be represented straightforwardly by "normal" reaction networks where each species is split up into multiple sub-species corresponding to each node of the graph, but doing this by hand is very tedious an error-prone. This issue has been brought up previously in #229 and I have attempted to augment Catalyst.jl by extending the reaction network language.

Implementation

In order not to break anything I implemented this as an experimental macro whose output can be sent to @reaction_network; this should prevent any code-breaking changes and retain full backwards compatibility.

Species are indexed by vertices in arbitrary graphs, and multiple graphs for one system are supported. Since I am a relative newcomer to Julia I am not sure if there are any standard tools for dealing with graphs, so I've hacked my own for demonstration purposes.

Each species can have an index into a graph, denoted using indexing: X becomes X[i], Y becomes Y[j], etc. Every index is declared to belong to one graph as macro parameters of the form G[i] H[j,k], which says that i is a vertex of G and j,k are vertices of H. The species X[i] is then expanded to a set of species X_i for every vertex i of G. Constants in reaction rates can be extended similarly, but they can depend on multiple indices, e.g. if a is a constant then a[i,j] mapping to a_(i,j) can be defined for each i in G, j in H. Every line using indexing is thus expanded into multiple lines for every valid combination of indices, e.g.

d, X[i] -> 0
c[i], 0 -> X[i]
D[j,k], Y[j] -> Y[k]
a, X[i] + Y[j] + Y[k] -> X[i]

become

d, X_i -> 0 (for each i in G)
c_i, 0 -> X_i (for each i in G)
D_(j,k), Y_j -> Y_k (for each edge j<->k in H)
a, X_i + Y_j + Y_k -> Y_j (for each i in G, edge j<->k in H)

Example:

# Cyclic graph on 4 elements
G = BasicGraph([ 1, 2, 3, 4 ], [ 1=>2, 2=>3, 3=>4, 4=>1 ]) 

@spatial_reaction_network begin
  sigma[i] --> X[i]
  D, X[i] --> X[j]
  delta, 2X[i]-->0
  a, 0 -> Z
  b, Z + X[i] -> 2X[i]
end G[i,j]

returns

quote                                               
    (sigma, 0 → X_1)                                
    (sigma, 0 → X_2)                                                                                     
    (sigma, 0 → X_3)                                
    (sigma, 0 → X_4)                                                                                     
    (D, X_2 → X_1)                                  
    (D, X_3 → X_2)                                  
    (D, X_4 → X_3)                                  
    (D, X_1 → X_4)                                  
    (D, X_1 → X_2)                                  
    (D, X_2 → X_3)                                                                                       
    (D, X_3 → X_4)                                                                                       
    (D, X_4 → X_1)                                  
    (delta, 2X_1 → 0)                                                                                    
    (delta, 2X_2 → 0)                               
    (delta, 2X_3 → 0)                               
    (delta, 2X_4 → 0)                               
    (a, 0 → Z)                                      
    (b, Z + X_1 → 2X_1)                             
    (b, Z + X_2 → 2X_2)                             
    (b, Z + X_3 → 2X_3)                             
    (b, Z + X_4 → 2X_4)                             
end

I haven't spent much time with either Julia or SciML, so I would greatly appreciate any feedback!

Remarks:

• If n indices belonging to a graph appear in a reaction, the reaction is expanded for every n-edge of the graph. If n==1 this corresponds to every node of the graph, if n==2 we get every edge and for n>=3 you can define your own higher-order edges to model more complex interactions.
• If two different graphs appear in a reaction we expand jointly over all valid permutations of indices of both graphs

Todos:

• UGLY: the macro needs the values of the graphs to generate its output, and for this I used Base.MainInclude.eval(). What's the proper way to do this?
• The generated variables have names of the form var"X_$i", where i ranges over the vertices of the graph. There is probably a better way to encode this information...
• Prevent one species from being indexed by multiple graphs to avoid name-conflicts?
• Do some testing
• Work on graph interface
• How best to combine this with arbitrary parameters?

[ChrisRackauckas]

We should really just build this on LightGraphs.jl instead of building out a whole new graph format.

[isaacsas]

@kaandocal looks cool. See also SciML/JumpProcesses.jl#131 where we started on a direct interface to the Gillespie methods in a similar way (giving a graph and flattening it to a giant set of reactions).

I'll take a look tonight or tomorrow and give you some feedback (sorry, need to finish up lectures now!).

[isaacsas]

@kaandocal haven't forgotten about this or your other PR; was just a busy week. Blocking out time tomorrow afternoon to go through them carefully.

[kaandocal]

No need to hurry, while I appreciate any feedback this is not an urgent issue!

[Vilin97]

Hi. This is currently being built. Check out https://vilin97.github.io/posts/post4/ where I describe how to set up a spatial jump problem. Or check out any of the spatial tests in DiffEqJump.

[yewalenikhil65]

Hi. This is currently being built. Check out https://vilin97.github.io/posts/post4/ where I describe how to set up a spatial jump problem. Or check out any of the spatial tests in DiffEqJump.

@Vilin97 i found this blog post really useful. I was curious whether we can define some say vector field(velocity) over the cartesian grid, and check out it's influence on the Jump process (where, say the rate laws is dependent on vector field)

[Vilin97]

I was curious whether we can define some say vector field(velocity) over the cartesian grid, and check out it's influence on the Jump process (where, say the rate laws is dependent on vector field)

You can define your own hopping rates. For example, your hopping rate to the nodes on the right might be 1, while the hopping rates to the nodes on the left might be 0.1. Thus your molecules will tend to drift right. Check out this test: /~https://github.com/SciML/DiffEqJump.jl/blob/master/test/spatial/diffusion.jl#L84. It sets up a similar thing: only hopping up and left is allowed.

[yewalenikhil65]

I was curious whether we can define some say vector field(velocity) over the cartesian grid, and check out it's influence on the Jump process (where, say the rate laws is dependent on vector field)

You can define your own hopping rates. For example, your hopping rate to the nodes on the right might be 1, while the hopping rates to the nodes on the left might be 0.1. Thus your molecules will tend to drift right. Check out this test: /~https://github.com/SciML/DiffEqJump.jl/blob/master/test/spatial/diffusion.jl#L84. It sets up a similar thing: only hopping up and left is allowed.

@Vilin97
This is awesome work! Hope advection of molecules effects are also incorporated in near future

[isaacsas]

Advection could be simulated right now. The user needs to specify the hopping rates currently, so they would need to calculate what the advective transition rates for their model are. Later we can hopefully move model specification to higher levels with something like this PR combined with automatic discretization of transport operators to transition rates. That's a bit off though.