Interfacing Problems, Integrators, and Solutions

When simulating a model, one begins with creating a problem. Next, a simulation is performed on the problem, during which the simulation's state is recorded through an integrator. Finally, the simulation output is returned as a solution. This tutorial describes how to access (or modify) the state (or parameter) values of problem, integrator, and solution structures.

Generally, when we have a structure simulation_struct and want to interface with the unknown (or parameter) x, we use simulation_struct[:x] to access the value. For situations where a value is accessed (or changed) a large number of times, it can improve performance to first create a specialised getter/setter function.

Interfacing problem objects

We begin by demonstrating how we can interface with problem objects. First, we create an ODEProblem representation of a chemical cross-coupling model (where a catalyst, $C$, couples two substrates, $S₁$ and $S₂$, to form a product, $P$).

using Catalyst
cc_system = @reaction_network begin
    k₁, S₁ + C --> S₁C
    k₂, S₁C + S₂ --> CP
    k₃, CP --> C + P
end

u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
tspan = (0., 10.0)
ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]
oprob = ODEProblem(cc_system, u0, tspan, ps)

We can find a species's (or variable's) initial condition value by simply indexing with the species of interest as input. Here we check the initial condition value of $C$:

oprob[:C]

0.05

An almost identical notation can be used for parameters, however, here we use oprob.ps (rather than oprob): with the notation being identical for parameters:

oprob.ps[:k₁]

5.0

To retrieve several species initial condition (or parameter) values, simply give a vector input. Here we check the values of the two substrates ($S₁$ and $S₂$):

oprob[[:S₁, :S₂]]

2-element Vector{Float64}:
 1.0
 1.2

A problem's time span can be accessed through the tspan field:

oprob.tspan

(0.0, 10.0)

Remaking problems using the remake function

To modify a problem, the remake function should be used. It takes an already created problem, and returns a new, updated, one (the input problem is unchanged). The remake function takes the following inputs:

• The problem that it remakes.
• (optionally) u0: A vector with initial conditions that should be updated. The vector takes the same form as normal initial condition vectors, but does not need to be complete (in which case only a subset of the initial conditions are updated).
• (optionally) tspan: An updated time span (using the same format as time spans normally are given in).
• (optionally) p: A vector with parameters that should be updated. The vector takes the same form as normal parameter vectors, but does not need to be complete (in which case only a subset of the parameters are updated).

Here we modify our problem to increase the initial condition concentrations of the two substrates ($S₁$ and $S₂$), and also confirm that the new problem is different from the old (unchanged) one:

using OrdinaryDiffEqDefault
oprob_new = remake(oprob; u0 = [:S₁ => 5.0, :S₂ => 2.5])
oprob_new != oprob

true

Here, we instead use remake to simultaneously update all three fields:

oprob_new_2 = remake(oprob; u0 = [:C => 0.2], tspan = (0.0, 20.0), p = [:k₁ => 2.0, :k₂ => 2.0])

Typically, when using remake to update a problem, the common workflow is to overwrite the old one with the output. E.g. to set the value of k₁ to 5.0 in oprob, you would do:

oprob = remake(oprob; p = [:k₁ => 5.0])

Interfacing integrator objects

During a simulation, the solution is stored in an integrator object. Here, we will describe how to interface with these. The almost exclusive circumstance when integrator-interfacing is relevant is when simulation events are implemented through callbacks. However, to demonstrate integrator indexing in this tutorial, we will create one through the init function (while circumstances where one might want to use init function exist, since integrators are automatically created during simulations, these are rare).

integrator = init(oprob)

We can interface with our integrator using an identical syntax as was used for problems. The primary exception is that there is no remake function for integrators. Instead, we can update species and parameter values using normal indexing. Here we update, and then check the values of, first the species $C$ and then the parameter $k₁$:

integrator[:C] = 0.0
integrator[:C]

0.0

or a parameter:

integrator.ps[:k₂] = 1.0
integrator.ps[:k₂]

1.0

Note that here, species-interfacing yields (or changes) a simulation's current value for a species, not its initial condition.

If you are interfacing with jump simulation integrators, you must always call reset_aggregated_jumps!(integrator) afterwards.

Interfacing solution objects

Finally, we consider solution objects. First, we simulate our problem:

sol = solve(oprob)

Next, we can access the simulation's values using the same notation as previously. When we access a species's, its values across the full simulation is returned as a vector:

sol[:P]

40-element Vector{Float64}:
 0.0
 6.767734672693565e-6
 3.695247961692878e-5
 0.0001460513767817545
 0.0003594540166565959
 0.0007906628934836733
 0.00150858206813723
 0.00273726393513229
 0.004746921986536083
 0.008678361994276148
 ⋮
 0.37712924701331096
 0.4376038651315273
 0.4973649041265619
 0.5525966629042953
 0.6046370737738512
 0.6515259949342557
 0.6944615150868407
 0.7326391998759442
 0.7534902945744667

