diff --git a/Project.toml b/Project.toml
index 03582b2f..9e420199 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,10 +32,12 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StyledStrings = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 
 [weakdeps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
@@ -43,6 +45,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [extensions]
 NetworkDynamicsCUDAExt = ["CUDA", "Adapt"]
+NetworkDynamicsDataFramesExt = ["DataFrames"]
 NetworkDynamicsMTKExt = ["ModelingToolkit", "SymbolicUtils", "Symbolics"]
 NetworkDynamicsSymbolicsExt = ["Symbolics", "MacroTools"]
 
@@ -81,4 +84,5 @@ SymbolicIndexingInterface = "0.3.27"
 SymbolicUtils = "3.24"
 Symbolics = "6.19.0"
 TimerOutputs = "0.5.23"
+YAML = "0.4.12"
 julia = "1.10"
diff --git a/docs/src/API.md b/docs/src/API.md
index b7be3b3c..5a92516a 100644
--- a/docs/src/API.md
+++ b/docs/src/API.md
@@ -19,7 +19,7 @@ EdgeModel()
 ```@docs
 VertexModel(::ModelingToolkit.ODESystem, ::Any, ::Any)
 EdgeModel(::ModelingToolkit.ODESystem, ::Any, ::Any, ::Any, ::Any)
-EdgeModel(::ModelingToolkit.ODESystem, ::Any, ::Any, ::Any)
+EdgeModel(::ModelingToolkit.ODESystem, ::Any, ::Any, ::NetworkDynamics.AnnotatedSym)
 ```
 
 ### Output Function Helpers/Wrappers 
diff --git a/ext/NetworkDynamicsDataFramesExt.jl b/ext/NetworkDynamicsDataFramesExt.jl
new file mode 100644
index 00000000..6e60fa5a
--- /dev/null
+++ b/ext/NetworkDynamicsDataFramesExt.jl
@@ -0,0 +1,96 @@
+module NetworkDynamicsDataFramesExt
+
+using NetworkDynamics: NetworkDynamics, Network
+using DataFrames: DataFrames, DataFrame
+using OrderedCollections: OrderedDict
+
+function _df_from_list(list, f)
+    df = DataFrame()
+    for (i, v) in list
+        row = OrderedDict{Symbol,Any}()
+        row[:idx] = i
+        for sym in f(v)
+            if NetworkDynamics.has_default_or_init(v, sym)
+                row[sym] = NetworkDynamics.get_default_or_init(v, sym)
+            else
+                row[sym] = missing
+            end
+        end
+        push!(df, NamedTuple(row); cols=:union)
+    end
+    df
+end
+function _df_from_pair(list, pair)
+    key, f = pair
+    df = DataFrame()
+    for (i, v) in list
+        row = OrderedDict{Symbol,Any}(:idx => i, key => f(v))
+        push!(df, NamedTuple(row); cols=:union)
+    end
+    df
+end
+
+function NetworkDynamics.describe_vertices(nw::Network, extras...; parameters=true, states=true, batch=nothing)
+    pairs = enumerate(nw.im.vertexm);
+    batches = map(idx -> findfirst(batch -> idx ∈ batch.indices, nw.vertexbatches), first.(pairs))
+
+    if !isnothing(batch)
+        idxs = findall(b -> b ∈ batch, batches)
+        pairs = collect(pairs)[idxs]
+        batch = collect(batches)[idxs]
+    end
+    isempty(pairs) && return DataFrame()
+
+    basedf = DataFrame(
+        idx = first.(pairs),
+        name = map(v->last(v).name, pairs),
+        batch = map(idx -> findfirst(batch -> idx ∈ batch.indices, nw.vertexbatches), first.(pairs)),
+    )
+
+    dfs = [basedf,]
+    if parameters
+        push!(dfs, _df_from_list(pairs, NetworkDynamics.psym))
+    end
+    if states
+        push!(dfs, _df_from_list(pairs, NetworkDynamics.sym))
+    end
+    for p in extras
+        push!(dfs, _df_from_pair(pairs, p))
+    end
+
+    foldl((a,b) -> DataFrames.leftjoin(a,b; on=:idx), dfs)
+end
+
+function NetworkDynamics.describe_edges(nw::Network, extras...; parameters=true, states=true, batch=nothing)
+    pairs = enumerate(nw.im.edgem);
+    batches = map(idx -> findfirst(batch -> idx ∈ batch.indices, nw.layer.edgebatches), first.(pairs))
+
+    if !isnothing(batch)
+        idxs = findall(b -> b ∈ batch, batches)
+        pairs = collect(pairs)[idxs]
+        batch = collect(batches)[idxs]
+    end
+    isempty(pairs) && return DataFrame()
+
+    basedf = DataFrame(
+        idx = first.(pairs),
+        srcdst = map(idx -> nw.im.edgevec[idx].src => nw.im.edgevec[idx].dst, first.(pairs)),
+        name = map(v->last(v).name, pairs),
+        batch = map(idx -> findfirst(batch -> idx ∈ batch.indices, nw.vertexbatches), first.(pairs)),
+    )
+
+    dfs = [basedf,]
+    if parameters
+        push!(dfs, _df_from_list(pairs, NetworkDynamics.psym))
+    end
+    if states
+        push!(dfs, _df_from_list(pairs, NetworkDynamics.sym))
+    end
+    for p in extras
+        push!(dfs, _df_from_pair(pairs, p))
+    end
+
+    foldl((a,b) -> DataFrames.leftjoin(a,b; on=:idx), dfs)
+end
+
+end
diff --git a/ext/NetworkDynamicsMTKExt.jl b/ext/NetworkDynamicsMTKExt.jl
index 14900fc0..7c09b25d 100644
--- a/ext/NetworkDynamicsMTKExt.jl
+++ b/ext/NetworkDynamicsMTKExt.jl
@@ -99,8 +99,26 @@ the symbol vector in
 - `Directed(dstout)`.
 
 Additional `kwargs` are the same as for the double-sided EdgeModel MTK constructor.
+
+Instead of wrapping, you can also provide the annotation as keyword argument like `annotation=:AntiSymmetric`.
 """
-EdgeModel(sys::ODESystem, srcin, dstin, dstout; kwargs...) = EdgeModel(sys, srcin, dstin, nothing, dstout; kwargs...)
+function EdgeModel(sys::ODESystem, srcin, dstin, dstout::AnnotatedSym; kwargs...)
+    EdgeModel(sys, srcin, dstin, nothing, dstout; kwargs...)
+end
+function EdgeModel(sys::ODESystem, srcin, dstin, dstout; annotation=nothing, kwargs...)
+    _dstout = if annotation==:AntiSymmetric
+        AntiSymmetric(dstout)
+    elseif annotation==:Symmetric
+        Symmetric(dstout)
+    elseif annotation==:Directed
+        Directed(dstout)
+    else
+        throw(ArgumentError("When constructing single sided EdgeModel, you must either wrap \
+            the output syms (e.g. AntiSymmetric(dstout)) orprovide the annotation as keyword \
+            argument."))
+    end
+    EdgeModel(sys, srcin, dstin, nothing, _dstout; kwargs...)
+end
 
 """
     EdgeModel(sys::ODESystem, srcin, srcout, dstin, dstout;
@@ -220,8 +238,12 @@ function _get_metadata(sys, name)
         if def isa Symbolic
             def = fixpoint_sub(def, alldefaults)
         end
-        def isa Symbolic && error("Could not resolve default $(ModelingToolkit.getdefault(sym)) for $name")
-        nt = (; nt..., default=def)
+
+        if def isa Symbolic
+            @warn "Could not resolve default term $(ModelingToolkit.getdefault(sym)) for $name, leave free."
+        else
+            nt = (; nt..., default=def)
+        end
     end
 
     # check for guess both in symbol metadata and in guesses of system
diff --git a/src/NetworkDynamics.jl b/src/NetworkDynamics.jl
index 65eb809d..5e171163 100644
--- a/src/NetworkDynamics.jl
+++ b/src/NetworkDynamics.jl
@@ -35,6 +35,8 @@ using InteractiveUtils: subtypes
 
 import SymbolicIndexingInterface as SII
 using StaticArrays: StaticArrays, SVector
+using YAML: YAML
+using OrderedCollections: OrderedDict
 
 include("utils.jl")
 
@@ -96,6 +98,12 @@ const CHECK_COMPONENT = Ref(true)
 export chk_component
 include("doctor.jl")
 
+export load_network
+include("importexport.jl")
+
+export describe_vertices, describe_edges
+function describe_vertices end
+function describe_edges end
 #=
 using StyledStrings
 s1 = styled"{bright_red:brred} {bright_green:brgreen} {bright_yellow:bryellow} {bright_blue:brblue} {bright_magenta:brmagenta} {bright_cyan:brcyan} {bright_black:brblack} {bright_white:brwhite}";
@@ -114,7 +122,20 @@ const ND_FACES = [
     :NetworkDynamics_fftype => StyledStrings.Face(foreground=:bright_blue),
 ]
 
-__init__() = foreach(StyledStrings.addface!, ND_FACES)
+function __init__()
+    if isdefined(Base.Experimental, :register_error_hint)
+        Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
+            if exc.f ∈ (describe_vertices, describe_edges)
+                ext = Base.get_extension(NetworkDynamics, :NetworkDynamicsDataFramesExt)
+                if isnothing(ext)
+                    printstyled(io, "\nLoad `DataFrames` in order to use `describe_vertices` and `describe_edges`."; bold=true)
+                end
+            end
+        end
+    end
+
+    foreach(StyledStrings.addface!, ND_FACES)
+end
 
 function reloadfaces!()
     StyledStrings.resetfaces!()
diff --git a/src/component_functions.jl b/src/component_functions.jl
index 6acb8667..086a098e 100644
--- a/src/component_functions.jl
+++ b/src/component_functions.jl
@@ -117,7 +117,7 @@ See also [`Symmetric`](@ref), [`Directed`](@ref), [`Fiducial`](@ref) and [`State
 struct AntiSymmetric{G} <: SingleSidedOutputWrapper
     g::G
 end
-AntiSymmetric(g::Union{AbstractVector,Number}) = AntiSymmetric(StateMask(g))
+AntiSymmetric(g::Union{AbstractVector{<:Number},Number}) = AntiSymmetric(StateMask(g))
 @inline function (c::AntiSymmetric)(osrc, odst, args...)
     @inline c.g(odst, args...)
     @inbounds for i in 1:length(osrc)
@@ -142,7 +142,7 @@ See also [`AntiSymmetric`](@ref), [`Directed`](@ref), [`Fiducial`](@ref) and [`S
 struct Symmetric{G} <: SingleSidedOutputWrapper
     g::G
 end
-Symmetric(g::Union{AbstractVector,Number}) = Symmetric(StateMask(g))
+Symmetric(g::Union{AbstractVector{<:Number},Number}) = Symmetric(StateMask(g))
 @inline function (c::Symmetric)(osrc, odst, args...)
     @inline c.g(odst, args...)
     @inbounds for i in 1:length(osrc)
@@ -167,7 +167,7 @@ See also [`AntiSymmetric`](@ref), [`Symmetric`](@ref), [`Fiducial`](@ref) and [`
 struct Directed{G} <: SingleSidedOutputWrapper
     g::G
 end
-Directed(g::Union{AbstractVector,Number}) = Directed(StateMask(g))
+Directed(g::Union{AbstractVector{<:Number},Number}) = Directed(StateMask(g))
 @inline function (c::Directed)(osrc, odst, args...)
     @inline c.g(odst, args...)
     nothing
@@ -226,7 +226,7 @@ Annotate a vector of output-symbols as `AntiSymmetric`, used when creating `Edge
 single-sided MTK models.
 """
 AntiSymmetric(s::Symbol) = AntiSymmetric([s])
-AntiSymmetric(s::AbstractVector{<:Symbol}) = AnnotatedSym(AntiSymmetric, s)
+AntiSymmetric(s::AbstractVector) = AnnotatedSym(AntiSymmetric, convert(Vector{Symbol}, s))
 """
     Symmetric(s::AbstractVector{<:Symbol})
 
@@ -234,7 +234,7 @@ Annotate a vector of output-symbols as `Symmetric`, used when creating `EdgeMode
 single-sided MTK models.
 """
 Symmetric(s::Symbol) = Symmetric([s])
-Symmetric(s::AbstractVector{<:Symbol}) = AnnotatedSym(Symmetric, s)
+Symmetric(s::AbstractVector) = AnnotatedSym(Symmetric, convert(Vector{Symbol}, s))
 """
     Directed(s::AbstractVector{<:Symbol})
 
@@ -242,7 +242,7 @@ Annotate a vector of output-symbols as `Directed`, used when creating `EdgeModel
 single-sided MTK models.
 """
 Directed(s::Symbol) = Directed([s])
-Directed(s::AbstractVector{<:Symbol}) = AnnotatedSym(Directed, s)
+Directed(s::AbstractVector) = AnnotatedSym(Directed, convert(Vector{Symbol}, s))
 
 
 abstract type ComponentModel end
diff --git a/src/construction.jl b/src/construction.jl
index 0d35a568..c49565ca 100644
--- a/src/construction.jl
+++ b/src/construction.jl
@@ -226,7 +226,7 @@ function _component_hash(c::ComponentModel)
     ))
 end
 
-function Network(vertexfs, edgefs; kwargs...)
+function Network(vertexfs, edgefs; warn_order=true, kwargs...)
     vertexfs = vertexfs isa VertexModel ? [vertexfs] : vertexfs
     edgefs   = edgefs isa EdgeModel     ? [edgefs]   : edgefs
     @argcheck all(has_graphelement, vertexfs) "All vertex models must have assigned `graphelement` to implicitly construct graph!"
@@ -269,6 +269,14 @@ function Network(vertexfs, edgefs; kwargs...)
 
     vfs_ordered = [vdict[k] for k in vertices(g)]
     efs_ordered = [edict[k] for k in edges(g)]
+    if warn_order && any(vfs_ordered .!== vertexfs)
+        @warn "Order of vertex models was changed to match the natural ordering of vertices in graph!\
+               Disable warning with kw `warn_order=false`."
+    end
+    if warn_order && any(efs_ordered .!== edgefs)
+        @warn "Order of edge models was changed to match the natural ordering edges in graph! \
+               Disable warning with kw `warn_order=false`."
+    end
 
     Network(g, vfs_ordered, efs_ordered; check_graphelement=false, kwargs...)
 end
diff --git a/src/importexport.jl b/src/importexport.jl
new file mode 100644
index 00000000..59bbe7a2
--- /dev/null
+++ b/src/importexport.jl
@@ -0,0 +1,236 @@
+"""
+For a given model, collect all referenced models
+"""
+collect_model_dependencies(model) = unique!(_collect_model_dependencies(String[], model))
+_collect_model_dependencies(deps, other) = deps
+function _collect_model_dependencies(deps, model::Union{OrderedDict,Vector})
+    for val in values(model)
+        if val isa String
+            m = match(r"^Models.(.*)$", val)
+            if !isnothing(m)
+                push!(deps, m[1])
+            end
+        else
+            _collect_model_dependencies(deps, val)
+        end
+    end
+    deps
+end
+
+"""
+    recursiv_resolve!(models, container)
+
+Go through container, and replace every occurence of "Models.<modelname>" with
+deepcopy of actual model dict.
+"""
+recursive_resolve!(models, _) = nothing
+function recursive_resolve!(models, container::Union{OrderedDict,Vector})
+    for (key, value) in pairs(container)
+        if value isa String
+            m = match(r"^Models.(.*)$", value)
+            if !isnothing(m)
+                container[key] = deepcopy(models[m[1]])
+            end
+        else
+            recursive_resolve!(models, value)
+        end
+    end
+    nothing
+end
+
+"""
+    update_parent_property!(parent, key, value)
+
+Update the property in the parent structure.
+- if key is nested (i.e. "a.b.c"), it will update parent["a"]["b"]["c"]
+- if key is indexed (i.e. "a[1]"), it will update parent["a"][1]
+"""
+function update_parent_property!(parent, key, value)
+    @debug "Update $key -> $value"
+    if contains(key, ".")
+        key_head, key_tail = split(key, ".", limit=2)
+        @debug "Found nested key $key_head -> $key_tail"
+        m = match(r"^(.*)\[(.*)\]$", key_head)
+        _parent = if !isnothing(m)
+            parent[m[1]][parse(Int, m[2])]
+        else
+            parent[key]
+        end
+        @debug m
+        update_parent_property!(_parent, key_tail, value)
+    else
+        m = match(r"^(.*)\[(.*)\]$", key)
+        if !isnothing(m)
+            if !(value isa eltype(parent[m[1]]))
+                # need to widen type
+                parent[m[1]] = collect(Any, parent[m[1]])
+            end
+            parent[m[1]][parse(Int, m[2])] = value
+        else
+            parent[key] = value
+        end
+    end
+end
+"""
+    merge_properties!(parent, child)
+
+Takes a child dict, and applies every key-value pair to the parent dict.
+keys will be resolved using update_parent_property!
+"""
+function merge_properties!(parent, child)
+    for (key, value) in pairs(child)
+        update_parent_property!(parent, key, value)
+    end
+    parent
+end
+
+"""
+    depth_first_flatten!(container)
+
+This goes through the container, and if it finds a dict with key "MODEL", it will
+merge the properties of the referenced parent model with the current dict.
+"""
+depth_first_flatten!(_) = nothing
+function depth_first_flatten!(container::Union{OrderedDict,Vector})
+    for (key, sub) in pairs(container)
+        depth_first_flatten!(sub)
+        if sub isa Union{OrderedDict} && haskey(sub, "MODEL")
+            parent_model = sub["MODEL"]
+            delete!(sub, "MODEL")
+            container[key] = merge_properties!(parent_model, sub)
+        end
+    end
+    nothing
+end
+
+widen_container_type(container::OrderedDict{K}) where {K} = convert(OrderedDict{K,Any}, container)
+widen_container_type(container::Vector) = convert(Vector{Any}, container)
+
+function replace_in_template(template, placeholder, model)
+    _replace_in_template!(deepcopy(template), placeholder, model)
+end
+_replace_in_template!(template, placeholder, model) = template
+function _replace_in_template!(container::Union{OrderedDict,Vector}, placeholder, model)
+    container = widen_container_type(container)
+    for (key, value) in pairs(container)
+        if value == placeholder
+            container[key] = deepcopy(model)
+        else
+            container[key] = _replace_in_template!(value, placeholder, model)
+        end
+    end
+    container
+end
+
+apply_wrappers!(::Nothing) = nothing
+function apply_wrappers!(container::Union{OrderedDict,Vector})
+    haskey(container, "WRAPPER") || return
+
+    wrapper = container["WRAPPER"]
+    delete!(container, "WRAPPER")
+    for (key, sub) in pairs(container)
+        container[key] = replace_in_template(wrapper, "MODEL", sub)
+    end
+    container
+end
+
+function parse_network(file)
+    data = YAML.load_file(file, dicttype=OrderedDict{Any,Any})
+    dependencies = OrderedDict(keys(data["Models"]) .=> map(collect_model_dependencies, values(data["Models"])))
+
+    # dependencies muss by subset of models
+    @assert reduce(vcat, values(dependencies)) ⊆ keys(data["Models"])
+    resolved_models = findall(deps -> isempty(deps), dependencies)
+    delete!.(Ref(dependencies), resolved_models) # delete resolve dependnecies
+
+    # resolve all models, i.e. replace parent/template models with full dict
+    while true
+        isempty(dependencies) && break
+        resolved_something = false
+        for (mod, deps) in dependencies
+            if deps ⊆ resolved_models
+                @debug "Resolving $mod -> $deps"
+                recursive_resolve!(data["Models"], data["Models"][mod])
+                delete!(dependencies, mod)
+                push!(resolved_models, mod)
+                resolved_something = true
+            else
+                @debug "Skip $mod -> $deps"
+            end
+        end
+        resolved_something && continue
+        error("Could not resolve all dependencies")
+    end
+
+    # apply wrappers to all models
+    apply_wrappers!(data["VertexModels"])
+    apply_wrappers!(data["EdgeModels"])
+
+    # with all deps resolve, we can replace the models edge and vertex list
+    recursive_resolve!(data["Models"], data["VertexModels"])
+    recursive_resolve!(data["Models"], data["EdgeModels"])
+
+    # traverses the tree and flattens the model properties, i.e. merges
+    # the local given properties with the parent/template model properties
+    depth_first_flatten!(data)
+    return data
+end
+
+recursive_construct(constructors, data) = _depth_first_construct!(constructors, deepcopy(data))
+_depth_first_construct!(_, data) = data
+_depth_first_construct!(_, s::String) = startswith(s, ":") ? Symbol(s[2:end]) : s
+function _depth_first_construct!(constructors, data::Union{OrderedDict,Vector})
+    data = widen_container_type(data)
+    for (key, value) in pairs(data)
+        data[key] = _depth_first_construct!(constructors, value)
+    end
+    if data isa OrderedDict && haskey(data, "CONSTRUCTOR")
+        if !haskey(constructors, data["CONSTRUCTOR"])
+            error("No constructor found for $(data["CONSTRUCTOR"])")
+        end
+        args = get(data, "ARGS", [])
+        kwargs = Dict{Symbol, Any}()
+        for (key, value) in data
+            key ∈ ("CONSTRUCTOR", "ARGS") && continue
+            kwargs[Symbol(key)] = value
+        end
+        @debug "call: $(data["CONSTRUCTOR"])(args...; $(kwargs))"
+        return constructors[data["CONSTRUCTOR"]](args...; kwargs...)
+    end
+    data
+end
+
+function build_network(data, constructors)
+    vertexm = VertexModel[]
+    for (k, v) in data["VertexModels"]
+        @debug "Constructing VertexModel $k"
+        vm = recursive_construct(constructors, v)
+        set_graphelement!(vm, k)
+        push!(vertexm, vm)
+    end
+
+    edgem = EdgeModel[]
+    for (k, v) in data["EdgeModels"]
+        @debug "Constructing EdgeModel $k"
+        em = recursive_construct(constructors, v)
+
+        m = match(r"^(.*)=>(.*)$", k)
+        isnothing(m) && error("Edge key must be of form 'src=>dst', got $k")
+        src = tryparse(Int, m[1])
+        dst = tryparse(Int, m[2])
+        if isnothing(src)
+            src = m[1]
+        end
+        if isnothing(dst)
+            dst = m[2]
+        end
+        set_graphelement!(em, src => dst)
+        push!(edgem, em)
+    end
+    Network(vertexm, edgem)
+end
+
+function load_network(constructors, file)
+    data = parse_network(file)
+    build_network(data, constructors)
+end
diff --git a/src/metadata.jl b/src/metadata.jl
index 9b3e13dc..551c3744 100644
--- a/src/metadata.jl
+++ b/src/metadata.jl
@@ -551,3 +551,7 @@ function _has_changed_hash(aliased_cfs)
     end
     changed
 end
+
+
+function describe_vertices end
+function describe_edges end
diff --git a/test/Project.toml b/test/Project.toml
index 514f06f1..0228ada6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -7,6 +7,7 @@ Bonito = "824d6782-a2ef-11e9-3a09-e5662e0c26f8"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
diff --git a/test/importexport_test.jl b/test/importexport_test.jl
new file mode 100644
index 00000000..fc47b522
--- /dev/null
+++ b/test/importexport_test.jl
@@ -0,0 +1,131 @@
+using NetworkDynamics
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as Dt
+using StableRNGs
+
+@mtkmodel Swing begin
+    @variables begin
+        θ(t) = 0.0, [description = "voltage angle", output=true]
+        ω(t) = 0.0, [description = "Rotor frequency"]
+        Pel(t), [description = "Electical power flowing from network into node", input=true]
+        Pdamping(t), [description = "Damping Power (observed)"]
+    end
+    @parameters begin
+        M = 1, [description = "Inertia"]
+        D = 0.1, [description = "Damping"]
+        Pmech=1.0, [description = "Mechanical Power"]
+    end
+    @equations begin
+        Dt(θ) ~ ω
+        Pdamping ~ - D * ω
+        Dt(ω) ~ 1/M * (Pmech + Pdamping + Pel)
+    end
+end
+
+@mtkmodel Load begin
+    @parameters begin
+        Pd = -1.0, [description = "Power Load"]
+    end
+    @variables begin
+        θ(t) = 0.0, [description = "voltage angle", output=true]
+        Pel(t), [description = "Electrical power flowing from network into node", input=true]
+    end
+    @equations begin
+        0 ~ Pel + Pd
+    end
+end
+
+@mtkmodel StaticPowerLine begin
+    @parameters begin
+        K = 1.63, [description = "line conductance parameter"]
+    end
+    @variables begin
+        srcθ(t), [description = "voltage angle at src end", input=true]
+        dstθ(t), [description = "voltage angle at dst end", input=true]
+        P(t), [description = "flow towards node at dst end", output=true]
+        Δθ(t), [description = "Voltage angle difference"]
+    end
+    @equations begin
+        Δθ ~ srcθ - dstθ
+        P ~ K*sin(Δθ)
+    end
+end
+
+function CombinedModel(subsystems...; kwargs...)
+    @variables begin
+        θ(t) = 0.0, [description = "voltage angle", output=true]
+        Pel(t), [description = "Electical power flowing from network into node", input=true]
+    end
+    eqs = [
+        (θ ~ sub.θ for sub in subsystems)...,
+        Pel ~ +((sub.Pel for sub in subsystems)...) # kirchoff constaint
+    ]
+    ODESystem(eqs, t, [θ, Pel], []; systems=collect(subsystems), kwargs...)
+end
+
+v1 = let
+    @named swing = Swing(; Pmech=1.0, M=1.0, D=1.0)
+    VertexModel(swing, [:Pel], [:θ], vidx=1)
+end
+
+v2 = let
+    @named load = Load(; Pd=-1.0)
+    VertexModel(load, [:Pel], [:θ], vidx=2)
+end
+
+v3 = let
+    @named swing = Swing(; Pmech=1.0, M=1.0, D=1.0)
+    @named load = Load(; Pd=-2.0)
+    @named swingload = CombinedModel(swing, load)
+    VertexModel(swingload, [:Pel], [:θ], vidx=3)
+end
+
+v4 = let
+    @named swing2 = Swing(; Pmech=2.0, M=1.0, D=1.0)
+    @named load = Load(; Pd=-1.0)
+    @named swingload2 = CombinedModel(swing2, load)
+    VertexModel(swingload2, [:Pel], [:θ], vidx=4)
+end
+
+e1 = let
+    @named line = StaticPowerLine(; K=0.5)
+    EdgeModel(line, [:srcθ], [:dstθ], [:P]; annotation=:AntiSymmetric, src=1, dst=2)
+end
+e2 = let
+    @named line = StaticPowerLine(; K=1.0)
+    EdgeModel(line, [:srcθ], [:dstθ], [:P]; annotation=:AntiSymmetric, src=2, dst=3)
+end
+e3 = let
+    @named line3 = StaticPowerLine()
+    EdgeModel(line3, [:srcθ], [:dstθ], [:P]; annotation=:AntiSymmetric, src=3, dst=4)
+end
+
+nw1 = Network([v1,v2,v3,v4],[e1,e2,e3])
+
+rng = StableRNG(1)
+u = randn(rng, dim(nw1))
+du1 = [NaN for _ in 1:dim(nw1)]
+nw1(du1, u, pflat(NWState(nw1)), NaN)
+
+constructors = Dict(
+    "Swing"=>Swing,
+    "VertexModel"=>VertexModel,
+    "Load"=>Load,
+    "CombinedModel"=>CombinedModel,
+    "StaticPowerLine"=>StaticPowerLine,
+    "EdgeModel"=>EdgeModel
+)
+nw2 = NetworkDynamics.load_network(constructors, "testgrid.yaml")
+
+du2 = [NaN for _ in 1:dim(nw2)]
+nw2(du2, u, pflat(NWState(nw2)), NaN)
+@test du1 ≈ du2
+
+_getname = x -> getproperty(x, :name)
+@test map(_getname, nw1.im.vertexm) == map(_getname, nw2.im.vertexm)
+@test map(_getname, nw1.im.edgem) == map(_getname, nw2.im.edgem)
+
+s1 = NWState(nw1)
+s2 = NWState(nw2)
+@test uflat(s1) == uflat(s2)
+@test pflat(s1) == pflat(s2)
diff --git a/test/testgrid.yaml b/test/testgrid.yaml
new file mode 100644
index 00000000..081139b6
--- /dev/null
+++ b/test/testgrid.yaml
@@ -0,0 +1,76 @@
+title: Power Grid Example
+version: "1.0"
+description: 4 Bus Powergrid Test Network
+
+# special words are "MODEL", "CONSTRUCTOR", "ARGS", "WRAPPER"
+Models:
+  Swing:
+    CONSTRUCTOR: Swing
+    Pmech: 1.0
+    M: 1.0
+    D: 1.0
+    name: :swing
+
+  Load:
+    CONSTRUCTOR: Load
+    Pd: -42.0
+    name: :load
+
+  SwingLoad:
+    CONSTRUCTOR: CombinedModel
+    ARGS:
+      - MODEL: Models.Swing
+      - MODEL: Models.Load
+        Pd: -2.0
+    name: :swingload
+
+  SwingLoad2:
+    CONSTRUCTOR: CombinedModel
+    ARGS:
+      - MODEL: Models.Swing
+        Pmech: 2.0
+      - MODEL: Models.Load
+        Pd: -1.0
+
+  StaticPowerLine:
+    CONSTRUCTOR: StaticPowerLine
+    K: 1.63
+    name: :line
+
+  StaticPowerLineK05:
+    MODEL: Models.StaticPowerLine
+    K: 0.5
+
+VertexModels:
+  WRAPPER:
+    CONSTRUCTOR: VertexModel
+    ARGS:
+      - MODEL
+      - [":Pel"]
+      - [":θ"]
+  1: Models.Swing
+  2:
+    MODEL: Models.Load
+    Pd: -1.0
+  3: Models.SwingLoad
+  4:
+    MODEL: Models.SwingLoad2
+    name: :swingload2
+    ARGS[1].name: :swing2
+
+EdgeModels:
+  WRAPPER:
+    CONSTRUCTOR: EdgeModel
+    ARGS:
+      - MODEL
+      - [":srcθ"] 
+      - [":dstθ"]
+      - [":P"]
+    annotation: :AntiSymmetric
+  "1=>2": Models.StaticPowerLineK05
+  "2=>3":
+    MODEL: Models.StaticPowerLine
+    K: 1.0
+  "3=>4":
+    MODEL: Models.StaticPowerLine
+    name: :line3