diff --git a/HISTORY.md b/HISTORY.md
index 2039c44398..aba2af7841 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -25,6 +25,22 @@
   StructuralIdentifiability has with Julia 1.10.5 and 1.11.
 - A tutorial on making interactive plot displays using Makie has been added.
 - The BifurcationKit extension has been updated to v.4.
+- There is a new DSL option `@require_declaration` that will turn off automatic inferring for species, parameters, and variables in the DSL. For example, the following will now error:
+  ```julia
+  rn = @reaction_network begin
+    @require_declaration
+    (k1, k2), A <--> B
+  end
+  ```
+  When this flag is set, all symbolics must be explicitly declared. 
+  ```julia
+  rn = @reaction_network begin
+    @species A(t) B(t)
+    @parameters k1 k2
+    @require_declaration
+    (k1, k2), A <--> B
+  end
+  ```
 
 ## Catalyst 14.4.1
 - Support for user-defined functions on the RHS when providing coupled equations 
diff --git a/Project.toml b/Project.toml
index cb7b3a5395..3bb0b0c1f1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -66,7 +66,7 @@ SciMLBase = "2.57.2"
 Setfield = "1"
 # StructuralIdentifiability = "0.5.8"
 SymbolicUtils = "2.1.2, 3.3.0"
-Symbolics = "5.30.1, 6"
+Symbolics = "6.22"
 Unitful = "1.12.4"
 julia = "1.10"
 
