diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml
index 3a84c647..df4414fc 100644
--- a/.github/workflows/Downgrade.yml
+++ b/.github/workflows/Downgrade.yml
@@ -25,7 +25,7 @@ jobs:
 #        if: ${{ matrix.version == '1.6' }}
         with:
         # skip standard libraries..
-          skip: Pkg,TOML,Logging,Random,Dates,LinearAlgebra
+          skip: Aqua,Compat,Dates,IterTools,LinearAlgebra,Logging,Pkg,Random,Test,TOML
           strict: 'false'
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
new file mode 100644
index 00000000..0fc32919
--- /dev/null
+++ b/.github/workflows/SpellCheck.yml
@@ -0,0 +1,13 @@
+name: Spell Check
+
+on: [pull_request]
+
+jobs:
+   typos-check:
+     name: Spell Check with Typos
+     runs-on: ubuntu-latest
+     steps:
+       - name: Checkout Actions Repository
+         uses: actions/checkout@v3
+       - name: Check spelling
+         uses: crate-ci/typos@v1.16.23
diff --git a/.typos.toml b/.typos.toml
new file mode 100644
index 00000000..aa9a6283
--- /dev/null
+++ b/.typos.toml
@@ -0,0 +1,2 @@
+[default.extend-words]
+ numer = "numer"
diff --git a/Project.toml b/Project.toml
index b8999467..9c061245 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ version = "0.5.2"
 
 [deps]
 AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -19,7 +18,6 @@ ParamPunPam = "3e851597-e36f-45a9-af0a-b7781937992f"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
@@ -32,8 +30,9 @@ ModelingToolkitExt = ["ModelingToolkit", "SymbolicUtils", "Symbolics"]
 
 [compat]
 AbstractAlgebra = "0.34.5, 0.35"
-BenchmarkTools = "1"
+Aqua = "0.8"
 Combinatorics = "1"
+CPUSummary = "0.2"
 DataStructures = "0.18"
 Dates = "1.6, 1.7"
 Groebner = "0.6.3"
@@ -41,27 +40,31 @@ IterTools = "1"
 LinearAlgebra = "1.6, 1.7"
 Logging = "1.6, 1.7"
 MacroTools = "0.5"
-ModelingToolkit = "8.74"
+ModelingToolkit = "8.75"
 Nemo = "0.38.3, 0.39"
 ParamPunPam = "0.3.1"
+Pkg = "1.6, 1.7"
 PrecompileTools = "1.2"
 Primes = "0.5"
 Random = "1.6, 1.7"
 SpecialFunctions = "2"
 SymbolicUtils = "1.4, 1.5"
 Symbolics = "5.16"
+Test = "1.6, 1.7"
 TestSetExtensions = "2"
 TimerOutputs = "0.5"
 julia = "1.6, 1.7"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CPUSummary = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 
 [targets]
-test = ["CPUSummary", "Pkg", "Test", "TestSetExtensions"]
+test = ["Aqua", "CPUSummary", "Pkg", "SpecialFunctions", "Test", "TestSetExtensions"]
diff --git a/README.md b/README.md
index 605339a9..0a10eb27 100644
--- a/README.md
+++ b/README.md
@@ -56,7 +56,7 @@ In this example:
   - `y(t)` is the **output variable** which is assumed to be observed in the experiments and, thus, known;
   - `a01, a21, a12, b` are unknown scalar **parameters**.
 
-Note that there may be mulitple inputs and outputs.
+Note that there may be multiple inputs and outputs.
 
 ### Assessing identifiability
 
diff --git a/benchmarking/IdentifiableFunctions/experiments.jl b/benchmarking/IdentifiableFunctions/experiments.jl
index f5b3077a..136421b5 100644
--- a/benchmarking/IdentifiableFunctions/experiments.jl
+++ b/benchmarking/IdentifiableFunctions/experiments.jl
@@ -1116,7 +1116,7 @@ end
 #! format: off
 
 new_rff = StructuralIdentifiability.RationalFunctionField(funcs1)
-cfs = StructuralIdentifiability.beautifuly_generators(new_rff)
+cfs = StructuralIdentifiability.beautiful_generators(new_rff)
 gb_rff = StructuralIdentifiability.RationalFunctionField(cfs)
 
 K = GF(2^31 - 1)
