diff --git a/Project.toml b/Project.toml
index b831249319..4448786224 100644
--- a/Project.toml
+++ b/Project.toml
@@ -62,7 +62,7 @@ LaTeXStrings = "1.3.0"
 Latexify = "0.16.6"
 MacroTools = "0.5.5"
 Makie = "0.22.1"
-ModelingToolkit = "9.73"
+ModelingToolkit = "10"
 NetworkLayout = "0.4.7"
 Parameters = "0.12"
 Reexport = "1.0"
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index a079816d87..cedebfe319 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -36,7 +36,7 @@ import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_v
 
 # internal but needed ModelingToolkit functions
 import ModelingToolkit: check_variables,
-                        check_parameters, _iszero, _merge, check_units,
+                        check_parameters, _iszero, check_units,
                         get_unit, check_equations, iscomplete
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index ad7f107961..cd15383c8e 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -276,7 +276,7 @@ Notes:
   units). Unit checking can be disabled by passing the keyword argument `checks=false`.
 """
 struct ReactionSystem{V <: NetworkProperties} <:
-       MT.AbstractTimeDependentSystem
+       MT.AbstractSystem
     """The equations (reactions and algebraic/differential) defining the system."""
     eqs::Vector{CatalystEqType}
     """The Reactions defining the system. """
@@ -398,7 +398,7 @@ function ReactionSystem(eqs, iv, unknowns, ps;
         name = nothing,
         default_u0 = Dict(),
         default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = merge(Dict(default_u0), Dict(default_p)),
         connection_type = nothing,
         checks = true,
         networkproperties = nothing,
@@ -481,13 +481,13 @@ function ReactionSystem(eqs, iv, unknowns, ps;
     end
 
     # Creates the continuous and discrete callbacks.
-    ccallbacks = MT.SymbolicContinuousCallbacks(continuous_events)
-    dcallbacks = MT.SymbolicDiscreteCallbacks(discrete_events)
+    ccallbacks = continuous_events === nothing ? MT.SymbolicContinuousCallback[] : continuous_events
+    dcallbacks = discrete_events === nothing ? MT.SymbolicDiscreteCallback[] : discrete_events
 
     ReactionSystem(
         eqs′, rxs, iv′, sivs′, unknowns′, spcs, ps′, var_to_name, observed, name,
         systems, defaults, connection_type, nps, combinatoric_ratelaws,
-        ccallbacks, dcallbacks, metadata; checks = checks)
+        ccallbacks, dcallbacks, metadata === nothing ? Base.ImmutableDict{Symbol,Any}() : metadata; checks = checks)
 end
 
 # Two-argument constructor (reactions/equations and time variable).
diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index 84d6060391..40674e3c96 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -507,7 +507,7 @@ COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be convert
 
 """
 ```julia
-Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
+Base.convert(::typeof(ODESystem),rs::ReactionSystem)
 ```
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.ODESystem`.
 
@@ -524,11 +524,11 @@ Keyword args and default values:
   with their rational function representation when converting to another system type. Set to
   `false`` to disable.
 """
-function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs),
+function Base.convert(::typeof(ODESystem), rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, remove_conserved = false, checks = false,
         default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true,
+        defaults = merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true,
         kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
@@ -544,7 +544,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
     ODESystem(eqs, get_iv(fullrs), us, ps;
         observed = obs,
         name,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         continuous_events = MT.get_continuous_events(fullrs),
         discrete_events = MT.get_discrete_events(fullrs),
@@ -569,7 +569,7 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
+Base.convert(::typeof(NonlinearSystem),rs::ReactionSystem)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.NonlinearSystem`.
@@ -592,11 +592,11 @@ Keyword args and default values:
   with their rational function representation when converting to another system type. Set to
   `false`` to disable.
 """
