Support for Spatial Models

[TorkelE]

I still have some old stuff to finish up in Cambridge, but the next step, which I think is relatively low hanging, would be to enable good support for spatial models.

Two main categories would be:

Discrete Space Models E.g. when modelling plants people often consider it a set of cells, each with a non-spatial CRN, and then add transportation/diffusion between neighbouring cells.

I imagine the best way to do this would be for the user to provide

1. A ReactionSystem (like those we already have).
2. A graph (each node would contain a CRN, (possibly weighted) edges would define neighbouring containers.
3. Some diffusion/transportation rules. These would take the state of two neighbouring nodes, the potential weight of their connecting edge, and define the rate of transportation between the two. From this, some big ODE could be constructed. I should be able to make the interface without problem, although you guys might have some opinions on how to make it actually performant.

Continuous Space Models E.g. typical Reaction/diffusion systems.

Again similar things would be provided, but not a graph but some domain. There might be some sensible way to define diffusion constants and possibly transportation rules for each species. The result would be a PDESystem. I think for this one I might want to discuss it through with you guys properly first.

[isaacsas]

This is exactly what I've been thinking about looking at now that most of the non-spatial work is finished, so that sounds great! (And spatial tooling / solvers should also be the focus of my supported efforts this Fall.)

For discrete spatial models, we could target a few grid types (similar to the spatial SSAs in JumpProcesses). These include a general graph, but also structured Cartesian grids. What is missing in JumpProcesses, and I'll work on it, is allowing a variety of boundary conditions, allowing non-mass action reactions, and support for non-local reactions (i.e. reactions with neighbors or other sites). The latter is a ways off though, so V1 could just target same-site reactions.

For PDEs see http://methodoflines.sciml.ai/dev/tutorials/brusselator/. A good first step would probably be seeing if we can get the DSL and ReactionSystem objects to work if we make the variables/parameters spatial functions too, and can automatically generate a function like brusselator_f there that just works with symbolically specified derivatives. It would be good if whatever approach we take remains compatible with specifying the symbolic, non-reaction, components as in that tutorial.

At the Catalyst modeling level, the main considerations seem to be specifying transport laws and functions (i.e. diffusivities and/or potential/advective terms), specifying the compartments and their geometries, and specifying boundary conditions on each compartment. I guess the first version could assume a single compartment / geometry, with one set of boundary conditions.

[TorkelE]

For discrete spatial models, would boundary conditions need to be defined? I guess one could e.g. define a series of containers with diffusion, with the two ends connected to the sea or something similar. I often see this in continuous PDE case, but is this something that occurs for discrete space?

Also, under the hood, would a general graph and a structured Cartesian grid be different? Having a function which conveniently generates a graph of e.g. Cartesian grid form would make sense in either case.

Once internal structures are all good, we could try and make a DSL for spatial stuff (when combined with a reaction system and a space generates a spatial problem).

[yewalenikhil65]

Out of curiosity, is it possible to club agent-based modeling approach with spatial chemical reactions ? I guess Agents.jl can facilitate this if such problems are relevant in comp-bio community ..

[isaacsas]

I think PDESystems -> MethodOfLines.jl are in good enough shape that we could tackle the reaction-diffusion PDE case today. It is really just making sure everything works when the species and rate constants are functions of space and time instead of just time (which might be a bit of work since I think lots of Catalyst implicitly assumes species are functions of time only and might need updating). Otherwise our current code should be able to generate the needed reaction terms for the PDE, and we can combine that with the derivative operators one can already define symbolically in ModelingToolkit. This would give us the needed symbolic backend.

As a second step we could think about adding a more DSL-based approach that let's users specify compartments, BCs, and transport information and then generates the symbolic derivatives from that, and updates the current DSL to support spatial species/parameters.

There isn't an analogous spatial jump process system in ModelingToolkit that we could target yet, but I'm hoping to get to that soon.

Cartesian grids, being structured, require much less memory to store than a large graph, and at least for spatial SSAs can lead to optimizations in not needing to store all the possible transition rates in many cases (for example, if the diffusivity is constant in space).

Discrete models can have analogous boundary conditions to continuous models, (absorbing, reflecting, partially reflecting, periodic, etc). For a general graph this is just encoded by the connectivity and the transition rates one specifies: reflecting are just leaf nodes in the graph, absorbing and partially absorbing can use leaf nodes you can't transition out of or a single global sink node (with possibly modified transition rates for partial absorbtion), periodic is just a choice of topology. But for Cartesian grids one might need to special case handling them since one generally doesn't store all transition rates.

[isaacsas]

@yewalenikhil65 Yes, one could create a particle-tracking version of spatial-stochastic models, i.e. a Brownian Dynamics model. That could potentially be built on top of Agents.jl (though there might be more optimized approaches that would make sense, but would need dedicated Julia solvers).

[isaacsas]

In terms of discrete, probably we should just target unstructured grids (i.e. graphs) first. Cartesian grids are really just a special case and could be added subsequently as an optimization.

[TorkelE]

(Simple) spatial models are currently possible. I have also created an issue for gathering future steps for improving spatial implementations at: https://github.com/SciML/Catalyst.jl/issues/726. Would it make sense to close this and gather future spatial discussions there?

[TorkelE]

Stuff from here is now intergrated to https://github.com/SciML/Catalyst.jl/issues/726