diff --git a/benchmarking/IdentifiableFunctions/homogenization.jl b/benchmarking/IdentifiableFunctions/homogenization.jl
index 94fafc0b..630b585d 100644
--- a/benchmarking/IdentifiableFunctions/homogenization.jl
+++ b/benchmarking/IdentifiableFunctions/homogenization.jl
@@ -15,7 +15,7 @@ Bilirubin2_io = @ODEmodel(
 funcs = find_identifiable_functions(Bilirubin2_io, with_states = true, strategy = (:gb,))
 
 rff = StructuralIdentifiability.RationalFunctionField(funcs)
-cfs = StructuralIdentifiability.beautifuly_generators(rff)
+cfs = StructuralIdentifiability.beautiful_generators(rff)
 rff = StructuralIdentifiability.RationalFunctionField(cfs)
 
 K = GF(2^31 - 1)
diff --git a/docs/src/tutorials/creating_ode.md b/docs/src/tutorials/creating_ode.md
index 00d36284..915a742a 100644
--- a/docs/src/tutorials/creating_ode.md
+++ b/docs/src/tutorials/creating_ode.md
@@ -12,7 +12,7 @@ which involves
 
   - a vector $\mathbf{x}(t)$ of the state variables of the system,
 
-  - a vector $\mathbf{u}(t)$ of extermal inputs,
+  - a vector $\mathbf{u}(t)$ of external inputs,
   - a vector $\mathbf{p}$ of scalar parameters,
   - a vector $\mathbf{y}(t)$ of outputs (i.e., observations),
   - and vectors of rational functions $\mathbf{f}$ and $\mathbf{g}$ (for discussion of the non-rational case, see this [issue](https://github.com/SciML/StructuralIdentifiability.jl/issues/144)).
@@ -54,7 +54,7 @@ assess_identifiability(ode)
 
 ## Defining using `ModelingToolkit`
 
-`StructuralIdentifiability` has an extension `ModelingToolkitExt` which allows to use `ODESystem` from `ModelingToolkit` to descibe
+`StructuralIdentifiability` has an extension `ModelingToolkitExt` which allows to use `ODESystem` from `ModelingToolkit` to describe
 a model. The extension is loaded automatically once `ModelingToolkit` is loaded via `using ModelingToolkit`.
 In this case, one should encode the equations for the states as `ODESystem` and specify the outputs separately.
 In order to do this, we first introduce all functions and scalars:
diff --git a/ext/ModelingToolkitExt.jl b/ext/ModelingToolkitExt.jl
index 4ca8de85..33447508 100644
--- a/ext/ModelingToolkitExt.jl
+++ b/ext/ModelingToolkitExt.jl
@@ -15,6 +15,7 @@ else
     using ..ModelingToolkit
 end
 
+export mtk_to_si
 export assess_local_identifiability, assess_identifiability, find_identifiable_functions
 
 # ------------------------------------------------------------------------------
diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl
index 100252e1..fbfea666 100644
--- a/src/RationalFunctionFields/RationalFunctionField.jl
+++ b/src/RationalFunctionFields/RationalFunctionField.jl
@@ -212,7 +212,7 @@ end
 # ------------------------------------------------------------------------------
 
 """
-    beautifuly_generators(rff::RationalFunctionField)
+    beautiful_generators(rff::RationalFunctionField)
 
 Given a field of rational functions `rff` returns a set of "simpler" and
 standardized generators for `rff`.
@@ -221,7 +221,7 @@ Applies the following passes:
 1. Filter constants,
 2. Remove redundant generators.
 """
-@timeit _to function beautifuly_generators(
+@timeit _to function beautiful_generators(
     rff::RationalFunctionField;
     discard_redundant = true,
     reversed_order = false,
@@ -402,7 +402,7 @@ Returns a set of Groebner bases for multiple different rankings of variables.
     # The first basis in some ordering
     ord = InputOrdering()
     new_rff = groebner_basis_coeffs(rff, seed = seed, ordering = ord)
-    cfs = beautifuly_generators(new_rff)
+    cfs = beautiful_generators(new_rff)
     ordering_to_generators[ord] = cfs
     if isempty(cfs)
         return ordering_to_generators
@@ -427,7 +427,7 @@ Returns a set of Groebner bases for multiple different rankings of variables.
                 ordering = ord,
                 up_to_degree = up_to_degree,
             )
-            cfs = beautifuly_generators(new_rff, discard_redundant = false)
+            cfs = beautiful_generators(new_rff, discard_redundant = false)
             ordering_to_generators[ord] = cfs
         end
     end
@@ -444,7 +444,7 @@ Returns a set of Groebner bases for multiple different rankings of variables.
                 ordering = ord,
                 up_to_degree = up_to_degree,
             )
-            cfs = beautifuly_generators(new_rff, discard_redundant = false)
+            cfs = beautiful_generators(new_rff, discard_redundant = false)
             ordering_to_generators[ord] = cfs
         end
     end
@@ -462,7 +462,7 @@ Returns a set of Groebner bases for multiple different rankings of variables.
                 ordering = ord,
                 up_to_degree = up_to_degree,
             )
-            cfs = beautifuly_generators(new_rff, discard_redundant = false)
+            cfs = beautiful_generators(new_rff, discard_redundant = false)
             ordering_to_generators[ord] = cfs
         end
     end
@@ -479,7 +479,7 @@ function monomial_generators_up_to_degree(
 ) where {T}
     @assert strategy in (:monte_carlo,)
     relations = linear_relations_between_normal_forms(
-        beautifuly_generators(rff),
+        beautiful_generators(rff),
         up_to_degree,
         seed = seed,
     )
@@ -544,7 +544,7 @@ Result is correct (in the Monte-Carlo sense) with probability at least `prob_thr
         seed = seed,
         rational_interpolator = rational_interpolator,
     )
-    new_fracs = beautifuly_generators(new_rff)
+    new_fracs = beautiful_generators(new_rff)
     if isempty(new_fracs)
         return new_fracs
     end
@@ -568,7 +568,7 @@ Result is correct (in the Monte-Carlo sense) with probability at least `prob_thr
 Final cleaning and simplification of generators. 
 Out of $(length(new_fracs)) fractions $(length(new_fracs_unique)) are syntactically unique."""
     runtime =
-        @elapsed new_fracs = beautifuly_generators(RationalFunctionField(new_fracs_unique))
+        @elapsed new_fracs = beautiful_generators(RationalFunctionField(new_fracs_unique))
     @debug "Checking inclusion with probability $prob_threshold"
     runtime =
         @elapsed result = issubfield(rff, RationalFunctionField(new_fracs), prob_threshold)
@@ -578,7 +578,7 @@ Out of $(length(new_fracs)) fractions $(length(new_fracs_unique)) are syntactica
         throw("The new subfield generators are not correct.")
     end
     @info "Inclusion checked with probability $prob_threshold in $(_runtime_logger[:id_inclusion_check]) seconds"
-    @debug "Out of $(length(rff.mqs.nums_qq)) initial generators there are $(length(new_fracs)) indepdendent"
+    @debug "Out of $(length(rff.mqs.nums_qq)) initial generators there are $(length(new_fracs)) independent"
     ranking = generating_set_rank(new_fracs)
     _runtime_logger[:id_ranking] = ranking
     @debug "The ranking of the new set of generators is $ranking"
diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl
index 9477fe7f..77b8dce5 100644
--- a/src/RationalFunctionFields/normalforms.jl
+++ b/src/RationalFunctionFields/normalforms.jl
@@ -1,5 +1,5 @@
 
-# Maintans a row echelon form of a set of vectors over the integrals.
+# Maintains a row echelon form of a set of vectors over the integrals.
 # Works well when the ambient dimension is small.
 mutable struct TinyRowEchelonForm{T}
     rows::Vector{Vector{T}}
diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl
index a18c9828..f5a8cb8e 100644
--- a/src/StructuralIdentifiability.jl
+++ b/src/StructuralIdentifiability.jl
@@ -21,7 +21,7 @@ using ParamPunPam: reduce_mod_p!, specialize_mod_p, AbstractBlackboxIdeal
 ParamPunPam.enable_progressbar(false)
 
 # defining a model
-export ODE, @ODEmodel, @DDSmodel, mtk_to_si
+export ODE, @ODEmodel, @DDSmodel
 
 # assessing identifiability
 export assess_local_identifiability, assess_identifiability
diff --git a/src/discrete.jl b/src/discrete.jl
index 3134fb5e..45bbe0c5 100644
--- a/src/discrete.jl
+++ b/src/discrete.jl
@@ -93,11 +93,11 @@ Input:
 - `param_values` - parameter values, must be a dictionary mapping parameter to a value
 - `initial_conditions` - initial conditions of `ode`, must be a dictionary mapping state variable to a value
 - `input_values` - input sequences in the form input => list of terms; length of the lists must be at least
-                   teh required number of terms in the result
+                   the required number of terms in the result
 - `num_terms` - number of terms to compute
 
 Output:
-- computes a sequence solution with teh required number of terms prec presented as a dictionary state_variable => corresponding sequence
+- computes a sequence solution with the required number of terms prec presented as a dictionary state_variable => corresponding sequence
 """
 function sequence_solution(
     dds::DDS{P},
@@ -301,8 +301,8 @@ function _assess_local_identifiability_discrete_aux(
     # Computing the bound from the Schwartz-Zippel-DeMilo-Lipton lemma
     deg_x = _degree_with_common_denom(values(x_equations(dds)))
     deg_y = _degree_with_common_denom(values(y_equations(dds)))
-    deg_known = reduce(+, map(total_degree, known_ic), init = 0)
-    deg_to_check = max(map(total_degree, funcs_to_check)...)
+    deg_known = reduce(+, map(total_degree_frac, known_ic), init = 0)
+    deg_to_check = max(map(total_degree_frac, funcs_to_check)...)
     Jac_degree = deg_to_check + deg_known
     if deg_x > 1
         Jac_degree += 2 * deg_y * ((deg_x^prec - 1) ÷ (deg_x - 1))
diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl
index da93e6d7..87060fd9 100644
--- a/src/global_identifiability.jl
+++ b/src/global_identifiability.jl
@@ -3,7 +3,7 @@
 
 Takes as input input-output equations, the corresponding ode, a list of functions assumed to be known
 and a flag `with_states`.
-Extracts generators of the field of identifiable functions (with or without states) withous
+Extracts generators of the field of identifiable functions (with or without states) without
 any simplifications.
 
 Returns a tuple consisting of
diff --git a/src/input_macro.jl b/src/input_macro.jl
index 3743117c..a4c8a6f1 100644
--- a/src/input_macro.jl
+++ b/src/input_macro.jl
@@ -263,7 +263,7 @@ Here,
 - `x1`, `x2` are state variables
 - `y` is an output variable
 - `u` is an input variable
-- `a`, `b`, `c` are time-indepdendent parameters
+- `a`, `b`, `c` are time-independent parameters
 
 """
 macro ODEmodel(ex::Expr...)
@@ -297,7 +297,7 @@ Here,
 - `x1`, `x2` are state variables
 - `y` is an output variable
 - `u` is an input variable
-- `a`, `b`, `c` are time-indepdendent parameters
+- `a`, `b`, `c` are time-independent parameters
 
 """
 macro DDSmodel(ex::Expr...)
diff --git a/src/states.jl b/src/states.jl
index 86e0269d..eee9a018 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -6,7 +6,7 @@ Input:
 - `vars` - list of variables
 
 The function considers `f` as `A / B`, where `A` and `B` are polynomials in `vars` with
-coefficients in rational fucntion field in the remaining variables such that at least one of the
+coefficients in rational function field in the remaining variables such that at least one of the
 coefficients is equal to one.
 
 Output:
diff --git a/src/submodels.jl b/src/submodels.jl
index 793e77ba..b492bc5f 100644
--- a/src/submodels.jl
+++ b/src/submodels.jl
@@ -90,7 +90,7 @@ end
 
 # ------------------------------------------------------------------------------
 
-# filters the models containin all states or no states
+# filters the models containing all states or no states
 function filter_max_empty(
     ode::ODE{P},
     submodels::Array{Set{QQMPolyRingElem}, 1},
diff --git a/src/util.jl b/src/util.jl
index 919c2bb6..7acc3bdf 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -1,11 +1,11 @@
 # ------------------------------------------------------------------------------
 
-function Nemo.vars(f::Generic.Frac{<:MPolyRingElem})
-    return collect(union(Set(vars(numerator(f))), Set(vars(denominator(f)))))
+function total_degree_frac(f::Generic.Frac{<:MPolyRingElem})
+    return sum(map(total_degree, unpack_fraction(f)))
 end
 
-function Nemo.total_degree(f::Generic.Frac{<:MPolyRingElem})
-    return sum(map(total_degree, unpack_fraction(f)))
+function total_degree_frac(f::MPolyRingElem)
+    return total_degree(f)
 end
 
 # ------------------------------------------------------------------------------
@@ -376,7 +376,7 @@ Input:
 - `poly` - multivariate polynomial
 - `variables` - a list of variables from the generators of the ring of p
 Output:
-- dictionary with keys being tuples of length `lenght(variables)` and values being polynomials in the variables other than those which are the coefficients at the corresponding monomials (in a smaller polynomial ring)
+- dictionary with keys being tuples of length `length(variables)` and values being polynomials in the variables other than those which are the coefficients at the corresponding monomials (in a smaller polynomial ring)
 """
 function extract_coefficients(poly::P, variables::Array{P, 1}) where {P <: MPolyRingElem}
     xs = gens(parent(poly))
diff --git a/test/extract_coefficients.jl b/test/extract_coefficients.jl
index 7a66fd8a..0d9302d3 100644
--- a/test/extract_coefficients.jl
+++ b/test/extract_coefficients.jl
@@ -1,4 +1,4 @@
-@testset "Coefficient extraction for rational fucntions" begin
+@testset "Coefficient extraction for rational functions" begin
     R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
     C = extract_coefficients_ratfunc(
         (x^2 + y * z - y^2 * z^3 + 3 * x * z^3) // (x + y + z + z^2 * (x^2 + 1)),
diff --git a/test/ode_ps_solution.jl b/test/ode_ps_solution.jl
index dcaa6073..35c59ccf 100644
--- a/test/ode_ps_solution.jl
+++ b/test/ode_ps_solution.jl
@@ -53,7 +53,7 @@
         end
     end
 
-    # Testing ps_ode_solution in conjuntion with the ODE class
+    # Testing ps_ode_solution in conjunction with the ODE class
     for i in 1:30
         # Setting up the ring
         NUMX = 3
@@ -68,7 +68,7 @@
         PType = fpMPolyRingElem
         TDict = Dict{PType, Union{PType, Generic.Frac{PType}}}
 
-        # Generating the intial conditions etc
+        # Generating the initial conditions etc
         ic = Dict(vars[i] => F(rand(-5:5)) for i in 1:NUMX)
         param_vals = Dict(vars[i + NUMX] => F(rand(-5:5)) for i in 1:NUMP)
         inputs =
diff --git a/test/qa.jl b/test/qa.jl
new file mode 100644
index 00000000..c0d2d861
--- /dev/null
+++ b/test/qa.jl
@@ -0,0 +1,11 @@
+using StructuralIdentifiability, Aqua
+@testset "Aqua" begin
+    Aqua.find_persistent_tasks_deps(StructuralIdentifiability)
+    Aqua.test_ambiguities(StructuralIdentifiability, recursive = false)
+    Aqua.test_deps_compat(StructuralIdentifiability)
+    Aqua.test_piracies(StructuralIdentifiability, treat_as_own = [])
+    Aqua.test_project_extras(StructuralIdentifiability)
+    Aqua.test_stale_deps(StructuralIdentifiability)
+    Aqua.test_unbound_args(StructuralIdentifiability)
+    Aqua.test_undefined_exports(StructuralIdentifiability)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index a2e3dc78..1c75f517 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@ using StructuralIdentifiability
 
 using Test
 using TestSetExtensions
+using SpecialFunctions
 
 using StructuralIdentifiability.DataStructures
 using StructuralIdentifiability.Nemo
@@ -136,9 +137,6 @@ end
 
 @info "Testing started"
 
-@test isempty(Test.detect_ambiguities(StructuralIdentifiability))
-@test isempty(Test.detect_unbound_args(StructuralIdentifiability))
-
 all_tests = get_test_files(GROUP)
 if !isempty(ARGS)
     all_tests = ARGS