diff --git a/docs/Project.toml b/docs/Project.toml
index 21e7b1811..e449e2a58 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -7,5 +7,5 @@ StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 [compat]
 BenchmarkTools = "1.3"
 Documenter = "0.27, 1"
-ModelingToolkit = "8.34"
+ModelingToolkit = "8.74"
 StructuralIdentifiability = "0.5"
diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml
index e4c519b6b..e449e2a58 100644
--- a/docs/src/assets/Project.toml
+++ b/docs/src/assets/Project.toml
@@ -7,5 +7,5 @@ StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 [compat]
 BenchmarkTools = "1.3"
 Documenter = "0.27, 1"
-ModelingToolkit = "8.34"
-StructuralIdentifiability = "0.4"
+ModelingToolkit = "8.74"
+StructuralIdentifiability = "0.5"
diff --git a/docs/src/input/input.md b/docs/src/input/input.md
index d41bf7374..5fa6ad68c 100644
--- a/docs/src/input/input.md
+++ b/docs/src/input/input.md
@@ -1,4 +1,4 @@
-# Parsing input ODE system
+# Parsing input system
 
 ```@docs
 @ODEmodel(ex::Expr...)
@@ -11,3 +11,9 @@ set_parameter_values
 ```@docs
 linear_compartment_model
 ```
+
+## Discrete-time systems
+
+```@docs
+@DDSmodel
+```
diff --git a/docs/src/tutorials/creating_ode.md b/docs/src/tutorials/creating_ode.md
index ec7d039ff..00d362846 100644
--- a/docs/src/tutorials/creating_ode.md
+++ b/docs/src/tutorials/creating_ode.md
@@ -54,7 +54,9 @@ assess_identifiability(ode)
 
 ## Defining using `ModelingToolkit`
 
-Alternatively, one can use `ModelingToolkit`: encode the equations for the states as `ODESystem` and specify the outputs separately.
+`StructuralIdentifiability` has an extension `ModelingToolkitExt` which allows to use `ODESystem` from `ModelingToolkit` to descibe
+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:
 
 ```@example 2; continued = true
diff --git a/docs/src/tutorials/discrete_time.md b/docs/src/tutorials/discrete_time.md
index 212f59db8..c8f0571f8 100644
--- a/docs/src/tutorials/discrete_time.md
+++ b/docs/src/tutorials/discrete_time.md
@@ -1,13 +1,20 @@
 # Identifiability of Discrete-Time Models (Local)
 
-Now we consider a discrete-time model in the state-space form
+Now we consider a discrete-time model in the state-space form. Such a model is typically written either in terms of **shift**:
+
+$\begin{cases}
+\mathbf{x}(t + 1) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}, \mathbf{u}(t)),\\
+\mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t), \mathbf{p}, \mathbf{u(t)}),
+\end{cases}$
+
+or in terms of **difference**
 
 $\begin{cases}
 \Delta\mathbf{x}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}, \mathbf{u}(t)),\\
 \mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t), \mathbf{p}, \mathbf{u(t)}),
-\end{cases} \quad \text{where} \quad \Delta \mathbf{x}(t) := \mathbf{x}(t + 1) - \mathbf{x}(t)$
+\end{cases} \quad \text{where} \quad \Delta \mathbf{x}(t) := \mathbf{x}(t + 1) - \mathbf{x}(t).$
 
-and $\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively,
+In both cases,$\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively,
 and $\mathbf{p}$ are scalar parameters.
 As in the ODE case, we will call that a parameter or a states (or a function of them) is **identifiable** if its value can be recovered from
 time series for inputs and outputs (in the generic case, see Definition 3 in [^1] for details).
@@ -21,50 +28,53 @@ Currently, `StructuralIdentifiability.jl` allows to assess only local identifiab
 and below we will describe how this can be done.
 As a running example, we will use the following discrete version of the [SIR](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) model:
 
-$\begin{cases}
-\Delta S(t) = S(t) - \beta S(t) I(t),\\
-\Delta I(t) = I(t) + \beta S(t) I(t) - \alpha I(t),\\
-\Delta R(t) = R(t) + \alpha I(t),\\
+$
+\begin{cases}
+S(t + 1) = S(t) - \beta S(t) I(t),\\
+I(t + 1) = I(t) + \beta S(t) I(t) - \alpha I(t),\\
+R(t + 1) = R(t) + \alpha I(t),\\
+y(t) = I(t),
+\end{cases}
+\quad \text{or}
+\quad
+\begin{cases}
+\Delta S(t) = -\beta S(t) I(t),\\
+\Delta I(t) = \beta S(t) I(t) - \alpha I(t),\\
+\Delta R(t) = \alpha I(t),\\
 y(t) = I(t),
 \end{cases}$
 
 where the observable is `I`, the number of infected people.