Parameter values can also be accessed (however, here we only get a single value):

sol.ps[:k₃]

100.0

Unlike for problems and integrators, species or parameter values of solutions cannot be changed.

A vector with the time values for all simulation time steps can be retrieved using

sol.t

40-element Vector{Float64}:
  0.0
  0.006877636354249612
  0.012703347894275383
  0.021456515180098498
  0.030894370672248542
  0.04334188107389527
  0.058215942527513456
  0.07776896100333044
  0.10349955221968685
  0.14515621694908853
  ⋮
  3.6234781654636397
  4.361840815227861
  5.1643597619354615
  5.984929309795551
  6.844006186652589
  7.7072436968947375
  8.591129091552913
  9.472113193562755
 10.0

To find simulation values at a specific time point, simply use this time point as input to your solution object (treating it as a function). I.e. here we get our simulation's values at time $t = 1.0$

sol(1.0)

6-element Vector{Float64}:
 0.8649253941887765
 0.027156866135933194
 0.02166469276091908
 1.0865900869496956
 0.0011784411031477487
 0.11223147194715666

This works whenever the simulations actually stopped at time $t = 1.0$ (if not, an interpolated value is returned). To get the simulation's values for a specific subset of species, we can use the idxs optional argument. I.e. here we get the value of $C$ at time $t = 1.0$

sol(1.0; idxs = [:C])

1-element Vector{Float64}:
 0.027156866135933194

Interfacing using specialised getter/setter functions

Internally, species and parameter values are stored in vectors. Whenever e.g. oprob[:C] is called, Julia must first find which index in the storage vector $C$ is stored in. Next, its value can be retrieved. If oprob[:C] is called a large number of times, this index must be found in each call. If a large number of such accesses are carried out, and performance is essential, it can be worthwhile to pre-compute a function to carry this out.

There exist four different functions, each returning a function for performing a specific type of interfacing:

• ModelingToolkit.getu: For accessing species values.
• ModelingToolkit.getp: For accessing parameter values.
• ModelingToolkit.setu: For changing species values.
• ModelingToolkit.setp: For changing parameter values.

For each species (or parameter) we wish to interface with, a new interfacing function must be created. Here we first creates a function for retrieving the value of $C$, and then use it for this purpose:

get_C = ModelingToolkit.getu(oprob, :C)
get_C(oprob)

0.05

Here, getu (as well as getp, setu, and setp) first takes the structure we wish to interface with, and then the target quantity. When using setu and setp, in the second step, we must also provide the update value:

set_C = ModelingToolkit.setu(oprob, :C)
set_C(oprob, 0.2)
get_C(oprob)

0.2

Like when indexing-based interfacing is used, these functions also work with vectors:

get_S = ModelingToolkit.getu(oprob, [:S₁, :S₂])
get_S(oprob)

2-element Vector{Float64}:
 1.0
 1.2

Interfacing using symbolic representations

When e.g. programmatic modelling is used, species and parameters can be represented as symbolic variables. These can be used to index a problem, just like symbol-based representations can. Here we create a simple two-state model programmatically, and use its symbolic variables to check, and update, an initial condition:

using Catalyst, OrdinaryDiffEqDefault
t = default_t()
@species X1(t) X2(t)
@parameters k1 k2
rxs = [
    Reaction(k1, [X1], [X2]),
    Reaction(k2, [X2], [X1])
]
@named two_state_model = ReactionSystem(rxs, t)
two_state_model = complete(two_state_model)

u0 = [X1 => 2.0, X2 => 0.0]
tspan = (0.0, 1.0)
ps = [k1 => 1.0, k2 => 2.0]
oprob = ODEProblem(two_state_model, u0, tspan, ps)

oprob = remake(oprob; u0 = [X1 => 5.0])
oprob[X1]

5.0

Symbolic variables can be used to access or update species or parameters for all the cases when Symbols can (including when using remake or e.g. getu).

An advantage when quantities are represented as symbolic variables is that symbolic expressions can be formed and used to index a structure. E.g. here we check the combined initial concentration of $X$ ($X1 + X2$) in our two-state problem:

oprob[X1 + X2]

5.0

Just like symbolic variables can be used to directly interface with a structure, symbolic variables stored in ReactionSystem models can be used:

oprob[two_state_model.X1 + two_state_model.X2]

5.0

This can be used to form symbolic expressions using model quantities when a model has been created using the DSL (as an alternative to @unpack). Alternatively, creating an observable, and then interface using its Symbol representation, is also possible.