Merge pull request #644 from SciML/lattice_reaction_system_may_2023

Spatial Reaction Network Implementation

Author: TorkelE

HISTORY.md

@@ -1,9 +1,32 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch)
+- Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases.
+- LatticeReactionSystem structure represents a spatial reaction network:
+  ```julia
+  rn = @reaction_network begin
+      (p,d), 0 <--> X
+  end
+  tr = @transport_reaction D X
+  lattice = Graphs.grid([5, 5])
+  lrs = LatticeReactionSystem(rn, [tr], lattice)
+```
+- Here, if a `u0` or `p` vector is given with scalar values:
+  ```julia
+  u0 = [:X => 1.0]
+  p = [:p => 1.0, :d => 0.5, :D => 0.1]
+  ```
+  this value will be used across the entire system. If their values are instead vectors, different values are used across the spatial system. Here
+  ```julia
+  X0 = zeros(25)
+  X0[1] = 1.0
+  u0 = [:X => X0]
+  ```
+  X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.)
+
 
 ## Catalyst 13.5
-- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reactin system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: 
+- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: 
 ```julia
 wilhelm_2009_model = @reaction_network begin
     k1, Y --> 2X

Project.toml

@@ -16,6 +16,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
@@ -40,6 +41,7 @@ ModelingToolkit = "8.66"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"
+RuntimeGeneratedFunctions = "0.5.12"
 SymbolicUtils = "1.0.3"
 Symbolics = "5.0.3"
 Unitful = "1.12.4"

src/Catalyst.jl

@@ -16,6 +16,10 @@ const MT = ModelingToolkit
 using Unitful
 @reexport using ModelingToolkit
 using Symbolics
+
+using RuntimeGeneratedFunctions
+RuntimeGeneratedFunctions.init(@__MODULE__)
+
 import Symbolics: BasicSymbolic
 import SymbolicUtils
 using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems,
@@ -33,6 +37,7 @@ import ModelingToolkit: check_variables,
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, Graphs
+import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne
 import DataStructures: OrderedDict, OrderedSet
 import Parameters: @with_kw_noshow
 
@@ -107,4 +112,22 @@ export balance_reaction
 function hc_steady_states end
 export hc_steady_states
 
+### Spatial Reaction Networks ###
+
+# spatial reactions
+include("spatial_reaction_systems/spatial_reactions.jl")
+export TransportReaction, TransportReactions, @transport_reaction
+export isedgeparameter
+
+# lattice reaction systems
+include("spatial_reaction_systems/lattice_reaction_systems.jl")
+export LatticeReactionSystem
+export spatial_species, vertex_parameters, edge_parameters
+
+# variosu utility functions
+include("spatial_reaction_systems/utility.jl")
+
+# spatial lattice ode systems.
+include("spatial_reaction_systems/spatial_ODE_systems.jl")
+
 end # module

src/spatial_reaction_systems/lattice_reaction_systems.jl

@@ -0,0 +1,91 @@
+### Lattice Reaction Network Structure ###
+# Describes a spatial reaction network over a graph.
+# Adding the "<: MT.AbstractTimeDependentSystem" part messes up show, disabling me from creating LRSs.
+struct LatticeReactionSystem{S,T} # <: MT.AbstractTimeDependentSystem 
+    # Input values.
+    """The reaction system within each compartment."""
+    rs::ReactionSystem{S}
+    """The spatial reactions defined between individual nodes."""
+    spatial_reactions::Vector{T}
+    """The graph on which the lattice is defined."""
+    lattice::SimpleDiGraph{Int64}
+
+    # Derived values.
+    """The number of compartments."""
+    num_verts::Int64
+    """The number of edges."""
+    num_edges::Int64
+    """The number of species."""
+    num_species::Int64
+    """Whenever the initial input was a digraph."""
+    init_digraph::Bool
+    """Species that may move spatially."""
+    spat_species::Vector{BasicSymbolic{Real}}
+    """
+    All parameters related to the lattice reaction system
+    (both with spatial and non-spatial effects).
+    """
+    parameters::Vector{BasicSymbolic{Real}}
+    """
+    Parameters which values are tied to vertexes (adjacencies), 
+    e.g. (possibly) have a unique value at each vertex of the system.
+    """
+    vertex_parameters::Vector{BasicSymbolic{Real}}
+    """
+    Parameters which values are tied to edges (adjacencies), 
+    e.g. (possibly) have a unique value at each edge of the system.
+    """
+    edge_parameters::Vector{BasicSymbolic{Real}}
+
+    function LatticeReactionSystem(rs::ReactionSystem{S}, spatial_reactions::Vector{T},
+                                   lattice::DiGraph; init_digraph = true) where {S, T}
+        # There probably some better way to ascertain that T has that type. Not sure how.
+        if !(T <: AbstractSpatialReaction) 
+            error("The second argument must be a vector of AbstractSpatialReaction subtypes.") 
+        end
+
+        if isempty(spatial_reactions)
+            spat_species = Vector{BasicSymbolic{Real}}[]
+        else
+            spat_species = unique(reduce(vcat, [spatial_species(sr) for sr in spatial_reactions]))
+        end
+        num_species = length(unique([species(rs); spat_species]))
+        rs_edge_parameters = filter(isedgeparameter, parameters(rs))
+        if isempty(spatial_reactions)
+            srs_edge_parameters = Vector{BasicSymbolic{Real}}[]
+        else
+            srs_edge_parameters = setdiff(reduce(vcat, [parameters(sr) for sr in spatial_reactions]), parameters(rs))
+        end
+        edge_parameters = unique([rs_edge_parameters; srs_edge_parameters])
+        vertex_parameters = filter(!isedgeparameter, parameters(rs))
+        # Ensures the parameter order begins similarly to in the non-spatial ReactionSystem.
+        ps = [parameters(rs); setdiff([edge_parameters; vertex_parameters], parameters(rs))]    
+
+        foreach(sr -> check_spatial_reaction_validity(rs, sr; edge_parameters=edge_parameters), spatial_reactions)   
+        return new{S,T}(rs, spatial_reactions, lattice, nv(lattice), ne(lattice), num_species, 
+                            init_digraph, spat_species, ps, vertex_parameters, edge_parameters)
+    end
+end
+function LatticeReactionSystem(rs, srs, lat::SimpleGraph)
+    return LatticeReactionSystem(rs, srs, DiGraph(lat); init_digraph = false)
+end
+
+### Lattice ReactionSystem Getters ###
+
+# Get all species.
+species(lrs::LatticeReactionSystem) = unique([species(lrs.rs); lrs.spat_species])
+# Get all species that may be transported.
+spatial_species(lrs::LatticeReactionSystem) = lrs.spat_species
+
+# Get all parameters.
+ModelingToolkit.parameters(lrs::LatticeReactionSystem) = lrs.parameters
+# Get all parameters which values are tied to vertexes (compartments).
+vertex_parameters(lrs::LatticeReactionSystem) = lrs.vertex_parameters
+# Get all parameters which values are tied to edges (adjacencies).
+edge_parameters(lrs::LatticeReactionSystem) = lrs.edge_parameters
+
+# Gets the lrs name (same as rs name).
+ModelingToolkit.nameof(lrs::LatticeReactionSystem) = nameof(lrs.rs)
+
+# Checks if a lattice reaction system is a pure (linear) transport reaction system.
+is_transport_system(lrs::LatticeReactionSystem) = all(sr -> sr isa TransportReaction, lrs.spatial_reactions)