-function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs),
+function Base.convert(::typeof(NonlinearSystem), rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         remove_conserved = false, conseqs_remake_warn = true, checks = false,
         default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = merge(Dict(default_u0), Dict(default_p)),
         all_differentials_permitted = false, expand_catalyst_funs = true, kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
@@ -625,7 +625,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     NonlinearSystem(eqs, us, ps;
         name,
         observed = obs, initialization_eqs = initeqs,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         kwargs...)
 end
@@ -661,7 +661,7 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
+Base.convert(::typeof(SDESystem),rs::ReactionSystem)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.SDESystem`.
@@ -679,11 +679,11 @@ Notes:
   with their rational function representation when converting to another system type. Set to
   `false`` to disable.
 """
-function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
+function Base.convert(::typeof(SDESystem), rs::ReactionSystem;
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, checks = false, remove_conserved = false,
         default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = merge(Dict(default_u0), Dict(default_p)),
         expand_catalyst_funs = true,
         kwargs...)
     # Error checks.
@@ -707,7 +707,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
     SDESystem(eqs, noiseeqs, get_iv(flatrs), us, ps;
         observed = obs,
         name,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         continuous_events = MT.get_continuous_events(flatrs),
         discrete_events = MT.get_discrete_events(flatrs),
@@ -747,7 +747,7 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
+Base.convert(::typeof(JumpSystem),rs::ReactionSystem; combinatoric_ratelaws=true)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.JumpSystem`.
@@ -769,10 +769,10 @@ Notes:
   `VariableRateJump` to save the solution before and/or after the jump occurs. Defaults to
   true for both.
 """
-function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs),
+function Base.convert(::typeof(JumpSystem), rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         remove_conserved = nothing, checks = false, default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true,
+        defaults = merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true,
         save_positions = (true, true), physical_scales = nothing, kwargs...)
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
@@ -814,7 +814,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
     JumpSystem(eqs, get_iv(flatrs), us, ps;
         observed = obs,
         name,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         discrete_events = MT.discrete_events(flatrs),
         continuous_events = MT.continuous_events(flatrs),
@@ -923,7 +923,7 @@ Inputs for a JumpProblem from a given `ReactionSystem`.
 # Fields
 $(FIELDS)
 """
-struct JumpInputs{S <: MT.JumpSystem, T <: SciMLBase.AbstractODEProblem}
+struct JumpInputs{S <: MT.AbstractSystem, T <: SciMLBase.AbstractODEProblem}
     """The `JumpSystem` to define the problem over"""
     sys::S
     """The problem the JumpProblem should be defined over, for example DiscreteProblem"""
diff --git a/src/spatial_reaction_systems/lattice_reaction_systems.jl b/src/spatial_reaction_systems/lattice_reaction_systems.jl
index e95c902176..33902af516 100644
--- a/src/spatial_reaction_systems/lattice_reaction_systems.jl
+++ b/src/spatial_reaction_systems/lattice_reaction_systems.jl
@@ -57,7 +57,7 @@ continuous space systems with them is possible, but requires the user to determi
 (the lattice). Better support for continuous space models is a work in progress.
 - Catalyst contains extensive documentation on spatial modelling, which can be found [here](https://docs.sciml.ai/Catalyst/stable/spatial_modelling/lattice_reaction_systems/).
 """
-struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractTimeDependentSystem
+struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractSystem
     # Input values.
     """The (non-spatial) reaction system within each vertex."""
     reactionsystem::ReactionSystem{Q}
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index d3d51823b2..7a381debbf 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -1247,7 +1247,7 @@ let
     # The `@discrete_events` option.
     rn61 = @reaction_network rn1 begin
         @species X(t)
-        @discrete_events [X > 3.0] => [X ~ X - 1]
+        @discrete_events [X > 3.0] => [X ~ Pre(X) - 1]
     end
     rn62 = @reaction_network rn1 begin
         @species X(t)
@@ -1258,7 +1258,7 @@ let
     @test isequal(rn61, rn62)
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @discrete_events [X > 3.0] => [X ~ X - 1] [X < 1.0] => [X ~ X + 1]
+        @discrete_events [X > 3.0] => [X ~ Pre(X) - 1] [X < 1.0] => [X ~ X + 1]
     end
 end
 
diff --git a/test/miscellaneous_tests/reactionsystem_serialisation.jl b/test/miscellaneous_tests/reactionsystem_serialisation.jl
index ff01784952..e0abe979b6 100644
--- a/test/miscellaneous_tests/reactionsystem_serialisation.jl
+++ b/test/miscellaneous_tests/reactionsystem_serialisation.jl
@@ -309,29 +309,29 @@ let
     ]
 
     # Prepares all events.
