From f94fe04f1f8d5695046699d693ec89dd80b59b28 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sat, 16 Nov 2024 20:43:41 +0000
Subject: [PATCH 01/17] save prorgess

---
 src/dsl.jl              | 48 +++++++++++++++++++++++------------------
 src/reactionsystem.jl   | 34 +++++++++++++++--------------
 test/dsl/dsl_options.jl | 32 +++++++++++++++++++++------
 3 files changed, 71 insertions(+), 43 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 91a2fa3851..4f5ea7c46d 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -315,11 +315,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads equations. 
-    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))
@@ -339,23 +334,29 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         combinatoric_ratelaws = true
     end
 
-    # 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. 
+    # Collect species and parameters.
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
-        variables)))
+        variables_declared)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(
         reactions, declared_syms)
 
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
+    # Reads equations.
+    designated_syms = [species; parameters; variables_declared]
+    vars_extracted, add_default_diff, equations = read_equations_options(
+        options, designated_syms)
+    variables = vcat(variables_declared, vars_extracted)
+
     # Create differential expression.
     diffexpr = create_differential_expr(
         options, add_default_diff, [species; parameters; variables], tiv)
 
+    # Reads observables.
+    observed_vars, observed_eqs, obs_syms = read_observed_options(
+        options, [species_declared; variables], all_ivs)
+
     # Checks for input errors.
     (sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
         error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.")
@@ -384,7 +385,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, equation)
     end
 
-    # Output code corresponding to the reaction system. 
+    # Output code corresponding to the reaction system.
     quote
         $ivexpr
         $ps
@@ -682,7 +683,7 @@ end
 # `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
 # `dtexpr`: If a differential equation is defined, the default derivative (D ~ Differential(t)) must be defined.
 # `equations`: a vector with the equations provided.
-function read_equations_options(options, variables_declared)
+function read_equations_options(options, syms_declared)
     # Prepares the equations. First, extracts equations from provided option (converting to block form if required).
     # 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] : :(begin end)
@@ -694,7 +695,8 @@ function read_equations_options(options, variables_declared)
     # Loops through all equations, checks for lhs of the form `D(X) ~ ...`.
     # When this is the case, the variable X and differential D are extracted (for automatic declaration).
     # Also performs simple error checks.
-    vars_extracted = Vector{Symbol}()
+    vars_extracted = OrderedSet{Union{Symbol, Expr}}()
+    excluded_syms = syms_declared
     add_default_diff = false
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)
@@ -702,22 +704,26 @@ function read_equations_options(options, variables_declared)
         end
 
         # Checks if the equation have the format D(X) ~ ... (where X is a symbol). This means that the
-        # default differential has been used. X is added as a declared variable to the system, and
-        # we make a note that a differential D = Differential(iv) should be made as well.
+        # default differential has been used and we make a note that it should be decalred in the DSL output.
         lhs = eq.args[2]
-        # if lhs: is an expression. Is a function call. The function's name is D. Calls a single symbol.
+        # If lhs: is an expression. Is a function call. The function's name is D. It has a single argument.
         if (lhs isa Expr) && (lhs.head == :call) && (lhs.args[1] == :D) &&
            (lhs.args[2] isa Symbol)
             diff_var = lhs.args[2]
             if in(diff_var, forbidden_symbols_error)
                 error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
             end
-            add_default_diff = true
-            in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
+            if !add_default_diff
+                add_default_diff = true
+                excluded_syms = [excluded_syms; :D]
+            end
         end
+
+        # Any undecalred symbolic variables encountered should be extracted as variables.
+        add_syms_from_expr!(vars_extracted, eq, excluded_syms)
     end
 
-    return vars_extracted, add_default_diff, equations
+    return collect(vars_extracted), add_default_diff, equations
 end
 
 # Creates an expression declaring differentials. Here, `tiv` is the time independent variables,
@@ -917,7 +923,7 @@ end
 
 ### Generic Expression Manipulation ###
 
-# Recursively traverses an expression and escapes all the user-defined functions. Special function calls like "hill(...)" are not expanded. 
+# 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_escape_functions!(expr.args[i]),
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 83ed069ea0..5536ed0c40 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -97,9 +97,9 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     terminallinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
 
-    checkedrobust::Bool = false 
+    checkedrobust::Bool = false
     robustspecies::Vector{Int} = Vector{Int}(undef, 0)
-    deficiency::Int = -1 
+    deficiency::Int = -1
 end
 #! format: on
 
@@ -215,11 +215,11 @@ end
 
 ### ReactionSystem Structure ###
 
-""" 
+"""
 WARNING!!!
 
