Interfacing with Lattice Problems, Integrators, and Solutions

We have previously described how to retrieve species and parameter values stored in non-spatial problems, integrators, and solutions. This section describes similar workflows for simulations based on DiscreteSpaceReactionSystems.

Generally, while for non-spatial systems these operations can typically be done by indexing a structure directly, e.g. through

sol[:X]

there are no equally straightforward interfaces for spatial simulations. Typically, helper functions have to be used, e.g

spat_getu(sol, :X, dsrs)

Furthermore, there are some cases of interfacing which are currently not supported (e.g. updating parameter values in JumpProblems). It is likely that these interfaces will be improved in the future (i.e. by introducing a similar syntax to the current non-spatial one). Finally, we note that many of the functions presented below have not been as extensively optimised for performance as other parts of the Catalyst package. Hence, you should take care when designing workflows which requires using them a large number of times.

Retrieving values from space simulations

Let us consider a simulation of a DiscreteSpaceReactionSystem:

using Catalyst, OrdinaryDiffEqDefault
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
diffusion_rx = @transport_reaction D X1
space = CartesianGrid((2,3))
dsrs = DiscreteSpaceReactionSystem(two_state_model, [diffusion_rx], space)

u0 = [:X1 => [0.0 0.0 0.0; 2.0 2.0 2.0], :X2 => 0.0]
tspan = (0.0, 1.0)
ps = [:k1 => 2.0, :k2 => 1.0, :D => 0.1]
oprob = ODEProblem(dsrs, u0, tspan, ps)
sol = solve(oprob)

To retrieve the values of $X1$ across the simulation we use the spat_getu function. It takes three arguments:

• The solution objects from which we wish to retrieve values.
• The species which values we wish to retrieve.
• The DiscreteSpaceReactionSystem which was simulated.

spat_getu(sol, :X1, dsrs)

13-element Vector{Matrix{Float64}}:
 [0.0 0.0 0.0; 2.0 2.0 2.0]
 [4.988650432250024e-5 4.988650432250024e-5 4.988650432250024e-5; 1.998952234018206 1.998952234018206 1.998952234018206]
 [0.0005458858864143923 0.0005458858864143923 0.0005458858864143923; 1.9885184230183552 1.9885184230183552 1.9885184230183552]
 [0.0041835786833196815 0.0041835786833196815 0.004183578683319682; 1.9110504808166604 1.9110504808166604 1.9110504808166604]
 [0.010555124017922137 0.010555124017922135 0.010555124017922137; 1.7708834128823219 1.770883412882322 1.7708834128823219]
 [0.01751130416769682 0.017511304167696813 0.017511304167696817; 1.6099324747184667 1.609932474718467 1.6099324747184667]
 [0.02475469979790144 0.024754699797901424 0.024754699797901434; 1.4307707484082977 1.4307707484082983 1.4307707484082977]
 [0.03151854077808042 0.0315185407780804 0.0315185407780804; 1.2486002831094976 1.2486002831094982 1.2486002831094976]
 [0.0372977162591279 0.03729771625912789 0.03729771625912789; 1.0766136067128713 1.076613606712872 1.0766136067128715]
 [0.04182477748198844 0.04182477748198845 0.04182477748198844; 0.928433562106622 0.9284335621066225 0.9284335621066222]
 [0.04528546428561481 0.04528546428561482 0.045285464285614835; 0.8115561816724685 0.8115561816724687 0.8115561816724685]
 [0.0481822071541792 0.048182207154179216 0.048182207154179216; 0.7274839359461986 0.7274839359461988 0.7274839359461985]
 [0.05042586336225185 0.05042586336225186 0.050425863362251855; 0.6826244818849131 0.6826244818849131 0.6826244818849131]

Here, the output is a vector with $X1$'s value at each simulation time step. How the species's value is represented at each time step depends on the space which was originally used to create the DiscreteSpaceReactionSystem:

• For Cartesian spaces, an array of the same size as the Cartesian space is used. Each array element corresponds to the species's value in the corresponding compartment.
• For masked spaces, a sparse array of the same size as the masked space is used. Each filled array element corresponds to the species's value in the corresponding compartment. Unfilled array elements correspond to positions without compartments.
• For unstructured (graph) spaces, vectors are used. The i'th element in the vectors corresponds to the species's value in the i'th compartment.

Unlike for non-spatial simulations, spat_getu does not take vector (e.g. spat_getu(sol, [:X1, :X2], dsrs)) or symbolic expression (e.g. spat_getu(sol, [X1 + X2], dsrs)) inputs. However, it is possible to use symbolic variables as input (e.g. spat_getu(sol, two_state_model.X1, dsrs)).

Retrieving space simulations values at specific time points

