diff --git a/Project.toml b/Project.toml
index 6d94040ef2..35122c7efd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -52,7 +52,7 @@ RuntimeGeneratedFunctions = "0.5.12"
 Setfield = "1"
 StructuralIdentifiability = "0.5.1"
 SymbolicUtils = "1.0.3"
-Symbolics = "5.14"
+Symbolics = "5.27"
 Unitful = "1.12.4"
 julia = "1.9"
 
diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index 2a29e3e918..1cfd8abf60 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -157,10 +157,17 @@ species
 nonspecies
 reactionparams
 reactions
-has_diff_equations
+is_reaction
+is_alg_equation
+is_diff_equation
+alg_equations
 diff_equations
 has_alg_equations
-alg_equations
+has_diff_equations
+get_alg_eqs
+get_diff_eqs
+has_alg_eqs
+has_diff_eqs
 numspecies
 numparams
 numreactions
diff --git a/src/networkapi.jl b/src/networkapi.jl
index 7cd1852db2..72f385d46d 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -97,54 +97,6 @@ function reactions(network)
     [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
 end
 
-"""
-    diff_equations(network)
-Given a [`ReactionSystem`](@ref), return a vector of all `Equations` in the system that are differential equations (contains a derivative with respect to any variable).
-Notes:
-- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
-"""
-function diff_equations(network)
-    eqs = equations(network)
-    filter!(!isreaction, eqs)
-    systems = filter_nonrxsys(network)
-    isempty(systems) && (return rxs)
-    [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
-end
-
-"""
-    has_diff_equations(network)
-Given a [`ReactionSystem`](@ref), check whether it contain any differential equations (i.e. in addition to those generated through reactions).
-Notes:
-- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
-"""
-function has_diff_equations(network)
-    return !isempty(diff_equations(network))
-end
-
-"""
-    alg_equations(network)
-Given a [`ReactionSystem`](@ref), return a vector of all `Equations` in the system that are algebraic equations (does not contain any derivatives).
-Notes:
-- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
-"""
-function alg_equations(network)
-    eqs = equations(network)
-    filter!(!isreaction, eqs)
-    filter!(!isreaction, eqs)
-    systems = filter_nonrxsys(network)
-    isempty(systems) && (return rxs)
-    [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
-end
-
-"""
-    has_alg_equations(network)
-Given a [`ReactionSystem`](@ref), check whether it contain any algebraic equations.
-Notes:
-- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
-"""
-function has_alg_equations(network)
-    return !isempty(alg_equations(network))
-end
 
 """
     speciesmap(network)
@@ -633,6 +585,94 @@ end
 symmap_to_varmap(sys, symmap) = symmap
 #error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.")
 
+
+### Equation Handling Accessors ###
+
+"""
+    is_reaction(eq::Equation)
+
+Returns `true` if the input is a `Reaction`, else false.
+"""
+function is_reaction(eq::Equation)
+    return eq isa Reaction
+end
+
+"""
+    is_alg_equation(eq::Equation)
+
+Returns `true` if the input is an algebraic equation that does not contain any differentials.
+"""
+function is_alg_equation(eq)
+    return isdefined(eq, :lhs) && isdefined(eq, :rhs) && !is_diff_equation(eq)
+end
+
+"""
+    is_diff_equation(eq::Equation)
+
+Returns `true` if the input is an that contains at least one differential.
+"""
+function is_diff_equation(eq)
+    isdefined(eq, :lhs) && occursin(is_derivative, wrap(eq.lhs)) && (return true)
+    isdefined(eq, :rhs) && occursin(is_derivative, wrap(eq.rhs)) && (return true)
+    return false
+end
+
+"""
+    alg_equations(rs)
+
+Returns all the algebraic equations (that does not contain differnetials) in `rs` and its subsystems.
+"""
+alg_equations(rs) = filter(is_alg_equation, equations(rs))
+
+"""
+    diff_equations(rs)
+
+Returns all the differnetials equations (equations that contain differnetials) in `rs` and its subsystems.
+"""
+diff_equations(rs) = filter(is_diff_equation, equations(rs))
+
+"""
+    has_alg_equations(rs)
+
+Returns true if `rs` (or any of its subsystems) has an algebraic equation (that does not contain a differential).
+"""
+has_alg_equations(rs) = any(is_alg_equation, equations(rs))
+
+"""
+    has_diff_equations(rs)
+
+Returns true if `rs` (or any of its subsystems) has a differnetials equations (equations that contain differnetials) .
+"""
+has_diff_equations(rs) = any(is_diff_equation, equations(rs))
+
+"""
+    get_alg_eqs(rs)
+
+Returns all the algebraic equations (that does not contain differnetials) in `rs` and (but not its subsystems).
+"""
+get_alg_eqs(rs) = filter(is_alg_equation, get_eqs(rs))
+
+"""
+    diff_equations(rs)
+
+Returns all the differnetial equations (equations that contain differnetials) in `rs` (but not its subsystems).
+"""
+get_diff_eqs(rs) = filter(is_diff_equation, get_eqs(rs))
+
+"""
+    has_alg_equations(rs)
+
+Returns true if `rs` has an algebraic equation (that does not contain a differential). Does not consider subsystems.
+"""
+has_alg_eqs(rs) = any(is_alg_equation, get_eqs(rs))
+
+"""
+    has_diff_equations(rs)
+
+Returns true if `rs` has a differnetial equations (equations that contain differnetials). Does not consider subsystems.
+"""
+has_diff_eqs(rs) =  any(is_diff_equation, get_eqs(rs))
+
 ######################## reaction complexes and reaction rates ###############################
 
 """
diff --git a/src/reaction_network.jl b/src/reaction_network.jl
index f314ad44f7..23ddb7ea2a 100644
--- a/src/reaction_network.jl
+++ b/src/reaction_network.jl
@@ -816,15 +816,35 @@ function make_observed_eqs(observables_expr)
     end 
 end
 
-# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a
-# default metadata value to the `default_reaction_metadata` vector.
-function check_default_noise_scaling!(default_reaction_metadata, options)
-    if haskey(options, :default_noise_scaling)
-        if (length(options[:default_noise_scaling].args) != 3) # Becasue of how expressions are, 3 is used.
-            error("@default_noise_scaling should only have a single input, this appears not to be the case: \"$(options[:default_noise_scaling])\"")
+# Reads the variables options. Outputs:
+# `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
+# `dtexpr`: If a differentialequation is defined, the default derrivative (D ~ Differential(t)) must be defined.
+# `equations`: a vector with the equations provided.
+function read_equations_options(options, variables_declared)
+    # Prepares the equations. First, extracts equations from provided option (converting to block form if requried).
+    # Next, uses MTK's `parse_equations!` function to split input into a vector with the equations.
+    eqs_input = haskey(options, :equations) ? options[:equations].args[3] : :(begin end)
+    eqs_input = option_block_form(eqs_input)
+    equations = Expr[]
+    ModelingToolkit.parse_equations!(Expr(:block), equations, Dict{Symbol, Any}(), eqs_input)
+
+    # 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 = []
+    add_default_diff = false
+    for eq in equations
+        ((eq.head != :call) || (eq.args[1] != :~)) && error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")     
+        (eq.args[2] isa Symbol || eq.args[2].head != :call) && continue
+        if (eq.args[2].args[1] == :D) && (eq.args[2].args[2] isa Symbol) && (length(eq.args[2].args) == 2)
+            diff_var = eq.args[2].args[2]
+            in(diff_var, forbidden_symbols_error) && error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
+            add_default_diff = true
+            in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
         end
-        push!(default_reaction_metadata.args, :(:noise_scaling => $(options[:default_noise_scaling].args[3])))
     end
+
+    return vars_extracted, add_default_diff, equations
 end
 
 # Creates an expression declaring differentials.
@@ -850,8 +870,8 @@ end
 # Read the events (continious or discrete) provided as options to the DSL. Returns an expression which evalutes to these.
 function read_events_option(options, event_type::Symbol)    
     # Prepares the events, if required to, converts them to block form.
-    (event_type in [:continuous_events]) || error("Trying to read an unsupported event type.")
-    events_input = haskey(options, event_type) ? options[event_type].args[3] : :(begin end)
+    (event_type in [:continuous_events, :discrete_events]) || error("Trying to read an unsupported event type.")
+    events_input = haskey(options, event_type) ? options[event_type].args[3] : MacroTools.striplines(:(begin end))
     events_input = option_block_form(events_input)
 
     # Goes throgh the events, checks for errors, and adds them to the output vector.
@@ -874,6 +894,17 @@ function read_events_option(options, event_type::Symbol)
     return events_expr
 end
 
+# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a
+# default metadata value to the `default_reaction_metadata` vector.
+function check_default_noise_scaling!(default_reaction_metadata, options)
+    if haskey(options, :default_noise_scaling)
+        if (length(options[:default_noise_scaling].args) != 3) # Becasue of how expressions are, 3 is used.
+            error("@default_noise_scaling should only have a single input, this appears not to be the case: \"$(options[:default_noise_scaling])\"")
+        end
+        push!(default_reaction_metadata.args, :(:noise_scaling => $(options[:default_noise_scaling].args[3])))
+    end
+end
+
 ### Functionality for expanding function call to custom and specific functions ###
 
 #Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 7eab89973e..088eb0cde0 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -546,9 +546,9 @@ struct ReactionSystem{V <: NetworkProperties} <:
                             name, systems, defaults, connection_type, nps, cls, cevs, devs,
                             metadata, complete; checks::Bool = true)
                             
-        # Filters away any potential obervables from `states` and `spcs`.
+        # Filters away any potential obervables from `unknowns` and `spcs`.
         obs_vars = [obs_eq.lhs for obs_eq in observed]
-        states = filter(state -> !any(isequal(state, obs_var) for obs_var in obs_vars), states)
+        unknowns = filter(state -> !any(isequal(state, obs_var) for obs_var in obs_vars), unknowns)
         spcs = filter(spc -> !any(isequal(spc, obs_var) for obs_var in obs_vars), spcs)
 
         # unit checks are for ODEs and Reactions only currently
@@ -558,7 +558,7 @@ struct ReactionSystem{V <: NetworkProperties} <:
             check_parameters(ps, iv)
             nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
             !isempty(nonrx_eqs) && check_equations(nonrx_eqs, iv)
-            check_equations(equations(cevs), iv)
+            !isempty(cevs) && check_equations(equations(cevs), iv)
         end
 
         if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0)
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 5ea24980c7..2c0e7b0dbe 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -657,7 +657,7 @@ let
 
     # Checks that the internal structures have the correct lengths.
     @test length(species(rn)) == 1
-    @test length(states(rn)) == 3
+    @test length(unknowns(rn)) == 3
     @test length(reactions(rn)) == 2
     @test length(equations(rn)) == 4
     @test !has_diff_equations(rn)
@@ -666,9 +666,9 @@ let
     @test isequal(alg_equations(rn), [X + 5 ~ k*S, 3Y + X  ~ S + X*d])
 
     # Checks that the internal structures contain the correct stuff, and are correctly sorted.
-    @test isspecies(states(rn)[1])
-    @test !isspecies(states(rn)[2])
-    @test !isspecies(states(rn)[3])
+    @test isspecies(unknowns(rn)[1])
+    @test !isspecies(unknowns(rn)[2])
+    @test !isspecies(unknowns(rn)[3])
     @test equations(rn)[1] isa Reaction
     @test equations(rn)[2] isa Reaction
     @test equations(rn)[3] isa Equation
@@ -769,7 +769,7 @@ let
 
     # Checks that the internal structures have the correct lengths.
     @test length(species(rn)) == 1
-    @test length(states(rn)) == 3
+    @test length(unknowns(rn)) == 3
     @test length(reactions(rn)) == 2
     @test length(equations(rn)) == 4
     @test has_diff_equations(rn)
@@ -778,9 +778,9 @@ let
     @test length(alg_equations(rn)) == 1
 
     # Checks that the internal structures contain the correct stuff, and are correctly sorted.
-    @test isspecies(states(rn)[1])
-    @test !isspecies(states(rn)[2])
-    @test !isspecies(states(rn)[3])
+    @test isspecies(unknowns(rn)[1])
+    @test !isspecies(unknowns(rn)[2])
+    @test !isspecies(unknowns(rn)[3])
     @test equations(rn)[1] isa Reaction
     @test equations(rn)[2] isa Reaction
     @test equations(rn)[3] isa Equation