@@ -77,7 +77,11 @@ Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
+OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -93,7 +97,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["DiffEqCallbacks", "DomainSets", "Graphviz_jll", "Logging", "NonlinearSolve", 
-        "OrdinaryDiffEq", "Pkg", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", 
-        "StableRNGs", "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", 
-        "Test", "Unitful"]
+test = ["DiffEqCallbacks", "DomainSets", "Graphviz_jll", "Logging", "NonlinearSolve", "OrdinaryDiffEqBDF", "OrdinaryDiffEqDefault", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Pkg", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 40471c42c9..086dcf6e38 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -87,4 +87,4 @@ SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 SteadyStateDiffEq = "2.2"
 StochasticDiffEq = "6.65"
-Symbolics = "5.30.1, 6"
+Symbolics = "6.22"
diff --git a/docs/src/faqs.md b/docs/src/faqs.md
index 095031c6c2..25f09ca330 100644
--- a/docs/src/faqs.md
+++ b/docs/src/faqs.md
@@ -309,3 +309,16 @@ end
 In some cases, it may be necessary or desirable to register functions with
 Symbolics.jl before their use in Catalyst, see the discussion
 [here](https://symbolics.juliasymbolics.org/stable/manual/functions/).
+
+## How can I turn off automatic inferring of species and parameters when using the DSL?
+This option can be set using the `@require_declaration` option inside `@reaction_network`. In this case all the species, parameters, and variables in the system must be pre-declared using one of the `@species`, `@parameters`, or `@variables` macros. For more information about what is inferred automatically and not, please see the section on [`@require_declaration`](@ref dsl_advanced_options_require_dec).
+
+```@example faq9
+using Catalyst
+rn = @reaction_network begin
+    @require_declaration
+    @species A(t) B(t)
+    @parameters k1 k2
+    (k1, k2), A <--> B
+end
+```
diff --git a/docs/src/model_creation/dsl_advanced.md b/docs/src/model_creation/dsl_advanced.md
index c59c4ba21e..af806864b6 100644
--- a/docs/src/model_creation/dsl_advanced.md
+++ b/docs/src/model_creation/dsl_advanced.md
@@ -205,7 +205,7 @@ ModelingToolkit.getdescription(two_state_system.kA)
 
 ### [Designating constant-valued/fixed species parameters](@id dsl_advanced_options_constant_species)
 
-Catalyst enables the designation of parameters as `constantspecies`. These parameters can be used as species in reactions, however, their values are not changed by the reaction and remain constant throughout the simulation (unless changed by e.g. the [occurrence of an event]@ref constraint_equations_events). Practically, this is done by setting the parameter's `isconstantspecies` metadata to `true`. Here, we create a simple reaction where the species `X` is converted to `Xᴾ` at rate `k`. By designating `X` as a constant species parameter, we ensure that its quantity is unchanged by the occurrence of the reaction.
+Catalyst enables the designation of parameters as `constantspecies`. These parameters can be used as species in reactions, however, their values are not changed by the reaction and remain constant throughout the simulation (unless changed by e.g. the [occurrence of an event](@ref constraint_equations_events). Practically, this is done by setting the parameter's `isconstantspecies` metadata to `true`. Here, we create a simple reaction where the species `X` is converted to `Xᴾ` at rate `k`. By designating `X` as a constant species parameter, we ensure that its quantity is unchanged by the occurrence of the reaction.
 ```@example dsl_advanced_constant_species
 using Catalyst # hide
 rn = @reaction_network begin
@@ -272,6 +272,40 @@ sol = solve(oprob)
 plot(sol)
 ```
 
+### [Turning off species, parameter, and variable inference](@id dsl_advanced_options_require_dec)
+In some cases it may be desirable for Catalyst to not infer species and parameters from the DSL, as in the case of reaction networks with very many variables, or as a sanity check that variable names are written correctly. To turn off inference, simply use the `@require_declaration` option when using the `@reaction_network` DSL. This will require every single variable, species, or parameter used within the DSL to be explicitly declared using the `@variable`, `@species`, or `@parameter` options. In the case that the DSL parser encounters an undeclared symbolic, it will error with an `UndeclaredSymbolicError` and print the reaction or equation that the undeclared symbolic was found in. 
+
+```julia
+using Catalyst
+rn = @reaction_network begin
+    @require_declaration
+    (k1, k2), A <--> B
+end
+```
+Running the code above will yield the following error: 
+```
+LoadError: UndeclaredSymbolicError: Unrecognized variables A detected in reaction expression: "((k1, k2), A <--> B)". Since the flag @require_declaration is declared, all species must be explicitly declared with the @species macro.
+```
+In order to avoid the error in this case all the relevant species and parameters will have to be declared.
+```@example dsl_advanced_require_dec
+# The following case will not error. 
+using Catalyst
+t = default_t()
+rn = @reaction_network begin
+    @require_declaration
+    @species A(t) B(t)
+    @parameters k1 k2
+    (k1, k2), A <--> B
+end
+```
+
+The following cases in which the DSL would normally infer variables will all throw errors if `@require_declaration` is set and the variables are not explicitly declared.
+- Occurrence of an undeclared species in a reaction, as in the example above.
+- Occurrence of an undeclared parameter in a reaction rate expression, as in the reaction line `k*n, A --> B`.
+- Occurrence of an undeclared parameter in the stoichiometry of a species, as in the reaction line `k, n*A --> B`.
+- Occurrence of an undeclared differential variable on the LHS of a coupled differential equation, as in `A` in `@equations D(A) ~ A^2`.
+- Occurrence of an undeclared [observable](@ref dsl_advanced_options_observables) in an `@observables` expression, such as `@observables X1 ~ A + B`.
+
 ## [Naming reaction networks](@id dsl_advanced_options_naming)
 Each reaction network model has a name. It can be accessed using the `nameof` function. By default, some generic name is used:
 ```@example dsl_advanced_names
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index ad12f0fa4b..73fbc7586f 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -41,6 +41,7 @@ import ModelingToolkit: check_variables,
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, Graphs
+using MacroTools: striplines
 import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne
 import DataStructures: OrderedDict, OrderedSet
 import Parameters: @with_kw_noshow
diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl
index 1dd993ff16..7fbf43e43b 100644
--- a/src/chemistry_functionality.jl
+++ b/src/chemistry_functionality.jl
@@ -63,11 +63,11 @@ t = default_t()
 @compound CO2(t) ~ C + 2O
 ```
 
-Notes: 
+Notes:
 - The component species must be defined before using the `@compound` macro.
 """
 macro compound(expr)
-    make_compound(MacroTools.striplines(expr))
+    make_compound(striplines(expr))
 end
 
 # Declares compound error messages:
@@ -83,12 +83,12 @@ function make_compound(expr)
         error(COMPOUND_CREATION_ERROR_BAD_SEPARATOR)
 
     # Loops through all components, add the component and the coefficients to the corresponding vectors
-    # Cannot extract directly using e.g. "getfield.(composition, :reactant)" because then 
+    # Cannot extract directly using e.g. "getfield.(composition, :reactant)" because then
     # we get something like :([:C, :O]), rather than :([C, O]).
     composition = Catalyst.recursive_find_reactants!(expr.args[3], 1,
         Vector{ReactantStruct}(undef, 0))
-    components = :([])                                      # Becomes something like :([C, O]).                                         
-    coefficients = :([])                                    # Becomes something like :([1, 2]). 
+    components = :([])                                      # Becomes something like :([C, O]).
+    coefficients = :([])                                    # Becomes something like :([1, 2]).
     for comp in composition
         push!(components.args, comp.reactant)
         push!(coefficients.args, comp.stoichiometry)
@@ -110,13 +110,13 @@ function make_compound(expr)
                for comp in $components])))
 
     # Creates the found expressions that will create the compound species.
-    # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in 
-    # the scope the users use the macro in (to e.g. detect already exiting species).     
+    # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in
+    # the scope the users use the macro in (to e.g. detect already exiting species).
     # Creates something like (where `compound_ivs` and `component_ivs` evaluates to all the compound's and components' ivs):
     #   `@species CO2(..)`
     #   `isempty([])` && length(component_ivs) && error("When ...)
-    #   `CO2 = CO2(component_ivs..)` 
-    #   `issetequal(compound_ivs, component_ivs) || error("The ...)` 
+    #   `CO2 = CO2(component_ivs..)`
+    #   `issetequal(compound_ivs, component_ivs) || error("The ...)`
     #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)`
     #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])`
     #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])`
@@ -159,7 +159,7 @@ Macro that creates several compound species, which each is composed of smaller c
 Example:
 ```julia
 t = default_t()
-@species C(t) H(t) O(t) 
+@species C(t) H(t) O(t)
 @compounds
     CH4(t) = C + 4H
     O2(t) = 2O
@@ -168,11 +168,11 @@ t = default_t()
 end
 ```
 
-Notes: 
+Notes:
 - The component species must be defined before using the `@compound` macro.
 """
 macro compounds(expr)
-    make_compounds(MacroTools.striplines(expr))
+    make_compounds(striplines(expr))
 end
 
 # Function managing the @compound macro.
@@ -183,7 +183,7 @@ function make_compounds(expr)
     # For each compound in `expr`, creates the set of 7 compound creation lines (using `make_compound`).
     # Next, loops through all 7*[Number of compounds] lines and add them to compound_declarations.
     compound_calls = [Catalyst.make_compound(line) for line in expr.args]
-    for compound_call in compound_calls, line in MacroTools.striplines(compound_call).args
+    for compound_call in compound_calls, line in striplines(compound_call).args
         push!(compound_declarations.args, line)
     end
 
@@ -249,7 +249,7 @@ brxs = balance_reaction(rx) # No solution.
 
 Notes:
 - Balancing reactions that contain compounds of compounds is currently not supported.
-- A reaction may not always yield a single solution; it could have an infinite number of solutions or none at all. When there are multiple solutions, a vector of all possible `Reaction` objects is returned. However, substrates and products may be interchanged as we currently do not solve for a linear combination that maintains the set of substrates and products. 
+- A reaction may not always yield a single solution; it could have an infinite number of solutions or none at all. When there are multiple solutions, a vector of all possible `Reaction` objects is returned. However, substrates and products may be interchanged as we currently do not solve for a linear combination that maintains the set of substrates and products.
 - If the reaction cannot be balanced, an empty `Reaction` vector is returned.
 """
 function balance_reaction(reaction::Reaction)
@@ -369,16 +369,16 @@ From a system, creates a new system where each reaction is a balanced version of
 reaction of the original system. For more information, consider the `balance_reaction` function
 (which is internally applied to each system reaction).
 
-Arguments 
+Arguments
 - `rs`: The reaction system that should be balanced.
 
 Notes:
 - If any reaction in the system cannot be balanced, throws an error.
-- If any reaction in the system have an infinite number of potential reactions, throws an error. 
+- If any reaction in the system have an infinite number of potential reactions, throws an error.
 Here, it would be possible to generate a valid reaction, however, no such routine is currently
 implemented in `balance_system`.
 - `balance_system` will not modify reactions of subsystems to the input system. It is recommended
-not to apply `balance_system` to non-flattened systems. 
+not to apply `balance_system` to non-flattened systems.
 """
 function balance_system(rs::ReactionSystem)
     @set! rs.eqs = CatalystEqType[get_balanced_reaction(eq) for eq in get_eqs(rs)]
@@ -391,7 +391,7 @@ end
 function get_balanced_reaction(rx::Reaction)
     brxs = balance_reaction(rx)
 
-    # In case there are no, or multiple, solutions to the balancing problem. 
+    # In case there are no, or multiple, solutions to the balancing problem.
     if isempty(brxs)
         error("Could not balance reaction `$rx`, unable to create a balanced `ReactionSystem`.")
     end
diff --git a/src/dsl.jl b/src/dsl.jl
index f17c17522c..a2bf01ec41 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -141,17 +141,17 @@ ReactionSystems generated through `@reaction_network` are complete.
 """
 macro reaction_network(name::Symbol, ex::Expr)
     :(complete($(make_reaction_system(
-        MacroTools.striplines(ex); name = :($(QuoteNode(name)))))))
+        striplines(ex); name = :($(QuoteNode(name)))))))
 end
 
 # Allows @reaction_network $name begin ... to interpolate variables storing a name.
 macro reaction_network(name::Expr, ex::Expr)
     :(complete($(make_reaction_system(
-        MacroTools.striplines(ex); name = :($(esc(name.args[1])))))))
+        striplines(ex); name = :($(esc(name.args[1])))))))
 end
 
 macro reaction_network(ex::Expr)
-    ex = MacroTools.striplines(ex)
+    ex = striplines(ex)
 
     # no name but equations: @reaction_network begin ... end ...
     if ex.head == :block
@@ -179,16 +179,16 @@ Equivalent to `@reaction_network` except the generated `ReactionSystem` is not m
 complete.
 """
 macro network_component(name::Symbol, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name))))
