diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index b4d5c7ae8..ed4e5a4db 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -161,7 +161,6 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/julia-buildpkg@v1
       - name: Load Julia packages from cache
         id: julia-cache
         uses: julia-actions/cache@v2
@@ -202,6 +201,9 @@ jobs:
             using Pkg
             using Test
             pkg"dev ."
+            # pkg"update"
+            pkg"instantiate"
+            pkg"build"
             @testset "PowerDynamics Downstream Tests" begin
                 @testset "Normal Package tests" begin
                     Pkg.test("PowerDynamics"; coverage=false)
diff --git a/NEWS.md b/NEWS.md
index f5ffb4bc1..f51335213 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,8 +1,16 @@
 # NetworkDynamics Release Notes
 
+## v0.10.4 Changelog
+- [#303](https://github.com/JuliaDynamics/NetworkDynamics.jl/pull/303) update for ModelingToolkit.jl v10 compatibility:
+  - rename all `ODESystem` -> `System` (follows MTK v10 API)
+  - MTK extension now uses `mtkcompile` instead of `structural_simplify` internally
+  - Add new `implicit_output` function to handle fully implicit output variables in MTK models
+  - Add documentation for handling fully implicit outputs in MTK integration
+  - Update minimum ModelingToolkit.jl requirement from v9.67 to v10
+
 ## v0.10.3 Changelog
 - [#301](https://github.com/JuliaDynamics/NetworkDynamics.jl/pull/301) improve callback system performance and flexibility:
-  - Add callback batching for better DiscreteComponentCallback performance 
+  - Add callback batching for better DiscreteComponentCallback performance
   - Allow `EIndex(1=>2)` as standalone edge index with relaxed type constraints
   - Optimize CallbackSet construction to prevent performance bottlenecks
   - Add important documentation warning about parameter array copying in callbacks
diff --git a/Project.toml b/Project.toml
index c3167f913..834219b14 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "NetworkDynamics"
 uuid = "22e9dc34-2a0d-11e9-0de0-8588d035468b"
 authors = ["Frank Hellmann <hellmann@pik-potsdam.de>, Michael Lindner <michaellindner@pik-potsdam.de>, Hans Würfel <git@wuerfel.io"]
-version = "0.10.3"
+version = "0.10.4"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
@@ -67,7 +67,7 @@ KernelAbstractions = "0.9.18"
 LinearAlgebra = "1"
 MacroTools = "0.5.15"
 Mixers = "0.1.2"
-ModelingToolkit = "9.82"
+ModelingToolkit = "10"
 NNlib = "0.9.13"
 NonlinearSolve = "4"
 OrderedCollections = "1.8.0"
diff --git a/docs/examples/gas_network.jl b/docs/examples/gas_network.jl
index f2af2c40b..cc789f839 100644
--- a/docs/examples/gas_network.jl
+++ b/docs/examples/gas_network.jl
@@ -233,7 +233,7 @@ As a workaround we had to explicitly define `LinearInterpolations` as unitless,
 
 To build the network, we first need to define the components. This is a two-step process:
 
-- first create the symbolic `ODESystem` using ModelingToolkit
+- first create the symbolic `System` using ModelingToolkit
 - secondly build a NetworkDynamics component model ([`VertexModel`](@ref)/[`EdgeModel`](@ref)) based on the symbolic system.
 
 In the first step we can use the keyword arguments to pass "default" values for our parameters and states.
diff --git a/docs/examples/init_tutorial.jl b/docs/examples/init_tutorial.jl
index e9cced6f5..9d0b65bb7 100644
--- a/docs/examples/init_tutorial.jl
+++ b/docs/examples/init_tutorial.jl
@@ -61,7 +61,11 @@ end
 nothing #hide
 
 #=
-**B) A static prosumer node which forces a certain flow (pressure is fully implicit)**
+**B) A static prosumer node which forces a certain flow**
+
+!!! note "Fully Implicit Output"
+    We need to use `implicit_output(p)` to handle the fully implicit pressure
+    output. See [fully implicit outputs](@ref Fully-Implicit-Outputs) for details.
 =#
 @mtkmodel StaticProsumerNode begin
     @extend GasNode()
@@ -69,7 +73,7 @@ nothing #hide
         q̃_prosumer, [description="flow injected by prosumer"]
     end
     @equations begin
-        -q̃_nw ~ q̃_prosumer
+        -q̃_nw ~ q̃_prosumer + implicit_output(p)
     end
 end
 nothing #hide