-    continuous_events_1 = [(A ~ t_1) => [A ~ A + 2.0, X ~ X/2]]
-    continuous_events_2 = [(A ~ t_2) => [A ~ A + 2.0, X ~ X/2]]
-    continuous_events_3 = [(A ~ t_3) => [A ~ A + 2.0, X ~ X/2]]
-    continuous_events_4 = [(A ~ t_4) => [A ~ A + 2.0, X ~ X/2]]
+    continuous_events_1 = [(A ~ t_1) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
+    continuous_events_2 = [(A ~ t_2) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
+    continuous_events_3 = [(A ~ t_3) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
+    continuous_events_4 = [(A ~ t_4) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
     discrete_events_1 = [
-        10.0 => [X2_1 ~ X2_1 + 1.0]
+        10.0 => [X2_1 ~ Pre(X2_1) + 1.0]
         [5.0, 10.0] => [b_1 ~ 2 * b_1]
-        (X > 5.0) => [X2_1 ~ X2_1 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_1 ~ Pre(X2_1) + 1.0, X ~ Pre(X) - 1]
     ]
     discrete_events_2 = [
-        10.0 => [X2_2 ~ X2_2 + 1.0]
+        10.0 => [X2_2 ~ Pre(X2_2) + 1.0]
         [5.0, 10.0] => [b_2 ~ 2 * b_2]
-        (X > 5.0) => [X2_2 ~ X2_2 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_2 ~ Pre(X2_2) + 1.0, X ~ Pre(X) - 1]
     ]
     discrete_events_3 = [
-        10.0 => [X2_3 ~ X2_3 + 1.0]
+        10.0 => [X2_3 ~ Pre(X2_3) + 1.0]
         [5.0, 10.0] => [b_3 ~ 2 * b_3]
-        (X > 5.0) => [X2_3 ~ X2_3 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_3 ~ Pre(X2_3) + 1.0, X ~ Pre(X) - 1]
     ]
     discrete_events_4 = [
-        10.0 => [X2_4 ~ X2_4 + 1.0]
+        10.0 => [X2_4 ~ Pre(X2_4) + 1.0]
         [5.0, 10.0] => [b_4 ~ 2 * b_4]
-        (X > 5.0) => [X2_4 ~ X2_4 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_4 ~ Pre(X2_4) + 1.0, X ~ Pre(X) - 1]
     ]
 
     # Creates the systems.
@@ -367,8 +367,8 @@ let
     rs = @reaction_network begin
         @equations D(V) ~ 1 - V
         @continuous_events begin
-            [X ~ 5.0] => [X ~ X + 1.0]
-            [X ~ 20.0] => [X ~ X - 1.0]
+            [X ~ 5.0] => [X ~ Pre(X) + 1.0]
+            [X ~ 20.0] => [X ~ Pre(X) - 1.0]
         end
         @discrete_events 5.0 => [d ~ d/2]
         d, X --> 0
diff --git a/test/reactionsystem_core/custom_crn_functions.jl b/test/reactionsystem_core/custom_crn_functions.jl
index cfa34a7374..1e72a29e28 100644
--- a/test/reactionsystem_core/custom_crn_functions.jl
+++ b/test/reactionsystem_core/custom_crn_functions.jl
@@ -305,7 +305,7 @@ let
     rates = getfield.(jsyseqs, :rate)
     affects = getfield.(jsyseqs, :affect!)
     reqs = [ Y*X*hill(X, v, K, n), Y*X*mm(X, v, K), hillr(X, v, K, n)*Y*X, Y*X*mmr(X, v, K)]
-    affeqs = [Z ~ 1 + Z, Y ~ -1 + Y, X ~ -1 + X]
+    affeqs = [Z ~ 1 + Pre(Z), Y ~ -1 + Pre(Y), X ~ -1 + Pre(X)]
     @test all(iszero, simplify(rates .- reqs))
     @test all(aff -> isequal(aff, affeqs), affects)
 
@@ -314,7 +314,7 @@ let
     rates = getfield.(jsyseqs, :rate)
     affects = getfield.(jsyseqs, :affect!)
     reqs = [ Y*X*hill2(X, v, K, n), Y*X*mm2(X, v, K), hillr2(X, v, K, n)*Y*X, Y*X*mmr2(X, v, K)]
-    affeqs = [Z ~ 1 + Z, Y ~ -1 + Y, X ~ -1 + X]
+    affeqs = [Z ~ 1 + Pre(Z), Y ~ -1 + Pre(Y), X ~ -1 + Pre(X)]
     @test all(iszero, simplify(rates .- reqs))
     @test all(aff -> isequal(aff, affeqs), affects)
 end
\ No newline at end of file
diff --git a/test/reactionsystem_core/events.jl b/test/reactionsystem_core/events.jl
index 6821cdd21d..12967e99e3 100644
--- a/test/reactionsystem_core/events.jl
+++ b/test/reactionsystem_core/events.jl
@@ -73,11 +73,11 @@ let
         Reaction(p, nothing, [X]),
         Reaction(d, [X], nothing)
     ]
-    continuous_events = [α ~ t] => [A ~ A + a]
+    continuous_events = [α ~ t] => [A ~ Pre(A) + a]
     discrete_events = [2.0 => [A ~ α + a]]
     @named rs_ce = ReactionSystem(rxs, t; continuous_events)
     @named rs_de = ReactionSystem(rxs, t; discrete_events)
-    continuous_events = [[α ~ t] => [A ~ A + α]]
+    continuous_events = [[α ~ t] => [A ~ Pre(A) + α]]
     discrete_events = [2.0 => [A ~ a]]
     @named rs_ce_de = ReactionSystem(rxs, t; continuous_events, discrete_events)
     rs_ce = complete(rs_ce)
@@ -114,7 +114,7 @@ let
     end
 
     # Handles `JumpInput`s and `JumpProblem`s (these cannot contain continuous events or variables).
-    discrete_events = [2.0 => [A ~ A + α]]
+    discrete_events = [2.0 => [A ~ Pre(A) + α]]
     @named rs_de_2 = ReactionSystem(rxs, t; discrete_events)
     rs_de_2 = complete(rs_de_2)
     jin = JumpInputs(rs_de_2, u0, (0.0, 10.0), ps)
@@ -196,13 +196,13 @@ let
         @parameters thres=7.0 dY_up
         @variables Z(t)
         @continuous_events begin
-            [t ~ 2.5] => [p ~ p + 0.2]
-            [X ~ thres, Y ~ X] => [X ~ X - 0.5, Z ~ Z + 0.1]
+            [t ~ 2.5] => [p ~ Pre(p) + 0.2]
+            [X ~ thres, Y ~ X] => [X ~ Pre(X) - 0.5, Z ~ Pre(Z) + 0.1]
         end
         @discrete_events begin
-            2.0 => [dX ~ dX + 0.01, dY ~ dY + dY_up]
-            [1.0, 5.0] => [p ~ p - 0.1]
-            (Z > Y) => [Z ~ Z - 0.1]
+            2.0 => [dX ~ Pre(dX) + 0.01, dY ~ Pre(dY) + dY_up]
+            [1.0, 5.0] => [p ~ Pre(p) - 0.1]
+            (Z > Y) => [Z ~ Pre(Z) - 0.1]
         end
 
         (p, dX), 0 <--> X
@@ -221,13 +221,13 @@ let
         Reaction(dY, [Y], nothing, [1], nothing)
     ]
     continuous_events = [
-        [t ~ 2.5] => [p ~ p + 0.2]
-        [X ~ thres, Y ~ X] => [X ~ X - 0.5, Z ~ Z + 0.1]
+        [t ~ 2.5] => [p ~ Pre(p) + 0.2]
+        [X ~ thres, Y ~ X] => [X ~ Pre(X) - 0.5, Z ~ Pre(Z) + 0.1]
     ]
     discrete_events = [
-        2.0 => [dX ~ dX + 0.01, dY ~ dY + dY_up]
-        [1.0, 5.0] => [p ~ p - 0.1]
-        (Z > Y) => [Z ~ Z - 0.1]
+        2.0 => [dX ~ Pre(dX) + 0.01, dY ~ Pre(dY) + dY_up]
+        [1.0, 5.0] => [p ~ Pre(p) - 0.1]
+        (Z > Y) => [Z ~ Pre(Z) - 0.1]
     ]
     rn_prog = ReactionSystem(rxs, t; continuous_events, discrete_events, name = :rn)
     rn_prog = complete(rn_prog)
@@ -248,60 +248,60 @@ end
 let
     # Quantity in event not declared elsewhere (continuous events).
     @test_throws Exception @eval @reaction_network begin
-        @continuous_events X ~ 2.0 => [X ~ X + 1]
+        @continuous_events X ~ 2.0 => [X ~ Pre(X) + 1]
     end
 
     # Scalar condition (continuous events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @continuous_events X ~ 2.0 => [X ~ X + 1]
+        @continuous_events X ~ 2.0 => [X ~ Pre(X) + 1]
     end
 
     # Scalar affect (continuous events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @continuous_events [X ~ 2.0] => X ~ X + 1
+        @continuous_events [X ~ 2.0] => X ~ Pre(X) + 1
     end
 
     # Tuple condition (continuous events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @continuous_events (X ~ 2.0,) => [X ~ X + 1]
+        @continuous_events (X ~ 2.0,) => [X ~ Pre(X) + 1]
     end
 
     # Tuple affect (continuous events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @continuous_events [X ~ 2.0] => (X ~ X + 1,)
+        @continuous_events [X ~ 2.0] => (X ~ Pre(X) + 1,)
     end
 
     # Non-equation condition (continuous events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @continuous_events [X - 2.0] => [X ~ X + 1]
+        @continuous_events [X - 2.0] => [X ~ Pre(X) + 1]
     end
 
     # Quantity in event not declared elsewhere (discrete events).
     @test_throws Exception @eval @reaction_network begin
-        @discrete_events X ~ 2.0 => [X ~ X + 1]
+        @discrete_events X ~ 2.0 => [X ~ Pre(X) + 1]
     end
 
     # Scalar affect (discrete events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @discrete_events 1.0 => X ~ X + 1
+        @discrete_events 1.0 => X ~ Pre(X) + 1
     end
 
     # Tuple affect (discrete events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @discrete_events 1.0 => (X ~ X + 1, )
+        @discrete_events 1.0 => (X ~ Pre(X) + 1, )
     end
 
     # Equation condition (discrete events).
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @discrete_events X ~ 1.0 => [X ~ X + 1]
+        @discrete_events X ~ 1.0 => [X ~ Pre(X) + 1]
     end
 end
 
@@ -382,12 +382,12 @@ let
         @default_noise_scaling 0.0
         @parameters add::Int64
         @continuous_events begin
-            [X ~ 90.0] => [X ~ X + 10.0]
+            [X ~ 90.0] => [X ~ Pre(X) + 10.0]
         end
         @discrete_events begin
-            [5.0, 10.0] => [X ~ X + add, Y ~ Y + add]
-            20.0 => [X ~ X + add]
-            (Y < X) => [Y ~ Y + add]
+            [5.0, 10.0] => [X ~ Pre(X) + add, Y ~ Pre(Y) + add]
+            20.0 => [X ~ Pre(X) + add]
+            (Y < X) => [Y ~ Pre(Y) + add]
         end
         (p,d), 0 <--> X
         (p,d), 0 <--> Y
@@ -395,9 +395,9 @@ let
     rn_dics_events = @reaction_network begin
         @parameters add::Int64
         @discrete_events begin
-            [5.0, 10.0] => [X ~ X + add, Y ~ Y + add]
-            20.0 => [X ~ X + add]
-            (Y < X) => [Y ~ Y + add]
+            [5.0, 10.0] => [X ~ Pre(X) + add, Y ~ Pre(Y) + add]
+            20.0 => [X ~ Pre(X) + add]
+            (Y < X) => [Y ~ Pre(Y) + add]
         end
         (p,d), 0 <--> X
         (p,d), 0 <--> Y
diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl
index 31529f82d6..419530f7e3 100644
--- a/test/reactionsystem_core/reactionsystem.jl
+++ b/test/reactionsystem_core/reactionsystem.jl
@@ -152,7 +152,7 @@ let
     du = oderhs(last.(u), last.(p), 0.0)
     G = sdenoise(last.(u), last.(p), 0.0)
     sdesys = complete(convert(SDESystem, rs))
-    sf = SDEFunction{false}(sdesys, unknowns(rs), parameters(rs))
+    sf = SDEFunction{false}(sdesys; u0=unknowns(rs), p=parameters(rs))
     sprob = SDEProblem(rs, u, (0.0, 0.0), p)
     du2 = sf.f(sprob.u0, sprob.p, 0.0)
 
@@ -193,9 +193,9 @@ let
     midxs = 1:14
     cidxs = 15:18
     vidxs = 19:20
-    @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.MassActionJump, midxs))
-    @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.ConstantRateJump, cidxs))
-    @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.VariableRateJump, vidxs))
+    @test all(map(i -> typeof(ModelingToolkit.get_jumps(js)[i]) <: JumpProcesses.MassActionJump, midxs))
+    @test all(map(i -> typeof(ModelingToolkit.get_jumps(js)[i]) <: JumpProcesses.ConstantRateJump, cidxs))
+    @test all(map(i -> typeof(ModelingToolkit.get_jumps(js)[i]) <: JumpProcesses.VariableRateJump, vidxs))
 
     p = rand(rng, length(k))
     pmap = [k => p]
