diff --git a/Project.toml b/Project.toml
index 5d9e1ae58d..66cddb75ea 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,6 +23,7 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
@@ -65,6 +66,7 @@ SciMLBase = "2.46"
 Setfield = "1"
 StructuralIdentifiability = "0.5.8"
 Symbolics = "5.30.1"
+SymbolicUtils = "2.1.2"
 Unitful = "1.12.4"
 julia = "1.10"
 
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 242323a0e2..2c5c6003a2 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -45,6 +45,7 @@ import DataStructures: OrderedDict, OrderedSet
 import Parameters: @with_kw_noshow
 import Symbolics: occursin, wrap
 import Symbolics.RewriteHelpers: hasnode, replacenode
+import SymbolicUtils: getmetadata, hasmetadata, setmetadata
 
 # globals for the modulate
 function default_time_deriv()
diff --git a/src/reaction.jl b/src/reaction.jl
index 8f85d77343..e6202a73ac 100644
--- a/src/reaction.jl
+++ b/src/reaction.jl
@@ -461,9 +461,10 @@ function getmetadata_dict(reaction::Reaction)
 end
 
 """
-hasmetadata(reaction::Reaction, md_key::Symbol)
+    SymbolicUtils.hasmetadata(reaction::Reaction, md_key::Symbol)
 
-Checks if a `Reaction` have a certain metadata field. If it does, returns `true` (else returns `false`).
+Checks if a `Reaction` have a certain metadata field. If it does, returns `true` (else
+returns `false`).
 
 Arguments:
 - `reaction`: The reaction for which we wish to check if a specific metadata field exist.
@@ -475,14 +476,15 @@ reaction = @reaction k, 0 --> X, [description="Production reaction"]
 hasmetadata(reaction, :description)
 ```
 """
-function hasmetadata(reaction::Reaction, md_key::Symbol)
+function SymbolicUtils.hasmetadata(reaction::Reaction, md_key::Symbol)
     return any(isequal(md_key, entry[1]) for entry in getmetadata_dict(reaction))
 end
 
 """
-getmetadata(reaction::Reaction, md_key::Symbol)
+    SymbolicUtils.getmetadata(reaction::Reaction, md_key::Symbol)
 
-Retrieves a certain metadata value from a `Reaction`. If the metadata does not exist, throws an error.
+Retrieves a certain metadata value from a `Reaction`. If the metadata does not exist, throws
+an error.
 
 Arguments:
 - `reaction`: The reaction for which we wish to retrieve a specific metadata value.
@@ -494,7 +496,7 @@ reaction = @reaction k, 0 --> X, [description="Production reaction"]
 getmetadata(reaction, :description)
 ```
 """
-function getmetadata(reaction::Reaction, md_key::Symbol)
+function SymbolicUtils.getmetadata(reaction::Reaction, md_key::Symbol)
     metadata = getmetadata_dict(reaction)
     idx = findfirst(isequal(md_key, entry[1]) for entry in metadata)
     (idx === nothing) &&
@@ -502,7 +504,29 @@ function getmetadata(reaction::Reaction, md_key::Symbol)
     return metadata[idx][2]
 end
 