-The following variable is used to check that code that should be updated when the `ReactionSystem` 
-fields are updated has in fact been updated. Do not just blindly update this without first checking 
+The following variable is used to check that code that should be updated when the `ReactionSystem`
+fields are updated has in fact been updated. Do not just blindly update this without first checking
 all such code and updating it appropriately (e.g. serialization). Please use a search for
 `reactionsystem_fields` throughout the package to ensure all places which should be updated, are updated.
 """
@@ -318,7 +318,7 @@ struct ReactionSystem{V <: NetworkProperties} <:
     """
     discrete_events::Vector{MT.SymbolicDiscreteCallback}
     """
-    Metadata for the system, to be used by downstream packages. 
+    Metadata for the system, to be used by downstream packages.
     """
     metadata::Any
     """
@@ -480,10 +480,10 @@ function ReactionSystem(iv; kwargs...)
     ReactionSystem(Reaction[], iv, [], []; kwargs...)
 end
 
-# Called internally (whether DSL-based or programmatic model creation is used). 
+# Called internally (whether DSL-based or programmatic model creation is used).
 # Creates a sorted reactions + equations vector, also ensuring reaction is first in this vector.
-# Extracts potential species, variables, and parameters from the input (if not provided as part of 
-# the model creation) and creates the corresponding vectors. 
+# Extracts potential species, variables, and parameters from the input (if not provided as part of
+# the model creation) and creates the corresponding vectors.
 # While species are ordered before variables in the unknowns vector, this ordering is not imposed here,
 # but carried out at a later stage.
 function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
@@ -495,7 +495,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     any(in(obs_vars), us_in) &&
         error("Found an observable in the list of unknowns. This is not allowed.")
 
-    # Creates a combined iv vector (iv and sivs). This is used later in the function (so that 
+    # Creates a combined iv vector (iv and sivs). This is used later in the function (so that
     # independent variables can be excluded when encountered quantities are added to `us` and `ps`).
     t = value(iv)
     ivs = Set([t])
@@ -524,6 +524,8 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     # Extracts any species, variables, and parameters that occur in (non-reaction) equations.
     # Creates the new reactions + equations vector, `fulleqs` (sorted reactions first, equations next).
     if !isempty(eqs)
+        println(eqs)
+        println(iv)
         osys = ODESystem(eqs, iv; name = gensym())
         fulleqs = CatalystEqType[rxs; equations(osys)]
         union!(us, unknowns(osys))
@@ -560,7 +562,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     end
     psv = collect(new_ps)
 
-    # Passes the processed input into the next `ReactionSystem` call.    
+    # Passes the processed input into the next `ReactionSystem` call.
     ReactionSystem(fulleqs, t, usv, psv; spatial_ivs, continuous_events,
         discrete_events, observed, kwargs...)
 end
@@ -1062,8 +1064,8 @@ end
 
 ### General `ReactionSystem`-specific Functions ###
 
-# Checks if the `ReactionSystem` structure have been updated without also updating the 
-# `reactionsystem_fields` constant. If this is the case, returns `false`. This is used in 
+# Checks if the `ReactionSystem` structure have been updated without also updating the
+# `reactionsystem_fields` constant. If this is the case, returns `false`. This is used in
 # certain functionalities which would break if the `ReactionSystem` structure is updated without
 # also updating these functionalities.
 function reactionsystem_uptodate_check()
@@ -1241,7 +1243,7 @@ end
 ### `ReactionSystem` Remaking ###
 
 """
