init

Author: TorkelE

docs/src/catalyst_functionality/dsl_description.md

@@ -567,3 +567,58 @@ species(rn)
 !!! note
     When using interpolation, expressions like `2$spec` won't work; the
     multiplication symbol must be explicitly included like `2*$spec`.
+
+## Including observables
+Sometimes, one might want to include observable variables. These are variables that can be computed directly from the other system variables (rather than having their values implicitly given through some differential equation). These can be introduced through the `@observables` option.
+
+Let us consider a simple example where two species ($X$ and $Y$) are produced and degraded at constant rates. They can also bind, forming a complex ($XY$). If we want to access the total amount of $X$ in the system we can create an observable that denotes this quantity ($Xtot = X + XY$). Here, we create observables for the total amount of $X$ and $Y$:
+```@example obs1
+using Catalyst # hide
+rn = @reaction_network begin
+  @observables begin
+    Xtot ~ X + XY
+    Ytot ~ Y + XY
+  end
+  (pX,dX), 0 <--> X
+  (pY,dY), 0 <--> Y
+  (kB,kD), X + Y <--> XY
+end
+```
+The `@observables` option is followed by one line for each observable formula (enclosed by a `begin ... end` block). The left-hand sides indicate the observables' names, and the right-hand sides how their values are computed. The two sides are separated by a `~`. 
+
+If we now simulate our model:
+```@example obs1
+using DifferentialEquations # hide
+u0 = [:X => 0.0, :Y => 0.0, :XY => 0.0]
+tspan = (0.0, 10.0)
+ps = [:pX => 1.0, :dX => 0.2, :pY => 1.0, :dY => 0.5, :kB => 1.0, :kD => 0.2]
+oprob = ODEProblem(rn, u0, tspan, ps)
+sol = solve(oprob)
+nothing # hide
+```
+we can index the solution using our observables (just like for [other variables](@ref simulation_structure_interfacing_solutions)). E.g. we can receive a vector with all $Xtot$ values using
+```@example obs1
+sol[:Xtot]
+```
+similarly, we can plot the values of $Xtot$ and $Ytot$ using
+```@example obs1
+plot(sol; idxs=[:Xtot, :Ytot])
+```
+
+If we only wish to provide a single observable, the `begin ... end` block is note required. E.g., to record only the total amount of $X$ we can use:
+```@example obs1
+using Catalyst # hide
+rn = @reaction_network begin
+  @observables  Xtot ~ X + XY
+  (pX,dX), 0 <--> X
+  (pY,dY), 0 <--> Y
+  (kB,kD), X + Y <--> XY
+end
+```
+
+Finally, some general rules for creating observables:
+- Observables can depend on any species, parameters, or variables, but not on other observables.
+- All observables components must be declared somewhere (i.e., they cannot only appear as a part of the observables formula).
+- Only a single `@observables` option block can be used in each `@reaction_network` call.
+- The left-hand side of the observables expression must be a single symbol, indicating the observable's name.
+- The right-hand side of the observables expression can be any valid algebraic expression.

src/reaction_network.jl

@@ -123,7 +123,7 @@ const forbidden_variables_error = let
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs, :compounds)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -358,8 +358,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
                        option_lines))
 
-    # Reads compounds options.
+    # Reads options.
     compound_expr, compound_species = read_compound_options(options)
+    observed_vars, observed_eqs = read_observed_options(options)
 
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
@@ -414,11 +415,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         $ivexpr
         $vars
         $sps
+        $observed_vars
         $comps
 
         Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym, $compssym),
-                                              $pssym; name = $name,
-                                              spatial_ivs = $sivs)
+                                              $pssym; name = $name, spatial_ivs = $sivs,
+                                              observed = $observed_eqs)
     end
 end
 