Just like for non-spatial solutions, it is possible to access the simulation's values at designated time points. This is possible even if the simulation did not stop at those specific time points (in which case an interpolated value is returned). To do this, the desired time points to sample are provided as a vector to spat_getu using the optional argument t. E.g. here we retrieve the simulation's (interpolated) values at time points 0.5 and 0.75:

spat_getu(sol, :X1, dsrs; t = [0.5, 0.75])

2-element Vector{Matrix{Float64}}:
 [0.042009527416151625 0.04200952741615163 0.042009527416151625; 0.9221641543195104 0.9221641543195109 0.9221641543195106]
 [0.04694938927951838 0.046949389279518396 0.0469493892795184; 0.7602471207337441 0.7602471207337443 0.7602471207337441]

Retrieving and updating species values in problems and integrators

Let us consider a spatial ODEProblem

using Catalyst, OrdinaryDiffEqDefault
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
diffusion_rx = @transport_reaction D X1
space = CartesianGrid((2,3))
dsrs = DiscreteSpaceReactionSystem(two_state_model, [diffusion_rx], space)

u0 = [:X1 => [0.0 0.0 0.0; 2.0 2.0 2.0], :X2 => 0.0]
tspan = (0.0, 1.0)
ps = [:k1 => 2.0, :k2 => 1.0, :D => 0.1]
oprob = ODEProblem(dsrs, u0, tspan, ps)

We can retrieve the species values stored in oprob using the spat_getu function. It uses identical syntax as for simulations (except that you cannot specify a time point). However, it returns a single set of species values (while for simulations it returns a vector across different time steps):

spat_getu(oprob, :X1, dsrs)

2×3 Matrix{Float64}:
 0.0  0.0  0.0
 2.0  2.0  2.0

Again, the format used corresponds to the space used to create the original DiscreteSpaceReactionSystem. Here, even if a species has homogeneous values, the full format is used.

spat_getu(oprob, :X2, dsrs)

2×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0

For both problems and integrators, species values can be updated using the spat_setu! function. It uses a similar syntax as spat_getu, but takes a fourth argument which is the new values to use for the designated species:

spat_setu!(oprob, :X1, dsrs, [1.0 2.0 3.0; 4.0 5.0 6.0])

Here, the same format (which depends on the used space) is used for the species's new values, as which is used when initially designating their initial conditions. I.e. to make $X1$'s initial condition values uniform we can call

spat_setu!(oprob, :X1, dsrs, 1.0)

Species values in integrators can be interfaced with using identical syntax as for problems.

Retrieving and updating parameter values in problems and integrators

Retrieval and updating of parameter values for problems and integrators works similarly as for species, but with the following differences:

• The spat_getp and spat_setp! functions are used.
• It is currently not possible to interface with parameter values of JumpProblems and their integrators.
• After parameter values are modified, the rebuild_spat_internals! function must be applied before the problem/integrator can be used for further analysis.
• Updating of edge parameters is limited and uses a different interface.

Let us consider the spatial ODEProblem we previously declared. We can check the value of $k1$ by using spat_getp

spat_getp(oprob, :k1, dsrs)

2×3 Matrix{Float64}:
 2.0  2.0  2.0
 2.0  2.0  2.0

Next, we can update it using spat_setp! (here we also confirm that it now has the updated values):

spat_setp!(oprob, :k1, dsrs, [1.0 2.0 3.0; 4.0 5.0 6.0])
spat_getp(oprob, :k1, dsrs)

2×3 Matrix{Float64}:
 1.0  2.0  3.0
 4.0  5.0  6.0

If we now were to simulate oprob, the simulation would not take the updated value of $k1$ into account. For our changes to take effect we might first need to call rebuild_spat_internals! with oprob as an input

rebuild_spat_internals!(oprob)

There are two different circumstances when rebuild_spat_internals! must be called:

• When modifying the value of an edge parameter.
• When changing a parameter from having spatially uniform values to spatially non-uniform values, or the other way around.

Parameter values of integrators can be interfaced with just like for problems (this is primarily relevant when using callbacks). Again, after doing so, the rebuild_spat_internals! function might need to be applied to the integrator.

Retrieving and updatingedge parameter values in problems and integrators

The spat_getp and spat_setp! functions cannot currently be applied to edge parameters. Instead, to access the value of an edge parameter, use

oprob.ps[:D]

1×1 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
 0.1

To update an edge parameter's value, use

oprob.ps[:D] = [0.2]

This interface is somewhat limited, and the following aspects should be noted:

• Edge parameter values can only be interfaced with if the edge parameter's value is spatially uniform.
• When accessing an (spatially uniform) edge parameter's value, its single value will be encapsulated in a vector.
• When setting an (spatially uniform) edge parameter's value, you must encapsulate the new value in a vector.