-    remake_ReactionSystem_internal(rs::ReactionSystem; 
+    remake_ReactionSystem_internal(rs::ReactionSystem;
         default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T}
 
 Takes a `ReactionSystem` and remakes it, returning a modified `ReactionSystem`. Modifications depend
@@ -1274,7 +1276,7 @@ function set_default_metadata(rs::ReactionSystem; default_reaction_metadata = []
     # Currently, `noise_scaling` is the only relevant metadata supported this way.
     drm_dict = Dict(default_reaction_metadata)
     if haskey(drm_dict, :noise_scaling)
-        # Finds parameters, species, and variables in the noise scaling term.       
+        # Finds parameters, species, and variables in the noise scaling term.
         ns_expr = drm_dict[:noise_scaling]
         ns_syms = [Symbolics.unwrap(sym) for sym in get_variables(ns_expr)]
         ns_ps = Iterators.filter(ModelingToolkit.isparameter, ns_syms)
@@ -1414,7 +1416,7 @@ function ModelingToolkit.compose(sys::ReactionSystem, systems::AbstractArray; na
         MT.collect_scoped_vars!(newunknowns, newparams, ssys, iv)
     end
 
-    if !isempty(newunknowns) 
+    if !isempty(newunknowns)
         @set! sys.unknowns = union(get_unknowns(sys), newunknowns)
         sort!(get_unknowns(sys), by = !isspecies)
         @set! sys.species = filter(isspecies, get_unknowns(sys))
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 8174598fd4..9adecbb030 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -901,6 +901,26 @@ let
     @test 5*sol[:Y][end] ≈ sol[:S][end] + sol[:X][end]
 end
 
+# Tests that the correct symbolic variables are infered as species, variables, and paraemters.
+let
+    rn = @reaction_network begin
+        @parameters p1 p2
+        @species X1(t) X2(t)
+        @variables W(t)
+        @equations begin
+            D(V1) ~ p1 * X1 - k1 * V1 + W
+            k2 * V2 ~ D(V2) + p2 * X2
+            V3 + X3 ~ V1^2 + X2^2
+        end
+        (k1, k2), X1 <--> X2
+        (k3, k4), X2 <--> X3
+    end
+
+    @test issetequal(species(rn), @species X1(t) X2(t) X3(t))
+    @test issetequal(parameters(rn), @parameters p1 p2 k1 k2 k3 k4)
+    @test issetequal(nonspecies(rn), @variables V1(t) V2(t) V3(t) W(t))
+end
+
 # Tests that various erroneous declarations throw errors.
 let
     # Using = instead of ~ (for equation).
@@ -951,7 +971,7 @@ let
     @test isequal(rl, k1*A^2)
 end
 
-# Test whether user-defined functions are properly expanded in equations. 
+# Test whether user-defined functions are properly expanded in equations.
 let
     f(A, t) = 2*A*t
 
@@ -965,7 +985,7 @@ let
     @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
 
 
-    # Test whether expansion happens properly for unregistered/registered functions. 
+    # Test whether expansion happens properly for unregistered/registered functions.
     hill_unregistered(A, v, K, n) = v*(A^n) / (A^n + K^n)
     rn2 = @reaction_network begin
         @parameters v K n
@@ -978,7 +998,7 @@ let
 
     hill2(A, v, K, n) = v*(A^n) / (A^n + K^n)
     @register_symbolic hill2(A, v, K, n)
-    # Registered symbolic function should not expand. 
+    # Registered symbolic function should not expand.
     rn2r = @reaction_network begin
         @parameters v K n
         @equations D(A) ~ hill2(A, v, K, n)
@@ -1009,9 +1029,9 @@ let
     @named rn3_sym = ReactionSystem(eq, t)
     rn3_sym = complete(rn3_sym)
     @test isequivalent(rn3, rn3_sym)
-    
-    
-    # Test more complicated expression involving both registered function and a user-defined function. 
+
+
+    # Test more complicated expression involving both registered function and a user-defined function.
     g(A, K, n) = A^n + K^n
     rn4 = @reaction_network begin
         @parameters v K n

From a1f916b08a922524771e3559f01558546303ecb0 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 20:26:38 +0000
Subject: [PATCH 02/17] new inference of D and tests and History file

---
 HISTORY.md                                    |  23 +++
 src/reactionsystem.jl                         |   2 -
 test/dsl/dsl_options.jl                       | 194 +++++++++++++++++-
 .../coupled_equation_crn_systems.jl           |  40 +---
 4 files changed, 208 insertions(+), 51 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 78d89beb1e..456f9a6657 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -7,6 +7,29 @@
   (at the time the release is made). If you need a dependency version increased,
   please open an issue and we can update it and make a new Catalyst release once
   testing against the newer dependency version is complete. 
+- New formula for inferring variables from equations (declared using the `@equations` options) in the DSL. The order of inference of species/variables/parameters is now:
+    (1) Every symbol explicitly declared using `@species`, `@variables`, and `@parameters` are assigned to the correct category.
+    (2) Every symbol used as a reaction reactant is inferred as a species.
+    (3) Every symbol not declared in (1) or (2) that occurs in an expression provided after `@equations` is inferred as a variable.
+    (4) Every symbol not declared in (1), (2), or (3) that occurs either as a reaction rate or stoichiometric coefficient is inferred to be a parameter.
+E.g. in
+```julia
+@reaction_network begin
+    @equations V1 + S ~ V2^2
+    (p + S + V1), S --> 0
+end
+```
+`S` is inferred as a species, `V1` and `V2` as variables, and `p` as a parameter. The previous special cases for the `@observables`, `@compounds`, and `@differentials` options still hold. Finally, the `@require_declaration` options (described in more detail below) can now be used to require everything to be explicitly declared.
+- New formula for determining whether the default differentials have been used within an `@equations` option. Right now, if any expression `D(...)` is encountered (where `...`) can be anything, this is inferred as usage of the default differential D. E.g. in the following equations `D` is inferred as a differential with respect to the default independent variable:
+```julia
+@reaction_network begin
+  @equations D(V) + V ~ 1
+end
+@reaction_network begin
+  @equations D(D(V)) ~ 1
+end
+```
+Please note that this cannot be used at the same time as `D` is used to represent a species, variable, or parameter. 
 - Array symbolics support is more consistent with ModelingToolkit v9. Parameter
   arrays are no longer scalarized by Catalyst, while species and variables
   arrays still are (as in ModelingToolkit). As such, parameter arrays should now
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 5536ed0c40..a2cd0ecf2e 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -524,8 +524,6 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     # Extracts any species, variables, and parameters that occur in (non-reaction) equations.
     # Creates the new reactions + equations vector, `fulleqs` (sorted reactions first, equations next).
     if !isempty(eqs)
-        println(eqs)
-        println(iv)
         osys = ODESystem(eqs, iv; name = gensym())
         fulleqs = CatalystEqType[rxs; equations(osys)]
         union!(us, unknowns(osys))
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 6a706b0044..51aead7a44 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -424,6 +424,74 @@ let
     end
 end
 
+# Tests that explicitly declaring a single symbol as several things does not work.
+# Several of these are broken, but note sure how to test broken-ness on `@test_throws false Exception @eval`.
+# Relevant issue: https://github.com/SciML/Catalyst.jl/issues/1173
+let
+    # Species + parameter.
+    @test_broken false #@test_throws Exception @eval @reaction_network begin
+        #@species X(t)
+        #@parameters X
+    #end
+
+    # Species + variable.
+    @test_broken false #@test_throws Exception @eval @reaction_network begin
+        #@species X(t)
+        #@variables X(t)
+    #end
+
+    # Variable + parameter.
+    @test_broken false #@test_throws Exception @eval @reaction_network begin
+        #@variables X(t)
+        #@parameters X
+    #end
+
+    # Species + differential.
+    @test_throws false Exception @eval @reaction_network begin
+        @species X(t)
+        @differentials X = Differential(t)
+    end
+
+    # Parameter + differential.
+    @test_throws false Exception @eval @reaction_network begin
+        @parameters X
+        @differentials X = Differential(t)
+    end
+
+    # Variable + differential.
+    @test_throws false Exception @eval @reaction_network begin
+        @variables X(t)
+        @differentials X = Differential(t)
+    end
+
+    # Parameter + observable (species/variable + observable is OK, as this e.g. provide additional observables information).
+    @test_broken false #@test_throws Exception @eval @reaction_network begin
+        #@species Y(t)
+        #@parameters X
+        #@observables X ~ Y
+    #end
+
+    # Species + compound.
+    @test_broken false #@test_throws Exception @eval @reaction_network begin
+        #@species X(t) O(t)
+        #@compounds begin X(t) ~ 2O end
+    #end
+
+    # Parameter + compound.
+    @test_broken false #@test_throws Exception @eval @reaction_network begin
+        #@species O(t)
+        #@parameters X
+        #@compounds begin X(t) ~ 2O end
+    #end
+
+    # Variable + compound.
+    @test_broken false #@test_throws Exception @eval @reaction_network begin
+        #@species O(t)
+        #@variables X(t)
+        #@compounds begin X(t) ~ 2O end
+    #end
+end
+
 ### Test Independent Variable Designations ###
 
 # Test ivs in DSL.
@@ -450,6 +518,118 @@ let
     @test issetequal(Catalyst.get_sivs(rn), [x])
 end
 
+### Test Symbolic Variable Inference ###
+
+# Basic checks that that symbolic variables not explicitly declared are correctly inferred.
+let
+    # Case 1 (a reaction only).
+    rn1 = @reaction_network begin
+        (p1/(S1+p2) + S2), S1 --> S2
+    end
+    @test issetequal(species(rn1), [rn1.S1, rn1.S2])
+    @test issetequal(parameters(rn1), [rn1.p1, rn1.p2])
+
+    # Case 2 (reactions and equations).
+    rn2 = @reaction_network begin
+        @equations V1 + log(V2 + S1) ~ V2^2
+        (p1/V1 + S1 + log(S2 + V2 + p2)), S1 --> S2
+    end
+    @test issetequal(species(rn2), [rn2.S1, rn2.S2])
+    @test issetequal(nonspecies(rn2), [rn2.V1, rn2.V2])
+    @test issetequal(parameters(rn2), [rn2.p1, rn2.p2])
+
+    # Case 3 (reaction and equations with a differential).
+    rn3 = @reaction_network begin
+        @equations begin
+            D(V1) ~ S1 + V1
+            V2 + S2 ~ V1^2 + V2^2
+        end
+        (p1/V1 + S1 + log(S2 + V2 + p2)), S1 --> S2
+    end
+    @test issetequal(species(rn3), [rn3.S1, rn3.S2])
+    @test issetequal(nonspecies(rn3), [rn3.V1, rn3.V2])
+    @test issetequal(parameters(rn3), [rn3.p1, rn3.p2])
+
+    # Case 4 (reactions and equations with a pre-declared parameter).
+    rn4 = @reaction_network begin
+        @parameters p1
+        @equations V1 + sin(p1 + S1) ~ S2*V2
+        (p1+p2+V1+V2+S1+S2), S1 --> S2
+    end
+    @test issetequal(species(rn4), [rn2.S1, rn2.S2])
+    @test issetequal(nonspecies(rn4), [rn2.V1, rn2.V2])
+    @test issetequal(parameters(rn4), [rn2.p1, rn2.p2])
+
+    # Case 5 (algebraic equation containing D, which is pre-declared as a species).
+    rn5 = @reaction_network begin
+        @species D(t)
+        @equations D * (S1 + V1 + V2) ~ S2
+        (p1 + p2*(D + V1 + V2 + S2 + S2)), S1 --> S2 + D
+    end
+    @test issetequal(species(rn5), [rn5.S1, rn5.S2, rn5.D])
+    @test issetequal(nonspecies(rn5), [rn5.V1, rn5.V2])
+    @test issetequal(parameters(rn5), [rn5.p1, rn5.p2])
+
+    # Case 6 (algebraic equation containing D, which is pre-declared as a parameter).
+    rn6 = @reaction_network begin
+        @parameters D
+        @equations D * (S1 + V1 + V2) ~ S2
+        (p1 + p2*(D + V1 + V2 + S2 + S2)), S1 --> S2
+    end
+    @test issetequal(species(rn6), [rn6.S1, rn6.S2])
+    @test issetequal(nonspecies(rn6), [rn6.V1, rn6.V2])
+    @test issetequal(parameters(rn6), [rn6.p1, rn6.p2, rn6.D])
+
+    # Case 7 (algebraic equation containing D, which is pre-declared as a variable).
+    rn7 = @reaction_network begin
+        @variables D(t)
+        @equations D * (S1 + V1 + V2) ~ S2
+        (p1 + p2*(D + V1 + V2 + S2 + S2)), S1 --> S2
+    end
+    @test issetequal(species(rn7), [rn7.S1, rn7.S2])
+    @test issetequal(nonspecies(rn7), [rn7.V1, rn7.V2, rn7.D])
+    @test issetequal(parameters(rn7), [rn7.p1, rn7.p2])
+
+    # Case 8 (reactions, equations, and a custom differential).
+    rn8 = @reaction_network begin
+        @differentials Δ = Differential(t)
+        @equations Δ(V1) + Δ(V2) + log(V2 + S1) ~  S2
+        (p1/V1 + S1 + log(S2 + V2 + p2)), S1 --> S2
+    end
+    @test issetequal(species(rn8), [rn8.S1, rn8.S2])
+    @test issetequal(nonspecies(rn8), [rn8.V1, rn8.V2])
+    @test issetequal(parameters(rn8), [rn8.p1, rn8.p2])
+end
+
+# Checks that various cases where symbolic variables cannot (or shouldn't) be inferred generate errors.
+let
+    # Species/variables/parameter named after default differential used as function call.
+    # In the future, species/variables should be usable this way (designating a time delay).
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(V) ~ 1 - V
+        d, D --> 0
+    end
+    @test_throws Exception @eval @reaction_network begin
+        @variables D(t)
+        @equations D(V) ~ 1 - V
+        d, X --> 0
+    end
+    @test_throws Exception @eval @reaction_network begin
+        @parameters D
+        @equations D(V) ~ 1 - V
+        d, X --> 0
+    end
+
+    # Symbol only occurring in events.
+    @test_throws Exception @eval @reaction_network begin
+        @discrete_event (X > 1.0) => [V => V/2]
+        d, X --> 0
+    end
+    @test_throws Exception @eval @reaction_network begin
+        @continuous_event [X > 1.0] => [V => V/2]
+        d, X --> 0
+    end
+end
 
 ### Observables ###
 
@@ -790,7 +970,7 @@ end
 # Check that DAE is solved correctly.
 let
     rn = @reaction_network rn begin
-        @parameters k
+        @parameters k d
         @variables X(t) Y(t)
         @equations begin
             X + 5 ~ k*S
@@ -940,10 +1120,10 @@ let
     @test 5*sol[:Y][end] ≈ sol[:S][end] + sol[:X][end]
 end
 
-# Tests that the correct symbolic variables are infered as species, variables, and paraemters.
+# Tests that the correct symbolic variables are inferred as species, variables, and parameters.
 let
     rn = @reaction_network begin
-        @parameters p1 p2
+        @parameters p1 p2 k1 k2 k3 k4
         @species X1(t) X2(t)
         @variables W(t)
         @equations begin
@@ -968,12 +1148,6 @@ let
         @equations X = 1 - S
         (p,d), 0 <--> S
     end
-
-    # Equation with component undeclared elsewhere.
-    @test_throws Exception @eval @reaction_network begin
-        @equations X ~ p - S
-        (P,D), 0 <--> S
-    end
 end
 
 # test combinatoric_ratelaws DSL option
@@ -1082,7 +1256,7 @@ let
     @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n))
 end
 
-### test that @no_infer properly throws errors when undeclared variables are written
+### test that @no_infer properly throws errors when undeclared variables are written ###
 
 import Catalyst: UndeclaredSymbolicError
 let
diff --git a/test/reactionsystem_core/coupled_equation_crn_systems.jl b/test/reactionsystem_core/coupled_equation_crn_systems.jl
index 3b4921ff96..8f631d58d3 100644
--- a/test/reactionsystem_core/coupled_equation_crn_systems.jl
+++ b/test/reactionsystem_core/coupled_equation_crn_systems.jl
@@ -899,38 +899,6 @@ end
 
 # Checks that various misformatted declarations yield errors.
 let
-    # Symbol in equation not appearing elsewhere (1).
-    @test_throws Exception @eval @reaction_network begin
-        @equations D(V) ~ -X
-    end
-
-    # Symbol in equation not appearing elsewhere (2).
-    @test_throws Exception @eval @reaction_network begin
-        @equations 1 + log(x) ~ 2X
-    end
-
-    # Attempting to infer differential variable not isolated on lhs (1).
-    @test_throws Exception @eval @reaction_network begin
-        @equations D(V) + 1 ~ 0
-    end
-
-    # Attempting to infer differential variable not isolated on lhs (2).
-    @test_throws Exception @eval @reaction_network begin
-        @equations -1.0 ~ D(V)
-    end
-
-    # Attempting to infer differential operator not isolated on lhs (1).
-    @test_throws Exception @eval @reaction_network begin
-        @variables V(t)
-        @equations D(V) + 1 ~ 0
-    end
-
-    # Attempting to infer a variable when using a non-default differential.
-    @test_throws Exception @eval @reaction_network begin
-        @differentials Δ = Differential(t)
-        @equations Δ(V) ~ -1,0
-    end
-
     # Attempting to create a new differential from an unknown iv.
     @test_throws Exception @eval @reaction_network begin
         @differentials D = Differential(τ)
@@ -944,14 +912,8 @@ let
 
     # Several equations without `begin ... end` block.
     @test_throws Exception @eval @reaction_network begin
-        @variables V(t)
         @equations D(V) + 1 ~ - 1.0
-    end
-
-    # Undeclared differential.
-    @test_throws Exception @eval @reaction_network begin
-        @species V
-        @equations Δ(V) ~ -1.0
+        @equations D(W) + 1 ~ - 1.0
     end
 
     # System using multiple ivs.

From 006ae00108ed6a0668302a2f08ec7996ed6258d1 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 20:40:18 +0000
Subject: [PATCH 03/17] update history file

---
 HISTORY.md              | 2 +-
 src/dsl.jl              | 7 +------
 test/dsl/dsl_options.jl | 6 +++---
 3 files changed, 5 insertions(+), 10 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 456f9a6657..d6d381585b 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -20,7 +20,7 @@ E.g. in
 end
 ```
 `S` is inferred as a species, `V1` and `V2` as variables, and `p` as a parameter. The previous special cases for the `@observables`, `@compounds`, and `@differentials` options still hold. Finally, the `@require_declaration` options (described in more detail below) can now be used to require everything to be explicitly declared.
-- New formula for determining whether the default differentials have been used within an `@equations` option. Right now, if any expression `D(...)` is encountered (where `...`) can be anything, this is inferred as usage of the default differential D. E.g. in the following equations `D` is inferred as a differential with respect to the default independent variable:
+- New formula for determining whether the default differentials have been used within an `@equations` option. Now, if any expression `D(...)` is encountered (where `...` can be anything), this is inferred as usage of the default differential D. E.g. in the following equations `D` is inferred as a differential with respect to the default independent variable:
 ```julia
 @reaction_network begin
   @equations D(V) + V ~ 1
diff --git a/src/dsl.jl b/src/dsl.jl
index 8d9f851224..886bbecaa6 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -290,7 +290,7 @@ struct UndeclaredSymbolicError <: Exception
     msg::String
 end
 
-function Base.showerror(io::IO, err::UndeclaredSymbolicError) 
+function Base.showerror(io::IO, err::UndeclaredSymbolicError)
     print(io, "UndeclaredSymbolicError: ")
     print(io, err.msg)
 end
@@ -328,11 +328,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads equations.
-    vars_extracted, add_default_diff, equations = read_equations_options(
-        options, variables_declared; requiredec)
-    variables = vcat(variables_declared, vars_extracted)
-
     # Handle independent variables
     if haskey(options, :ivs)
         ivs = Tuple(extract_syms(options, :ivs))
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 51aead7a44..f1e52586b3 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -447,19 +447,19 @@ let
     #end
 
     # Species + differential.
-    @test_throws false Exception @eval @reaction_network begin
+    @test_throws Exception @eval @reaction_network begin
         @species X(t)
         @differentials X = Differential(t)
     end
 
     # Parameter + differential.
-    @test_throws false Exception @eval @reaction_network begin
+    @test_throws Exception @eval @reaction_network begin
         @parameters X
         @differentials X = Differential(t)
     end
 
     # Variable + differential.
-    @test_throws false Exception @eval @reaction_network begin
+    @test_throws Exception @eval @reaction_network begin
         @variables X(t)
         @differentials X = Differential(t)
     end

From 7d6aa95e73c194bb456650f007672b801c57d5b3 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 20:41:36 +0000
Subject: [PATCH 04/17] dsl file fix

---
 src/dsl.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 886bbecaa6..ff9217f5f9 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -347,10 +347,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         combinatoric_ratelaws = true
     end
 
-    # Reads observables.
-    observed_vars, observed_eqs, obs_syms = read_observed_options(
-        options, [species_declared; variables], all_ivs; requiredec)
-
     # Collect species and parameters, including ones inferred from the reactions.
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
         variables_declared)))

