GitHub

Introduction to Spatial Modelling with Catalyst

Catalyst supports the expansion of non-spatial ReactionSystems (created using e.g. the @reaction_network DSL) to spatial domains. Spatial simulation of Catalyst models is a work in progress. Currently, the following is supported:

  • Spatial ODE and Jump simulations.
  • Discrete spatial domains.
  • Constant-rate transportation reactions (species moving spatially at constant rates).

Features for which support is planned in future updates include:

  • Models on continuous domains with automatic discretisation (these models can already be simulated if the user provides a discretisation).
  • SDE simulations.
  • Transport reactions with non-constant rates as well as more general spatial reactions.

This tutorial introduces spatial modelling on discrete domains, here called lattices. It describes the basics of creating and simulating such models. To do so, it uses ODE simulations as examples. Additional tutorials provide further details on how to interact with spatial simulation structures and plot spatial simulations, and also provide further details on ODE and jump simulations, respectively.

Basic example of a spatial simulation on a discrete domain

To perform discrete-space spatial simulations, the user must first define a LatticeReactionSystem. These combine:

  • A (non-spatial) ReactionSystem(@ref) model (created using standard Catalyst syntax).
  • A vector of spatial reactions, describing how species can move spatially across the domain.
  • A lattice defining the spatial domain's compartments and how they are connected.

Here, as an example, we will simulate a spatial two-state model. To do so, we first define our (non-spatial) model, the spatial reactions, and the lattice. These are then bundled into a LatticeReactionSystem.

using Catalyst
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
diffusion_rx = @transport_reaction D X1
lattice = CartesianGrid((5,5))
lrs = LatticeReactionSystem(two_state_model, [diffusion_rx], lattice)

This model contains:

  • A single spatial reaction, a transport reaction where $X1$ moves at constant rate $D$ between adjacent compartments.
  • A 2d Cartesian grid of 5x5 compartments to simulate our model on.

More details on spatial reactions are available here. In addition to Cartesian grid lattices (in 1, 2, and 3 dimensions), masked and unstructured (graph) lattices are also supported. The different lattice types described in more detail here.

Once created, LatticeReactionSystems can be used as input to various problem types, which then can be simulated using the same syntax as non-spatial models. Here, we prepare an ODE simulation by creating an ODEProblem:

u0 = [:X1 => rand(5, 5), :X2 => 2.0]
tspan = (0.0, 10.0)
ps = [:k1 => 1.0, :k2 => 2.0, :D => 0.2]
oprob = ODEProblem(lrs, u0, tspan, ps)

In this example we used non-uniform values for $X1$'s initial condition, but uniform values for the remaining initial condition and parameter values. More details of uniform and non-uniform initial conditions and parameter values are provided here. We also note that the diffusion reaction introduces a new parameter, $D$ (determining $X1$'s diffusion rate), whose value must be designated in the parameter vector.

We can now simulate our model:

using OrdinaryDiffEqDefault
sol = solve(oprob)

We note that simulations of spatial models are often computationally expensive. Advice on the performance of spatial ODE simulations is provided here.

Finally, we can access "$X1$'s value across the simulation using