@@ -682,6 +684,46 @@ function recursive_find_reactants!(ex::ExprValues, mult::ExprValues,
     reactants
 end
 
+### DSL Options Handling ###
+# Most options handled in previous sections, when code re-organised, these should ideally be moved to the same place.
+
+# Reads the observables options. Outputs an expression ofr creating the obervable variables, and a vector of observable equations.
+function read_observed_options(options)
+    if haskey(options, :observables)
+        # Gets list of observable equations and prepares variable declaration expression.
+        # (`options[:observables]` inlucdes `@observables`, `.args[3]` removes this part)
+        observed_eqs = make_observed_eqs(options[:observables].args[3])
+        observed_vars = :(@variables)
+
+        # For each observable, checks for errors and adds the variable to `observed_vars`.
+        for obs_eq in observed_eqs.args
+            (obs_eq.args[1] != :~) && error("Malformed observable formula :\"$(obs_eq)\". Formula should contain two expressions separated by a '~'.")
+            (obs_eq.args[2] isa Symbol) || error("Malformed observable formula :\"$(obs_eq)\". Left hand side should be a single symbol.")
+            push!(observed_vars.args, obs_eq.args[2])
+        end
+    else  
+        # If option is not used, return empty expression and vector.
+        observed_vars = :()
+        observed_eqs = :([])
+    end
+    return observed_vars, observed_eqs
+end
+
+# From the input to the @observables options, creates a vector containing one equation for each observable.
+# Checks separate cases for "@obervables O ~ ..." and "@obervables begin ... end". Other cases errors.
+function make_observed_eqs(observables_expr)
+    if observables_expr.head == :call
+        return :([$(observables_expr)])
+    elseif observables_expr.head == :block
+        observed_eqs = :([])
+        for arg in observables_expr.args
+            push!(observed_eqs.args, arg)
+        end
+        return observed_eqs
+    else
+        error("Malformed observables option usage: $(observables_expr).")
+    end 
+end
 ### Functionality for expanding function call to custom and specific functions ###
 
 #Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.

test/dsl/dsl_options.jl

@@ -1,7 +1,7 @@
 #! format: off
 
 ### Fetch Packages and Set Global Variables ###
-using Catalyst, ModelingToolkit
+using Catalyst, ModelingToolkit, Plots
 @variables t
 
 ### Run Tests ###
@@ -339,3 +339,112 @@ let
     @test ModelingToolkit.defaults(rn27) == defs29
     @test merge(ModelingToolkit.defaults(rn28), defs28) == ModelingToolkit.defaults(rn27)
 end
+
+### Observables ###
+
+# Test basic functionality.
+# Tests various types of indexing.
+let 
+    rn = @reaction_network begin
+        @observables begin
+            X ~ Xi + Xa
+            Y ~ Y1 + Y2
+        end
+        (p,d), 0 <--> Xi
+        (k1,k2), Xi <--> Xa
+        (k3,k4), Y1 <--> Y2
+    end
+    @unpack X, Xi, Xa, Y, Y1, Y2, p, d, k1, k2, k3, k4 = rn
+
+    # Test that ReactionSystem have the correct properties.
+    @test length(species(rn)) == 4
+    @test length(states(rn)) == 4
+    @test length(observed(rn)) == 2
+    @test length(equations(rn)) == 6
+
+    @test isequal(observed(rn)[1], X ~ Xi + Xa)
+    @test isequal(observed(rn)[2], Y ~ Y1 + Y2)
+
+    # Tests correct indexing of solution.
+    u0 = [Xi => 0.0, Xa => 0.0, Y1 => 1.0, Y2 => 2.0]
+    ps = [p => 1.0, d => 0.2, k1 => 1.5, k2 => 1.5, k3 => 5.0, k4 => 5.0]
+
+    oprob = ODEProblem(rn, u0, (0.0, 1000.0), ps)
+    sol = solve(oprob, Tsit5())
+    @test sol[X][end] ≈ 10.0
+    @test sol[Y][end] ≈ 3.0
+    @test sol[rn.X][end] ≈ 10.0
+    @test sol[rn.Y][end] ≈ 3.0
+    @test sol[:X][end] ≈ 10.0
+    @test sol[:Y][end] ≈ 3.0
+
+    # Tests that observables can be used for plot indexing.
+    @test plot(sol; idxs=X).series_list[1].plotattributes[:y][end] ≈ 2.5
+    @test plot(sol; idxs=rn.X).series_list[1].plotattributes[:y][end] ≈ 2.5
+    @test plot(sol; idxs=:X).series_list[1].plotattributes[:y][end] ≈ 2.5
+    @test plot(sol; idxs=[X, Y]).series_list[2].plotattributes[:y][end] ≈ 3.0
+    @test plot(sol; idxs=[rn.X, rn.Y]).series_list[2].plotattributes[:y][end] ≈ 3.0
+    @test plot(sol; idxs=[:X, :Y]).series_list[2].plotattributes[:y][end] ≈ 3.0
+end
+
+# Tests for complicated observable formula.
+# Tests using a single observable (without begin/end statement).
+# Tests using observable component not part of reaction.
+# Tests using parameters in observables formula.
+let 
+    rn = @reaction_network begin
+        @parameters op_1 op_2
+        @species X4(t)
+        @observables X ~ X1^2 + op_1*(X2 + 2X3) + X1*X4/op_2 + p        
+        (p,d), 0 <--> X1
+        (k1,k2), X1 <--> X2
+        (k3,k4), X2 <--> X3
+    end
+    
+    u0 = Dict([:X1 => 1.0, :X2 => 2.0, :X3 => 3.0, :X4 => 4.0])
+    ps = Dict([:p => 1.0, :d => 0.2, :k1 => 1.5, :k2 => 1.5, :k3 => 5.0, :k4 => 5.0, :op_1 => 1.5, :op_2 => 1.5])
+
+    oprob = ODEProblem(rn, u0, (0.0, 1000.0), ps)
+    sol = solve(oprob, Tsit5())
+
+    @test sol[:X][1] == u0[:X1]^2 + ps[:op_1]*(u0[:X2] + 2*u0[:X3]) + u0[:X1]*u0[:X4]/ps[:op_2] + ps[:p]  
+end
+
+# Tests various erroneous declarations throw errors.
+let 
+    # System with undeclared component as observable.
+    @test_throws Exception @eval @reaction_network begin
+        @observables begin
+            X ~ X1 + X2
+        end
+        (p,d), 0 <--> X1
+    end
+
+    # System with observable in observable formula.
+    @test_throws Exception @eval @reaction_network begin
+        @observables begin
+            X ~ X1 + X2
+            X2 ~ 2X
+        end
+        (p,d), 0 <--> X1 + X2
+    end
+
+    # System with multiple observables blocks.
+    @test_throws Exception @eval @reaction_network begin
+        @observables begin
+            X ~ X1 + X2
+        end
+        @observables begin
+            X2 ~ 2(X1 + X2)
+        end
+        (p,d), 0 <--> X1 + X2
+    end
+
+    # System with no trivial observable left-hand side.
+    @test_throws Exception @eval @reaction_network begin
+        @observables begin
+            X + X2 ~ 2X1
+        end
+        (p,d), 0 <--> X1 + X2
+    end
+end