From da50bf8fa30842868fef2b97a53bf27d8362f5ac Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:02:42 +0000
Subject: [PATCH 05/17] add update to read equations function back in

---
 src/dsl.jl | 46 +++++++++++++++++++++++++---------------------
 1 file changed, 25 insertions(+), 21 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index ff9217f5f9..a9a259d946 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -702,7 +702,7 @@ end
 # `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
 # `dtexpr`: If a differential equation is defined, the default derivative (D ~ Differential(t)) must be defined.
 # `equations`: a vector with the equations provided.
-function read_equations_options(options, variables_declared; requiredec = false)
+function read_equations_options(options, syms_declared; requiredec = false)
     # Prepares the equations. First, extracts equations from provided option (converting to block form if required).
     # 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] : :(begin end)
@@ -722,35 +722,39 @@ function read_equations_options(options, variables_declared; requiredec = false)
             error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
         end
 
-        # Checks if the equation have the format D(X) ~ ... (where X is a symbol). This means that the
-        # default differential has been used and we make a note that it should be decalred in the DSL output.
-        lhs = eq.args[2]
-        # If lhs: is an expression. Is a function call. The function's name is D. It has a single argument.
-        if (lhs isa Expr) && (lhs.head == :call) && (lhs.args[1] == :D) &&
-           (lhs.args[2] isa Symbol)
-            diff_var = lhs.args[2]
-            if in(diff_var, forbidden_symbols_error)
-                error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
-            elseif (!in(diff_var, variables_declared)) && requiredec
-                throw(UndeclaredSymbolicError(
-                                              "Unrecognized symbol $(diff_var) was used as a variable in an equation: \"$eq\". Since the @require_declaration flag is set, all variables in equations must be explicitly declared via @variables, @species, or @parameters."))
-            else
-                add_default_diff = true
-                in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
-            end
-            if !add_default_diff
-                add_default_diff = true
-                excluded_syms = [excluded_syms; :D]
-            end
+        # If the default differential (`D`) is used, record that it should be decalred later on.
+        if !in(eq, excluded_syms) && find_D_call(eq)
+            requiredec && throw(UndeclaredSymbolicError(
+                "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option."))
+            add_default_diff = true
+            excluded_syms = [excluded_syms; :D]
         end
 
         # Any undecalred symbolic variables encountered should be extracted as variables.