lat_getu(sol, :X1, lrs)
21-element Vector{Matrix{Float64}}:
 [0.5211411238735933 0.08265198170584165 … 0.12380555805699667 0.8475203901009595; 0.5481703303476589 0.001017352378020031 … 0.11195344563549503 0.8704358896787583; … ; 0.8709749218177816 0.8695528842481064 … 0.06913098768378156 0.8622860226771659; 0.30986705063176767 0.11133967669621114 … 0.09386267557193417 0.8216652566054319]
 [0.5324011320726877 0.09631116263457853 … 0.13753653479629427 0.8575062045594158; 0.5593962014790875 0.015048083380178176 … 0.12578312375417036 0.8798925651610332; … ; 0.8809236884704853 0.8786382592296224 … 0.0831959877928885 0.8717404061472732; 0.322341169750587 0.12504607953293195 … 0.10747085235449638 0.8317457388280489]
 [0.6231642412020931 0.20644282323995872 … 0.24825739278975745 0.9379695077323111; 0.649885620371928 0.1282039052316013 … 0.23726757281146385 0.9561185985111984; … ; 0.9610669641531648 0.9518738251107386 … 0.19658112236643613 0.9479489842038699; 0.42293095392901814 0.23554478880527824 … 0.21721063872759125 0.9129665303732556]
 [0.8064217831663169 0.4290433315831617 … 0.47212044222546273 1.1002378490309588; 0.8326027333844734 0.3570975188841845 … 0.4624816675166776 1.109981074243309; … ; 1.1225576781646647 1.0996953527944122 … 0.4256611042107394 1.1017810399761152; 0.6263065678764326 0.45880974883453374 … 0.43913813470918167 1.076730031291637]
 [1.011788797319829 0.6791957161224171 … 0.7238697419847683 1.2816285671160768; 1.0374009704291938 0.6147143919357021 … 0.715378395885817 1.2821714575152188; … ; 1.3028079275695588 1.2650396104261727 … 0.6829662600722733 1.2739415524922153; 0.8549056440887073 0.7095766703302432 … 0.6887777597276536 1.2597267167717103]
 [1.2007681733240227 0.9107733336325672 … 0.957233634098885 1.4477903406165331; 1.2259313136231034 0.853794286071939 … 0.9493325970423994 1.4400079258369485; … ; 1.46754337741848 1.4163339581662495 … 0.9211190972214953 1.431745389431645; 1.066456319964537 0.941589676580724 … 0.920220273005265 1.427270178199781]
 [1.3601320603678755 1.108484202618954 … 1.1569523487474571 1.5867404416422937; 1.3850462986573087 1.0587492824313673 … 1.149007375882964 1.5718885631606219; … ; 1.6047908168740699 1.5421876188913113 … 1.1245722382868861 1.5635757935515717; 1.246762040320754 1.139569730854355 … 1.1182597060707735 1.567260615783056]
 [1.4770281447593454 1.2571851223785955 … 1.3078297385368833 1.6869910319168056; 1.5019539346867523 1.214000809190364 … 1.299264457266209 1.666650187589749; … ; 1.7031593174412196 1.631695554365876 … 1.277950773524856 1.6582549849371404; 1.3817654947157572 1.2884228182759396 … 1.2677527893169367 1.6681207032897105]
 [1.5538717157994788 1.3600245703347167 … 1.413006061637878 1.7506175762150094; 1.5790669740103982 1.3227580082238681 … 1.403397756956604 1.7261219940748842; … ; 1.7647437813661409 1.68654670003153 … 1.3846091620963963 1.7176024448290295; 1.474213466651746 1.391362393729838 … 1.3717957905185165 1.731959504613163]
 [1.5979677435892137 1.4256073198529267 … 1.4810089716359751 1.7841442019538791; 1.6236422317691468 1.3937788571438043 … 1.4701072958524768 1.7565569562916177; … ; 1.7960668285113044 1.7128686106349293 … 1.4533651458969397 1.7478732636857124; 1.5319952523828093 1.4570176253008404 … 1.438915324425763 1.7653762724244424]
 ⋮
 [1.627374800359404 1.5074987344042412 … 1.5693504697985072 1.7873896816860724; 1.6549569776230746 1.4914639070357054 … 1.554349369237827 1.7552648583109334; … ; 1.7911256555189543 1.7017206229481303 … 1.5422609757773358 1.7459744994102793; 1.5980601441869215 1.538582751264939 … 1.5268718223398432 1.7672070014951833]
 [1.624377579502656 1.5208556596688194 … 1.5839180520598592 1.7750323227235347; 1.6524420238836517 1.5103080924377033 … 1.5678936558862304 1.743029337045385; … ; 1.7753436885084355 1.6873341446642527 … 1.5571066408697871 1.7335505797163835; 1.6069822687649136 1.5514408294909086 … 1.5425511210714067 1.754255293049952]
 [1.6202327515983097 1.5329766109942078 … 1.5963310906594328 1.7604322019390715; 1.6485079734464334 1.5279677089129169 … 1.5797348725972764 1.7295555031905783; … ; 1.756966056149856 1.6726166657673942 … 1.570155008807058 1.7199199879949425; 1.6145799437044743 1.562775783807136 … 1.556993399798323 1.7391258576181299]
 [1.6159582519360463 1.5448258549318392 … 1.6071845922097356 1.7442765804775304; 1.6439905537757562 1.5451862716055005 … 1.5907111628112736 1.7156620369595128; … ; 1.7367724479434485 1.6587002534007287 … 1.582237449856552 1.7059217515769918; 1.621512717090523 1.5735099644541644 … 1.5709807425546505 1.7226287210556592]
 [1.6120529485617252 1.5567290365060344 … 1.616341718996768 1.7267258440223285; 1.6391149734549162 1.5618996494502222 … 1.6009701508502594 1.7016057381325023; … ; 1.715137508442656 1.6462553167716691 … 1.593503712661951 1.691851850599939; 1.6276470229507691 1.5839222396425232 … 1.5845389104137446 1.7051249163042181]
 [1.6089397169079678 1.5685998895041364 … 1.6232758147684045 1.7081063882038097; 1.6339456759770847 1.577367417730447 … 1.6101938894654735 1.6876285382077447; … ; 1.6928555146163557 1.6359816687494135 … 1.6036595459871723 1.6780341327628618; 1.632304527127893 1.5938795882248578 … 1.5971090220627808 1.6872482283865922]
 [1.6071288801966452 1.5793206445191454 … 1.6273517596606022 1.6905302763135066; 1.628987022768824 1.5897026321507446 … 1.6173430737038188 1.6750285941495064; … ; 1.6728629703656712 1.6289359320661025 … 1.611676032578608 1.665864902627683; 1.6345388584060114 1.6023861083085258 … 1.6071610855250926 1.671292917844304]
 [1.60654514418106 1.5872752670819312 … 1.6288862938824875 1.6769045634470228; 1.6250310241807606 1.5975688983428764 … 1.6214680074283319 1.6655600080930866; … ; 1.6583290415780927 1.6252926977689413 … 1.6165702883822468 1.6570464602319128; 1.6344983170661866 1.6081612729654646 … 1.6135481984375242 1.6597184391124953]
 [1.6067513607242119 1.5917849132762867 … 1.6294937694046707 1.6697098441254405; 1.6230667919743458 1.601503827976721 … 1.6234707954079288 1.6605686092607441; … ; 1.6511172168874109 1.6239833155952748 … 1.6190733274906757 1.6525711670914922; 1.6339885904811435 1.611281131580729 … 1.61675088779101 1.653990740291758]