+    make_reaction_system(striplines(ex); name = :($(QuoteNode(name))))
 end
 
 # Allows @network_component $name begin ... to interpolate variables storing a name.
 macro network_component(name::Expr, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(esc(name.args[1]))))
+    make_reaction_system(striplines(ex); name = :($(esc(name.args[1]))))
 end
 
 macro network_component(ex::Expr)
-    ex = MacroTools.striplines(ex)
+    ex = striplines(ex)
 
     # no name but equations: @network_component begin ... end ...
     if ex.head == :block
@@ -328,7 +328,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads equations. 
+    # Reads equations.
     vars_extracted, add_default_diff, equations = read_equations_options(
         options, variables_declared; requiredec)
     variables = vcat(variables_declared, vars_extracted)
@@ -356,7 +356,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     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. 
+    # 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!(
@@ -397,7 +397,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
@@ -448,7 +448,7 @@ function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(u
             push_reactions!(reactions, reaction.args[3], reaction.args[2],
                 rate, metadata, arrow, line)
         else
-            throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))")
+            throw("Malformed reaction, invalid arrow type used in: $(striplines(line))")
         end
     end
     return reactions
@@ -670,7 +670,7 @@ function read_events_option(options, event_type::Symbol)
         error("Trying to read an unsupported event type.")
     end
     events_input = haskey(options, event_type) ? options[event_type].args[3] :