diff --git a/docs/src/API.md b/docs/src/API.md
index 1837b5076..43f5c012e 100644
--- a/docs/src/API.md
+++ b/docs/src/API.md
@@ -18,9 +18,9 @@ EdgeModel()
 
 ## Component Models with MTK
 ```@docs
-VertexModel(::ModelingToolkit.ODESystem, ::Any, ::Any)
-EdgeModel(::ModelingToolkit.ODESystem, ::Any, ::Any, ::Any, ::Any)
-EdgeModel(::ModelingToolkit.ODESystem, ::Any, ::Any, ::Any)
+VertexModel(::ModelingToolkit.System, ::Any, ::Any)
+EdgeModel(::ModelingToolkit.System, ::Any, ::Any, ::Any, ::Any)
+EdgeModel(::ModelingToolkit.System, ::Any, ::Any, ::Any)
 ```
 
 ### Output Function Helpers/Wrappers
@@ -245,6 +245,7 @@ save_parameters!
 ff_to_constraint
 Base.copy(::NetworkDynamics.ComponentModel)
 extract_nw
+implicit_output
 ```
 
 ## NetworkDynamicsInspector API
diff --git a/docs/src/metadata.md b/docs/src/metadata.md
index ed2cf3f7a..21c9f32c8 100644
--- a/docs/src/metadata.md
+++ b/docs/src/metadata.md
@@ -33,7 +33,7 @@ Special cases for symbol metadata are:
 
 For those, there are special functions `has_*`, `get_*`, `set_*!`, `delete_*!` and `strip_*!`. The `strip_*!` functions remove all metadata of a specific type from all symbols in a component. See [Per Symbol Metadata API](@ref).
 
-These are closely aligned with the [metadata use in ModelingToolkit](@extref ModelingToolkit symbolic_metadata). They are automatically copied from the `ODESystem` if you use MTK models to create NetworkDynamics models.
+These are closely aligned with the [metadata use in ModelingToolkit](@extref ModelingToolkit symbolic_metadata). They are automatically copied from the `System` if you use MTK models to create NetworkDynamics models.
 
 ## Metadata Utils
 Accessing metadata (especially defaults) of states and parameters is a very
diff --git a/docs/src/mtk_integration.md b/docs/src/mtk_integration.md
index 77fc09f80..e2a547483 100644
--- a/docs/src/mtk_integration.md
+++ b/docs/src/mtk_integration.md
@@ -6,13 +6,13 @@ which are then connected on network scale using NetworkDynamics.
 
 The main entry point for this interop are the constructors
 ```julia
-VertexModel(::ODESystem, inputs, outputs)
-EdgeModel(::ODESystem, srcin, dstin, [srscout], dstout)
+VertexModel(::System, inputs, outputs)
+EdgeModel(::System, srcin, dstin, [srscout], dstout)
 ```
 whose docstrings can be found in the [Component Models with MTK](@ref) section in the API.
 
 These constructors will:
-- transform the states marked as input to parameters and `structural_simplify`ing the system,
+- transform the states marked as input to parameters and `mtkcompile`ing the system,
 - generate the `f` and `g` functions,
 - generate code for observables,
 - port all supported [Metadata](@ref) from MTK symbols to component symbols and
@@ -43,7 +43,7 @@ The system to model is 2 node, 1 edge network. The node output states are the vo
 
 ideal v source      Resistor     Capacitor
             v1 o─←────MMM────→─o v2
-               │               ┴ 
+               │               ┴
               (↗)              ┬
                │               │
                ⏚               ⏚
@@ -76,7 +76,7 @@ An ideal voltage source is just a model which pins its output voltage to a fixed
 The source ejects whatever current is necessary. We introduce another variable `i(t)`
 to "capture" this current. This variable will be removed during structural simplify, but will
 be available for plotting through the [Observables](@ref) mechanism.
-The `VertexModel` can be generated from an `ODESystem` by providing names of the input and output states:
+The `VertexModel` can be generated from an `System` by providing names of the input and output states:
 
 ```@example mtk
 @mtkmodel VoltageSource begin