+        # Additional step required to handle `requiredec = true` (to be improved later).
+        prev_vars_extracted = deepcopy(vars_extracted)
         add_syms_from_expr!(vars_extracted, eq, excluded_syms)
+        if requiredec && length(prev_vars_extracted) < length(vars_extracted)
+            throw(UndeclaredSymbolicError(
+                "Unrecognized symbols $(setdiff(vars_extracted, prev_vars_extracted)) was used in an equation: \"$eq\".  Since the flag @require_declaration is set, all variables must be declared with the @species, @parameters, or @variables macros."))
+        end
     end
 
     return collect(vars_extracted), add_default_diff, equations
 end
 
+# Searches an expresion `expr` and returns true if it have any subexpression `D(...)` (where `...` can be anything).
+# Used to determine whether the default differential D has been used in any equation provided to `@equations`.
+function find_D_call(expr)
+    return if Base.isexpr(expr, :call) && expr.args[1] == :D
+        true
+    elseif expr isa Expr
+        any(find_D_call, expr.args)
+    else
+        false
+    end
+end
+
 # Creates an expression declaring differentials. Here, `tiv` is the time independent variables,
 # which is used by the default differential (if it is used).
 function create_differential_expr(options, add_default_diff, used_syms, tiv)

From d4cdcd5a581edc58464b62e4b9de7ddbca984550 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:14:07 +0000
Subject: [PATCH 06/17] merge fixes

