init

docs/src/catalyst_functionality/dsl_description.md

@@ -622,3 +622,67 @@ Finally, some general rules for creating observables:
 - 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.
+
+
+## Incorporating (differential) equations into reaction network models
+Some models cannot be purely described as a reaction network. E.g. consider the growth of a cell, where the rate of change in cell's volume depends on some growth factor. Here, the cell's volume would be described by a normal equation. Such equations can be incorporated into a model using the `@equations` option. Here, we create a model where a growth factor ($G$) is produced and degraded at a linear rates, and the rate of change in cell volume ($V$) os linear to the amount of growth factor:
+```@example eqs1
+using Catalyst #hide
+rn = @reaction_network begin
+  @equations begin
+    D(V) ~ G
+  end
+  (p,d), 0 <--> G
+end
+```
+Here, `D(V)` indicates the (time) derivative with respect to `D`. The differential equation left and right hand sides are separated by a `~`. The left-hand side should contain differential only, the right hand side can contain any algebraic expression.
+
+We can check the differential equation corresponding to this reaction network using latexify:
+```@example eqs1
+using Latexify
+latexify(rn; form=:ode)
+```
+We can also simulate it using the normal syntax
+```@example eqs1
+using DifferentialEquations, Plots # hide
+u0 = [:G => 0.0, :V => 0.1]
+ps = [:p => 1.0, :d => 0.5]
+oprob = ODEProblem(rn, u0, (0.0, 1.0), ps)
+sol = solve(oprob)
+plot(sol)
+```
+Here, growth is indefinite. To improve the model, [a callback](@ref advanced_simulations_callbacks) can be used to half the volume (cell division) once some threshold is reached.
+
+When creating differential equations this way, the subject of the differential is automatically inferred to be a variable, however, any component on the right-hand side must be declare somewhere in the macro. E.g. to add a scaling parameter ($k$), we must declare that $k$ is a parmaeter using the `@paraemters` option:
+```@example eqs1
+rn = @reaction_network begin
+  @parameters k
+  @equations begin
+    D(V) ~ k*G
+  end
+  (p,d), 0 <--> G
+end
+nothing #hide
+```
+
+It is possible to add several equations to the model. In this case, each have a separate line. E.g. to keep track of a supply of nutrition ($N$) in the growth media, we can use:
+```@example eqs1
+rn = @reaction_network begin
+  @parameters k
+  @equations begin
+    D(V) ~ G
+    D(N) ~ -G
+  end
+  (p,d), 0 <--> G
+end
+nothing #hide
+```
+
+When only a single equation is added, the `begin ... end` statement can be omitted. E.g., the first model can be declared equivalently using:
+```@example eqs1
+rn = @reaction_network begin
+  @equations D(V) ~ G
+  (p,d), 0 <--> G
+end
+nothing # hide
+```

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, :observables, :equations)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :differentials, :equations, :observables)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -361,14 +361,17 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     # Reads options.
     compound_expr, compound_species = read_compound_options(options)
     observed_vars, observed_eqs = read_observed_options(options)
-    vars_extracted, add_default_diff, equations = read_equations_options(options)
 
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
     species_declared = [extract_syms(options, :species); compound_species]
     parameters_declared = extract_syms(options, :parameters)
-    variables = [extract_syms(options, :variables); vars_extracted]
-
+    variables_declared = extract_syms(options, :variables)
+
+    # Reads more options.
+    vars_extracted, add_default_diff, equations = read_equations_options(options, variables_declared)
+    variables = vcat(variables_declared, vars_extracted)
+
     # handle independent variables
     if haskey(options, :ivs)
         ivs = Tuple(extract_syms(options, :ivs))
@@ -389,7 +392,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters = vcat(parameters_declared, parameters_extracted)
 
     # Create differential expression.