-We start with creating a system as a `DiscreteSystem` from `ModelingToolkit`:
+The native way to define such a model in `StructuralIdentifiability` is to use `@DDSmodel` macro which
+uses the shift notation:
 
-```@example discrete
-using ModelingToolkit
+```@example discrete_dds
 using StructuralIdentifiability
 
-@parameters α β
-@variables t S(t) I(t) R(t) y(t)
-D = Difference(t; dt = 1.0)
-
-eqs = [D(S) ~ S - β * S * I, D(I) ~ I + β * S * I - α * I, D(R) ~ R + α * I]
-@named sir = DiscreteSystem(eqs)
+dds = @DDSmodel(
+    S(t + 1) = S(t) - β * S(t) * I(t),
+    I(t + 1) = I(t) + β * S(t) * I(t) - α * I(t),
+    R(t + 1) = R(t) + α * I(t),
+    y(t) = I(t)
+)
 ```
 
-Once the model is defined, we can assess identifiability by providing the formula for the observable:
+Then local identifiability can be assessed using `assess_local_identifiability` function:
 
-```@example discrete
-assess_local_identifiability(sir; measured_quantities = [y ~ I])
+```@example discrete_dds
+assess_local_identifiability(dds)
 ```
 
 For each parameter or state, the value in the returned dictionary is `true` (`1`) if the parameter is locally identifiable and `false` (`0`) otherwise.
 We see that `R(t)` is not identifiable, which makes sense: it does not affect the dynamics of the observable in any way.
 
-In principle, it is not required to give a name to the observable, so one can write this shorter
-
-```@example discrete
-assess_local_identifiability(sir; measured_quantities = [I])
-```
-
 The `assess_local_identifiability` function has three important keyword arguments:
 
   - `funcs_to_check` is a list of functions for which one want to assess identifiability, for example, the following code
     will check if `β * S` is locally identifiable.
 
-```@example discrete
-assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [β * S])
+```@example discrete_dds
+assess_local_identifiability(dds; funcs_to_check = [β * S])
 ```
 
   - `prob_threshold` is the probability of correctness (default value `0.99`, i.e., 99%). The underlying algorithm is a Monte-Carlo algorithm, so in
@@ -77,6 +87,39 @@ assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [
 As other main functions in the package, `assess_local_identifiability` accepts an optional parameter `loglevel` (default: `Logging.Info`)
 to adjust the verbosity of logging.
 
+If one loads `ModelingToolkit` (and thus the `ModelingToolkitExt` extension), one can use `DiscreteSystem` from `ModelingToolkit` to
+describe the input model (now in terms of difference!):
+
+```@example discrete_mtk
+using ModelingToolkit
+using StructuralIdentifiability
+
+@parameters α β
+@variables t S(t) I(t) R(t) y(t)
+D = Difference(t; dt = 1.0)
+
+eqs = [D(S) ~ S - β * S * I, D(I) ~ I + β * S * I - α * I, D(R) ~ R + α * I]
+@named sir = DiscreteSystem(eqs)
+```
+
+Then the same computation can be carried out with the models defined this way:
+
+```@example discrete_mtk
+assess_local_identifiability(sir; measured_quantities = [y ~ I])
+```
+
+In principle, it is not required to give a name to the observable, so one can write this shorter
+
+```@example discrete_mtk
+assess_local_identifiability(sir; measured_quantities = [I])
+```
+
+The same example but with specified functions to check
+
+```@example discrete_mtk
+assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [β * S])
+```
+
 The implementation is based on a version of the observability rank criterion and will be described in a forthcoming paper.
 
 [^1]: > S. Nõmm, C. Moog, [*Identifiability of discrete-time nonlinear systems*](https://doi.org/10.1016/S1474-6670(17)31245-4), IFAC Proceedings Volumes, 2004.
diff --git a/src/discrete.jl b/src/discrete.jl
index 4b96e6123..a25004994 100644
--- a/src/discrete.jl
+++ b/src/discrete.jl
@@ -391,6 +391,9 @@ function assess_local_identifiability(
     restart_logging(loglevel = loglevel)
     reset_timings()
     with_logger(_si_logger[]) do
+        if isempty(funcs_to_check)
+            funcs_to_check = vcat(parameters(dds), x_vars(dds))
+        end
         return _assess_local_identifiability_discrete_aux(
             dds,
             funcs_to_check,