@@ -236,10 +236,11 @@ let
                                  integrator -> (integrator.u[4] -= 2; integrator.u[5] -= 1; integrator.u[6] += 2))
 
     unknownoid = Dict(unknown => i for (i, unknown) in enumerate(unknowns(js)))
-    dprob = DiscreteProblem(js, u0map, (0.0, 10.0), pmap)
+    dprob = DiscreteProblem(js, u0map, (0.0, 10.0), pmap; check_compatibility = false, check_length = false)
     mtkpars = dprob.p
     jspmapper = ModelingToolkit.JumpSysMajParamMapper(js, mtkpars)
-    symmaj = ModelingToolkit.assemble_maj(equations(js).x[1], unknownoid, jspmapper)
+    maj_jumps = [ModelingToolkit.get_jumps(js)[i] for i in midxs]
+    symmaj = ModelingToolkit.assemble_maj(maj_jumps, unknownoid, jspmapper)
     maj = MassActionJump(symmaj.param_mapper(mtkpars), symmaj.reactant_stoch, symmaj.net_stoch,
                          symmaj.param_mapper, scale_rates = false)
     for i in midxs
@@ -248,7 +249,7 @@ let
         @test jumps[i].net_stoch == maj.net_stoch[i]
     end
     for i in cidxs
-        crj = ModelingToolkit.assemble_crj(js, equations(js)[i], unknownoid)
+        crj = ModelingToolkit.assemble_crj(js, ModelingToolkit.get_jumps(js)[i], unknownoid)
         @test isapprox(crj.rate(u0, mtkpars, ttt), jumps[i].rate(u0, p, ttt))
         fake_integrator1 = (u = zeros(6), p = mtkpars, t = 0.0)
         fake_integrator2 = (u = zeros(6), p, t = 0.0)