---
 src/dsl.jl | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index a9a259d946..2e9076a245 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -353,14 +353,17 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(
         reactions, declared_syms; requiredec)
 
-    species = vcat(species_declared, species_extracted)
-    parameters = vcat(parameters_declared, parameters_extracted)
-
-    # Reads equations.
-    designated_syms = [species; parameters; variables_declared]
+    # Reads equations (and infers potential variables).
+    # Exlucdes any parameters already extracted (if they also was a variable).
+    designated_syms = [declared_syms; species_extracted]
     vars_extracted, add_default_diff, equations = read_equations_options(
         options, designated_syms)
     variables = vcat(variables_declared, vars_extracted)
+    parameters_extracted = setdiff(parameters_extracted, vars_extracted)
+
+    # Creates the finalised parameter and species lists.
+    species = vcat(species_declared, species_extracted)
+    parameters = vcat(parameters_declared, parameters_extracted)
 
     # Create differential expression.
     diffexpr = create_differential_expr(
@@ -723,6 +726,7 @@ function read_equations_options(options, syms_declared; requiredec = false)
         end
 
         # If the default differential (`D`) is used, record that it should be decalred later on.
+
         if !in(eq, excluded_syms) && find_D_call(eq)
             requiredec && throw(UndeclaredSymbolicError(
                 "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option."))

From 018fcf6693a04196d51b816c80c2445804eb6b50 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:28:38 +0000
Subject: [PATCH 07/17] another merrge fix

---
 src/dsl.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 2e9076a245..c7c4383f29 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -354,10 +354,10 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         reactions, declared_syms; requiredec)
 
     # Reads equations (and infers potential variables).
-    # Exlucdes any parameters already extracted (if they also was a variable).
-    designated_syms = [declared_syms; species_extracted]
+    # Excludes any parameters already extracted (if they also was a variable).
+    declared_syms = union(declared_syms, species_extracted)
     vars_extracted, add_default_diff, equations = read_equations_options(
-        options, designated_syms)
+        options, declared_syms)
     variables = vcat(variables_declared, vars_extracted)
     parameters_extracted = setdiff(parameters_extracted, vars_extracted)
 
@@ -731,7 +731,7 @@ function read_equations_options(options, syms_declared; requiredec = false)
             requiredec && throw(UndeclaredSymbolicError(
                 "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option."))
             add_default_diff = true
-            excluded_syms = [excluded_syms; :D]
+            excluded_syms = push!(excluded_syms, :D)
         end
 
         # Any undecalred symbolic variables encountered should be extracted as variables.

From f2ed0a35b53fc4f99c18df4c985aac0472ea12a4 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:33:57 +0000
Subject: [PATCH 08/17] imporve @require_declaration handling

---
 src/dsl.jl | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index c7c4383f29..501fba5377 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -735,13 +735,9 @@ function read_equations_options(options, syms_declared; requiredec = false)
         end
 
         # Any undecalred symbolic variables encountered should be extracted as variables.
-        # Additional step required to handle `requiredec = true` (to be improved later).
-        prev_vars_extracted = deepcopy(vars_extracted)
         add_syms_from_expr!(vars_extracted, eq, excluded_syms)
