diff --git a/docs/pages.jl b/docs/pages.jl
index caef31720..2b809b7d2 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -2,11 +2,8 @@ pages = [
     "Home" => "index.md",
     "Tutorials" => Any[
         "tutorials/creating_ode.md",
-        "tutorials/local_identifiability.md",
-        "tutorials/global_identifiability.md",
+        "tutorials/identifiability.md",
         "tutorials/identifiable_functions.md",
-        "tutorials/reparametrization.md",
-        "tutorials/using_modeling_toolkit.md",
     ],
     "Basics" => Any["input/input.md", "identifiability/identifiability.md"],
     "Library" => Any[
diff --git a/docs/src/tutorials/global_identifiability.md b/docs/src/tutorials/global_identifiability.md
deleted file mode 100644
index 0f96d52ba..000000000
--- a/docs/src/tutorials/global_identifiability.md
+++ /dev/null
@@ -1,44 +0,0 @@
-# Global Identifiability of Differential Models
-
-In this tutorial, let us cover an example problem of querying the ODE for globally identifiable parameters.
-
-## Input System
-
-Let us consider the following four-dimensional model with two outputs:
-
-$\begin{cases}x'(t) = lm - d \, x(t) - \beta \, x(t) \, v(t),\\
-y'(t) = \beta \, x(t) \, v(t) - a \, y(t),\\
-v'(t) = k \, y(t) - u \, v(t),\\
-w'(t) = c \, x(t) \, y(t) \, w(t) - c \, q \, y(t) \, w(t) - b \, w(t),\\
-z'(t) = c \, q \, y(t) \, w(t) - h \, z(t),\\
-y_1(t) = w(t),\\
-y_2(t) = z(t)\end{cases}$
-
-This model describes HIV dynamics[^1]. Let us run a global identifiability check on this model to get the result with probability of correctness being `p=0.99`. To do this, we will use `assess_identifiability(ode, p)` function.
-
-Global identifiability needs information about local identifiability first, hence the function we chose here will take care of that extra step for us.
-
-```@repl
-using StructuralIdentifiability
-
-ode = @ODEmodel(
-    x'(t) = lm - d * x(t) - beta * x(t) * v(t),
-    y'(t) = beta * x(t) * v(t) - a * y(t),
-    v'(t) = k * y(t) - u * v(t),
-    w'(t) = c * x(t) * y(t) * w(t) - c * q * y(t) * w(t) - b * w(t),
-    z'(t) = c * q * y(t) * w(t) - h * z(t),
-    y1(t) = w(t),
-    y2(t) = z(t)
-)
-global_id = assess_identifiability(ode, 0.99)
-```
-
-We also note that it's usually inexpensive to obtain the result with higher probability of correctness. For example, taking `p=0.9999` in the system above will result only in a slight slowdown.
-
-## Note on the probability of correctness
-
-Currently, the probability of correctness does not include the probability of correctness of the modular reconstruction for Groebner bases.
-This probability is ensured by an additional check modulo a large prime, and can be neglected for practical purposes. However, in the future versions, we plan to
-eliminate this possible error.
-
-[^1]: > D. Wodarz, M. Nowak, [*Specific therapy regimes could lead to long-term immunological control of HIV*](https://doi.org/10.1073/pnas.96.25.14464), PNAS December 7, 1999 96 (25) 14464-14469;
diff --git a/docs/src/tutorials/identifiability.md b/docs/src/tutorials/identifiability.md
new file mode 100644
index 000000000..4d6339adf
--- /dev/null
+++ b/docs/src/tutorials/identifiability.md
@@ -0,0 +1,125 @@
+# Identifiability of Differential Models (Local and Global)
+
+Recall that we consider ODE models in the state-space form
+
+$\begin{cases}
+\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}$
+
+where $\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively,
+and $\mathbf{p}$ are scalar parameters.
+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.
+Typically, two types of identifiability are distinguished
+
+  - **local** identifiability: the value can be recovered up to finitely many options;
+
+  - **global** identifiability: the value can be recoevered uniquely.
+
+Note that in the case of states, term **observability** it typically used. We will use **identifiability** for both
+states and parameters for brevity and uniformity.
+While the notion of global identifiability is much more precise, assessing local identifiability is typically much faster,
+and can be performed for the models whose global identifiability analysis is out of reach.
+
+## Local identifiability
+
+We consider the Lotka-Volterra model:
+
+$\begin{cases}
+x_1'(t) = a x_1(t) - b x_1(t) x_2(t) + u(t),\\
+x_2'(t) = -c x_2(t) + d x_1(t) x_2(t),\\
+y(t) = x_1(t)
+\end{cases}$
+
+The local identifiability of all parameters and states in this model can be assessed as follows
+
+```@example local
+using StructuralIdentifiability
+
+ode = @ODEmodel(
+    x1'(t) = a * x1(t) - b * x1(t) * x2(t) + u(t),
+    x2'(t) = -c * x2(t) + d * x1(t) * x2(t),
+    y(t) = x1(t)
+)
+
+assess_local_identifiability(ode)
+```
+
+We see that $x_1(t)$ is locally identifiable (no surprises, this is an output), $a, c,$ and $d$ are identifiable as well.
+On the other hand, $x_2(t)$ and $b$ are nonidentifiable. This can be explained by the following scaling symmetry
+
+$x_2(t) \to \lambda x_2(t), \quad b \to \frac{b}{\lambda}$
+
+which preserves input and output for every nonzero $\lambda$.
+The algorithm behind this call is the one due to Sedoglavic[^1].
+
+Function `assess_local_identifiability` has several optional parameters
+
+  - `funcs_to_check` a list of specific functions of parameters and states to check identifiability for (see an example below).
+    If not provided, the identifiability is assessed for all parameters and states.
+
+  - `p` (default $0.99$) is the probability of correctness. The algorithm can, in theory, produce wrong result, but the probability that it is correct
+    is guaranteed to be at least `p`. However, the probability bounds we use are quite conservative, so the actual probability of correctness is
+    likely to be much higher.
+  - `type` (default `:SE`). By default, the algorithm checks the standard single-experiment identifiability. If one sets `type = :ME`, then the algorithm
+    checks multi-experiment identifiability, that is, identifiability from several experiments with independent initial conditions (the algorithm from [^2] is used).
+
+Note that the scaling symmetry given above suggests that $b x_2(t)$ may in fact be identifiable. This can be checked using `funcs_to_check` parameter:
+
+```@example local
+assess_local_identifiability(ode, funcs_to_check = [b * x2])
+```
+
+Indeed!
+
+## Global identifiability
+
+One can obtain more refined information about a model using `assess_identifiability` function.
+We will showcase it using the Goodwin oscillator model [^3].
+
+```@example global
+using StructuralIdentifiability
+
+ode = @ODEmodel(
+    x1'(t) = -b * x1(t) + 1 / (c + x4(t)),
+    x2'(t) = alpha * x1(t) - beta * x2(t),
+    x3'(t) = gama * x2(t) - delta * x3(t),
+    x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t),
+    y(t) = x1(t)
+)
+
+assess_identifiability(ode)
+```
+
+As a result, each parameter/state is assigned one of the labels `:globally` (globally identifiable), `:locally` (locally but not globally identifiable),
+or `:nonidentifiable` (not identifiable, even locally).
+The algorithm behind this computation follows [^4].
+
+Similarly to `assess_local_identifiability`, this function has optional parameters:
+
+  - `funcs_to_check` a list of specific functions of parameters and states to check identifiability for (see an example below).
+    If not provided, the identifiability is assessed for all parameters and states. Note that the computations for states may be
+    more involved than for the parameters, so one may want to call the function with `funcs_to_check = ode.parameters` if the
+    call `assess_identifiability(ode)` takes too long.
+
+  - `p` (default $0.99$) is the probability of correctness. Same story as above: the probability estimates are very conservative, so the actual
+    error probability is much lower than 1%.
+    Also, currently, the probability of correctness does not include the probability of correctness of the modular reconstruction for Groebner bases.
+    This probability is ensured by an additional check modulo a large prime, and can be neglected for practical purposes.
+
+Using `funcs_to_check` parameter, one can further inverstigate the nature of the lack of identifiability in the model at hand.
+For example, for the Goodwin oscillator, we can check if `beta + delta` and `beta * delta` are identifiable:
+
+```@example global
+assess_identifiability(ode, funcs_to_check = [beta + delta, beta * delta])
+```
+
+And we see that they indeed are. This means, in particular, that the reason why `beta` and `delta` are not identifiable is because their values
+can be exchanged. One may wonder how could we guess these functions `beta + delta, beta * delta`. In fact, they can be just computed using
+`find_identifiable_functions` function as we will explain in the next tutorial. Stay tuned!
+
+[^1]: > A. Sedoglavic, [*A probabilistic algorithm to test local algebraic observability in polynomial time*](https://doi.org/10.1006/jsco.2002.0532), Journal of Symbolic Computation, 2002.
+[^2]: > A. Ovchinnikov, A. Pillay, G. Pogudin, T. Scanlon, [*Multi-experiment Parameter Identifiability of ODEs and Model Theory*](https://doi.org/10.1137/21M1389845), SIAM Journal on Applied Algebra and Geometry, 2022.
+[^3]: > D. Gonze, P. Ruoff, [*The Goodwin Oscillator and its Legacy*](https://doi.org/10.1007/s10441-020-09379-8), Acta Biotheoretica, 2020.
+[^4]: > R. Dong, C. Goodbrake, H. Harrington, G. Pogudin, [*Differential elimination for dynamical models via projections with applications to structural identifiability*](https://doi.org/10.1137/22M1469067), SIAM Journal on Applied Algebra and Geometry, 2023.
diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md
index fcac01d6b..db98e6043 100644
--- a/docs/src/tutorials/identifiable_functions.md
+++ b/docs/src/tutorials/identifiable_functions.md
@@ -1,10 +1,13 @@
 # Globally Identifiable Functions
 
-StructuralIdentifiability.jl provides the function `find_identifiable_functions` that returns the identifiable combinations of the ODE variables.
+In addition to assessing identifiabuility of given functions of parameters or states, `StructuralIdentifiability.jl`
+provides the function `find_identifiable_functions` which finds a set of identifiable functions such that any other
+identifiable function can be expressed using them.
+This allows to find out what actually is identifiable and what does non-identifiability in the model at hand looks like.
 
 For example, consider the following model[^1].
 
-```@example
+```@example funcs
 using StructuralIdentifiability # hide
 LLW1987 = @ODEmodel(
     x1'(t) = -p1 * x1(t) + p2 * u(t),
@@ -14,19 +17,27 @@ LLW1987 = @ODEmodel(
 )
 ```
 
-Several decades ago, this model was introduced to demonstrate the presence of nontrivial **un**identifiability in nonlinear systems of ODEs. Nowadays, we can automatically find the identifiable combinations of parameters:
+Several decades ago, this model was introduced to demonstrate the presence of nontrivial **un**identifiability in nonlinear systems of ODEs.
+Nowadays, we can automatically find the identifiable combinations of parameters:
 
-```@example
+```@example funcs
 using StructuralIdentifiability # hide
 find_identifiable_functions(LLW1987)
 ```
 
-And even of parameters and states:
+From these expressions, we see that the values of `p1` and `p3` are not identifiable but an unordered pair
+of numbers `{p1, p3}` is uniquely determined since `p1 + p3` and `p1 * p3` are known.
+Furthermore, we see that, for fixed input and output, `p2` and `p4` can take infinitely many values but
+knowing one of them, we would also be able to determine the other.
 
-```@example
+Moreover, we can find generators of all identifiable functions in parameters and states:
+
+```@example funcs
 find_identifiable_functions(LLW1987, with_states = true)
 ```
 
-By default, `find_identifiable_functions` tries to simplify the output functions as much as possible. This feature is useful, for example, in [Model Reparametrization](@ref).
+By default, `find_identifiable_functions` tries to simplify the output functions as much as possible, and it has `simplify` keyword responsible for
+the degree of simplification. The default value is `:standard` but one could use `:string` to try to simplify further
+(at the expense of heavier computation) or use `:weak` to simplify less (but compute faster).
 
 [^1]: > Y. Lecourtier, F. Lamnabhi-Lagarrigue, and E. Walter, [*A method to prove that nonlinear models can be unidentifiable*](https://doi.org/10.1109/CDC.1987.272467), Proceedings of the 26th Conference on Decision and Control, December 1987, 2144-2145;
diff --git a/docs/src/tutorials/local_identifiability.md b/docs/src/tutorials/local_identifiability.md
deleted file mode 100644
index a8511c046..000000000
--- a/docs/src/tutorials/local_identifiability.md
+++ /dev/null
@@ -1,60 +0,0 @@
-# Local Identifiability of Differential Models
-
-In this tutorial, we will go over an example of solving a local identifiability problem for a simple system of ordinary differential equations.
-
-We will introduce how to use the input parsing in `StructuralIdentifiability.jl` and the local identifiability assessment functionality.
-
-## Input System
-
-We will consider a simple two-species competition model
-
-$x'_1 = k \,(1 - x_1 - x_2)\\ x'_2=r\,(1-x_1-x_2).$
-
-To make it a proper input for the algorithm, we add an output function $y=x_1$ that equals to the population density of species 1 at any time $t$.
-
-### Using the `@ODEmodel` macro
-
-To parse the system of ordinary differential equations as above, we will use `@ODEmodel` macro. This is the easiest way to do so.
-
-We have two state variables `x1, x2` (population densities), two parameters `k, r` (intrinsic growth rates), and one output function `y`. Note that there must be `(t)` to indicate time-dependent functions.
-
-After using the macro, we use the `assess_local_identifiability` function for that. This function accepts the ODE model, the probability of correctness, and the type of identifiability we would like to inquire about.
-
-```@example
-using StructuralIdentifiability
-
-ode = @ODEmodel(
-    x1'(t) = k * (1 - x1(t) - x2(t)),
-    x2'(t) = r * (1 - x1(t) - x2(t)),
-    y(t) = x1(t)
-)
-
-local_id = assess_local_identifiability(ode, 0.99)
-```
-
-The result shows that only the state variable's initial value $x'_1(0)$ is locally identifiable.
-
-Let us now add another output function `y2(t)`:
-
-```@example
-using StructuralIdentifiability
-
-ode = @ODEmodel(
-    x1'(t) = k * (1 - x1(t) - x2(t)),
-    x2'(t) = r * (1 - x1(t) - x2(t)),
-    y1(t) = x1(t),
-    y2(t) = x2(t) # new output function!
-)
-
-local_id = assess_local_identifiability(ode, 0.99) # this is a different result!
-```
-
-As you can see, for this new model with an additional output, all parameters are reported as locally identifiable with probability 0.99.
-
-## Note on Probability of Correctness
-
-We set the probability of correctness $p$ to be `0.99`. Why would we ever want a lower value? Great question! The underlying algorithm relies on operations being modulo a large enough prime characteristic $\mathcal{P}\geq \kappa p$ where $\kappa$ is determined by the algorithm internally.
-
-The algorithm's complexity is proportional to the size of operands (see proposition 3.1 in the main paper[^1]) and high probability of correctness may thus lead to higher size of coefficients during computation for some systems. Hence, one may wish to lower $p$ to save on runtime (though in practice this is _very_ rare).
-
-[^1]: > A. Sedoglavic, [*A probabilistic algorithm to test local algebraic observability in polynomial time*](https://doi.org/10.1006/jsco.2002.0532), Journal of Symbolic Computation, 2002.
diff --git a/docs/src/tutorials/reparametrization.md b/docs/src/tutorials/reparametrization.md
deleted file mode 100644
index 769f5d6bc..000000000
--- a/docs/src/tutorials/reparametrization.md
+++ /dev/null
@@ -1,5 +0,0 @@
-# Model Reparametrization
-
-*This functionality is experimental*
-
-TODO: a couple of  examples
diff --git a/docs/src/tutorials/using_modeling_toolkit.md b/docs/src/tutorials/using_modeling_toolkit.md
deleted file mode 100644
index 2c3dc663c..000000000
--- a/docs/src/tutorials/using_modeling_toolkit.md
+++ /dev/null
@@ -1,88 +0,0 @@
-# Using ModelingToolkit.jl With StructuralIdentifiability.jl
-
-In this tutorial, we will cover examples of solving identifiability problems for models defined with the syntax of `ModelingToolkit.jl`.
-
-## Input System
-
-Let us consider the following ODE model with two outputs:
-
-$\begin{cases}
-\dot{S} = -b \, S \, (I + J + q \, A) \, N_{inv},\\
-\dot{E} = b \, S \, (I + J + q \, A) \, N_{inv} - k \, E,\\
-\dot{A} = k \, (1 - r) \, E - g_1 \, A,\\
-\dot{I} = k \, r \, E - (\alpha + g_1) \, I,\\
-\dot{J} = \alpha \, I - g_2 \, J,\\
-\dot{C} = \alpha \, I,\\
-y_1 = C,\\
-y_2 = N_{inv}
-\end{cases}$
-
-This is an infectious disease model defined in [^1].
-
-The main difference between the input formats in `ModelingToolkit.jl` and `StructuralIdentifiability.jl` is that the output (measured values/functions) must be specified separately in `ModelingToolkit.jl`. In this example, measured quantities are presented by $y_1$, $y_2$.
-
-First, let us define the ODE. We will use `@parameters` and `@variables` macro to define parameters and time-depended functions in the ODE.
-
-```julia
-using StructuralIdentifiability, ModelingToolkit
-
-@parameters b q N_inv k r alpha g1 g2
-@variables t S(t) E(t) A(t) I(t) J(t) C(t) y1(t) y2(t)
-```
-
-The actual ODE will be defined using `ODESystem` structure from `ModelingToolkit.jl`:
-
-```julia
-D = Differential(t)
-
-eqs = [
-    D(S) ~ -b * S * (I + J + q * A) * N_inv,
-    D(E) ~ b * S * (I + J + q * A) * N_inv - k * E,
-    D(A) ~ k * (1 - r) * E - g1 * A,
-    D(I) ~ k * r * E - (alpha + g1) * I,
-    D(J) ~ alpha * I - g2 * J,
-    D(C) ~ alpha * I,
-]
-
-ode = ODESystem(eqs, t, name = :SEIAJRCmodel)
-```
-
-Finally, let us define the array of measured quantities and call the `assess_identifiability` function. This is the main function that determines local/global identifiability properties of each parameter and state. We will use the probability of correctness $p=0.99$.
-
-For `ModelingToolkit.jl`, both `assess_identifiability` and `assess_local_identifiability` functions accept keyword arguments:
-
-  - `measured_quantities`, also called “output functions” in identifiability literature; these are crucial for answering identifiability questions.
-  - `p`, probability of correctness. This value equals 0.99 by default.
-  - `funcs_to_check`, functions of parameters of which we wish to check identifiability.
-
-```julia
-measured_quantities = [y1 ~ C, y2 ~ N_inv]
-@time global_id = assess_identifiability(ode, measured_quantities = measured_quantities)
-```
-
-Let us put all the code above together:
-
-```@repl
-using StructuralIdentifiability, ModelingToolkit
-
-@parameters b q N_inv k r alpha g1 g2
-@variables t S(t) E(t) A(t) I(t) J(t) C(t) y1(t) y2(t)
-
-D = Differential(t)
-
-eqs = [
-    D(S) ~ -b * S * (I + J + q * A) * N_inv,
-    D(E) ~ b * S * (I + J + q * A) * N_inv - k * E,
-    D(A) ~ k * (1 - r) * E - g1 * A,
-    D(I) ~ k * r * E - (alpha + g1) * I,
-    D(J) ~ alpha * I - g2 * J,
-    D(C) ~ alpha * I,
-]
-
-ode = ODESystem(eqs, t, name = :SEIAJRCmodel)
-
-measured_quantities = [y1 ~ C, y2 ~ N_inv]
-@time global_id = assess_identifiability(ode, measured_quantities = measured_quantities)
-```
-
-[^1]: > K. Roosa and G. Chowell. [*Assessing parameter identifiability in compartmental dynamic models using a computational approach: application to infectious disease transmission models*](https://doi.org/10.1186/s12976-018-0097-6), Theor Biol Med Model 16, 1 (2019)
diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl
index d0b49995b..8632868b7 100644
--- a/src/StructuralIdentifiability.jl
+++ b/src/StructuralIdentifiability.jl
@@ -82,47 +82,40 @@ include("submodels.jl")
 include("discrete.jl")
 
 """
-    assess_identifiability(ode::ODE{P}, p::Float64=0.99) where P <: MPolyElem{fmpq}
-
-Input:
-- `ode` - the ODE model
-- `p` - probability of correctness.
-
-Assesses identifiability (both local and global) of a given ODE model (parameters detected automatically). The result is guaranteed to be correct with the probability
-at least `p`.
-"""
-function assess_identifiability(ode::ODE{P}, p::Float64 = 0.99) where {P <: MPolyElem{fmpq}}
-    result = assess_identifiability(ode, vcat(ode.parameters, ode.x_vars), p)
-    return result #Dict(param => res for (param, res) in zip(ode.parameters, result_list))
-end
-
-"""
-    assess_identifiability(ode, [funcs_to_check, p=0.99])
+    assess_identifiability(ode; funcs_to_check = [], p=0.99)
 
 Input:
 - `ode` - the ODE model
+- `funcs_to_check` - list of functions to check identifiability for; if empty, all parameters
+   and states are taken
 - `p` - probability of correctness.
     
 Assesses identifiability of a given ODE model. The result is guaranteed to be correct with the probability
 at least `p`.
-
-If `funcs_to_check` are given, then the function will assess the identifiability of the provided functions.
-Otherwise, it will assess identifiability of all the parameters and states.
 The function returns a dictionary from the functions to check to their identifiability properties 
 (one of `:nonidentifiable`, `:locally`, `:globally`).
 """
 function assess_identifiability(
-    ode::ODE{P},
-    funcs_to_check::Array{<:RingElem, 1},
+    ode::ODE{P};
+    funcs_to_check::Array{<:RingElem, 1} = Array{RingElem, 1}(),
     p::Float64 = 0.99,
 ) where {P <: MPolyElem{fmpq}}
     p_glob = 1 - (1 - p) * 0.9
     p_loc = 1 - (1 - p) * 0.1
 
+    if isempty(funcs_to_check)
+        funcs_to_check = vcat(ode.parameters, ode.x_vars)
+    end
+
     @info "Assessing local identifiability"
     trbasis = Array{fmpq_mpoly, 1}()
-    runtime = @elapsed local_result =
-        assess_local_identifiability(ode, funcs_to_check, p_loc, :SE, trbasis)
+    runtime = @elapsed local_result = assess_local_identifiability(
+        ode,
+        funcs_to_check = funcs_to_check,
+        p = p_loc,
+        type = :SE,
+        trbasis = trbasis,
+    )
     @info "Local identifiability assessed in $runtime seconds"
     @debug "Trasncendence basis to be specialized is $trbasis"
     _runtime_logger[:loc_time] = runtime
@@ -188,7 +181,7 @@ function assess_identifiability(
     end
     funcs_to_check_ = [eval_at_nemo(each, conversion) for each in funcs_to_check]
 
-    result = assess_identifiability(ode, funcs_to_check_, p)
+    result = assess_identifiability(ode, funcs_to_check = funcs_to_check_, p = p)
     nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check)
     out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_)
     return out_dict
diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl
index 91f4a60c2..c99baf6ba 100644
--- a/src/local_identifiability.jl
+++ b/src/local_identifiability.jl
@@ -195,55 +195,31 @@ function assess_local_identifiability(
     funcs_to_check_ = [eval_at_nemo(x, conversion) for x in funcs_to_check]
 
     if isequal(type, :SE)
-        result = assess_local_identifiability(ode, funcs_to_check_, p, type)
+        result = assess_local_identifiability(
+            ode,
+            funcs_to_check = funcs_to_check_,
+            p = p,
+            type = type,
+        )
         nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check)
         out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_)
         return out_dict
     elseif isequal(type, :ME)
-        result, bd = assess_local_identifiability(ode, funcs_to_check_, p, type)
+        result, bd = assess_local_identifiability(
+            ode,
+            funcs_to_check = funcs_to_check_,
+            p = p,
+            type = type,
+        )
         nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check)
         out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_)
         return (out_dict, bd)
     end
 end
 # ------------------------------------------------------------------------------
-"""
-    assess_local_identifiability(ode::ODE{P}, p::Float64 = 0.99, type=:SE) where P <: MPolyElem{Nemo.fmpq}
-
-Input:
-- `ode` - the ODE model
-- `p` - probability of correctness
-- `type` - identifiability type (`:SE` for single-experiment, `:ME` for multi-experiment)
-
-Output:
-- for `type=:SE`, the result is a dictionary from each parameter to boolean;
-- for `type=:ME`, the result is a tuple with the dictionary as in `:SE` case and array of number of experiments.
-
-The main entry point for local identifiability checks.
-Call this function to automatically take care of local identifiability of all parameters and initial conditions.
-The result is correct with probability at least `p`.
 
-`type` can be either `:SE` (single-experiment identifiability) or `:ME` (multi-experiment identifiability).
-The return value is a tuple consisting of the array of bools and the number of experiments to be performed.
 """
-function assess_local_identifiability(
-    ode::ODE{P},
-    p::Float64 = 0.99,
-    type = :SE,
-) where {P <: MPolyElem{Nemo.fmpq}}
-    funcs_to_check = ode.parameters
-    if type == :SE
-        funcs_to_check = vcat(funcs_to_check, ode.x_vars)
-    end
-    result = assess_local_identifiability(ode, funcs_to_check, p, type)
-    if type == :SE
-        return Dict(a => result[a] for a in funcs_to_check)
-    end
-    return (Dict(a => result[1][a] for a in funcs_to_check), result[2])
-end
-
-"""
-    assess_local_identifiability(ode::ODE{P}, funcs_to_check::Array{<: Any, 1}, p::Float64=0.99, type=:SE) where P <: MPolyElem{Nemo.fmpq}
+    assess_local_identifiability(ode::ODE{P}; funcs_to_check::Array{<: Any, 1}, p::Float64=0.99, type=:SE) where P <: MPolyElem{Nemo.fmpq}
 
 Checks the local identifiability/observability of the functions in `funcs_to_check`. The result is correct with probability at least `p`.
 
@@ -253,12 +229,18 @@ Call this function if you have a specific collection of parameters of which you
 If the type is `:ME`, states are not allowed to appear in the `funcs_to_check`.
 """
 function assess_local_identifiability(
-    ode::ODE{P},
-    funcs_to_check::Array{<:Any, 1},
+    ode::ODE{P};
+    funcs_to_check::Array{<:Any, 1} = Array{Any, 1}(),
     p::Float64 = 0.99,
     type = :SE,
     trbasis = nothing,
 ) where {P <: MPolyElem{Nemo.fmpq}}
+    if isempty(funcs_to_check)
+        funcs_to_check = ode.parameters
+        if type == :SE
+            funcs_to_check = vcat(funcs_to_check, ode.x_vars)
+        end
+    end
 
     # Checking whether the states appear in the ME case
     if type == :ME
diff --git a/test/identifiability.jl b/test/identifiability.jl
index 2dda61295..aa992907e 100644
--- a/test/identifiability.jl
+++ b/test/identifiability.jl
@@ -209,7 +209,7 @@
     #-------------------------------------------------------------------------- 
 
     for case in test_cases
-        result = assess_identifiability(case[:ode], case[:funcs])
+        result = assess_identifiability(case[:ode], funcs_to_check = case[:funcs])
         @test result == case[:correct]
     end
 end
diff --git a/test/linear_compartment.jl b/test/linear_compartment.jl
index 01e2cfe34..f52fad656 100644
--- a/test/linear_compartment.jl
+++ b/test/linear_compartment.jl
@@ -161,7 +161,7 @@
         for (e, id) in case[:result]
             correct[str_to_var("a_$(e[2])_$(e[1])", bring)] = id
         end
-        result = assess_identifiability(ode, collect(keys(correct)))
+        result = assess_identifiability(ode, funcs_to_check = collect(keys(correct)))
         @test correct == result
     end
 end
diff --git a/test/local_identifiability.jl b/test/local_identifiability.jl
index b2e90602f..976583097 100644
--- a/test/local_identifiability.jl
+++ b/test/local_identifiability.jl
@@ -117,10 +117,14 @@
     for case in test_cases
         trbasis = []
         ode = case[:ode]
-        result = assess_local_identifiability(ode, case[:funcs], 0.99, :SE, trbasis)
+        result = assess_local_identifiability(
+            ode,
+            funcs_to_check = case[:funcs],
+            trbasis = trbasis,
+        )
         @test result == case[:correct]
         for (i, p) in enumerate(trbasis)
-            res_for_p = assess_local_identifiability(ode, [p])
+            res_for_p = assess_local_identifiability(ode, funcs_to_check = [p])
             @test !first(values(res_for_p))
             ode = add_outputs(ode, Dict("YYY$i" => p))
         end
diff --git a/test/local_identifiability_me.jl b/test/local_identifiability_me.jl
index 1ca946683..6c6541473 100644
--- a/test/local_identifiability_me.jl
+++ b/test/local_identifiability_me.jl
@@ -136,7 +136,7 @@ end
         for n in n_min:n_max
             model = _linear_compartment_model(case[:graph](n), [1])
             println(case[:name] * ", n = $n")
-            @time result = assess_local_identifiability(model, 0.97, :ME)
+            @time result = assess_local_identifiability(model, type = :ME)
             correct = undef
             if n - n_min + 1 > length(case[:bound])
                 correct = case[:bound][end]
@@ -210,7 +210,12 @@ end
     #--------------------------------------------------------------------------
 
     for case in test_cases
-        result = assess_local_identifiability(case[:ode], case[:funcs], 0.932, :ME)
+        result = assess_local_identifiability(
+            case[:ode],
+            funcs_to_check = case[:funcs],
+            p = 0.932,
+            type = :ME,
+        )
         @test result == case[:correct]
     end
 end