@@ -257,7 +258,7 @@ let
         @test fake_integrator1.u == fake_integrator2.u
     end
     for i in vidxs
-        crj = ModelingToolkit.assemble_vrj(js, equations(js)[i], unknownoid)
+        crj = ModelingToolkit.assemble_vrj(js, ModelingToolkit.get_jumps(js)[i], unknownoid)
         @test isapprox(crj.rate(u0, mtkpars, ttt), jumps[i].rate(u0, p, ttt))
         fake_integrator1 = (u = zeros(6), p = mtkpars, t = 0.0)
         fake_integrator2 = (u = zeros(6), p, t = 0.0)
@@ -501,7 +502,7 @@ let
     @test issetequal(_ps, [A, k1, k2])
     du1 = zeros(4)
     du2 = zeros(4)
-    sprob = SDEProblem(ssys, u0map, tspan, pmap; check_length = false)
+    sprob = SDEProblem(ssys, u0map, tspan, pmap; check_length = false, check_compatibility = false)
     sprob.f(du1, sprob.u0, sprob.p, 1.0)
     fs!(du2, u0, p, 1.0)
     @test isapprox(du1, du2)
diff --git a/test/simulation_and_solving/hybrid_models.jl b/test/simulation_and_solving/hybrid_models.jl
index a628e4a7a8..61a1b885f2 100644
--- a/test/simulation_and_solving/hybrid_models.jl
+++ b/test/simulation_and_solving/hybrid_models.jl
@@ -254,7 +254,7 @@ let
         @parameters X0 Y0
         @species X(τ) = X0 [description = "A jump only species"] Y(τ) = Y0 [description = "An ODE only species"]
         @observables Ztot ~ Z1 + Z2