and plot the simulation using

import CairoMakie
lattice_animation(sol, :X1, lrs, "lattice_simulation.mp4")

More information on how to retrieve values from spatial simulations can be found here, and for plotting them, here. Finally, a list of functions for querying LatticeReactionSystems for various properties can be found here.

Spatial reactions

Spatial reactions describe reaction events which involve species across two connected compartments. Currently, only so-called transportation reactions are supported. These consist of:

  • A rate at which the reaction occurs. As for non-spatial reactions, this can be any expression. However, currently, it may only consist of parameters and other constants.
  • A single species which is transported from one compartment to an adjacent one.

At the occurrence of a transport reaction, the specific species moves to the adjacent compartment. Many common spatial models can be represented using transport reactions only. These can model phenomena such as diffusion or constant flux. A transportation reaction can be created using the @transportation_reaction macro. E.g. above we used

diffusion_rx = @transport_reaction D X1

to create a reaction where species $X$ moves at a constant rate $D$ between adjacent compartments (in the ODE this creates terms $D\cdot X1_i$, where $X1_i$ is the concentration of $X1$ in compartment $i$). Transport reactions may have rates depending on several parameters. E.g. to model a system with two species $X1$ and $X2$, where both species are transported at a rate which depends both on the species, but also on some non-uniform parameter which is unique to each connection (e.g. representing the area connecting two cells in a tissue) we could do:

dr_X1 = @transport_reaction D1*a X1
dr_X2 = @transport_reaction D2*a X2
Note

Any species which occurs is occurs in a transport reaction that is used to construct a LatticeReactionSystem must also occur in the corresponding non-spatial ReactionSystem.

Creating transport reactions programmatically

If models are created programmatically it is also possible to create transportation reactions programmatically. To do so, use the TransportReaction constructor, providing first the rate and then the transported species:

@variables t
@species X1(t) X2(t)
@parameters k1 k2 D [edgeparameter=true]
tr_X = TransportReaction(D, X1)

Note that in this example, we specifically designate $D$ as an edge parameter.

Defining discrete spatial domains (lattices)

Discrete spatial domains can represent:

  1. Systems which are composed of a (finite number of) compartments, where each compartment can be considered well-mixed (e.g. can be modelled non-spatially) and where (potentially) species can move between adjacent compartments. Tissues, where each compartment corresponds to a biological cell, are examples of such systems.
  2. Systems that are continuous in nature, but have been approximated as a discrete domain. Future Catalyst updates will include the ability for the definition, and automatic discretisation, of continuous domains. Currently, however, the user has to perform this discretisation themselves.