-                   MacroTools.striplines(:(begin end))
+                   striplines(:(begin end))
     events_input = option_block_form(events_input)
 
     # Goes through the events, checks for errors, and adds them to the output vector.
@@ -750,7 +750,7 @@ function create_differential_expr(options, add_default_diff, used_syms, tiv)
     # 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].args[3] :
-                MacroTools.striplines(:(begin end)))
+                striplines(:(begin end)))
     diffexpr = option_block_form(diffexpr)
 
     # Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere.
@@ -810,14 +810,8 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted; req
             # For Observables that have already been declared using @species/@variables,
             # or are interpolated, this parts should not be carried out.
             if !((obs_name in species_n_vars_declared) || is_escaped_expr(obs_eq.args[2]))
-                # Appends (..) to the observable (which is later replaced with the extracted ivs).
-                # Adds the observable to the first line of the output expression (starting with `@variables`).
-                obs_expr = insert_independent_variable(obs_eq.args[2], :(..))
-                push!(observed_vars.args[1].args, obs_expr)
-
-                # Adds a line to the `observed_vars` expression, setting the ivs for this observable.
-                # Cannot extract directly using e.g. "getfield.(dependants_structs, :reactant)" because
-                # then we get something like :([:X1, :X2]), rather than :([X1, X2]).
+                # Creates an expression which extracts the ivs of the species & variables the
+                # observable depends on, and splats them out in the correct order.
                 dep_var_expr = :(filter(!MT.isparameter,
                     Symbolics.get_variables($(obs_eq.args[3]))))
                 ivs_get_expr = :(unique(reduce(
@@ -825,8 +819,9 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted; req
                            for dep in $dep_var_expr])))
                 ivs_get_expr_sorted = :(sort($(ivs_get_expr);
                     by = iv -> findfirst(MT.getname(iv) == ivs for ivs in $ivs_sorted)))
-                push!(observed_vars.args,
-                    :($obs_name = $(obs_name)($(ivs_get_expr_sorted)...)))
+
+                obs_expr = insert_independent_variable(obs_eq.args[2], :($ivs_get_expr_sorted...))
+                push!(observed_vars.args[1].args, obs_expr)
             end
 
             # In case metadata was given, this must be cleared from `observed_eqs`.
@@ -840,7 +835,8 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted; req
         end
 
         # If nothing was added to `observed_vars`, it has to be modified not to throw an error.
-        (length(observed_vars.args) == 1) && (observed_vars = :())
+        (striplines(observed_vars) == striplines(Expr(:block, :(@variables)))) &&
+            (observed_vars = :())
     else
         # If option is not used, return empty expression and vector.
         observed_vars = :()
@@ -944,7 +940,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/spatial_reaction_systems/spatial_reactions.jl b/src/spatial_reaction_systems/spatial_reactions.jl
index 18c2a8d8e8..3b846d3049 100644
--- a/src/spatial_reaction_systems/spatial_reactions.jl
+++ b/src/spatial_reaction_systems/spatial_reactions.jl
@@ -43,7 +43,7 @@ end
 
 # Macro for creating a `TransportReaction`.
 macro transport_reaction(rateex::ExprValues, species::ExprValues)
-    make_transport_reaction(MacroTools.striplines(rateex), species)
+    make_transport_reaction(striplines(rateex), species)
 end
 function make_transport_reaction(rateex, species)
     # Handle interpolation of variables