-        if requiredec && length(prev_vars_extracted) < length(vars_extracted)
-            throw(UndeclaredSymbolicError(
-                "Unrecognized symbols $(setdiff(vars_extracted, prev_vars_extracted)) was used in an equation: \"$eq\".  Since the flag @require_declaration is set, all variables must be declared with the @species, @parameters, or @variables macros."))
-        end
+        (!isempty(vars_extracted) && requiredec) && throw(UndeclaredSymbolicError(
+            "Unrecognized symbolic variables $(join(vars_extracted, ", ")) detected in equation expression: \"$(string(eq))\". Since the flag @require_declaration is declared, all symbolic variables must be explicitly declared with the @species, @variables, and @parameters options."))
     end
 
     return collect(vars_extracted), add_default_diff, equations

From e4c3e4e2ffcb429e48a4133ff9a6509273c1de2d Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:42:17 +0000
Subject: [PATCH 09/17] add @require_declaration test

---
 src/dsl.jl              |  2 +-
 test/dsl/dsl_options.jl | 32 ++++++++++++++++++++++----------
 2 files changed, 23 insertions(+), 11 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 501fba5377..1bc597bbb4 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -357,7 +357,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     # Excludes any parameters already extracted (if they also was a variable).
     declared_syms = union(declared_syms, species_extracted)
     vars_extracted, add_default_diff, equations = read_equations_options(
-        options, declared_syms)
+        options, declared_syms; requiredec)
     variables = vcat(variables_declared, vars_extracted)
     parameters_extracted = setdiff(parameters_extracted, vars_extracted)
 
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index f1e52586b3..7d6752b5c4 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -609,16 +609,16 @@ let
         @equations D(V) ~ 1 - V
         d, D --> 0
     end
-    @test_throws Exception @eval @reaction_network begin
-        @variables D(t)
-        @equations D(V) ~ 1 - V
-        d, X --> 0
-    end
-    @test_throws Exception @eval @reaction_network begin
-        @parameters D
-        @equations D(V) ~ 1 - V
-        d, X --> 0
-    end
+    @test_broken false # @test_throws Exception @eval @reaction_network begin
+        #@variables D(t)
+        #@equations D(V) ~ 1 - V
+        #d, X --> 0
+    #end
+    @test_broken false # @test_throws Exception @eval @reaction_network begin
+        #@parameters D
+        #@equations D(V) ~ 1 - V
+        #d, X --> 0
+    #end
 
     # Symbol only occurring in events.
     @test_throws Exception @eval @reaction_network begin
@@ -1323,4 +1323,16 @@ let
         @variables X1(t) X2(t)
         @observables X2 ~ X1
     end
+
+    # Test when the default differential D is inferred
+    @test_throws UndeclaredSymbolicError @macroexpand @reaction_network begin
+        @require_declaration
+        @variables V(t)
+        @equations D(V) ~ 1 - V
+    end
+    @test_nowarn @macroexpand @reaction_network begin
+        @differentials D = Differential(t)
+        @variables X1(t) X2(t)
+        @observables X2 ~ X1
+    end
 end

From 6e57098137e3fff65a9d6a5469b2d0d82f560b79 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:44:40 +0000
Subject: [PATCH 10/17] up

---
 src/dsl.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 1bc597bbb4..b0be2cde3c 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -718,7 +718,6 @@ function read_equations_options(options, syms_declared; requiredec = false)
     # When this is the case, the variable X and differential D are extracted (for automatic declaration).
     # Also performs simple error checks.
     vars_extracted = OrderedSet{Union{Symbol, Expr}}()
-    excluded_syms = syms_declared
     add_default_diff = false
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)

From 26b0fb2501e8f05aadd3a8ec7e96cdce6e16d48c Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:46:59 +0000
Subject: [PATCH 11/17] test improvement

---
 test/dsl/dsl_options.jl | 20 --------------------
 1 file changed, 20 deletions(-)

diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 7d6752b5c4..9846cd2149 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -1120,26 +1120,6 @@ let
     @test 5*sol[:Y][end] ≈ sol[:S][end] + sol[:X][end]
 end
 
-# Tests that the correct symbolic variables are inferred as species, variables, and parameters.
-let
-    rn = @reaction_network begin
-        @parameters p1 p2 k1 k2 k3 k4
-        @species X1(t) X2(t)
-        @variables W(t)
-        @equations begin
-            D(V1) ~ p1 * X1 - k1 * V1 + W
-            k2 * V2 ~ D(V2) + p2 * X2
-            V3 + X3 ~ V1^2 + X2^2
-        end
-        (k1, k2), X1 <--> X2
-        (k3, k4), X2 <--> X3
-    end
-
-    @test issetequal(species(rn), @species X1(t) X2(t) X3(t))
-    @test issetequal(parameters(rn), @parameters p1 p2 k1 k2 k3 k4)
-    @test issetequal(nonspecies(rn), @variables V1(t) V2(t) V3(t) W(t))
-end
-
 # Tests that various erroneous declarations throw errors.
 let
     # Using = instead of ~ (for equation).

From fc4d3876ae551690175f0a24a4d9bef2f3662369 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:48:17 +0000
Subject: [PATCH 12/17] more mrege fixing

---
 src/dsl.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index b0be2cde3c..6d7ed08666 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -730,11 +730,11 @@ function read_equations_options(options, syms_declared; requiredec = false)
             requiredec && throw(UndeclaredSymbolicError(
                 "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option."))
             add_default_diff = true
-            excluded_syms = push!(excluded_syms, :D)
+            push!(syms_declared, :D)
         end
 
         # Any undecalred symbolic variables encountered should be extracted as variables.
-        add_syms_from_expr!(vars_extracted, eq, excluded_syms)
+        add_syms_from_expr!(vars_extracted, eq, syms_declared)
         (!isempty(vars_extracted) && requiredec) && throw(UndeclaredSymbolicError(
             "Unrecognized symbolic variables $(join(vars_extracted, ", ")) detected in equation expression: \"$(string(eq))\". Since the flag @require_declaration is declared, all symbolic variables must be explicitly declared with the @species, @variables, and @parameters options."))
     end

From 734756d23d84b98b4b00aad1309374e3b673db6b Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 21:58:05 +0000
Subject: [PATCH 13/17] minor bugfix

