# Efficient simulation of spatially coupled ODEs on arbitrary graphs

**URL:** <https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation, catalyst\
**Created:** [June 26, 2023, 5:51pm UTC](https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852 "2023-06-26T17:51:20Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [June 26, 2023, 5:51pm UTC](https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852/1 "2023-06-26T17:51:20Z")

</div>

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](https://docs.sciml.ai/Catalyst/stable/), [JumpProcesses.jl](https://docs.sciml.ai/JumpProcesses/stable/) and [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/)), I found a few prospects for solving such systems.

The main prospect that has explicit documentation on spatial processes on graphs using [Graphs.jl](https://juliagraphs.org/Graphs.jl/dev/) is `JumpProcesses`: [Spatial SSAs with JumpProcesses.jl · JumpProcesses.jl](https://docs.sciml.ai/JumpProcesses/stable/tutorials/spatial/). 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](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/#Specifying-alternative-time-variables-and/or-extra-independent-variables)), 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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [June 26, 2023, 8:08pm UTC](https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852/2 "2023-06-26T20:08:43Z")

</div>

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

> <https://github.com/SciML/Catalyst.jl/pull/644>
>
> A revamped version of https://github.com/SciML/Catalyst.jl/pull/571. While we mi…ght be able to tune some stuff, this should work well. Jacobian can be built, so stiff systems can be implemented properly. 
> 
> Syntax is similar to previously, see this example:
> 
> When making a spatial simulation:
> \`u0\` is a NxM matrix, where N is the number of species in the (non-spatial) reaction network, and M is the number of compartments.
> \`p\` is a tupple \`(pV,pE)\`, 
> \`pV\` is a LxM matrix, where L is the number of parameters in the (non-spatial) reaction network, and M is the number of compartments.
> \`pE\` is a RxS matrix, where R is the number of parameters tied to spatial reactions, and S is the number of connections.
> 
> Sample code:
> \`\`\`julia
> using Catalyst, OrdinaryDiffEq, Graphs
> 
> rs = @reaction\_network begin
> A, ∅ → X
> 1, 2X + Y → 3X
> B, X → Y
> 1, X → ∅
> end
> srs = \[DiffusionReaction(:D, :X)\]
> lattice = grid(\[20, 20\])
> lrs = LatticeReactionSystem(rs, srs, lattice);
> 
> u0\_in = \[:X =\> 10 \* rand(nv(lattice)), :Y =\> 10 \* rand(nv(lattice))\]
> tspan = (0.0, 100.0)
> p\_in = \[:A =\> 1.0, :B =\> 4.0, :D =\> 0.2\]
> 
> @time oprob = ODEProblem(lrs, u0\_in, tspan, p\_in)
> @time sol = solve(oprob, Tsit5());
> \`\`\`
> 
> Plotting is not added to Catalyst right now, but running this code:
> \`\`\`julia
> using GLMakie
> using Makie.Colors
> using GraphMakie
> function animate\_spatial\_sol(sol, x\_pos, y\_pos; animation\_name="spatial\_animation.mp4", framerate=100, min\_val=nothing, max\_val=nothing, var=1, node\_size=2.0, node\_marker=:circle, edge\_width=0.0, timetitle=true)
> trajectories = map(vals -\> vals\[var, :\], sol.u)
> mini, maxi = extrema(vcat(trajectories...))
> !isnothing(min\_val) && (mini = min\_val)
> !isnothing(max\_val) && (maxi = min\_val)
> fixed\_layout(\_) = map(i -\> (x\_pos\[i\], y\_pos\[i\]), 1:length(trajectories\[1\]))
> 
> idx = Observable(1)
> colors = @lift map(val -\> RGB{Float64}(val, val, 0.0), (trajectories\[$idx\] .- mini) ./ maxi)
> title = timetitle ? @lift("t = $(sol.t\[($idx)\])") : ""
> fig = graphplot(lattice, layout=fixed\_layout, node\_size=node\_size, node\_marker=node\_marker, edge\_width=edge\_width, node\_color=colors,
> axis=(title=title,))
> 
> record(fig, animation\_name, 1:length(sol.t);
> framerate=framerate) do i
> idx\[\] = i
> end
> end
> \`\`\`
> enable us to create a plot through:
> \`\`\`julia
> x\_pos = vcat(fill(1:20, 20)...)
> y\_pos = vcat(map(i -\> fill(i, 20), 1:20)...)
> @time animate\_spatial\_sol(sol, x\_pos, y\_pos; node\_size=40.0, node\_marker=:rect)
> \`\`\`

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

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [June 27, 2023, 1:55am UTC](https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852/3 "2023-06-27T01:55:24Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [June 27, 2023, 8:56am UTC](https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852/4 "2023-06-27T08:56:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [June 27, 2023, 9:13am UTC](https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852/5 "2023-06-27T09:13:29Z")

</div>

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)](https://doi.org/10.1016/j.cpc.2017.06.007)) 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)](https://doi.org/10.1038/s41567-023-02020-8) 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.

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [June 29, 2023, 9:00pm UTC](https://discourse.julialang.org/t/efficient-simulation-of-spatially-coupled-odes-on-arbitrary-graphs/100852/6 "2023-06-29T21:00:37Z")

</div>

If you check the latest update of [Spatial Reaction Network Implementation by TorkelE · Pull Request #644 · SciML/Catalyst.jl · GitHub](https://github.com/SciML/Catalyst.jl/pull/644) 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.