Catalyst supports three distinct types of lattices:

  • Cartesian lattices. These are grids where each grid point corresponds to a compartment. Spatial transportation is permitted between adjacent compartments.
  • Masked lattices. In these grids, only a subset of the grid point actually corresponds to compartments. Spatial transportation is permitted between adjacent compartments.
  • Unstructured (or graph) lattices. These are defined by graphs, where vertices correspond to compartments and edges connect adjacent compartments.

Here, Cartesian lattices are a subset of the masked lattices, which are a subset of the unstructured lattices. If possible, it is advantageous to use as narrow a lattice definition as possible (this may both improve simulation performance and simplify syntax). Cartesian and masked lattices can be defined as one, two, and three-dimensional. By default, these lattices assume that diagonally neighbouring compartments are non-adjacent (do not permit direct movement of species in between themselves). To change this, provide the diagonally_adjacent = true argument to your LatticeReactionSystem when it is created.

Defining Cartesian lattices

A Cartesian lattice is defined using the CartesianGrid function, which takes a single argument. For a 1d grid, simply provide the length of the grid as a single argument:

cgrid_1d = CartesianGrid(5)

For 2d and 3d grids, we instead provide a Tuple with the length of the grid in each dimension:

cgrid_2d = CartesianGrid((3, 9))
cgrid_3d = CartesianGrid((2, 4, 8))

Defining masked lattices

Masked lattices are defined through 1d, 2d, or 3d Boolean arrays. Each position in the array is true if it corresponds to a compartment, and false if it does not. E.g. to define a 1d grid corresponding to two disjoint sets, each with 3 compartments, use:

rgrid_1d = [true, true, true, false, true, true, true]

To define a 2d grid corresponding to the shape of an (laying) "8", we can use:

rgrid_2d = [
    true  true  true  true  true;
    true  false true  false true;
    true  true  true  true  true
]

Finally, a 4x5x6 3d grid of randomly distributed compartments can be created using:

rgrid_3d = rand([true, false], 4, 5, 6)

Defining unstructured lattices

To define unstructured lattices, we must first import the Graphs.jl package. Next, we can either use some pre-defined formula for building graphs, or build a graph from scratch. Here we create a cyclic graph (where each compartment is connected to exactly two other compartments):

using Graphs
cycle_graph(7)

Since graphs can represent any network of connected compartments, they do not have dimensions (like Cartesian or masked lattices). Another feature of graph lattices is that they can have non-symmetric connections, i.e. pairs of compartments where spatial movement of species is only permitted in one direction (in practice, this can be done for Cartesian and masked lattices as well, by defining non-uniform spatial rates and setting them to zero in one direction). This can be done by using a directed graph as input. E.g. here we define a directed cyclic graph, where movement is only allowed in one direction of the cycle:

cycle_digraph(7)

Non-uniform initial conditions and parameter values

For spatial models, initial conditions and parameter values are provided similarly as for non-spatial models. Wherever a single value is provided, it is used uniformly across the lattice. E.g. if we, for our previous two-state model, set

u0 = [:X1 => 1.0, :X2 => 2.0]
ps = [:k1 => 1.0, :k2 => 2.0, :D => 0.2]

The initial condition will be $1.0$ for $X1$ across compartments, and $2.0$ for $X2$. Furthermore, for each simulation, in each compartment, the value of $1.0$ will be used for $k1$ and $2.0$ for $k2$. Finally, the transportation rate of $X1$ (set by the parameter $D$) will be $0.2$ across all connections. To set non-uniform values, non-scalar values must be provided in the map. How to do this depends on how the lattice was defined. Furthermore, some parameters that are part of spatial reactions may have their value tied to connections between compartments, rather than compartments (we call these edge parameters). These are handled slightly differently. How to designate parameters as either edge parameters or compartment parameters is described here.

Below we describe how to set non-uniform values in the various cases.

Non-uniform compartment values for Cartesian lattices

To provide non-uniform values across a Cartesian lattice, simply provide the values in an array of the same dimension and size as the Cartesian lattice. E.g. for a 5x10 Cartesian lattice:

ccart_lattices = CartesianGrid((5, 10))

random values (uniformly distributed between $0$ and $1$) can be provided using

[:X1 => rand(5, 10), :X2 => 10.0]