diff --git a/test/compositional_modelling/component_based_model_creation.jl b/test/compositional_modelling/component_based_model_creation.jl
index 5efd20967f..bcbcf393db 100644
--- a/test/compositional_modelling/component_based_model_creation.jl
+++ b/test/compositional_modelling/component_based_model_creation.jl
@@ -3,7 +3,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, LinearAlgebra, OrdinaryDiffEq, SciMLNLSolve, Test
+using Catalyst, LinearAlgebra, OrdinaryDiffEqTsit5, SciMLNLSolve, Test
 using ModelingToolkit: nameof
 
 # Sets the default `t` to use.
@@ -543,4 +543,4 @@ let
     @test length(unknowns(sys3)) == 4
     @test length(parameters(sys4)) == 5
     @test any(isequal(p5), parameters(sys4))
-end
\ No newline at end of file
+end
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 872ba4f0b2..e07f0cdbc3 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -3,7 +3,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, ModelingToolkit, OrdinaryDiffEq, StochasticDiffEq, Plots, Test
+using Catalyst, ModelingToolkit, OrdinaryDiffEqTsit5, OrdinaryDiffEqRosenbrock, StochasticDiffEq, Plots, Test
 using Symbolics: unwrap
 
 # Sets stable rng number.
@@ -565,6 +565,45 @@ 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
 
+# Checks that models created w/o specifying `@variables` for observables are identical.
+# Compares both to model with explicit declaration, and programmatically created model.
+let
+    # With default ivs.
+    rn1 = @reaction_network rn begin
+        @variables X(t) X1(t) X2(t)
+        @observables X ~ X1 + X2
+    end
+    rn2 = @reaction_network rn begin
+        @variables X1(t) X2(t)
+        @observables X ~ X1 + X2
+    end
+    @variables X(t) X1(t) X2(t)
+    rn3 = complete(ReactionSystem([], t, [X1, X2], []; name = :rn, observed = [X ~ X1 + X2]))
+    @test isequal(rn1, rn2)
+    @test isequal(rn1, rn3)
+    @test isequal(rn1.X, rn2.X)
+    @test isequal(rn1.X, rn3.X)
+
+    # With non-default ivs.
+    rn4 = @reaction_network rn begin
+        @ivs τ x
+        @variables X(τ,x) X1(τ,x) X2(τ,x)
+        @observables X ~ X1 + X2
+    end
+    rn5 = @reaction_network rn begin
+        @ivs τ x
+        @variables X1(τ,x) X2(τ,x)
+        @observables X ~ X1 + X2
+    end
+    @parameters τ x
+    @variables X(τ,x) X1(τ,x) X2(τ,x)
+    rn6 = complete(ReactionSystem([], τ, [X1, X2], []; name = :rn, observed = [X ~ X1 + X2], spatial_ivs = [x]))
+    @test isequal(rn4, rn5)
+    @test isequal(rn4, rn6)
+    @test isequal(rn4.X, rn5.X)
+    @test isequal(rn4.X, rn6.X)
+end
+
 # Checks that ivs are correctly found.
 let
     rn = @reaction_network begin
@@ -604,7 +643,7 @@ let
         k, 0 --> X1 + X2
     end
     @test isequal(observed(rn1)[1].rhs, observed(rn2)[1].rhs)
-    @test_broken isequal(observed(rn1)[1].lhs.metadata, observed(rn2)[1].lhs.metadata)
+    @test isequal(observed(rn1)[1].lhs.metadata, observed(rn2)[1].lhs.metadata)
     @test isequal(unknowns(rn1), unknowns(rn2))
 
     # Case with metadata.
@@ -618,7 +657,7 @@ let
         k, 0 --> X1 + X2
     end
     @test isequal(observed(rn3)[1].rhs, observed(rn4)[1].rhs)
-    @test_broken isequal(observed(rn3)[1].lhs.metadata, observed(rn4)[1].lhs.metadata)
+    @test isequal(observed(rn3)[1].lhs.metadata, observed(rn4)[1].lhs.metadata)
     @test isequal(unknowns(rn3), unknowns(rn4))
 end
 
@@ -951,7 +990,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 +1004,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 +1017,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 +1048,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
@@ -1030,7 +1069,7 @@ let
     # Test error when species are inferred
     @test_throws UndeclaredSymbolicError @macroexpand @reaction_network begin
         @require_declaration
-        @parameters k 
+        @parameters k
         k, A --> B
     end
     @test_nowarn @macroexpand @reaction_network begin
@@ -1091,4 +1130,3 @@ let
         @observables X2 ~ X1
     end
 end
-
diff --git a/test/extensions/lattice_simulation_plotting.jl b/test/extensions/lattice_simulation_plotting.jl
index 0092916748..808acc877f 100644
--- a/test/extensions/lattice_simulation_plotting.jl
+++ b/test/extensions/lattice_simulation_plotting.jl
@@ -1,7 +1,7 @@
 ### Preparations ###
 
 # Fetch packages.
