Efficient simulation of spatially coupled ODEs on arbitrary graphs

I know this question is not very specific, but I would like to check whether there are existing approaches for simulating coupled ODEs on nodes of (arbitrarily large) graphs/networks with diffusion between the nodes. Looking into the documentation of the main packages related to solving such problems (Catalyst.jl, JumpProcesses.jl and DifferentialEquations.jl), I found a few prospects for solving such systems.

The main prospect that has explicit documentation on spatial processes on graphs using Graphs.jl is JumpProcesses: Spatial SSAs with JumpProcesses.jl · JumpProcesses.jl. I have no reason to believe that the implementation is well-rounded and efficient and would for sure suit much of my needs. However, it remains a stochastic simulation algorithm and does not provide deterministic solutions to the system of equations, correct?

Further, in the documentation of Catalyst.jl I noticed that one can introduce spatial variables using the @ivs macro (as seen here), but the documentation does not mention the introduction spatial variables as nodes on a graph. Of course I could just define a flattened version using the ReactionNetwork function, e.g. something like:

function generate_rs(S::Int, G::SimpleGraph)
    M = size(G)[begin]

    @variables t
    @parameters r[1:S], hop_rate[1:M,1:M]
    @species (x(t))[1:S,1:M]

    adj = Graphs.LinAlg.adjacency_matrix(G)

    growth_rxs = []
    hop_rxs = []
    for m in 1:M
        for n in 1:M, i in 1:S
            # Hops between nodes
            if adj[m,n] != 0
                _hop_rxn = Reaction(hop_rate[m,n], [x[i,m]], [x[i,n]], [1], [1])
                push!(hop_rxs, _hop_rxn)
        for i in 1:S
            # Reaction on node
            _growth_rxn = Reaction(r[i], [x[i,m]], [x[i,m]], [1], [2])
            push!(growth_rxs, _growth_rxn)
    return @named rs = ReactionSystem(vcat(growth_rxs, hop_rxs))

Which appears to create the equations that I am after. However some piece of text in the documentation of JumpProcesses.jl warns me about this:

Additionally, all standard solvers are supported as well, although they are expected to use more memory and be slower. They “flatten” the problem, i.e., turn all hops into reactions, resulting in a much larger system.

This can create issues for me down the line as my systems of interest are relatively large. Does creating an ODEProblem with such a structure result in fast and (memory) efficient solutions to the ODE, even when the number of species and the size of the graph greatly increase? Or is there perhaps another way one can define spatial ODEs on graphs/networks?

Thank you for any insight on this topic you can give me.

Just to clarify, do you want to simulate a system of ODEs on a graph or a system of jump processes on a graph? JumpProcesses is only for the latter, and should work if you are ok with only mass action reactions at each node of the graph. There is a Catalyst PR for generating more efficient representations of chemical reaction network ODE models on graphs, see

but I’m not sure its current state. Maybe @Gaussia can comment on how well it is working right now.

Yes, I have recently developed something that should pretty much exactly do this. You provide a reaction system, a graph, and a list of species that are difussing it, and then an ODEProblem is created.

I have been out for a bit, but am currently back. In principle the PR is ready, but need to do some benchmarking tests. Hopefully there’s a working version end of this week. Unsure when it will get merged though.

How big is your graph and ODE?

Initially I will most likely have only mass-action reactions, so JumpProcesses should suffice – for now. But as I am interested in general reaction-diffusion systems that will most likely do not only contain mass-action reactions. Nevertheless, I will try out spatial JumpProcesses and report back. I will also compare the Catalyst PR and provide feedback if applicable.

Good to hear that there is more interest in these types of systems. I have also encountered work on solvers for Markovian (epidemic) processes on graphs (see e.g. Cota, Ferreira (2017)) so I can also implement efficient Gillespie algorithms on graphs for the general reaction-diffusion equations I am interested in.

The graphs need not be very large, but a few hundred or a few thousand nodes is not out of the question. An example of the systems I am interested epidemics models with local dynamics on the nodes (e.g. an SIR model within a city) and dispersal between the nodes (movement of individuals between cities, for example). Real-world networks of such systems can grow quite rapidly (see e.g. the work mentioned above, which has \sim10^4 nodes. Other relevant works, such as Meena et al. (2022) have networks with \sim 10^3 nodes.

I am also looking into how other authors that study dynamical processes on top of networks numerically solve their systems.

If you check the latest update of Spatial Reaction Network Implementation by TorkelE · Pull Request #644 · SciML/Catalyst.jl · GitHub I have added simulations of the SIR model with diffusion to the tests (test/spatial_reaction_systems/lattice_reaction_system). It can simulate a 100x100 grid in 0.1 seconds and a 100x100x100 in 20 seconds. Currently we do not have any specific support for spatial Gillespie simulations (although that is very much in the pipeline).

Please ask if you have any further questions.