-### Implemented Reaction Metadata ###
+"""
+    SymbolicUtils.setmetadata(rx::Reaction, key::Symbol, val)
+
+Sets the metadata with key `key` to the value `val`, overwriting if already present or
+adding if not present.
+
+Arguments:
+- `rx`: The reaction to add the metadata too.
+- `key`: `Symbol` representing the metadata's key (i.e. name).
+- `val`: value for the metadata.
+"""
+function SymbolicUtils.setmetadata(rx::Reaction, key::Symbol, val)
+    mdvec = getmetadata_dict(rx)
+    idx = findfirst(isequal(key, first(md)) for md in mdvec)
+    if idx === nothing
+        push!(mdvec, key => val)
+    else
+        mdvec[idx] = key => val
+    end
+    nothing
+end
+
+### Catalyst Defined Reaction Metadata ###
 
 # Noise scaling.
 """
diff --git a/test/dsl/dsl_advanced_model_construction.jl b/test/dsl/dsl_advanced_model_construction.jl
index 5638c9173c..c4e46bb521 100644
--- a/test/dsl/dsl_advanced_model_construction.jl
+++ b/test/dsl/dsl_advanced_model_construction.jl
@@ -189,15 +189,15 @@ let
     @test isequal(Catalyst.getmetadata_dict(r3), Catalyst.getmetadata_dict(rxs[3]))
 
     # Checks that accessor functions works on the DSL.
-    @test Catalyst.hasmetadata(rxs[1], :noise_scaling)
-    @test !Catalyst.hasmetadata(rxs[1], :md_1)
-    @test !Catalyst.hasmetadata(rxs[2], :noise_scaling)
-    @test Catalyst.hasmetadata(rxs[2], :md_1)
-    @test !Catalyst.hasmetadata(rxs[3], :noise_scaling)
-    @test !Catalyst.hasmetadata(rxs[3], :md_1)
-
-    @test isequal(Catalyst.getmetadata(rxs[1], :noise_scaling), η)
-    @test isequal(Catalyst.getmetadata(rxs[2], :md_1), 1.0)
+    @test hasmetadata(rxs[1], :noise_scaling)
+    @test !hasmetadata(rxs[1], :md_1)
+    @test !hasmetadata(rxs[2], :noise_scaling)
+    @test hasmetadata(rxs[2], :md_1)
+    @test !hasmetadata(rxs[3], :noise_scaling)
+    @test !hasmetadata(rxs[3], :md_1)
+
+    @test isequal(getmetadata(rxs[1], :noise_scaling), η)
+    @test isequal(getmetadata(rxs[2], :md_1), 1.0)
 
     # Test that metadata works for @reaction macro.
     rx1 = @reaction k, 2X --> X2, [noise_scaling=$η]
diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl
index 36d5d946d0..3d3b5396f7 100644
--- a/test/reactionsystem_core/reaction.jl
+++ b/test/reactionsystem_core/reaction.jl
@@ -133,10 +133,14 @@ let
     r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata)
 
     @test Catalyst.getmetadata_dict(r) == [:noise_scaling => 0.0]
-    @test Catalyst.hasmetadata(r, :noise_scaling)
-    @test !Catalyst.hasmetadata(r, :nonexisting_metadata)
-    @test Catalyst.getmetadata(r, :noise_scaling) == 0.0
-    @test_throws Exception Catalyst.getmetadata(r, :misc)
+    @test hasmetadata(r, :noise_scaling)
+    @test !hasmetadata(r, :nonexisting_metadata)
+    @test getmetadata(r, :noise_scaling) == 0.0
+    @test_throws Exception getmetadata(r, :misc)
+    setmetadata(r, :test_metadata, 1234)
+    @test getmetadata(r, :test_metadata) == 1234
+    setmetadata(r, :test_metadata, 1111)
+    @test getmetadata(r, :test_metadata) == 1111
 
     metadata_repeated = [:noise_scaling => 0.0, :noise_scaling => 1.0, :metadata_entry => "unused"]
     @test_throws Exception Reaction(k, [X], [X2], [2], [1]; metadata=metadata_repeated)
@@ -153,7 +157,7 @@ let
 
     @test isequal(r1, r2)
     @test Catalyst.getmetadata_dict(r1) == Pair{Symbol,Any}[]
-    @test !Catalyst.hasmetadata(r1, :md)
+    @test !hasmetadata(r1, :md)
 end
 
 # Tests creation.
@@ -173,20 +177,20 @@ let
     r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata)
 
     @test Catalyst.getmetadata_dict(r) isa Vector{Pair{Symbol,Any}}
-    @test Catalyst.hasmetadata(r, :md_1)
-    @test Catalyst.hasmetadata(r, :md_2)
-    @test Catalyst.hasmetadata(r, :md_3)
-    @test Catalyst.hasmetadata(r, :md_4)
-    @test Catalyst.hasmetadata(r, :md_5)
-    @test Catalyst.hasmetadata(r, :md_6)
-    @test !Catalyst.hasmetadata(r, :md_8)
-
-    @test isequal(Catalyst.getmetadata(r, :md_1), 1.0)
-    @test isequal(Catalyst.getmetadata(r, :md_2), false)
-    @test isequal(Catalyst.getmetadata(r, :md_3), "Hello world")
-    @test isequal(Catalyst.getmetadata(r, :md_4), :sym)
-    @test isequal(Catalyst.getmetadata(r, :md_5), X + X2^k -1)
-    @test isequal(Catalyst.getmetadata(r, :md_6), (0.1, 2.0))
+    @test hasmetadata(r, :md_1)
+    @test hasmetadata(r, :md_2)
+    @test hasmetadata(r, :md_3)
+    @test hasmetadata(r, :md_4)
+    @test hasmetadata(r, :md_5)
+    @test hasmetadata(r, :md_6)
+    @test !hasmetadata(r, :md_8)
+
+    @test isequal(getmetadata(r, :md_1), 1.0)
+    @test isequal(getmetadata(r, :md_2), false)
+    @test isequal(getmetadata(r, :md_3), "Hello world")
+    @test isequal(getmetadata(r, :md_4), :sym)
+    @test isequal(getmetadata(r, :md_5), X + X2^k -1)
+    @test isequal(getmetadata(r, :md_6), (0.1, 2.0))
 end
 
 # Tests the noise scaling metadata.