-using Catalyst, CairoMakie, GraphMakie, Graphs, JumpProcesses, OrdinaryDiffEq, Test
+using Catalyst, CairoMakie, GraphMakie, Graphs, JumpProcesses, OrdinaryDiffEqDefault, OrdinaryDiffEq, Test
 
 
 ### Checks Basic Plot Cases ###
@@ -130,4 +130,4 @@ let
         @test isfile("animation_tmp.mp4")
         rm("animation_tmp.mp4")
     end
-end
\ No newline at end of file
+end
diff --git a/test/extensions/stability_computation.jl b/test/extensions/stability_computation.jl
index 24ddfe7a29..eeeba72131 100644
--- a/test/extensions/stability_computation.jl
+++ b/test/extensions/stability_computation.jl
@@ -1,7 +1,7 @@
 ### Fetch Packages ###
 
 # Fetch packages.
-using Catalyst, OrdinaryDiffEq, SteadyStateDiffEq, Test
+using Catalyst, OrdinaryDiffEqVerner, SteadyStateDiffEq, Test
 import HomotopyContinuation
 
 # Sets rnd number.
@@ -105,4 +105,4 @@ let
     p = [:p => 1.0, :d => 1.0]
     u = [1.0]
     @test_throws Exception steady_state_stability(u, rn, p)
-end
\ No newline at end of file
+end
diff --git a/test/miscellaneous_tests/api.jl b/test/miscellaneous_tests/api.jl
index 9f2cddbb20..85b51ff43b 100644
--- a/test/miscellaneous_tests/api.jl
+++ b/test/miscellaneous_tests/api.jl
@@ -3,7 +3,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, NonlinearSolve, OrdinaryDiffEq, SparseArrays, StochasticDiffEq, Test
+using Catalyst, NonlinearSolve, OrdinaryDiffEqTsit5, SparseArrays, StochasticDiffEq, Test
 using LinearAlgebra: norm
 using ModelingToolkit: value
 
@@ -505,4 +505,4 @@ let
         f(d,t), X --> 0
     end
     @test !isautonomous(rn8)
-end
\ No newline at end of file
+end
diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index fe79920a02..31736d8fa2 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, JumpProcesses, LinearAlgebra, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq, StochasticDiffEq, Test
+using Catalyst, JumpProcesses, LinearAlgebra, NonlinearSolve, OrdinaryDiffEqTsit5, OrdinaryDiffEqDefault, OrdinaryDiffEqVerner, SteadyStateDiffEq, StochasticDiffEq, Test
 
 # Sets stable rng number.
 using StableRNGs
diff --git a/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl b/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl
index 804dc4bc05..ca044d87e1 100644
--- a/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl
+++ b/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl
@@ -3,7 +3,7 @@
 ### Preparations ###
 
 # Fetch packages.
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 using Random, Statistics, SparseArrays, Test
 
 # Fetch test networks.
diff --git a/test/reactionsystem_core/coupled_equation_crn_systems.jl b/test/reactionsystem_core/coupled_equation_crn_systems.jl
index d6c8bd5dc6..9968eb05f3 100644
--- a/test/reactionsystem_core/coupled_equation_crn_systems.jl
+++ b/test/reactionsystem_core/coupled_equation_crn_systems.jl
@@ -1,5 +1,5 @@
 # Fetch packages.
-using Catalyst, NonlinearSolve, OrdinaryDiffEq, Statistics, SteadyStateDiffEq, StochasticDiffEq, Test
+using Catalyst, NonlinearSolve, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqRosenbrock, Statistics, SteadyStateDiffEq, StochasticDiffEq, Test
 using ModelingToolkit: getdefault, getdescription, getdefault
 using Symbolics: BasicSymbolic, unwrap
 
diff --git a/test/reactionsystem_core/events.jl b/test/reactionsystem_core/events.jl
index 5264114d83..f98d77c6e9 100644
--- a/test/reactionsystem_core/events.jl
+++ b/test/reactionsystem_core/events.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, DiffEqCallbacks, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Test
+using Catalyst, DiffEqCallbacks, JumpProcesses, OrdinaryDiffEqTsit5, StochasticDiffEq, Test
 
 # Sets stable rng number.
 using StableRNGs
diff --git a/test/reactionsystem_core/parameter_type_designation.jl b/test/reactionsystem_core/parameter_type_designation.jl
index 2a95bfa312..9b535e352c 100644
--- a/test/reactionsystem_core/parameter_type_designation.jl
+++ b/test/reactionsystem_core/parameter_type_designation.jl
@@ -1,7 +1,7 @@
 ### Fetch Packages and Set Global Variables ###
 
 # Fetch packages.
-using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, StochasticDiffEq, Test
+using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEqTsit5, StochasticDiffEq, Test
 using Symbolics: BasicSymbolic, unwrap
 
 # Sets stable rng number.
diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl
index 3acdf008e3..505781305a 100644
--- a/test/reactionsystem_core/reactionsystem.jl
+++ b/test/reactionsystem_core/reactionsystem.jl
@@ -1,7 +1,7 @@
 ### Fetch Packages and Set Global Variables ###
 
 # Fetch packages.
-using Catalyst, LinearAlgebra, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Test
+using Catalyst, LinearAlgebra, JumpProcesses, OrdinaryDiffEqTsit5, OrdinaryDiffEqVerner, StochasticDiffEq, Test
 const MT = ModelingToolkit
 
 # Sets stable rng number.
@@ -1021,4 +1021,4 @@ let
     jprob = JumpProblem(jinput; rng)
     sol = solve(jprob)
     @test sol(10.0; idxs = :A) > 0
-end
\ No newline at end of file
+end
diff --git a/test/reactionsystem_core/symbolic_stoichiometry.jl b/test/reactionsystem_core/symbolic_stoichiometry.jl
index c2ec784919..2b51bba630 100644
--- a/test/reactionsystem_core/symbolic_stoichiometry.jl
+++ b/test/reactionsystem_core/symbolic_stoichiometry.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Statistics, Test
+using Catalyst, JumpProcesses, OrdinaryDiffEqTsit5, StochasticDiffEq, Statistics, Test
 using Symbolics: BasicSymbolic, unwrap
 
 # Sets stable rng number.
@@ -297,4 +297,4 @@ let
     end_vals_ref = [[sol[s][end] for sol in sols_ref.u] for s in species(sir_ref)]
     @test mean.(end_vals_ref) ≈ mean.(end_vals) atol=1e-2 rtol = 1e-2
     @test var.(end_vals_ref) ≈ var.(end_vals) atol=1e-1 rtol = 1e-1
-end
\ No newline at end of file
+end
diff --git a/test/simulation_and_solving/simulate_ODEs.jl b/test/simulation_and_solving/simulate_ODEs.jl
index ecb2a468ce..3bff83eae5 100644
--- a/test/simulation_and_solving/simulate_ODEs.jl
+++ b/test/simulation_and_solving/simulate_ODEs.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, OrdinaryDiffEq, Test
+using Catalyst, OrdinaryDiffEqVerner, OrdinaryDiffEqRosenbrock, Test
 
 # Sets stable rng number.
 using StableRNGs
@@ -200,4 +200,4 @@ let
     oprob.f(du1, oprob.u0, oprob.p, 90.0)
     oderhs(du2, u_2, p_2, 0.0)
     @test du1 ≈ du2
-end
\ No newline at end of file
+end
diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl
index a5d91d758b..c0b19a81da 100644
--- a/test/simulation_and_solving/simulate_jumps.jl
+++ b/test/simulation_and_solving/simulate_jumps.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, JumpProcesses, OrdinaryDiffEq, Statistics, Test
+using Catalyst, JumpProcesses, OrdinaryDiffEqTsit5, Statistics, Test
 
 # Sets stable rng number.
 using StableRNGs
diff --git a/test/simulation_and_solving/solve_nonlinear.jl b/test/simulation_and_solving/solve_nonlinear.jl
index 0315a65ac5..4f5b19a8a6 100644
--- a/test/simulation_and_solving/solve_nonlinear.jl
+++ b/test/simulation_and_solving/solve_nonlinear.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq
+using Catalyst, NonlinearSolve, OrdinaryDiffEqRosenbrock, SteadyStateDiffEq
 using Random, Test
 
 # Sets stable rng number.