-    diffexpr = create_differential_expr(options, add_default_diff, [species; species; variables])
+    diffexpr = create_differential_expr(options, add_default_diff, [species; parameters; variables])
 
     # Checks for input errors.
     (sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
@@ -738,10 +741,10 @@ end
 # `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
 # `dtexpr`: If a differentialequation is defined, the default derrivative (D ~ Differential(t)) must be defined.
 # `equations`: a vector with the equations provided.
-function read_equations_options(options)
+function read_equations_options(options, variables_declared)
     # Prepares the equations. First, extracts equations from provided option (converting to block form if requried).
     # Next, uses MTK's `parse_equations!` function to split input into a vector with the equations.
-    eqs_input = haskey(options, :equations) ? options[:equations].args[3] : :()
+    eqs_input = haskey(options, :equations) ? options[:equations].args[3] : :(begin end)
     eqs_input = option_block_form(eqs_input)
     equations = Expr[]
     ModelingToolkit.parse_equations!(Expr(:block), equations, Dict{Symbol, Any}(), eqs_input)
@@ -752,11 +755,13 @@ function read_equations_options(options)
     vars_extracted = []
     add_default_diff = false
     for eq in equations
-        ((eq.head != :call) || (eq.args[1] != :~)) && error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
-        (eq.args[2].head != :call) && continue
+        ((eq.head != :call) || (eq.args[1] != :~)) && error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")     
+        (eq.args[2] isa Symbol || eq.args[2].head != :call) && continue
         if (eq.args[2].args[1] == :D) && (eq.args[2].args[2] isa Symbol) && (length(eq.args[2].args) == 2)
+            diff_var = eq.args[2].args[2]
+            in(diff_var, forbidden_symbols_error) && error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
             add_default_diff = true
-            push!(vars_extracted, eq.args[2].args[2])
+            in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
         end
     end
 
@@ -768,14 +773,16 @@ function create_differential_expr(options, add_default_diff, used_syms)
     # Creates the differential expression.
     # If differentials was provided as options, this is used as the initial expression.
     # If the default differential (D(...)) was used in equations, this is added to the expression.
-    diffexpr = (haskey(options, :differentials) ? options[:differentials] : MacroTools.striplines(:(begin end)))
+    diffexpr = (haskey(options, :differentials) ? options[:differentials].args[3] : MacroTools.striplines(:(begin end)))
+    diffexpr = option_block_form(diffexpr)
     add_default_diff && push!(diffexpr.args, :(D = Differential($(DEFAULT_IV_SYM))))
 
     # Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere.
     for dexpr in diffexpr.args
         (dexpr.head != :(=)) && error("Differential declaration must have form like D = Differential(t), instead \"$(dexpr)\" was given.")
         (dexpr.args[1] isa Symbol) || error("Differential left-hand side must be a single symbol, instead \"$(dexpr.args[1])\" was given.")
         in(dexpr.args[1], used_syms) && error("Differential name ($(dexpr.args[1])) is also a species, variable, or parameter. This is ambigious and not allowed.")
+        in(dexpr.args[1], forbidden_symbols_error) && error("A forbidden symbol ($(dexpr.args[1])) was used as a differential name.")
     end
 
     return diffexpr

test/dsl/dsl_options.jl

@@ -410,7 +410,7 @@ let
     @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.
+# Tests that various erroneous declarations throw errors.
 let 
     # System with undeclared component as observable.
     @test_throws Exception @eval @reaction_network begin
@@ -455,4 +455,198 @@ let
         end
         (p,d), 0 <--> X1 + X2
     end
+end
+
+### Equations ###
+
+# Basic checks on simple case with additional differential equations.
+# Checks indexing works.
+# Checks that non-block form for single equation works.
+let
+    # Creates model.
+    rn = @reaction_network rn begin
+        @parameters k
+        @equations begin
+            D(V) ~ X - k*V
+        end 
+        (p,d), 0 <--> X
+    end
+
+    @unpack k, p, d, X, V = rn
+    @variables t
+    D = Differential(t)
+
+    # Checks that the internal structures have the correct lengths.
+    @test length(species(rn)) == 1
+    @test length(states(rn)) == 2
+    @test length(reactions(rn)) == 2
+    @test length(equations(rn)) == 3
+
+    # Checks that the internal structures contain the correct stuff, and are correctly sorted.
+    @test isspecies(states(rn)[1])
+    @test !isspecies(states(rn)[2])
+    @test equations(rn)[1] isa Reaction
+    @test equations(rn)[2] isa Reaction
+    @test equations(rn)[3] isa Equation
+    @test isequal(equations(rn)[3], D(V) ~ X - k*V)
+
+    # Checks that simulations has the correct output
+    u0 = Dict([X => 1 + rand(rng), V => 1 + rand(rng)])
+    ps = Dict([p => 1 + rand(rng), d => 1 + rand(rng), k => 1 + rand(rng)])
+    oprob = ODEProblem(rn, u0, (0.0, 10000.0), ps)
+    sol = solve(oprob, Tsit5(); abstol=1e-9, reltol=1e-9)
+    @test sol[X][end] ≈ ps[p]/ps[d]
+    @test sol[V][end] ≈ ps[p]/(ps[d]*ps[k])
+
+    # Checks that set and get index works for variables.
+    @test oprob[V] == u0[V]
+    oprob[V] = 2.0
+    @test_broken oprob[V] == 2.0
+    integrator = init(oprob, Tsit5())
+    @test_broken integrator[V] == 2.0
+    integrator[V] = 5.0
+    @test integrator[V] == 5.0
+
+    # Checks that block form is not required when only a single equation is used.
+    rn2 = @reaction_network rn begin
+        @parameters k
+        @equations D(V) ~ X - k*V
+        (p,d), 0 <--> X
+    end
+    @test rn == rn2
+end
+
+# Tries complicated set of equations.
+# Tries with pre-declaring some variables (and not others).
+# Tries using default values.
+# Tries using mixing of parameters, variables, and species in different places.
+let 
+    rn = @reaction_network begin
+        @species S(t)=2.0
+        @variables Y(t)=0.5 Z(t)
+        @parameters p1 p2 p3
+        @equations begin
+            D(X) ~ p1*S - X
+            D(Y) ~ p2 + X - Y 
+            D(Z) ~ p + Y - Z 
+        end
+        (p*T,d), 0 <--> S
+        (p*Z,d), 0 <--> T
+    end
+
+    u0 = [:X => 1.0, :Z => 10.0, :S => 1.0, :T => 1.0]
+    ps = [:p1 => 0.5, :p2 => 1.0, :p3 => 1.0, :p => 1.0, :d => 1.0]
+    oprob = ODEProblem(rn, u0, (0.0, 1000.0), ps)
+    sol = solve(oprob, Rosenbrock23(); abstol=1e-9, reltol=1e-9)
+    @test sol[:S][end] ≈ 4
+    @test sol[:T][end] ≈ 4
+    @test sol[:X][end] ≈ 2
+    @test sol[:Y][end] ≈ 3
+    @test sol[:Z][end] ≈ 4
+end
+
+# Tries for reaction system without any reactions (just an equation).
+# Tries with interpolating a value into an equation.
+# Tries using rn.X notation for designating variables.
+# Tries for empty parameter vector.
+let 
+    c = 4.0
+    rn = complete(@reaction_network begin
+        @equations D(X) ~ $c - X
+    end)
+
+    u0 = [rn.X => 0.0]
+    ps = []
+    oprob = ODEProblem(rn, u0, (0.0, 100.0), ps)
+    sol = solve(oprob, Tsit5(); abstol=1e-9, reltol=1e-9)
+    @test sol[rn.X][end] ≈ 4.0
+end
+
+# Checks hierarchical model.
+let 
+    base_rn = @reaction_network begin
+        @equations begin
+            D(V1) ~ X - 2V1
+        end 
+        (p,d), 0 <--> X
+    end
+    @unpack X, V1, p, d = base_rn
+
+    internal_rn = @reaction_network begin
+        @equations begin
+            D(V2) ~ X - 3V2
+        end 
+        (p,d), 0 <--> X
+    end
+
+    rn = compose(base_rn, [internal_rn])
+
+    u0 = [V1 => 1.0, X => 3.0, internal_rn.V2 => 2.0, internal_rn.X => 4.0]
+    ps = [p => 1.0, d => 0.2, internal_rn.p => 2.0, internal_rn.d => 0.5]
+    oprob = ODEProblem(rn, u0, (0.0, 1000.0), ps)
+    sol = solve(oprob, Tsit5(); abstol=1e-9, reltol=1e-9)
+
+    @test sol[X][end] ≈ 5.0
+    @test sol[V1][end] ≈ 2.5
+    @test sol[internal_rn.X][end] ≈ 4.0
+    @test sol[internal_rn.V2][end] ≈ 4/3
+end
+
+# Tests that various erroneous declarations throw errors.
+let 
+    # Using = instead of ~ (for equation).
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(pi) = 1 - pi
+        (p,d), 0 <--> S
+    end
+
+    # Using ~ instead of = (for differential).
+    @test_throws Exception @eval @reaction_network begin
+        @differentials D ~ Differential(t)
+        (p,d), 0 <--> S
+    end
+
+    # Equation with component undeclared elsewhere.
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(X) ~ p - X
+        (P,D), 0 <--> S
+    end
+
+    # Using default differential D and a symbol D.
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(X) ~ -X
+        (P,D), 0 <--> S
+    end
+
+    # Declaring a symbol as a differential when it is used elsewhere.
+    @test_throws Exception @eval @reaction_network begin
+        @differentials d = Differential(t)
+        (p,d), 0 <--> S
+    end
+
+    # Declaring differential equation using a forbidden variable.
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(pi) ~ 1 - pi
+        (p,d), 0 <--> S
+    end
+
+    # Declaring forbidden symbol as differential.
+    @test_throws Exception @eval @reaction_network begin
+        @differentials pi = Differential(t)
+        (p,d), 0 <--> S
+    end
+
+    # System with derivatives with respect to several independent variables.
+    @test_throws Exception @eval @reaction_network begin
+        @ivs s t
+        @variables X(s) Y(t)
+        @differentials begin
+            Ds = Differential(s)
+            Dt = Differential(t)
+        end
+        @equations begin
+            Ds(X) ~ 1 - X
+            Dt(Y) ~ 1 - Y
+        end
+    end
 end