expanding equations when passed to @reaction_network

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "14.4"
+version = "14.4.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/dsl.jl

@@ -297,7 +297,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
 
     # Get macro options.
     if length(unique(arg.args[1] for arg in option_lines)) < length(option_lines)
-        error("Some options where given multiple times.")
+        error("Some options were given multiple times.")
     end
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
         option_lines))
@@ -315,12 +315,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads more options.
+    # Reads equations. 
     vars_extracted, add_default_diff, equations = read_equations_options(
         options, variables_declared)
     variables = vcat(variables_declared, vars_extracted)
 
-    # handle independent variables
+    # Handle independent variables
     if haskey(options, :ivs)
         ivs = Tuple(extract_syms(options, :ivs))
         ivexpr = copy(options[:ivs])
@@ -339,14 +339,16 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         combinatoric_ratelaws = true
     end
 
-    # Reads more options.
+    # Reads observables.
     observed_vars, observed_eqs, obs_syms = read_observed_options(
         options, [species_declared; variables], all_ivs)
 
+    # Collect species and parameters, including ones inferred from the reactions. 
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
         variables)))
-    species_extracted, parameters_extracted = extract_species_and_parameters!(reactions,
-        declared_syms)
+    species_extracted, parameters_extracted = extract_species_and_parameters!(
+        reactions, declared_syms)
+
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
@@ -376,9 +378,11 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, get_rxexprs(reaction))
     end
     for equation in equations
+        equation = escape_equation_RHS!(equation)
         push!(rxexprs.args, equation)
     end
 
+    # Output code corresponding to the reaction system. 
     quote
         $ivexpr
         $ps
@@ -572,7 +576,7 @@ function get_rxexprs(rxstruct)
     subs_stoich_init = deepcopy(subs_init)
     prod_init = isempty(rxstruct.products) ? nothing : :([])
     prod_stoich_init = deepcopy(prod_init)
-    reaction_func = :(Reaction($(recursive_expand_functions!(rxstruct.rate)), $subs_init,
+    reaction_func = :(Reaction($(recursive_escape_functions!(rxstruct.rate)), $subs_init,
         $prod_init, $subs_stoich_init, $prod_stoich_init,
         metadata = $(rxstruct.metadata)))
     for sub in rxstruct.substrates
@@ -904,17 +908,28 @@ end
 
 ### Generic Expression Manipulation ###
 
-# Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
-function recursive_expand_functions!(expr::ExprValues)
+# Recursively traverses an expression and escapes all the user-defined functions. Special function calls like "hill(...)" are not expanded. 
+function recursive_escape_functions!(expr::ExprValues)
     (typeof(expr) != Expr) && (return expr)
-    foreach(i -> expr.args[i] = recursive_expand_functions!(expr.args[i]),
+    foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i]),
         1:length(expr.args))
     if expr.head == :call
         !isdefined(Catalyst, expr.args[1]) && (expr.args[1] = esc(expr.args[1]))
     end
     expr
 end
 
+# Recursively escape functions in the right-hand-side of an equation written using user-defined functions. Special function calls like "hill(...)" are not expanded.
+function escape_equation_RHS!(eqexpr::Expr)
+    lhs = recursive_escape_functions!(eqexpr.args[2])
+    eqexpr.args[2] = lhs
+
+    rhs = recursive_escape_functions!(eqexpr.args[3])
+    eqexpr.args[3] = rhs
+
+    eqexpr
+end
+
 # Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical).
 function tup_leng(ex::ExprValues)
     (typeof(ex) == Expr && ex.head == :tuple) && (return length(ex.args))

test/dsl/dsl_options.jl

@@ -13,6 +13,7 @@ seed = rand(rng, 1:100)
 
 # Sets the default `t` to use.
 t = default_t()
+D = default_time_deriv()
 
 ### Tests `@parameters`, `@species`, and `@variables` Options ###
 
@@ -152,6 +153,68 @@ let
     @test issetequal(parameters(rn5), @parameters k B X)
 end
 
+# Test whether user-defined functions are properly expanded in equations. 
+let
+    f(A, t) = 2*A*t
+
+    # Test user-defined function
+    rn = @reaction_network begin
+        @equations D(A) ~ f(A, t)
+    end
+    @test length(equations(rn)) == 1
+    @test equations(rn)[1] isa Equation
+    @species A(t)
+    @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
+
+
+    hill_unregistered(A, v, K, n) = v*(A^n) / (A^n + K^n)
+    rn2 = @reaction_network begin
+        @parameters v K n
+        @equations D(A) ~ hill_unregistered(A, v, K, n)
+    end
+    @test length(equations(rn2)) == 1
+    @test equations(rn2)[1] isa Equation
+    @parameters v K n
+    @test isequal(equations(rn2)[1], D(A) ~ v*(A^n) / (A^n + K^n))
+
+    @register_symbolic hill2(A, v, K, n)
+    # Test Catalyst function
+    rn2r = @reaction_network begin
+        @parameters v K n
+        @equations D(A) ~ hill2(A, v, K, n)
+    end
+    @test length(equations(rn2r)) == 1
+    @test equations(rn2r)[1] isa Equation
+    @parameters v K n
+    @test isequal(equations(rn2r)[1], D(A) ~ hill2(A, v, K, n))
+
+
+    rn3 = @reaction_network begin
+        @species Iapp(t)
+        @equations begin
+            D(A) ~ Iapp
+            Iapp ~ f(A,t)
+        end
+    end
+    @test length(equations(rn3)) == 2
+    @test equations(rn3)[1] isa Equation
+    @test equations(rn3)[2] isa Equation
+    @variables Iapp(t)
+    @test isequal(equations(rn3)[1], D(A) ~ Iapp)
+    @test isequal(equations(rn3)[2], Iapp ~ 2*A*t)
+
+
+    g(A, K, n) = A^n + K^n
+    rn4 = @reaction_network begin
+        @parameters v K n
+        @equations D(A) ~ hill(A, v, K, n)*g(A, K, n)
+    end
+    @test length(equations(rn4)) == 1
+    @test equations(rn4)[1] isa Equation
+    @parameters v n
+    @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n))
+end
+
 # Test inferring with stoichiometry symbols and interpolation.
 let
     @parameters k g h gg X y [isconstantspecies = true]