-        @discrete_events [1.0] => [Z1 ~ Z1 + 1.0]
+        @discrete_events [1.0] => [Z1 ~ Pre(Z1) + 1.0]
         @continuous_events [Y ~ 1.0] => [Y ~ 5.0]
         @equations Δ(V) ~ Z1 + X^2 - V
         (p,d), 0 <--> X
diff --git a/test/spatial_modelling/lattice_reaction_systems.jl b/test/spatial_modelling/lattice_reaction_systems.jl
index 5c6534eacf..f24aa172f6 100644
--- a/test/spatial_modelling/lattice_reaction_systems.jl
+++ b/test/spatial_modelling/lattice_reaction_systems.jl
@@ -247,7 +247,7 @@ end
 
     # Events.
     rs3 = @reaction_network begin
-        @discrete_events [1.0] => [p ~ p + 1]
+        @discrete_events [1.0] => [p ~ Pre(p) + 1]
         (p,d), 0 <--> X
     end
     @test_throws ArgumentError LatticeReactionSystem(rs3, [tr], short_path)
diff --git a/test_pre.jl b/test_pre.jl
new file mode 100644
index 0000000000..23426c292c
--- /dev/null
+++ b/test_pre.jl
@@ -0,0 +1,28 @@
+using Catalyst
+using ModelingToolkit: Pre
+
+# Set up variables
+t = default_t()
+
+# Test if Pre() operator works with discrete events
+@species X(t)
+@parameters p
+
+# Test basic discrete event with Pre() operator
+try
+    rn = @reaction_network begin
+        @discrete_events [1.0] => [X ~ Pre(X) + 1]
+        p, 0 --> X
+    end
+    println("✓ Pre() operator works with discrete events")
+    
+    # Test continuous event with Pre() operator  
+    rn2 = @reaction_network begin
+        @continuous_events [X ~ 2.0] => [X ~ Pre(X) - 0.5]
+        p, 0 --> X
+    end
+    println("✓ Pre() operator works with continuous events")
+
+catch e
+    println("✗ Error with Pre() operator: ", e)
+end
\ No newline at end of file