Non-uniform values for parameters (which values are tied to compartments) are provided similarly.

Non-uniform compartment values for masked lattices

Non-uniform values for masked lattices are provided in the same manner as for Cartesian lattices (however, values at coordinates that do not hold compartments are ignored). E.g. To provide random values for a masked lattice contained within a 5x10 Cartesian lattices we can again set:

[:X1 => rand(5, 10), :X2 => 10.0]

If we want, it is also possible to provide the values as a sparse array with values only in the coordinates that corresponds to compartments.

Non-uniform compartment values for unstructured lattices

In graphs (which are used to represent unstructured lattices) each vertex (i.e. compartment) has a specific index. To set non-uniform values for unstructured lattices, provide a vector where the $i$'th value corresponds to the value in the compartment with index $i$ in the graph. E.g. for a graph with 5 vertices, where we want $X$ to be zero in all compartments bar one (where it is $1.0$) we use:

[:X1 => [0.0, 0.0, 0.0, 0.0, 1.0], :X2 => 10.0]

Non-uniform values for edge-parameters

Adjacent compartments are connected by edges (with which compartments are connected by edges being defined by the lattice). For unstructured lattices, it is possible (if a directed graph was used) to have edges from one compartment to another, but not in the opposite direction. For a lattice with $N$ compartments, edge values are set by a $NxN$ matrix, where value $(i,j)$ corresponds to the parameter's values in the edge going from compartment $i$to compartment $j$. This matrix can be either sparse or non-sparse. In the latter cases, values corresponding to non-existing edges are ignored.

Let's consider a 1d Cartesian lattice with 4 compartments. Here, an edge parameter's values are provided in a 4x4 matrix. For the Brusselator model described previously, $D$'s value was tied to edges. If we wish to set the value of $D$ to various values between $0.1$ and $0.4$ we can do:

ps = [:k1 => 1.0, :k2 => 2.0,
      :D => [
        0.0 0.1 0.0 0.0;
        0.1 0.0 0.2 0.0;
        0.0 0.2 0.0 0.3;
        0.0 0.0 0.3 0.0]
]

Here, the value at index $i,j$ corresponds to $D$'s value in the edge from compartment $i$ to compartment $j$. 0.0 is used for elements that do not correspond to an edge. The make_edge_p_values and make_directed_edge_values provide convenient interfaces for generating non-uniform edge parameter values.

Edge parameters and compartment parameters

Parameters can be divided into edge parameters and compartment parameters (initial condition values are always tied to compartments). Here, edge parameters have their values tied to edges, while compartment parameters have their values tied to compartments. All parameters that are part of the rates (or stoichiometries) of non-spatial reactions must be compartment parameters. Parameters that are part of spatial reactions can be either compartment parameters or edge parameters. When a spatial reaction's rate is computed, edge parameters fetch their values for from the edge of the transition, and compartment parameters from the compartment from which the edge originates.

When a LatticeReactionSystem is created, its parameters is the union of all parameters occurring in the (non-spatial) ReactionSystem and in all spatial reactions. By default, parameters occurring only in spatial reactions are considered edge parameters (and if they occur in the non-spatial ReactionSystem they are considered compartment parameters). It is, however, possible to designate a parameter specifically as an edge parameter (or not), by using the edgeparametermetadata. E.g. to designate that D (when declared in a non-spatial ReactionSystem using the DSL) is an edge parameter, not a compartment parameter, we use:

two_state_model = @reaction_network begin
    @parameters D [edgeparameter=true]
    (k1,k2), X1 <--> X2
end

To learn the compartment and edge parameters of a LatticeReaction, the vertex_parameters and edge_parameters functions can be used:

two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
diffusion_rx = @transport_reaction D X1
lattice = CartesianGrid((20,20))
lrs = LatticeReactionSystem(two_state_model, [diffusion_rx], lattice)
edge_parameters(lrs)
1-element Vector{Any}:
 D

Spatial modelling limitations

Many features which are supported for non-spatial ReactionSystems are currently unsupported for LatticeReactionSystems. This includes observables, algebraic and differential equations, hierarchical models, and events. It is possible that these features will be supported in the future. Furthermore, removal of conserved quantities is not supported when creating spatial ODEProblems.

If you are using Catalyst's features for spatial modelling, please give us feedback on how we can improve these features. Additionally, just letting us know that you use these features is useful, as it helps inform us whether continued development of spatial modelling features is worthwhile.