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)
end
end
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)
end
end
return @named rs = ReactionSystem(vcat(growth_rxs, hop_rxs))
end
```

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.