@@ -100,4 +100,4 @@ let
     # Checks output using the ODE's drift function
     @test f_eval(steady_state_network_3, [:X => sol1[X], :Y => sol1[Y], :Y2 => sol1[Y2], :XY2 => sol1[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10
     @test f_eval(steady_state_network_3, [:X => sol2[X], :Y => sol2[Y], :Y2 => sol2[Y2], :XY2 => sol2[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10
-end
\ No newline at end of file
+end
diff --git a/test/spatial_modelling/lattice_reaction_systems.jl b/test/spatial_modelling/lattice_reaction_systems.jl
index 875f9a45ee..b7fcb607bb 100644
--- a/test/spatial_modelling/lattice_reaction_systems.jl
+++ b/test/spatial_modelling/lattice_reaction_systems.jl
@@ -1,7 +1,7 @@
 ### Preparations ###
 
 # Fetch packages.
-using Catalyst, Graphs, OrdinaryDiffEq, Test
+using Catalyst, Graphs, OrdinaryDiffEqTsit5, Test
 
 # Fetch test networks.
 include("../spatial_test_networks.jl")
@@ -392,4 +392,4 @@ let
     @test_throws Exception make_directed_edge_values(lrs_2d, (1.0, 0.0), (1.0, 0.0), (1.0, 0.0))
     @test_throws Exception make_directed_edge_values(lrs_3d, (1.0, 0.0))
     @test_throws Exception make_directed_edge_values(lrs_3d, (1.0, 0.0), (1.0, 0.0))
-end
\ No newline at end of file
+end
diff --git a/test/spatial_modelling/lattice_reaction_systems_ODEs.jl b/test/spatial_modelling/lattice_reaction_systems_ODEs.jl
index f2a247a6b6..3bce253820 100644
--- a/test/spatial_modelling/lattice_reaction_systems_ODEs.jl
+++ b/test/spatial_modelling/lattice_reaction_systems_ODEs.jl
@@ -1,7 +1,7 @@
 ### Preparations ###
 
 # Fetch packages.
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqBDF
 using Random, Statistics, SparseArrays, Test
 
 # Fetch test networks.
@@ -665,4 +665,4 @@ let
     @test_throws ArgumentError ODEProblem(lrs, u0, tspan, [d1 => 1.0])
     @test_throws ArgumentError ODEProblem(lrs, [X1 => [1.0, 2.0, 3.0]], tspan, ps)
     @test_throws ArgumentError ODEProblem(lrs, [X1 => ones(3,1)], tspan, ps)
-end
\ No newline at end of file
+end
diff --git a/test/spatial_modelling/lattice_reaction_systems_lattice_types.jl b/test/spatial_modelling/lattice_reaction_systems_lattice_types.jl
index dbdc233b89..bf37344d30 100644
--- a/test/spatial_modelling/lattice_reaction_systems_lattice_types.jl
+++ b/test/spatial_modelling/lattice_reaction_systems_lattice_types.jl
@@ -1,7 +1,7 @@
 ### Preparations ###
 
 # Fetch packages.
-using Catalyst, Graphs, OrdinaryDiffEq, Test
+using Catalyst, Graphs, OrdinaryDiffEqTsit5, OrdinaryDiffEqBDF, Test
 
 # Fetch test networks.
 include("../spatial_test_networks.jl")
@@ -257,4 +257,4 @@ let
         sol = solve(oprob, Tsit5(); saveat = 1.0, abstol = 1e-9, reltol = 1e-9)
         @test hcat(sol.u...) ≈ sol_base[13:end,:]
     end
-end
\ No newline at end of file
+end
diff --git a/test/spatial_modelling/lattice_simulation_struct_interfacing.jl b/test/spatial_modelling/lattice_simulation_struct_interfacing.jl
index df49fe0743..6a31d21030 100644
--- a/test/spatial_modelling/lattice_simulation_struct_interfacing.jl
+++ b/test/spatial_modelling/lattice_simulation_struct_interfacing.jl
@@ -1,7 +1,7 @@
 ### Preparations ###
 
 # Fetch packages.
-using Catalyst, Graphs, JumpProcesses, OrdinaryDiffEq, SparseArrays, Test
+using Catalyst, Graphs, JumpProcesses, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqRosenbrock, SparseArrays, Test
 
 # Fetch test networks.
 include("../spatial_test_networks.jl")
diff --git a/test/spatial_modelling/simulate_PDEs.jl b/test/spatial_modelling/simulate_PDEs.jl
index f0051e69fd..3ccf17f0e9 100644
--- a/test/spatial_modelling/simulate_PDEs.jl
+++ b/test/spatial_modelling/simulate_PDEs.jl
@@ -1,7 +1,7 @@
 ### Fetch Packages and Set Global Variables ###
 
 # Fetch packages.
-using Catalyst, OrdinaryDiffEq, Test
+using Catalyst, Test
 using ModelingToolkit, DomainSets
 const MT = ModelingToolkit
 
diff --git a/test/upstream/mtk_problem_inputs.jl b/test/upstream/mtk_problem_inputs.jl
index c973596f03..4997a36dfc 100644
--- a/test/upstream/mtk_problem_inputs.jl
+++ b/test/upstream/mtk_problem_inputs.jl
@@ -3,7 +3,7 @@
 ### Prepares Tests ###
 
 # Fetch packages
-using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, StaticArrays, SteadyStateDiffEq, 
+using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEqTsit5, StaticArrays, SteadyStateDiffEq, 
     StochasticDiffEq, Test
 
 # Sets rnd number.
diff --git a/test/upstream/mtk_structure_indexing.jl b/test/upstream/mtk_structure_indexing.jl
index f8768b8ed0..ab2b369bba 100644
--- a/test/upstream/mtk_structure_indexing.jl
+++ b/test/upstream/mtk_structure_indexing.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages
-using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, Plots, SteadyStateDiffEq, StochasticDiffEq, Test
+using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEqTsit5, Plots, SteadyStateDiffEq, StochasticDiffEq, Test
 import ModelingToolkit: getp, getu, setp, setu
 
 # Sets rnd number.