@@ -140,7 +140,7 @@ resistor_edge = EdgeModel(resistor, [:src₊v], [:dst₊v], [:src₊i], [:dst₊
 
 Having all those components defined, we can build the network. We don't need to provide a graph
 object here because we specified the placement in the graph on a per component basis.
- 
+
 ```@example mtk
 nw = Network([vs_vertex, cap_vertex], [resistor_edge])
 ```
@@ -149,7 +149,7 @@ We can see, that NetworkDynamics internally is able to reduce all of the "output
 
 Now we can simulate the system. For that we generate the `u0` object. Since the metadata (such as default values) was automatically transferred, we can straight away construct the `ODEProblem`
 and solve the system.
- 
+
 ```@example mtk
 u0 = NWState(nw) # generate state based on default values
 prob = ODEProblem(nw, uflat(u0), (0, 10.0), pflat(u0))
@@ -164,3 +164,67 @@ axislegend(ax2)
 fig # hide
 ```
 
+## Fully Implicit Outputs
+When working with MTK systems in NetworkDynamics, you may encounter situations where
+your desired output variables don't explicitly appear in the equations. This creates **fully
+implicit outputs** - variables that are determined by the system's constraints but aren't
+directly computed.
+
+!!! tip "tl;dr"
+    Introduce "fake" dependencies to your input-forcing equations `0 ~ in + implicit_output(y)`.
+    Which is mathematically equivalent to `0 ~ in` but helps MTK to reason about dependencies.
+
+Consider a system with a fully implicit output:
+```
+   u┌───────┐y
+  ─→┤ 0 ~ u ├→─
+    └───────┘
+```
+Here, $y$ does not appear in the equations at all. In general, that doesn't make too much sense.
+During simplification, MTK will potentially get rid of the equation as it does not contribute to the system's state.
+
+However, in NetworkDynamics, we're always dealing with **open loop models** on the equation level, which is not exactly what MTK was made for.
+If you build a closed loop between a subsystem A which **has input forcing** and a subsystem
+B which has **input feed forward**, the resulting system can be solved:
+```
+    (system with input forcing)
+          ua┌─────────┐ya
+        ┌──→┤  0 ~ ua ├→──┐
+        │   └─────────┘   │
+        │ yb┌─────────┐ub │
+        └──←┤ yb ~ ub ├←──┘
+            └─────────┘
+(system with input feed forward)
+```
+
+Since MTK does not know about the closed loop (which is only introduced on the NetworkDynamics level once we leave the equation based domain) we need to help MTK to figure out those dependencies.
+We can do so by introducing "fake" dependencies using [`implicit_output`](@ref).
+This function is defined as
+```julia
+implicit_output(x) = 0
+ModelingToolkit.@register_symbolic implicit_output(x)
+```
+which makes it numerically equivalent to zero (no effect on the simulation) but is
+opaque to the Symbolic Simplification.
+
+### Example
+
+Consider a "Kirchhoff Node" between multiple resistors:
+- the currents through the resistors directly depend on the voltage output of the node (input feed forward) and
+- the Kirchhoff node requires the sum of all inflowing currents to be zero (input forcing).
+
+We can model this type of node like this:
+```@example mtk
+@mtkmodel KirchhoffNode begin
+    @variables begin
+        v(t), [description="Node voltage", output=true]
+        i_sum(t), [description="Sum of incoming currents", input=true]
+    end
+    @equations begin
+        0 ~ i_sum + implicit_output(v)  # Kirchhoff's current law
+    end
+end
+@named kirchhoff = KirchhoffNode()
+VertexModel(kirchhoff, [:i_sum], [:v])
+```
+where we "trick" MTK into believing that the input forcing equation depends on the output too.
diff --git a/ext/MTKExt_utils.jl b/ext/MTKExt_utils.jl
index 491d72892..0bea0a12e 100644
--- a/ext/MTKExt_utils.jl
+++ b/ext/MTKExt_utils.jl
@@ -152,7 +152,7 @@ function warn_missing_features(sys)
     cev = ModelingToolkit.get_continuous_events(sys)
     dev = ModelingToolkit.get_discrete_events(sys)
     if !isempty(cev) || !isempty(dev)
-        @warn "Model has attached events, which is not supportet."
+        @warn "Model has attached events, which is not supported."
     end
 
     if !isempty(ModelingToolkit.initialization_equations(sys))
@@ -215,3 +215,16 @@ function fix_metadata!(invalid_eqs, sys)
     end
     invalid_eqs .= fixedeqs
 end
+
+function remove_implicit_output_fn!(eqs)
+    r = SymbolicUtils.@rule implicit_output(~~x) => 0
+    chain = SymbolicUtils.Chain([r])
+    rewriter = SymbolicUtils.Prewalk(chain)
+
+    for i in eachindex(eqs)
+        eq = eqs[i]
+        eqs[i] = rewriter(eq.lhs) ~ rewriter(eq.rhs)
+    end
+
+    eqs
+end
diff --git a/ext/NetworkDynamicsMTKExt.jl b/ext/NetworkDynamicsMTKExt.jl
index 54a5fc11f..1853c87dd 100644
--- a/ext/NetworkDynamicsMTKExt.jl
+++ b/ext/NetworkDynamicsMTKExt.jl
@@ -1,14 +1,15 @@
 module NetworkDynamicsMTKExt
 
 using ModelingToolkit: Symbolic, iscall, operation, arguments, build_function
-using ModelingToolkit: ModelingToolkit, Equation, ODESystem, Differential
-using ModelingToolkit: equations, full_equations, get_variables, structural_simplify, getname, unwrap
+using ModelingToolkit: ModelingToolkit, Equation, System, Differential
+using ModelingToolkit: equations, full_equations, get_variables, mtkcompile, getname, unwrap
 using ModelingToolkit: parameters, unknowns, independent_variables, observed, defaults
 using Symbolics: Symbolics, fixpoint_sub, substitute
 using RecursiveArrayTools: RecursiveArrayTools
 using ArgCheck: @argcheck
 using LinearAlgebra: Diagonal, I
 using SymbolicUtils.Code: Let, Assignment
+using SymbolicUtils: SymbolicUtils
 using OrderedCollections: OrderedDict
 
 using NetworkDynamics: NetworkDynamics, set_metadata!,
@@ -17,23 +18,26 @@ import NetworkDynamics: VertexModel, EdgeModel, AnnotatedSym
 
 include("MTKExt_utils.jl")
 
+import NetworkDynamics: implicit_output
+ModelingToolkit.@register_symbolic implicit_output(x)
+
 """
-    VertexModel(sys::ODESystem, inputs, outputs;
+    VertexModel(sys::System, inputs, outputs;
                 verbose=false, name=getname(sys), extin=nothing, ff_to_constraint=true, kwargs...)
 
-Create a `VertexModel` object from a given `ODESystem` created with ModelingToolkit.
+Create a `VertexModel` object from a given `System` created with ModelingToolkit.
 You need to provide 2 lists of symbolic names (`Symbol` or `Vector{Symbols}`):
 - `inputs`: names of variables in you equation representing the aggregated edge states
 - `outputs`: names of variables in you equation representing the node output
 
 Additional kw arguments:
-- `name`: Set name of the component model. Will be lifted from the ODESystem name.
+- `name`: Set name of the component model. Will be lifted from the System name.
 - `extin=nothing`: Provide external inputs as pairs, i.e. `extin=[:extvar => VIndex(1, :a)]`
    will bound the variable `extvar(t)` in the equations to the state `a` of the first vertex.
-- `ff_to_constraint=true`: Controlls, whether output transformations `g` which depend on inputs should be
+- `ff_to_constraint=true`: Controls, whether output transformations `g` which depend on inputs should be
   transformed into constraints. Defaults to true since ND.jl does not handle vertices with FF yet.
 """
-function VertexModel(sys::ODESystem, inputs, outputs; verbose=false, name=getname(sys),
+function VertexModel(sys::System, inputs, outputs; verbose=false, name=getname(sys),
                      ff_to_constraint=true, extin=nothing, kwargs...)
     warn_missing_features(sys)
     inputs = inputs isa AbstractVector ? inputs : [inputs]
@@ -87,12 +91,12 @@ function VertexModel(sys::ODESystem, inputs, outputs; verbose=false, name=getnam
 end
 
 """
-    EdgeModel(sys::ODESystem, srcin, dstin, AntiSymmetric(dstout); kwargs...)
+    EdgeModel(sys::System, srcin, dstin, AntiSymmetric(dstout); kwargs...)
 
-Create a `EdgeModel` object from a given `ODESystem` created with ModelingToolkit for **single sided models**.
+Create a `EdgeModel` object from a given `System` created with ModelingToolkit for **single sided models**.
 
 Here you only need to provide one list of output symbols: `dstout`.
-To make it clear how to handle the single-sided output definiton, you musst wrap
+To make it clear how to handle the single-sided output definition, you must wrap
 the symbol vector in
 - `AntiSymmetric(dstout)`,
 - `Symmetric(dstout)`, or
@@ -100,13 +104,13 @@ the symbol vector in
 
 Additional `kwargs` are the same as for the double-sided EdgeModel MTK constructor.
 """
-EdgeModel(sys::ODESystem, srcin, dstin, dstout; kwargs...) = EdgeModel(sys, srcin, dstin, nothing, dstout; kwargs...)
+EdgeModel(sys::System, srcin, dstin, dstout; kwargs...) = EdgeModel(sys, srcin, dstin, nothing, dstout; kwargs...)
 
 """
-    EdgeModel(sys::ODESystem, srcin, dstin, srcout, dstout;
+    EdgeModel(sys::System, srcin, dstin, srcout, dstout;
               verbose=false, name=getname(sys), extin=nothing, ff_to_constraint=false, kwargs...)
 
-Create a `EdgeModel` object from a given `ODESystem` created with ModelingToolkit.
+Create a `EdgeModel` object from a given `System` created with ModelingToolkit.
 You need to provide 4 lists of symbolic names (`Symbol` or `Vector{Symbols}`):
 - `srcin`: names of variables in you equation representing the node state at the source
 - `dstin`: names of variables in you equation representing the node state at the destination
@@ -114,13 +118,13 @@ You need to provide 4 lists of symbolic names (`Symbol` or `Vector{Symbols}`):
 - `dstout`: names of variables in you equation representing the output at the destination
 
 Additional kw arguments:
-- `name`: Set name of the component model. Will be lifted from the ODESystem name.
+- `name`: Set name of the component model. Will be lifted from the System name.
 - `extin=nothing`: Provide external inputs as pairs, i.e. `extin=[:extvar => VIndex(1, :a)]`
    will bound the variable `extvar(t)` in the equations to the state `a` of the first vertex.
-- `ff_to_constraint=false`: Controlls, whether output transformations `g` which depend on inputs should be
+- `ff_to_constraint=false`: Controls, whether output transformations `g` which depend on inputs should be
   transformed into constraints.
 """
-function EdgeModel(sys::ODESystem, srcin, dstin, srcout, dstout; verbose=false, name=getname(sys),
+function EdgeModel(sys::System, srcin, dstin, srcout, dstout; verbose=false, name=getname(sys),
                    ff_to_constraint=false, extin=nothing, kwargs...)
     warn_missing_features(sys)
     srcin = srcin isa AbstractVector ? srcin : [srcin]
@@ -137,7 +141,7 @@ function EdgeModel(sys::ODESystem, srcin, dstin, srcout, dstout; verbose=false,
     singlesided = isnothing(srcout)
     if singlesided && !(dstout isa AnnotatedSym)
         throw(ArgumentError("If you only provide one output (single sided \
-        model), it musst be wrapped either in `AntiSymmetric`, `Symmetric` or \
+        model), it must be wrapped either in `AntiSymmetric`, `Symmetric` or \
         `Directed`!"))
     end
 
@@ -283,7 +287,6 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     alloutputs = reduce(union, outputss)
 
     missing_inputs = Set{Symbolic}()
-    implicit_outputs = Set{Symbolic}() # fully implicit outputs which do not appear in the equations
     sys = if ModelingToolkit.iscomplete(_sys)
         deepcopy(_sys)
     else
@@ -294,13 +297,28 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
             verbose && @warn "The specified inputs ($missing_inputs) do not appear in the equations of the system!"
             _openinputs = setdiff(_openinputs, missing_inputs)
         end
-        _definedoutputs = alloutputs ∩ all_eq_vars
-        if !(Set(_definedoutputs) == Set(alloutputs))
-            implicit_outputs = setdiff(alloutputs, _definedoutputs)
-            verbose && @warn "The specified outputs $implicit_outputs do not appear in the equations of the system!"
+
+        implicit_outputs = setdiff(alloutputs, all_eq_vars)
+        if !isempty(implicit_outputs)
+            throw(
+                ArgumentError("The outputs $(getname.(implicit_outputs)) do not appear in the equations of the system! \
+                    Try to to make them explicit using `implicit_output`\n" *
+                    NetworkDynamics.implicit_output_docstring)
+            )
+        end
+
+        verbose && @info "Simplifying system with inputs $_openinputs and outputs $alloutputs"
+        try
+            mtkcompile(_sys; inputs=_openinputs, outputs=alloutputs, simplify=false)
+        catch e
+            if e isa ModelingToolkit.ExtraEquationsSystemException
+                msg = "The system could not be compiled because of extra equations! \
+                       Sometimes, this can be related to fully implicit output equations. \
+                       Check `@doc implicit_output` for more information."
+                throw(ArgumentError(msg))
+            end
+            rethrow(e)
         end
-        verbose && @info "Simplifying system with inputs $_openinputs and outputs $_definedoutputs"
-        structural_simplify(_sys, (_openinputs, _definedoutputs); simplify=false)[1]
     end
 
     allparams = parameters(sys) # contains inputs!
@@ -308,14 +326,18 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     params = setdiff(allparams, Set(allinputs))
 
     # extract the main equations and observed equations
-    eqs::Vector{Equation} = ModelingToolkit.subs_constants(equations(sys))
-    obseqs_sorted::Vector{Equation} = ModelingToolkit.subs_constants(observed(sys))
+    eqs::Vector{Equation} = equations(sys)
+    obseqs_sorted::Vector{Equation} = observed(sys)
     fix_metadata!(eqs, sys);
     fix_metadata!(obseqs_sorted, sys);
 
+    # get rid of the implicit_output(⋅) terms
+    remove_implicit_output_fn!(eqs)
+    remove_implicit_output_fn!(obseqs_sorted)
+
     # assert the ordering of states and equations
     explicit_states = Symbolic[eq_type(eq)[2] for eq in eqs if !isnothing(eq_type(eq)[2])]
-    implicit_states = setdiff(unknowns(sys), explicit_states) ∪ implicit_outputs
+    implicit_states = setdiff(unknowns(sys), explicit_states)
 
     if length(explicit_states) + length(implicit_states) !== length(eqs)
         buf = IOBuffer()
@@ -420,7 +442,7 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
         end
         verbose && @info "Transformed algebraic eqs" eqs
 
-        # create massmatrix, we don't use the method provided by ODESystem because of reordering
+        # create massmatrix, we don't use the method provided by System because of reordering
         mm = generate_massmatrix(eqs)
         verbose && @info "Generated mass matrix" mm
         mm
diff --git a/src/NetworkDynamics.jl b/src/NetworkDynamics.jl
index 083ed86a8..4ee043e57 100644
--- a/src/NetworkDynamics.jl
+++ b/src/NetworkDynamics.jl
@@ -37,6 +37,7 @@ import SymbolicIndexingInterface as SII
 using SymbolicIndexingInterface: variable_symbols, parameter_symbols
 using StaticArrays: StaticArrays, SVector
 
+export implicit_output
 include("utils.jl")
 
 export VertexModel, EdgeModel
diff --git a/src/utils.jl b/src/utils.jl
index 9ef6bc640..11056b77c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -222,3 +222,32 @@ function _sym_to_int(x::SymbolicView, idx::Int)
     idx
 end
 _sym_to_int(x::SymbolicView, idx) = _sym_to_int.(Ref(x), idx)
+
+# temp variable to splice docstring into ArgumentError in MTK Ext
+const implicit_output_docstring = """
+This is a helper function to define MTK models with **fully implicit outputs**.
+It is sort of a barrier for `Symbolics` to not descent in to the equation. When added to
+an equation, it does nothing (defined as 0), but it tricks MTK/Symbolics into believing the
+equation depends on `x`. This can be necessary to define a model with fully implicit outputs.
+
+    @mtkmodel ImplicitForcing begin
+        @variables begin
+            u(t), [description = "Input Variable", input=true]
+            y(t), [description = "fully implicit output", output=true]
+        end
+        @equations begin
+            # 0 ~ u  # WRONG!
+            0 ~ u + implicit_output(y) # CORRECT!
+        end
+    end
+    VertexModel(ImplicitForcing(name=:implicit), [:u], [:y])
+"""
+"""
+    implicit_output(x) = 0
+    ModelingToolkit.@register_symbolic implicit_output(x)
+
+$implicit_output_docstring
+
+For more information see the NetworkDynamics docs on [fully implicit outputs](@ref Fully-Implicit-Outputs).
+"""
+implicit_output(x) = 0
diff --git a/test/MTK_test.jl b/test/MTK_test.jl
index 20d325353..f1399dab1 100644
--- a/test/MTK_test.jl
+++ b/test/MTK_test.jl
@@ -8,8 +8,6 @@ using Graphs
 using Chairmarks: @b
 using Test
 
-# TODO: Clean up MTK tests
-
 @mtkmodel Bus begin
     @variables begin
         θ(t), [description = "voltage angle", output=true]
@@ -235,30 +233,8 @@ b = @b $(NetworkDynamics.compfg(v))($data...)
 @test b.allocs == 0
 
 
-# test fully implicit outputs
-@mtkmodel FullyImplicit begin
-    @variables begin
-        u(t), [description = "Input Variable", input=true]
-        x(t), [description = "Explicit Variable"]
-        y(t), [description = "Implicit Variable, present in equations"]
-        z(t), [description = "fully implicit variable, not present but output", output=true]
-    end
-    @equations begin
-        Dt(x) ~ -x
-        0 ~ sqrt(y+x)
-        0 ~ u # implicitly forces z becaus u(z)
-    end
-end
-@named fullyimplicit = FullyImplicit()
-v = VertexModel(fullyimplicit, [:u], [:z])
-@test v.mass_matrix == LinearAlgebra.Diagonal([1,0,0])
-@test v.sym == [:x, :z, :y]
-
-data = NetworkDynamics.rand_inputs_fg(v)
-b = @b $(NetworkDynamics.compfg(v))($data...)
-@test b.allocs == 0
-
 @testset "Test constants in MTK models" begin
+    # from MTK@10 onwars, constants are just non-tunable parameters!
     @mtkmodel DQSwing_Constants begin
         @extend DQBus()
         @variables begin
@@ -358,8 +334,96 @@ end
 
     @named pv = PVConstraint(; P=2, V=1)
     @named busbar = BusBar()
-    mtkbus = ODESystem(connect(busbar.terminal, pv.terminal), t; systems=[busbar, pv], name=:pvbus)
+    mtkbus = System(connect(busbar.terminal, pv.terminal), t; systems=[busbar, pv], name=:pvbus)
     vf = VertexModel(mtkbus, [:busbar₊i_r, :busbar₊i_i], [:busbar₊u_r, :busbar₊u_i], verbose=false)
 
     @test vf.sym == [:busbar₊u_r,:busbar₊u_i]
 end
+
+@mtkmodel GasNode begin
+    @variables begin
+        p(t), [description="Pressure"] # node output
+        q̃_nw(t), [description="aggregated flow from pipes into node"] # node input
+        q̃_inj(t), [description="flow injected into the network"]
+    end
+    @equations begin
+        q̃_inj ~ -q̃_nw
+    end
+end
+@mtkmodel StaticProsumerNode begin
+    @extend GasNode()
+    @parameters begin
+        q̃_prosumer, [description="flow injected by prosumer"]
+    end
+    @equations begin
+        -q̃_nw ~ q̃_prosumer
+    end
+end
+@mtkmodel Wrapper begin
+    @components begin
+        prosumer = StaticProsumerNode()
+    end
+    @variables begin
+        p(t), [description="Pressure at prosumer"]
+        q̃_nw(t), [description="aggregated flow from pipes into node"] # node input
+    end
+    @equations begin
+        p ~ prosumer.p # connect the pressure output of the prosumer to the wrapper
+        q̃_nw ~ prosumer.q̃_nw # connect the flow input of the prosumer to the wrapper
+    end
+end
+@mtkmodel WrapperFixed begin
+    @components begin
+        prosumer = StaticProsumerNode()
+    end
+    @variables begin
+        p(t), [description="Pressure at prosumer"]
+        q̃_nw(t), [description="aggregated flow from pipes into node"] # node input
+    end
+    @equations begin
+        p ~ prosumer.p # connect the pressure output of the prosumer to the wrapper
+        q̃_nw + implicit_output(p) ~ prosumer.q̃_nw # connect the flow input of the prosumer to the wrapper
+    end
+end
+@testset "Test transformation of implicit outputs" begin
+    # first, lets test if the underlying MTK problem still exists
+    # see https://github.com/SciML/ModelingToolkit.jl/pull/3686
+    @mtkmodel ImplicitForcing begin
+        @variables begin
+            u(t), [description = "Input Variable", input=true]
+            y(t), [description = "fully implicit output", output=true]
+        end
+        @equations begin
+            0 ~ sqrt(u) # implicitly forces output y because u=f(y) in  closed loop
+        end
+    end
+    @named implicit = ImplicitForcing()
+    simp = mtkcompile(implicit; inputs = ModelingToolkit.unbound_inputs(implicit))
+    @test isempty(equations(simp)) # the equation was dropped!
+
+    # test fully implicit outputs
+    @mtkmodel FullyImplicit begin
+        @variables begin
+            u(t), [description = "Input Variable", input=true]
+            x(t), [description = "Explicit Variable"]
+            y(t), [description = "Implicit Variable, present in equations"]
+            z(t), [description = "fully implicit variable, not present but output", output=true]
+        end
+        @equations begin
+            Dt(x) ~ -x
+            0 ~ sqrt(y+x)
+            0 ~ u # implicitly forces z becaus u(z)
+        end
+    end
+    @named fullyimplicit = FullyImplicit()
+    @test_throws ArgumentError VertexModel(fullyimplicit, [:u], [:z])
+
+    # from init_tutorial
+    # dependent MTKModels need to be defined at top level, so they are in front of the testset
+    @named prosumer = StaticProsumerNode() # consumer
+    @test_throws ArgumentError VertexModel(prosumer, [:q̃_nw], [:p])
+    @named prosumer_wrapped = Wrapper()
+    @test_throws ArgumentError VertexModel(prosumer_wrapped, [:q̃_nw], [:p])
+    @named prosumer_fixed = WrapperFixed()
+    VertexModel(prosumer_fixed, [:q̃_nw], [:p]) # no throw
+end
diff --git a/test/initialization_test.jl b/test/initialization_test.jl
index 6492e46ee..08292071d 100644
--- a/test/initialization_test.jl
+++ b/test/initialization_test.jl
@@ -76,11 +76,12 @@ end
 
     @test_throws ErrorException NetworkDynamics.initialize_component!(vf, tol=0.0; verbose=true)
 
-    # make empty problem
+    # make empty problem (no free vars)
     vf_comp = copy(vf)
     set_default!(vf_comp, :Pm, get_init(vf_comp, :Pm))
     set_default!(vf_comp, :θ, get_init(vf_comp, :θ))
     set_default!(vf_comp, :ω, get_init(vf_comp, :ω))
+    set_default!(vf_comp, :V, get_init(vf_comp, :V))
 
     NetworkDynamics.initialize_component!(vf_comp; verbose=true)
     dump_initial_state(vf)
@@ -111,7 +112,7 @@ end
     # Verify the metadata was updated and values were applied correctly
     @test get_default(vf_conflict, :u_r) == 0.0
     @test get_default(vf_conflict, :u_i) == 1.0
-    @test get_init(vf_conflict, :θ) ≈ pi/2
+    @test get_init(vf_conflict, :θ) % 2pi ≈ pi/2
 
     # change metadtaa by providing custom input
     vf_sync = copy(vf)
@@ -438,7 +439,10 @@ end
     @test_throws ArgumentError InitConstraint(c10, c11)
 end
 
+_recursive_replace(el::Symbol, dict) = haskey(dict, el) ? dict[el] : NaN
+_recursive_replace(containter, dict) = map(el -> _recursive_replace(el, dict), containter)
 @testset "test input mapping" begin
+    using NetworkDynamics: insym_normalized, outsym_normalized
     vm = Lib.swing_mtk()
     c = @initconstraint begin
         :Pdamping # observable
@@ -447,9 +451,18 @@ end
         :ω # state
         :Pmech # parameter
     end
+    dict = Dict(:Pdamping => 1, :P=> 2, :θ => 3, :ω => 4, :Pmech => 5)
+
+    outdata = _recursive_replace(outsym_normalized(vm), dict)
+    udata   = _recursive_replace(sym(vm), dict)
+    indata  = _recursive_replace(insym_normalized(vm), dict)
+    pdata   = _recursive_replace(psym(vm), dict)
+    obsdata = _recursive_replace(obssym(vm), dict)
+
     mapping! = NetworkDynamics.generate_init_input_mapping(vm, c)
+
     syms = zeros(5)
-    mapping!(syms, ([3],), [NaN, 4], ([2],), [NaN, 5, NaN], [1])
+    mapping!(syms, outdata, udata, indata, pdata, obsdata)
     @test syms == 1:5
 
     em = Lib.line_mtk()
@@ -459,9 +472,17 @@ end
         :active #param
         :Δθ # obs
     end
+    dict = Dict(:P => 1, :dstθ => 2, :active => 3, :Δθ => 4)
+
+    outdata = _recursive_replace(outsym_normalized(em), dict)
+    udata   = _recursive_replace(sym(em), dict)
+    indata  = _recursive_replace(insym_normalized(em), dict)
+    pdata   = _recursive_replace(psym(em), dict)
+    obsdata = _recursive_replace(obssym(em), dict)
+
     mapping! = NetworkDynamics.generate_init_input_mapping(em, c)
     syms = zeros(4)
-    mapping!(syms, ([NaN],[1]), [], ([NaN],[2]), [NaN, NaN, 3], [4])
+    mapping!(syms, outdata, udata, indata, pdata, obsdata)
     @test syms == 1:4
 end