---
 src/dsl.jl       | 4 ++--
 test/runtests.jl | 7 -------
 2 files changed, 2 insertions(+), 9 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 6d7ed08666..5d3194d6e8 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -371,7 +371,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
 
     # Reads observables.
     observed_vars, observed_eqs, obs_syms = read_observed_options(
-        options, [species_declared; variables], all_ivs)
+        options, [species_declared; variables], all_ivs; requiredec)
 
     # Checks for input errors.
     (sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
@@ -726,7 +726,7 @@ function read_equations_options(options, syms_declared; requiredec = false)
 
         # If the default differential (`D`) is used, record that it should be decalred later on.
 
-        if !in(eq, excluded_syms) && find_D_call(eq)
+        if !in(eq, syms_declared) && find_D_call(eq)
             requiredec && throw(UndeclaredSymbolicError(
                 "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option."))
             add_default_diff = true
diff --git a/test/runtests.jl b/test/runtests.jl
index 29c2d0ca74..188bced121 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,13 +15,6 @@ end
 ### Run Tests ###
 @time begin
     if GROUP == "All" || GROUP == "Core"
-        # Tests the `ReactionSystem` structure and its properties.
-        @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end
-        @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end
-        @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
-        @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
-        @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
-        @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
         @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
         @time @safetestset "Events" begin include("reactionsystem_core/events.jl") end
 

From 6fbaff4cd62df9ac14056769ef8e0ede7e14f8eb Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 14 Jan 2025 23:32:13 +0100
Subject: [PATCH 14/17] Update dsl_options.jl

---
 test/dsl/dsl_options.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 9846cd2149..f3eef8c70f 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -1284,12 +1284,12 @@ let
     # Test error when a variable in an equation is inferred
     @test_throws UndeclaredSymbolicError @macroexpand @reaction_network begin
         @require_declaration
-        @equations D(V) ~ V^2
+        @equations V ~ V^2 + 2
     end
     @test_nowarn @macroexpand @reaction_network begin
         @require_declaration
         @variables V(t)
-        @equations D(V) ~ V^2
+        @equations V ~ V^2 + 2
     end
 
     # Test error when a variable in an observable is inferred

From 33858075ad5cca6cba1ed5014fa33d8dede5a000 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 15 Jan 2025 10:37:35 +0000
Subject: [PATCH 15/17] test fix

---
 test/runtests.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 188bced121..29c2d0ca74 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,13 @@ end
 ### Run Tests ###
 @time begin
     if GROUP == "All" || GROUP == "Core"
+        # Tests the `ReactionSystem` structure and its properties.
+        @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end
+        @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end
+        @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
+        @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
+        @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
+        @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
         @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
         @time @safetestset "Events" begin include("reactionsystem_core/events.jl") end
 

From 6f73003719301fc97ae2e3f45cf7813c2c267b7a Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Thu, 16 Jan 2025 14:20:56 +0000
Subject: [PATCH 16/17] up

---
 HISTORY.md              |  2 +-
 src/dsl.jl              |  7 +++----
 test/dsl/dsl_options.jl | 10 +++++-----
 3 files changed, 9 insertions(+), 10 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index d6d381585b..3214f9e566 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -29,7 +29,7 @@ end
   @equations D(D(V)) ~ 1
 end
 ```
-Please note that this cannot be used at the same time as `D` is used to represent a species, variable, or parameter. 
+Please note that this cannot be used at the same time as `D` is used to represent a species, variable, or parameter (including is these are implicitly designated as such by e.g. appearing as a reaction reactant). 
 - Array symbolics support is more consistent with ModelingToolkit v9. Parameter
   arrays are no longer scalarized by Catalyst, while species and variables
   arrays still are (as in ModelingToolkit). As such, parameter arrays should now
diff --git a/src/dsl.jl b/src/dsl.jl
index 5d3194d6e8..ac1bb88de9 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -357,7 +357,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     # Excludes any parameters already extracted (if they also was a variable).
     declared_syms = union(declared_syms, species_extracted)
     vars_extracted, add_default_diff, equations = read_equations_options(
-        options, declared_syms; requiredec)
+        options, declared_syms, parameters_extracted; requiredec)
     variables = vcat(variables_declared, vars_extracted)
     parameters_extracted = setdiff(parameters_extracted, vars_extracted)
 
@@ -705,7 +705,7 @@ end
 # `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
 # `dtexpr`: If a differential equation is defined, the default derivative (D ~ Differential(t)) must be defined.
 # `equations`: a vector with the equations provided.
-function read_equations_options(options, syms_declared; requiredec = false)
+function read_equations_options(options, syms_declared, parameters_extracted; requiredec = false)
     # Prepares the equations. First, extracts equations from provided option (converting to block form if required).
     # 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] : :(begin end)
@@ -725,8 +725,7 @@ function read_equations_options(options, syms_declared; requiredec = false)
         end
 
         # If the default differential (`D`) is used, record that it should be decalred later on.
-
-        if !in(eq, syms_declared) && find_D_call(eq)
+        if (:D ∉ union(syms_declared, parameters_extracted)) && find_D_call(eq)
             requiredec && throw(UndeclaredSymbolicError(
                 "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option."))
             add_default_diff = true
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index f3eef8c70f..16ef93d1f7 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -614,11 +614,11 @@ let
         #@equations D(V) ~ 1 - V
         #d, X --> 0
     #end
-    @test_broken false # @test_throws Exception @eval @reaction_network begin
-        #@parameters D
-        #@equations D(V) ~ 1 - V
-        #d, X --> 0
-    #end
+    @test_throws Exception @eval @reaction_network begin
+        @parameters D
+        @equations D(V) ~ 1 - V
+        d, X --> 0
+    end
 
     # Symbol only occurring in events.
     @test_throws Exception @eval @reaction_network begin

From 64f1e395778c2818b6c9e9cfb133d09fbbff497e Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 17 Jan 2025 13:45:15 +0000
Subject: [PATCH 17/17] up

---
 docs/src/inverse_problems/examples/ode_fitting_oscillation.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
index ff3fb63805..a72ebd6b44 100644
--- a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
+++ b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
@@ -56,7 +56,7 @@ function optimise_p(pinit, tend)
         newprob = remake(prob; tspan = (0.0, tend), p = p)
         sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes))
         loss = sum(abs2, sol .- sample_vals[:, 1:size(sol,2)])
-        return loss, sol
+        return loss
     end
 
     # optimize for the parameters that minimize the loss