diff --git a/examples/electrical_components.jl b/examples/electrical_components.jl
index 7e1f01380d..9a4ca63f7a 100644
--- a/examples/electrical_components.jl
+++ b/examples/electrical_components.jl
@@ -3,22 +3,10 @@ using ModelingToolkit, OrdinaryDiffEq
 
 @parameters t
 @connector function Pin(;name)
-    sts = @variables v(t)=1.0 i(t)=1.0
+    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
     ODESystem(Equation[], t, sts, []; name=name)
 end
 
-function ModelingToolkit.connect(::Type{Pin}, ps...)
-    eqs = [
-           0 ~ sum(p->p.i, ps) # KCL
-          ]
-    # KVL
-    for i in 1:length(ps)-1
-        push!(eqs, ps[i].v ~ ps[i+1].v)
-    end
-
-    return eqs
-end
-
 function Ground(;name)
     @named g = Pin()
     eqs = [g.v ~ 0]
diff --git a/examples/rc_model.jl b/examples/rc_model.jl
index c7bc4b901f..88b112bee8 100644
--- a/examples/rc_model.jl
+++ b/examples/rc_model.jl
@@ -11,7 +11,8 @@ V = 1.0
 rc_eqs = [
           connect(source.p, resistor.p)
           connect(resistor.n, capacitor.p)
-          connect(capacitor.n, source.n, ground.g)
+          connect(capacitor.n, source.n)
+          connect(capacitor.n, ground.g)
          ]
 
 @named rc_model = ODESystem(rc_eqs, t)
diff --git a/examples/serial_inductor.jl b/examples/serial_inductor.jl
index 63d3215a8e..feff486a5d 100644
--- a/examples/serial_inductor.jl
+++ b/examples/serial_inductor.jl
@@ -10,7 +10,8 @@ eqs = [
        connect(source.p, resistor.p)
        connect(resistor.n, inductor1.p)
        connect(inductor1.n, inductor2.p)
-       connect(source.n, inductor2.n, ground.g)
+       connect(source.n, inductor2.n)
+       connect(inductor2.n, ground.g)
       ]
 
 @named ll_model = ODESystem(eqs, t)
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 898a61ab0f..15444ee086 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -123,6 +123,7 @@ if !isdefined(Graphs, :IncrementalCycleTracker)
 end
 
 include("systems/abstractsystem.jl")
+include("systems/connectors.jl")
 
 include("systems/diffeqs/odesystem.jl")
 include("systems/diffeqs/sdesystem.jl")
@@ -157,8 +158,6 @@ for S in subtypes(ModelingToolkit.AbstractSystem)
     @eval convert_system(::Type{<:$S}, sys::$S) = sys
 end
 
-struct Flow end
-
 export AbstractTimeDependentSystem, AbstractTimeIndependentSystem, AbstractMultivariateSystem
 export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system
 export DAEFunctionExpr, DAEProblemExpr
@@ -173,7 +172,8 @@ export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem
 export ControlSystem
-export alias_elimination, flatten, connect, @connector
+export alias_elimination, flatten
+export connect, @connector, Connection, Flow, Stream, instream
 export ode_order_lowering, liouville_transform
 export runge_kutta_discretize
 export PDESystem
@@ -182,7 +182,7 @@ export Equation, ConstrainedEquation
 export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variables, independent_variable, states, parameters, equations, controls, observed, structure
-export structural_simplify
+export structural_simplify, expand_connections
 export DiscreteSystem, DiscreteProblem
 
 export calculate_jacobian, generate_jacobian, generate_function
@@ -204,7 +204,7 @@ export toexpr, get_variables
 export simplify, substitute
 export build_function
 export modelingtoolkitize
-export @variables, @parameters, Flow
+export @variables, @parameters
 export @named, @nonamespace, @namespace, extend, compose
 
 end # module
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 5c8bc9495f..e357a3bb7b 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -219,7 +219,8 @@ for prop in [
              :domain
              :ivs
              :dvs
-             :connection_type
+             :connector_type
+             :connections
              :preface
             ]
     fname1 = Symbol(:get_, prop)
@@ -282,7 +283,7 @@ function getvar(sys::AbstractSystem, name::Symbol; namespace=false)
     elseif !isempty(systems)
         i = findfirst(x->nameof(x)==name, systems)
         if i !== nothing
-            return namespace ? rename(systems[i], renamespace(sys, name)) : systems[i]
+            return namespace ? renamespace(sys, systems[i]) : systems[i]
         end
     end
 
@@ -355,6 +356,7 @@ GlobalScope(sym::Union{Num, Symbolic}) = setmetadata(sym, SymScope, GlobalScope(
 
 renamespace(sys, eq::Equation) = namespace_equation(eq, sys)
 
+renamespace(names::AbstractVector, x) = foldr(renamespace, names, init=x)
 function renamespace(sys, x)
     x = unwrap(x)
     if x isa Symbolic
@@ -367,6 +369,8 @@ function renamespace(sys, x)
                 x
             end
         end
+    elseif x isa AbstractSystem
+        rename(x, renamespace(sys, nameof(x)))
     else
         Symbol(getname(sys), :₊, x)
     end
@@ -562,7 +566,24 @@ function round_trip_expr(t, var2name)
     args = map(Base.Fix2(round_trip_expr, var2name), arguments(t))
     return :($f($(args...)))
 end
-round_trip_eq(eq, var2name) = Expr(:call, :~, round_trip_expr(eq.lhs, var2name), round_trip_expr(eq.rhs, var2name))
+
+function round_trip_eq(eq::Equation, var2name)
+    if eq.lhs isa Connection
+        syss = get_systems(eq.rhs)
+        call = Expr(:call, connect)
+        for sys in syss
+            strs = split(string(nameof(sys)), "₊")
+            s = Symbol(strs[1])
+            for st in strs[2:end]
+                s = Expr(:., s, Meta.quot(Symbol(st)))
+            end
+            push!(call.args, s)
+        end
+        call
+    else
+        Expr(:call, (~), round_trip_expr(eq.lhs, var2name), round_trip_expr(eq.rhs, var2name))
+    end
+end
 
 function push_eqs!(stmt, eqs, var2name)
     eqs_name = gensym(:eqs)
@@ -854,6 +875,7 @@ topological sort of the observed equations. When `simplify=true`, the `simplify`
 function will be applied during the tearing process.
 """
 function structural_simplify(sys::AbstractSystem; simplify=false)
+    sys = expand_connections(sys)
     sys = initialize_system_structure(alias_elimination(sys))
     check_consistency(sys)
     if sys isa ODESystem
@@ -905,71 +927,6 @@ function check_eqs_u0(eqs, dvs, u0; check_length=true, kwargs...)
     return nothing
 end
 
-###
-### Connectors
-###
-
-function with_connection_type(expr)
-    @assert expr isa Expr && (expr.head == :function || (expr.head == :(=) &&
-                                       expr.args[1] isa Expr &&
-                                       expr.args[1].head == :call))
-
-    sig = expr.args[1]
-    body = expr.args[2]
-
-    fname = sig.args[1]
-    args = sig.args[2:end]
-
-    quote
-        struct $fname
-            $(gensym()) -> 1 # this removes the default constructor
-        end
-        function $fname($(args...))
-            function f()
-                $body
-            end
-            res = f()
-            $isdefined(res, :connection_type) ? $Setfield.@set!(res.connection_type = $fname) : res
-        end
-    end
-end
-
-macro connector(expr)
-    esc(with_connection_type(expr))
-end
-
-promote_connect_rule(::Type{T}, ::Type{S}) where {T, S} = Union{}
-promote_connect_rule(::Type{T}, ::Type{T}) where {T} = T
-promote_connect_type(t1::Type, t2::Type, ts::Type...) = promote_connect_type(promote_connect_rule(t1, t2), ts...)
-@inline function promote_connect_type(::Type{T}, ::Type{S}) where {T,S}
-    promote_connect_result(
-        T,
-        S,
-        promote_connect_rule(T,S),
-        promote_connect_rule(S,T)
-    )
-end
-
-promote_connect_result(::Type, ::Type, ::Type{T}, ::Type{Union{}}) where {T} = T
-promote_connect_result(::Type, ::Type, ::Type{Union{}}, ::Type{S}) where {S} = S
-promote_connect_result(::Type, ::Type, ::Type{T}, ::Type{T}) where {T} = T
-function promote_connect_result(::Type{T}, ::Type{S}, ::Type{P1}, ::Type{P2}) where {T,S,P1,P2}
-    throw(ArgumentError("connection promotion for $T and $S resulted in $P1 and $P2. " *
-                        "Define promotion only in one direction."))
-end
-
-throw_connector_promotion(T, S) = throw(ArgumentError("Don't know how to connect systems of type $S and $T"))
-promote_connect_result(::Type{T},::Type{S},::Type{Union{}},::Type{Union{}}) where {T,S} = throw_connector_promotion(T,S)
-
-promote_connect_type(::Type{T}, ::Type{T}) where {T} = T
-function promote_connect_type(T, S)
-    error("Don't know how to connect systems of type $S and $T")
-end
-
-function connect(syss...)
-    connect(promote_connect_type(map(get_connection_type, syss)...), syss...)
-end
-
 ###
 ### Inheritance & composition
 ###
@@ -990,7 +947,7 @@ function Base.hash(sys::AbstractSystem, s::UInt)
 end
 
 """
-    $(TYPEDSIGNATURES)
+$(TYPEDSIGNATURES)
 
 entend the `basesys` with `sys`, the resulting system would inherit `sys`'s name
 by default.
@@ -1026,7 +983,7 @@ end
 Base.:(&)(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys)) = extend(sys, basesys; name=name)
 
 """
-    $(SIGNATURES)
+$(SIGNATURES)
 
 compose multiple systems together. The resulting system would inherit the first
 system's name.
diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 3043f3b5e7..d738d2abb4 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -179,7 +179,10 @@ struct AliasGraphKeySet <: AbstractSet{Int}
 end
 Base.keys(ag::AliasGraph) = AliasGraphKeySet(ag)
 Base.iterate(agk::AliasGraphKeySet, state...) = Base.iterate(agk.ag.eliminated, state...)
-Base.in(i::Int, agk::AliasGraphKeySet) = agk.ag.aliasto[i] !== nothing
+function Base.in(i::Int, agk::AliasGraphKeySet)
+    aliasto = agk.ag.aliasto
+    1 <= i <= length(aliasto) && aliasto[i] !== nothing
+end
 
 count_nonzeros(a::AbstractArray) = count(!iszero, a)
 
diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl
new file mode 100644
index 0000000000..77eb046384
--- /dev/null
+++ b/src/systems/connectors.jl
@@ -0,0 +1,425 @@
+function with_connector_type(expr)
+    @assert expr isa Expr && (expr.head == :function || (expr.head == :(=) &&
+                                       expr.args[1] isa Expr &&
+                                       expr.args[1].head == :call))
+
+    sig = expr.args[1]
+    body = expr.args[2]
+
+    fname = sig.args[1]
+    args = sig.args[2:end]
+
+    quote
+        function $fname($(args...))
+            function f()
+                $body
+            end
+            res = f()
+            $isdefined(res, :connector_type) && $getfield(res, :connector_type) === nothing ? $Setfield.@set!(res.connector_type = $connector_type(res)) : res
+        end
+    end
+end
+
+macro connector(expr)
+    esc(with_connector_type(expr))
+end
+
+abstract type AbstractConnectorType end
+struct StreamConnector <: AbstractConnectorType end
+struct RegularConnector <: AbstractConnectorType end
+
+function connector_type(sys::AbstractSystem)
+    sts = get_states(sys)
+    #TODO: check the criteria for stream connectors
+    n_stream = 0
+    n_flow = 0
+    for s in sts
+        vtype = getmetadata(s, ModelingToolkit.VariableConnectType, nothing)
+        vtype === Stream && (n_stream += 1)
+        vtype === Flow   && (n_flow += 1)
+    end
+    (n_stream > 0 && n_flow > 1) && error("There are multiple flow variables in $(nameof(sys))!")
+    n_stream > 0 ? StreamConnector() : RegularConnector()
+end
+
+Base.@kwdef struct Connection
+    inners = nothing
+    outers = nothing
+end
+
+# everything is inner by default until we expand the connections
+Connection(syss) = Connection(inners=syss)
+get_systems(c::Connection) = c.inners
+function Base.in(e::Symbol, c::Connection)
+    any(k->nameof(k) === e, c.inners) || any(k->nameof(k) === e, c.outers)
+end
+
+const EMPTY_VEC = []
+
+function Base.show(io::IO, ::MIME"text/plain", c::Connection)
+    # It is a bit unfortunate that the display of an array of `Equation`s won't
+    # call this.
+    @unpack outers, inners = c
+    if outers === nothing && inners === nothing
+        print(io, "<Connection>")
+    else
+        syss = Iterators.flatten((something(inners, EMPTY_VEC), something(outers, EMPTY_VEC)))
+        splitting_idx = length(inners)
+        sys_str = join((string(nameof(s)) * (i <= splitting_idx ? ("::inner") : ("::outer")) for (i, s) in enumerate(syss)), ", ")
+        print(io, "<", sys_str, ">")
+    end
+end
+
+# symbolic `connect`
+function connect(sys1::AbstractSystem, sys2::AbstractSystem, syss::AbstractSystem...)
+    syss = (sys1, sys2, syss...)
+    length(unique(nameof, syss)) == length(syss) || error("connect takes distinct systems!")
+    Equation(Connection(), Connection(syss)) # the RHS are connected systems
+end
+
+# the actual `connect`.
+function connect(c::Connection; check=true)
+    @unpack inners, outers = c
+
+    flow_eqs = Equation[]
+    other_eqs = Equation[]
+
+    cnts = Iterators.flatten((inners, outers))
+    fs, ss = Iterators.peel(cnts)
+    splitting_idx = length(inners) # anything after the splitting_idx is outer.
+    first_sts = get_states(fs)
+    first_sts_set = Set(getname.(first_sts))
+    for sys in ss
+        current_sts = getname.(get_states(sys))
+        Set(current_sts) == first_sts_set || error("$(nameof(sys)) ($current_sts) doesn't match the connection type of $(nameof(fs)) ($first_sts).")
+    end
+
+    ceqs = Equation[]
+    for s in first_sts
+        name = getname(s)
+        vtype = getmetadata(s, VariableConnectType, Equality)
+        vtype === Stream && continue
+        isflow = vtype === Flow
+        rhs = 0 # only used for flow variables
+        fix_val = getproperty(fs, name) # used for equality connections
+        for (i, c) in enumerate(cnts)
+            isinner = i <= splitting_idx
+            # https://specification.modelica.org/v3.4/Ch15.html
+            var = getproperty(c, name)
+            if isflow
+                rhs += isinner ? var : -var
+            else
+                i == 1 && continue # skip the first iteration
+                push!(ceqs, fix_val ~ getproperty(c, name))
+            end
+        end
+        isflow && push!(ceqs, 0 ~ rhs)
+    end
+
+    return ceqs
+end
+
+instream(a) = term(instream, unwrap(a), type=symtype(a))
+SymbolicUtils.promote_symtype(::typeof(instream), _) = Real
+
+isconnector(s::AbstractSystem) = has_connector_type(s) && get_connector_type(s) !== nothing
+isstreamconnector(s::AbstractSystem) = isconnector(s) && get_connector_type(s) isa StreamConnector
+isstreamconnection(c::Connection) = any(isstreamconnector, c.inners) || any(isstreamconnector, c.outers)
+
+function print_with_indent(n, x)
+    print(" " ^ n)
+    show(stdout, MIME"text/plain"(), x)
+    println()
+end
+
+function split_sys_var(var)
+    var_name = string(getname(var))
+    sidx = findlast(isequal('₊'), var_name)
+    sidx === nothing && error("$var is not a namespaced variable")
+    connector_name = Symbol(var_name[1:prevind(var_name, sidx)])
+    streamvar_name = Symbol(var_name[nextind(var_name, sidx):end])
+    connector_name, streamvar_name
+end
+
+function flowvar(sys::AbstractSystem)
+    sts = get_states(sys)
+    for s in sts
+        vtype = getmetadata(s, ModelingToolkit.VariableConnectType, nothing)
+        vtype === Flow && return s
+    end
+    error("There in no flow variable in $(nameof(sys))")
+end
+
+collect_instream!(set, eq::Equation) = collect_instream!(set, eq.lhs) | collect_instream!(set, eq.rhs)
+
+function collect_instream!(set, expr, occurs=false)
+    istree(expr) || return occurs
+    op = operation(expr)
+    op === instream && (push!(set, expr); occurs = true)
+    for a in SymbolicUtils.unsorted_arguments(expr)
+        occurs |= collect_instream!(set, a, occurs)
+    end
+    return occurs
+end
+
+positivemax(m, ::Any; tol=nothing)= max(m, something(tol, 1e-8))
+
+function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10)
+    subsys = get_systems(sys)
+    isempty(subsys) && return sys
+
+    # post order traversal
+    @set! sys.systems = map(s->expand_connections(s, debug=debug), subsys)
+
+    outer_connectors = Symbol[]
+    for s in subsys
+        n = nameof(s)
+        isconnector(s) && push!(outer_connectors, n)
+    end
+    isouter = let outer_connectors=outer_connectors
+        function isouter(sys)::Bool
+            s = string(nameof(sys))
+            isconnector(sys) || error("$s is not a connector!")
+            idx = findfirst(isequal('₊'), s)
+            parent_name = Symbol(idx === nothing ? s : s[1:prevind(s, idx)])
+            parent_name in outer_connectors
+        end
+    end
+
+    sys = flatten(sys)
+    eqs′ = get_eqs(sys)
+    eqs = Equation[]
+    instream_eqs = Equation[]
+    instream_exprs = []
+    cts = [] # connections
+    for eq in eqs′
+        if eq.lhs isa Connection
+            push!(cts, get_systems(eq.rhs))
+        elseif collect_instream!(instream_exprs, eq)
+            push!(instream_eqs, eq)
+        else
+            push!(eqs, eq) # split connections and equations
+        end
+    end
+
+    # if there are no connections, we are done
+    isempty(cts) && return sys
+
+    sys2idx = Dict{Symbol,Int}() # system (name) to n-th connect statement
+    narg_connects = Connection[]
+    for (i, syss) in enumerate(cts)
+        # find intersecting connections
+        exclude = 0 # exclude the intersecting system
+        idx = 0     # idx of narg_connects
+        for (j, s) in enumerate(syss)
+            idx′ = get(sys2idx, nameof(s), nothing)
+            idx′ === nothing && continue
+            idx = idx′
+            exclude = j
+        end
+        if exclude == 0
+            outers = []
+            inners = []
+            for s in syss
+                isouter(s) ? push!(outers, s) : push!(inners, s)
+            end
+            push!(narg_connects, Connection(outers=outers, inners=inners))
+            for s in syss
+                sys2idx[nameof(s)] = length(narg_connects)
+            end
+        else
+            # fuse intersecting connections
+            for (j, s) in enumerate(syss); j == exclude && continue
+                sys2idx[nameof(s)] = idx
+                c = narg_connects[idx]
+                isouter(s) ? push!(c.outers, s) : push!(c.inners, s)
+            end
+        end
+    end
+
+    isempty(narg_connects) && error("Unreachable reached. Please file an issue.")
+
+    @set! sys.connections = narg_connects
+
+    # Bad things happen when there are more than one intersections
+    for c in narg_connects
+        @unpack outers, inners = c
+        len = length(outers) + length(inners)
+        allconnectors = Iterators.flatten((outers, inners))
+        dups = find_duplicates(nameof(c) for c in allconnectors)
+        length(dups) == 0 || error("$(Connection(sys)) has duplicated connections: $(dups).")
+    end
+
+    if debug
+        println("============BEGIN================")
+        println("Connections for [$(nameof(sys))]:")
+        foreach(Base.Fix1(print_with_indent, 4), narg_connects)
+    end
+
+    connection_eqs = Equation[]
+    for c in narg_connects
+        ceqs = connect(c)
+        debug && append!(connection_eqs, ceqs)
+        append!(eqs, ceqs)
+    end
+
+    if debug
+        println("Connection equations:")
+        foreach(Base.Fix1(print_with_indent, 4), connection_eqs)
+        println("=============END=================")
+    end
+
+    # stream variables
+    stream_connects = filter(isstreamconnection, narg_connects)
+    instream_eqs, additional_eqs = expand_instream(instream_eqs, instream_exprs, stream_connects; debug=debug, tol=tol)
+
+    @set! sys.eqs = [eqs; instream_eqs; additional_eqs]
+    return sys
+end
+
+function expand_instream(instream_eqs, instream_exprs, connects; debug=false, tol)
+    sub = Dict()
+    seen = Set()
+    for ex in instream_exprs
+        var = only(arguments(ex))
+        connector_name, streamvar_name = split_sys_var(var)
+
+        # find the connect
+        cidx = findfirst(c->connector_name in c, connects)
+        cidx === nothing && error("$var is not a variable inside stream connectors")
+        connect = connects[cidx]
+        connectors = Iterators.flatten((connect.inners, connect.outers))
+        # stream variable
+        sv = getproperty(first(connectors), streamvar_name; namespace=false)
+        inner_sc = connect.inners
+        outer_sc = connect.outers
+
+        n_outers = length(outer_sc)
+        n_inners = length(inner_sc)
+        outer_names = (nameof(s) for s in outer_sc)
+        inner_names = (nameof(s) for s in inner_sc)
+        if debug
+            println("Expanding: $ex")
+            isempty(inner_names) || println("Inner connectors: $(collect(inner_names))")
+            isempty(outer_names) || println("Outer connectors: $(collect(outer_names))")
+        end
+
+        # expand `instream`s
+        # https://specification.modelica.org/v3.4/Ch15.html
+        # Based on the above requirements, the following implementation is
+        # recommended:
+        if n_inners == 1 && n_outers == 0
+            connector_name === only(inner_names) || error("$var is not in any stream connector of $(nameof(ogsys))")
+            sub[ex] = var
+        elseif n_inners == 2 && n_outers == 0
+            connector_name in inner_names || error("$var is not in any stream connector of $(nameof(ogsys))")
+            idx = findfirst(c->nameof(c) === connector_name, inner_sc)
+            other = idx == 1 ? 2 : 1
+            sub[ex] = states(inner_sc[other], sv)
+        elseif n_inners == 1 && n_outers == 1
+            isinner = connector_name === only(inner_names)
+            isouter = connector_name === only(outer_names)
+            (isinner || isouter) || error("$var is not in any stream connector of $(nameof(ogsys))")
+            if isinner
+                outerstream = states(only(outer_sc), sv) # c_1.h_outflow
+                sub[ex] = outerstream
+            end
+        else
+            fv = flowvar(first(connectors))
+            i = findfirst(c->nameof(c) === connector_name, inner_sc)
+            if i !== nothing
+                si = isempty(outer_sc) ? 0 : sum(s->max(states(s, fv), 0), outer_sc)
+                for j in 1:n_inners; j == i && continue
+                    f = states(inner_sc[j], fv)
+                    si += max(-f, 0)
+                end
+
+                num = 0
+                den = 0
+                for j in 1:n_inners; j == i && continue
+                    f = states(inner_sc[j], fv)
+                    tmp = positivemax(-f, si; tol=tol)
+                    den += tmp
+                    num += tmp * states(inner_sc[j], sv)
+                end
+                for k in 1:n_outers
+                    f = states(outer_sc[k], fv)
+                    tmp = positivemax(f, si; tol=tol)
+                    den += tmp
+                    num += tmp * instream(states(outer_sc[k], sv))
+                end
+                sub[ex] = num / den
+            end
+
+        end
+    end
+
+    # additional equations
+    additional_eqs = Equation[]
+    for c in connects
+        outer_sc = c.outers
+        isempty(outer_sc) && continue
+        inner_sc = c.inners
+        n_outers = length(outer_sc)
+        n_inners = length(inner_sc)
+        connector_representative = first(outer_sc)
+        fv = flowvar(connector_representative)
+        for sv in get_states(connector_representative)
+            vtype = getmetadata(sv, ModelingToolkit.VariableConnectType, nothing)
+            vtype === Stream || continue
+            if n_inners == 1 && n_outers == 1
+                innerstream = states(only(inner_sc), sv)
+                outerstream = states(only(outer_sc), sv)
+                push!(additional_eqs, outerstream ~ innerstream)
+            elseif n_inners == 0 && n_outers == 2
+                # we don't expand `instream` in this case.
+                v1 = states(outer_sc[1], sv)
+                v2 = states(outer_sc[2], sv)
+                push!(additional_eqs, v1 ~ instream(v2))
+                push!(additional_eqs, v2 ~ instream(v1))
+            else
+                sq = 0
+                for q in 1:n_outers
+                    sq += sum(s->max(-states(s, fv), 0), inner_sc)
+                    for k in 1:n_outers; k == q && continue
+                        f = states(outer_sc[k], fv)
+                        sq += max(f, 0)
+                    end
+
+                    num = 0
+                    den = 0
+                    for j in 1:n_inners
+                        f = states(inner_sc[j], fv)
+                        tmp = positivemax(-f, sq; tol=tol)
+                        den += tmp
+                        num += tmp * states(inner_sc[j], sv)
+                    end
+                    for k in 1:n_outers; k == q && continue
+                        f = states(outer_sc[k], fv)
+                        tmp = positivemax(f, sq; tol=tol)
+                        den += tmp
+                        num += tmp * instream(states(outer_sc[k], sv))
+                    end
+                    push!(additional_eqs, states(outer_sc[q], sv) ~ num / den)
+                end
+            end
+        end
+    end
+
+    instream_eqs = map(Base.Fix2(substitute, sub), instream_eqs)
+    if debug
+        println("===========BEGIN=============")
+        println("Expanded equations:")
+        for eq in instream_eqs
+            print_with_indent(4, eq)
+        end
+        if !isempty(additional_eqs)
+            println("Additional equations:")
+            for eq in additional_eqs
+                print_with_indent(4, eq)
+            end
+        end
+        println("============END==============")
+    end
+    return instream_eqs, additional_eqs
+end
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 45b56a66e4..e573230cae 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -80,9 +80,13 @@ struct ODESystem <: AbstractODESystem
     """
     structure::Any
     """
-    connection_type: type of the system
+    connector_type: type of the system
     """
-    connection_type::Any
+    connector_type::Any
+    """
+    connections: connections in a system
+    """
+    connections::Any
     """
     preface: inject assignment statements before the evaluation of the RHS function.
     """
@@ -93,7 +97,7 @@ struct ODESystem <: AbstractODESystem
     """
     continuous_events::Vector{SymbolicContinuousCallback}
 
-    function ODESystem(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connection_type, preface, events; checks::Bool = true)
+    function ODESystem(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connector_type, connections, preface, events; checks::Bool = true)
         if checks
             check_variables(dvs,iv)
             check_parameters(ps,iv)
@@ -101,7 +105,7 @@ struct ODESystem <: AbstractODESystem
             check_equations(equations(events),iv)
             all_dimensionless([dvs;ps;iv]) || check_units(deqs)
         end
-        new(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connection_type, preface, events)
+        new(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connector_type, connections, preface, events)
     end
 end
 
@@ -114,7 +118,7 @@ function ODESystem(
                    default_u0=Dict(),
                    default_p=Dict(),
                    defaults=_merge(Dict(default_u0), Dict(default_p)),
-                   connection_type=nothing,
+                   connector_type=nothing,
                    preface=nothing,
                    continuous_events=nothing,
                    checks = true,
@@ -148,7 +152,7 @@ function ODESystem(
         throw(ArgumentError("System names must be unique."))
     end
     cont_callbacks = SymbolicContinuousCallbacks(continuous_events)
-    ODESystem(deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing, connection_type, preface, cont_callbacks, checks = checks)
+    ODESystem(deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing, connector_type, nothing, preface, cont_callbacks, checks = checks)
 end
 
 function ODESystem(eqs, iv=nothing; kwargs...)
@@ -171,7 +175,9 @@ function ODESystem(eqs, iv=nothing; kwargs...)
     end
     iv = value(iv)
     iv === nothing && throw(ArgumentError("Please pass in independent variables."))
+    compressed_eqs = Equation[] # equations that need to be expanded later, like `connect(a, b)`
     for eq in eqs
+        eq.lhs isa Union{Symbolic,Number} || (push!(compressed_eqs, eq); continue)
         collect_vars!(allstates, ps, eq.lhs, iv)
         collect_vars!(allstates, ps, eq.rhs, iv)
         if isdiffeq(eq)
@@ -186,7 +192,7 @@ function ODESystem(eqs, iv=nothing; kwargs...)
     end
     algevars = setdiff(allstates, diffvars)
     # the orders here are very important!
-    return ODESystem(append!(diffeq, algeeq), iv, collect(Iterators.flatten((diffvars, algevars))), ps; kwargs...)
+    return ODESystem(Equation[diffeq; algeeq; compressed_eqs], iv, collect(Iterators.flatten((diffvars, algevars))), ps; kwargs...)
 end
 
 # NOTE: equality does not check cached Jacobian
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index a31b0cf3e5..6518416819 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -84,16 +84,16 @@ struct SDESystem <: AbstractODESystem
     """
     type: type of the system
     """
-    connection_type::Any
+    connector_type::Any
 
-    function SDESystem(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type; checks::Bool = true)
+    function SDESystem(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type; checks::Bool = true)
         if checks
             check_variables(dvs,iv)
             check_parameters(ps,iv)
             check_equations(deqs,iv)
             all_dimensionless([dvs;ps;iv]) || check_units(deqs,neqs)
         end
-        new(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type)
+        new(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type)
     end
 end
 
@@ -105,7 +105,7 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
                    default_p=Dict(),
                    defaults=_merge(Dict(default_u0), Dict(default_p)),
                    name=nothing,
-                   connection_type=nothing,
+                   connector_type=nothing,
                    checks = true,
                    )
     name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
@@ -134,7 +134,7 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
     ctrl_jac = RefValue{Any}(Matrix{Num}(undef, 0, 0))
     Wfact   = RefValue(Matrix{Num}(undef, 0, 0))
     Wfact_t = RefValue(Matrix{Num}(undef, 0, 0))
-    SDESystem(deqs, neqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type, checks = checks)
+    SDESystem(deqs, neqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type, checks = checks)
 end
 
 SDESystem(sys::ODESystem, neqs; kwargs...) = SDESystem(equations(sys), neqs, get_iv(sys), states(sys), parameters(sys); kwargs...)
diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl
index 60891eae0d..df694c59e3 100644
--- a/src/systems/discrete_system/discrete_system.jl
+++ b/src/systems/discrete_system/discrete_system.jl
@@ -54,14 +54,14 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
     """
     type: type of the system
     """
-    connection_type::Any
-    function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connection_type; checks::Bool = true)
+    connector_type::Any
+    function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connector_type; checks::Bool = true)
         if checks
             check_variables(dvs, iv)
             check_parameters(ps, iv)
             all_dimensionless([dvs;ps;iv;ctrls]) ||check_units(discreteEqs)
         end
-        new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connection_type)
+        new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connector_type)
     end
 end
 
@@ -79,7 +79,7 @@ function DiscreteSystem(
                    default_u0=Dict(),
                    default_p=Dict(),
                    defaults=_merge(Dict(default_u0), Dict(default_p)),
-                   connection_type=nothing,
+                   connector_type=nothing,
                    kwargs...,
                   )
     name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
@@ -103,7 +103,7 @@ function DiscreteSystem(
     if length(unique(sysnames)) != length(sysnames)
         throw(ArgumentError("System names must be unique."))
     end
-    DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, defaults, connection_type, kwargs...)
+    DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, defaults, connector_type, kwargs...)
 end
 
 
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 9dffa5e789..b0bc0d8f32 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -52,14 +52,14 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
     """
     type: type of the system
     """
-    connection_type::Any
-    function JumpSystem{U}(ap::U, iv, states, ps, var_to_name, observed, name, systems, defaults, connection_type; checks::Bool = true) where U <: ArrayPartition
+    connector_type::Any
+    function JumpSystem{U}(ap::U, iv, states, ps, var_to_name, observed, name, systems, defaults, connector_type; checks::Bool = true) where U <: ArrayPartition
         if checks
             check_variables(states, iv)
             check_parameters(ps, iv)
             all_dimensionless([states;ps;iv]) || check_units(ap,iv)
         end
-        new{U}(ap, iv, states, ps, var_to_name, observed, name, systems, defaults, connection_type)
+        new{U}(ap, iv, states, ps, var_to_name, observed, name, systems, defaults, connector_type)
     end
 end
 
@@ -70,7 +70,7 @@ function JumpSystem(eqs, iv, states, ps;
                     default_p=Dict(),
                     defaults=_merge(Dict(default_u0), Dict(default_p)),
                     name=nothing,
-                    connection_type=nothing,
+                    connector_type=nothing,
                     checks = true,
                     kwargs...)
     name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
@@ -102,7 +102,7 @@ function JumpSystem(eqs, iv, states, ps;
     process_variables!(var_to_name, defaults, states)
     process_variables!(var_to_name, defaults, ps)
 
-    JumpSystem{typeof(ap)}(ap, value(iv), states, ps, var_to_name, observed, name, systems, defaults, connection_type, checks = checks)
+    JumpSystem{typeof(ap)}(ap, value(iv), states, ps, var_to_name, observed, name, systems, defaults, connector_type, checks = checks)
 end
 
 function generate_rate_function(js::JumpSystem, rate)
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index f873de51ac..6b50f4d04e 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -53,12 +53,12 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
     """
     type: type of the system
     """
-    connection_type::Any
-    function NonlinearSystem(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults, structure, connection_type; checks::Bool = true)
+    connector_type::Any
+    function NonlinearSystem(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults, structure, connector_type; checks::Bool = true)
         if checks
             all_dimensionless([states;ps]) ||check_units(eqs)
         end
-        new(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults, structure, connection_type)
+        new(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults, structure, connector_type)
     end
 end
 
@@ -69,7 +69,7 @@ function NonlinearSystem(eqs, states, ps;
                          default_p=Dict(),
                          defaults=_merge(Dict(default_u0), Dict(default_p)),
                          systems=NonlinearSystem[],
-                         connection_type=nothing,
+                         connector_type=nothing,
                          continuous_events=nothing, # this argument is only required for ODESystems, but is added here for the constructor to accept it without error
                          checks = true,
                          )
@@ -99,7 +99,7 @@ function NonlinearSystem(eqs, states, ps;
     process_variables!(var_to_name, defaults, states)
     process_variables!(var_to_name, defaults, ps)
 
-    NonlinearSystem(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults, nothing, connection_type, checks = checks)
+    NonlinearSystem(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults, nothing, connector_type, checks = checks)
 end
 
 function calculate_jacobian(sys::NonlinearSystem; sparse=false, simplify=false)
diff --git a/src/systems/pde/pdesystem.jl b/src/systems/pde/pdesystem.jl
index 9c5a8623f8..56d168ea54 100644
--- a/src/systems/pde/pdesystem.jl
+++ b/src/systems/pde/pdesystem.jl
@@ -55,7 +55,7 @@ struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
     """
     type: type of the system
     """
-    connection_type::Any
+    connector_type::Any
     """
     name: the name of the system
     """
@@ -63,14 +63,14 @@ struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
     @add_kwonly function PDESystem(eqs, bcs, domain, ivs, dvs,
                                    ps=SciMLBase.NullParameters();
                                    defaults=Dict(),
-                                   connection_type = nothing,
+                                   connector_type = nothing,
                                    checks::Bool = true,
                                    name
                                   )
         if checks
             all_dimensionless([dvs;ivs;ps]) ||check_units(eqs)
         end
-        new(eqs, bcs, domain, ivs, dvs, ps, defaults, connection_type, name)
+        new(eqs, bcs, domain, ivs, dvs, ps, defaults, connector_type, name)
     end
 end
 
diff --git a/src/utils.jl b/src/utils.jl
index c34839858b..8a868889f5 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -353,3 +353,21 @@ function get_postprocess_fbody(sys)
     end
     return pre_
 end
+
+"""
+$(SIGNATURES)
+
+find duplicates in an iterable object.
+"""
+function find_duplicates(xs, ::Val{Ret}=Val(false)) where Ret
+    appeared = Set()
+    duplicates = Set()
+    for x in xs
+        if x in appeared
+            push!(duplicates, x)
+        else
+            push!(appeared, x)
+        end
+    end
+    return Ret ? (appeared, duplicates) : duplicates
+end
diff --git a/src/variables.jl b/src/variables.jl
index 5d9d1789f9..1c27d3b098 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -11,6 +11,11 @@ Symbolics.option_to_metadata_type(::Val{:description}) = VariableDescriptionType
 Symbolics.option_to_metadata_type(::Val{:input}) = VariableInput
 Symbolics.option_to_metadata_type(::Val{:output}) = VariableOutput
 
+abstract type AbstractConnectType end
+struct Equality <: AbstractConnectType end # Equality connection
+struct Flow <: AbstractConnectType end     # sum to 0
+struct Stream <: AbstractConnectType end   # special stream connector
+
 function isvarkind(m, x)
     p = getparent(x, nothing)
     p === nothing || (x = p)
diff --git a/test/components.jl b/test/components.jl
index 8a9a8a321e..da79b67c78 100644
--- a/test/components.jl
+++ b/test/components.jl
@@ -1,80 +1,119 @@
-using Test
-using ModelingToolkit, OrdinaryDiffEq
-
-include("../examples/rc_model.jl")
-
-sys = structural_simplify(rc_model)
-@test !isempty(ModelingToolkit.defaults(sys))
-u0 = [
-      capacitor.v => 0.0
-      capacitor.p.i => 0.0
-      resistor.v => 0.0
-     ]
-prob = ODEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Rodas4())
-
-@test sol[resistor.p.i] == sol[capacitor.p.i]
-@test sol[resistor.n.i] == -sol[capacitor.p.i]
-@test sol[capacitor.n.i] == -sol[capacitor.p.i]
-@test iszero(sol[ground.g.i])
-@test iszero(sol[ground.g.v])
-@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
-
-u0 = [
-      capacitor.v => 0.0
-     ]
-prob = ODAEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Tsit5())
-
-@test sol[resistor.p.i] == sol[capacitor.p.i]
-@test sol[resistor.n.i] == -sol[capacitor.p.i]
-@test sol[capacitor.n.i] == -sol[capacitor.p.i]
-@test iszero(sol[ground.g.i])
-@test iszero(sol[ground.g.v])
-@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
-#using Plots
-#plot(sol)
-
-include("../examples/serial_inductor.jl")
-sys = structural_simplify(ll_model)
-u0 = [
-      inductor1.i => 0.0
-      inductor2.i => 0.0
-      inductor2.v => 0.0
-     ]
-prob = ODEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Rodas4())
-
-prob = ODAEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Tsit5())
-
-@variables t x1(t) x2(t) x3(t) x4(t)
-D = Differential(t)
-@named sys1_inner = ODESystem([D(x1) ~ x1], t)
-@named sys1_partial = compose(ODESystem([D(x2) ~ x2], t; name=:foo), sys1_inner)
-@named sys1 = extend(ODESystem([D(x3) ~ x3], t; name=:foo), sys1_partial)
-@named sys2 = compose(ODESystem([D(x4) ~ x4], t; name=:foo), sys1)
-@test_nowarn sys2.sys1.sys1_inner.x1 # test the correct nesting
-
-
-# compose tests
-@parameters t
-
-function record_fun(;name)
-    pars = @parameters a=10 b=100
-    ODESystem(Equation[], t, [], pars; name)
-end
-
-function first_model(;name)
-    @named foo=record_fun()
-
-    defs = Dict()
-    defs[foo.a] = 3
-    defs[foo.b] = 300
-    pars = @parameters x=2 y=20
-    compose(ODESystem(Equation[], t, [], pars; name, defaults=defs), foo)
-end
-@named goo = first_model()
-@unpack foo = goo
-@test ModelingToolkit.defaults(goo)[foo.a] == 3
-@test ModelingToolkit.defaults(goo)[foo.b] == 300
+using Test
+using ModelingToolkit, OrdinaryDiffEq
+
+include("../examples/rc_model.jl")
+
+sys = structural_simplify(rc_model)
+@test !isempty(ModelingToolkit.defaults(sys))
+u0 = [
+      capacitor.v => 0.0
+      capacitor.p.i => 0.0
+      resistor.v => 0.0
+     ]
+prob = ODEProblem(sys, u0, (0, 10.0))
+sol = solve(prob, Rodas4())
+
+@test sol[resistor.p.i] == sol[capacitor.p.i]
+@test sol[resistor.n.i] == -sol[capacitor.p.i]
+@test sol[capacitor.n.i] == -sol[capacitor.p.i]
+@test iszero(sol[ground.g.i])
+@test iszero(sol[ground.g.v])
+@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
+
+# Outer/inner connections
+function rc_component(;name)
+    R = 1
+    C = 1
+    @named p = Pin()
+    @named n = Pin()
+    @named resistor = Resistor(R=R)
+    @named capacitor = Capacitor(C=C)
+    eqs = [
+           connect(p, resistor.p);
+           connect(resistor.n, capacitor.p);
+           connect(capacitor.n, n);
+          ]
+    @named sys = ODESystem(eqs, t)
+    compose(sys, [p, n, resistor, capacitor]; name=name)
+end
+
+@named ground = Ground()
+@named source = ConstantVoltage(V=1)
+@named rc_comp = rc_component()
+eqs = [
+       connect(source.p, rc_comp.p)
+       connect(source.n, rc_comp.n)
+       connect(source.n, ground.g)
+      ]
+@named sys′ = ODESystem(eqs, t)
+@named sys_inner_outer = compose(sys′, [ground, source, rc_comp])
+expand_connections(sys_inner_outer, debug=true)
+sys_inner_outer = structural_simplify(sys_inner_outer)
+@test !isempty(ModelingToolkit.defaults(sys_inner_outer))
+u0 = [
+      rc_comp.capacitor.v => 0.0
+      rc_comp.capacitor.p.i => 0.0
+      rc_comp.resistor.v => 0.0
+     ]
+prob = ODEProblem(sys_inner_outer, u0, (0, 10.0))
+sol_inner_outer = solve(prob, Rodas4())
+@test sol[capacitor.v] ≈ sol_inner_outer[rc_comp.capacitor.v]
+
+u0 = [
+      capacitor.v => 0.0
+     ]
+prob = ODAEProblem(sys, u0, (0, 10.0))
+sol = solve(prob, Tsit5())
+
+@test sol[resistor.p.i] == sol[capacitor.p.i]
+@test sol[resistor.n.i] == -sol[capacitor.p.i]
+@test sol[capacitor.n.i] == -sol[capacitor.p.i]
+@test iszero(sol[ground.g.i])
+@test iszero(sol[ground.g.v])
+@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
+#using Plots
+#plot(sol)
+
+include("../examples/serial_inductor.jl")
+sys = structural_simplify(ll_model)
+u0 = [
+      inductor1.i => 0.0
+      inductor2.i => 0.0
+      inductor2.v => 0.0
+     ]
+prob = ODEProblem(sys, u0, (0, 10.0))
+sol = solve(prob, Rodas4())
+
+prob = ODAEProblem(sys, u0, (0, 10.0))
+sol = solve(prob, Tsit5())
+
+@variables t x1(t) x2(t) x3(t) x4(t)
+D = Differential(t)
+@named sys1_inner = ODESystem([D(x1) ~ x1], t)
+@named sys1_partial = compose(ODESystem([D(x2) ~ x2], t; name=:foo), sys1_inner)
+@named sys1 = extend(ODESystem([D(x3) ~ x3], t; name=:foo), sys1_partial)
+@named sys2 = compose(ODESystem([D(x4) ~ x4], t; name=:foo), sys1)
+@test_nowarn sys2.sys1.sys1_inner.x1 # test the correct nesting
+
+
+# compose tests
+@parameters t
+
+function record_fun(;name)
+    pars = @parameters a=10 b=100
+    ODESystem(Equation[], t, [], pars; name)
+end
+
+function first_model(;name)
+    @named foo=record_fun()
+
+    defs = Dict()
+    defs[foo.a] = 3
+    defs[foo.b] = 300
+    pars = @parameters x=2 y=20
+    compose(ODESystem(Equation[], t, [], pars; name, defaults=defs), foo)
+end
+@named goo = first_model()
+@unpack foo = goo
+@test ModelingToolkit.defaults(goo)[foo.a] == 3
+@test ModelingToolkit.defaults(goo)[foo.b] == 300
diff --git a/test/connectors.jl b/test/connectors.jl
deleted file mode 100644
index f492c9f26c..0000000000
--- a/test/connectors.jl
+++ /dev/null
@@ -1,44 +0,0 @@
-using Test, ModelingToolkit
-
-@parameters t
-
-@connector function Foo(;name)
-    @variables x(t)
-    ODESystem(Equation[], t, [x], [], defaults=Dict(x=>1.0), name=name)
-end
-
-@connector function Goo(;name)
-    @variables x(t)
-    @parameters p
-    ODESystem(Equation[], t, [x], [p], defaults=Dict(x=>1.0, p=>1.0), name=name)
-end
-
-function ModelingToolkit.connect(::Type{<:Foo}, ss...)
-    n = length(ss)-1
-    eqs = Vector{Equation}(undef, n)
-    for i in 1:n
-        eqs[i] = ss[i].x ~ ss[i+1].x
-    end
-    eqs
-end
-
-@named f1 = Foo()
-@named f2 = Foo()
-@named f3 = Foo()
-@named f4 = Foo()
-@named g = Goo()
-
-@test isequal(connect(f1, f2), [f1.x ~ f2.x])
-@test_throws ArgumentError connect(f1, g)
-
-# Note that since there're overloadings, these tests are not re-runable.
-ModelingToolkit.promote_connect_rule(::Type{<:Foo}, ::Type{<:Goo}) = Foo
-@test isequal(connect(f1, g), [f1.x ~ g.x])
-@test isequal(connect(f1, f2, g), [f1.x ~ f2.x; f2.x ~ g.x])
-@test isequal(connect(f1, f2, g, f3), [f1.x ~ f2.x; f2.x ~ g.x; g.x ~ f3.x])
-@test isequal(connect(f1, f2, g, f3, f4), [f1.x ~ f2.x; f2.x ~ g.x; g.x ~ f3.x; f3.x ~ f4.x])
-ModelingToolkit.promote_connect_rule(::Type{<:Goo}, ::Type{<:Foo}) = Foo
-@test isequal(connect(f1, g), [f1.x ~ g.x])
-# test conflict
-ModelingToolkit.promote_connect_rule(::Type{<:Goo}, ::Type{<:Foo}) = Goo
-@test_throws ArgumentError connect(f1, g)
diff --git a/test/discretesystem.jl b/test/discretesystem.jl
index da227d9b94..fe993a8661 100644
--- a/test/discretesystem.jl
+++ b/test/discretesystem.jl
@@ -112,13 +112,3 @@ linearized_eqs = [
     y(t - 2.0) ~ y(t)
 ]
 @test all(eqs2 .== linearized_eqs)
-
-# Test connection_type
-@connector function DiscreteComponent(;name)
-    @variables v(t) i(t)
-    DiscreteSystem(Equation[], t, [v, i], [], name=name, defaults=Dict(v=>1.0, i=>1.0))
-end
-
-@named d1 = DiscreteComponent()
-
-@test ModelingToolkit.get_connection_type(d1) == DiscreteComponent
diff --git a/test/runtests.jl b/test/runtests.jl
index f7e510dcd3..c38be89307 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,6 +22,7 @@ using SafeTestsets, Test
 @safetestset "Constraints Test" begin include("constraints.jl") end
 @safetestset "Reduction Test" begin include("reduction.jl") end
 @safetestset "Components Test" begin include("components.jl") end
+@safetestset "Stream Connnect Test" begin include("stream_connectors.jl") end
 @safetestset "PDE Construction Test" begin include("pde.jl") end
 @safetestset "Lowering Integration Test" begin include("lowering_solving.jl") end
 @safetestset "Test Big System Usage" begin include("bigsystem.jl") end
@@ -37,6 +38,5 @@ println("Last test requires gcc available in the path!")
 @safetestset "StructuralTransformations" begin include("structural_transformation/runtests.jl") end
 @testset "Serialization" begin include("serialization.jl") end
 @safetestset "print_tree" begin include("print_tree.jl") end
-@safetestset "connectors" begin include("connectors.jl") end
 @safetestset "error_handling" begin include("error_handling.jl") end
 @safetestset "root_equations" begin include("root_equations.jl") end
diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl
new file mode 100644
index 0000000000..ad7f8b11b3
--- /dev/null
+++ b/test/stream_connectors.jl
@@ -0,0 +1,217 @@
+using Test
+using ModelingToolkit
+@variables t
+
+@connector function TwoPhaseFluidPort(;name, P=0.0, m_flow=0.0, h_outflow=0.0)
+    vars = @variables h_outflow(t)=h_outflow [connect=Stream] m_flow(t)=m_flow [connect=Flow] P(t)=P
+    ODESystem(Equation[], t, vars, []; name=name)
+end
+
+function MassFlowSource_h(;name,
+        h_in=420e3,
+        m_flow_in=-0.01,
+    )
+
+    pars = @parameters begin
+        h_in=h_in
+        m_flow_in=m_flow_in
+    end
+
+    vars = @variables begin
+        P(t)
+    end
+
+    @named port = TwoPhaseFluidPort()
+
+    subs = [port]
+
+    eqns = Equation[]
+
+    push!(eqns, port.P ~ P)
+    push!(eqns, port.m_flow ~ -m_flow_in)
+    push!(eqns, port.h_outflow ~ h_in)
+
+    compose(ODESystem(eqns, t, vars, pars; name=name), subs)
+end
+
+# Simplified components.
+function AdiabaticStraightPipe(;name,
+        kwargs...,
+    )
+
+    vars = []
+    pars = []
+
+    @named port_a = TwoPhaseFluidPort()
+    @named port_b = TwoPhaseFluidPort()
+
+    subs = [port_a; port_b]
+
+    eqns = Equation[]
+
+    #=
+    push!(eqns, port_a.P ~ port_b.P)
+    push!(eqns, 0 ~ port_a.m_flow + port_b.m_flow)
+    push!(eqns, port_b.h_outflow ~ instream(port_a.h_outflow))
+    push!(eqns, port_a.h_outflow ~ instream(port_b.h_outflow))
+    =#
+
+    push!(eqns, connect(port_a, port_b))
+    sys = ODESystem(eqns, t, vars, pars; name=name)
+    sys = compose(sys, subs)
+end
+
+
+function SmallBoundary_Ph(;name,
+        P_in=1e6,
+        h_in=400e3,
+    )
+
+    vars = []
+
+    pars = @parameters begin
+        P=P_in
+        h=h_in
+    end
+
+    @named port1 = TwoPhaseFluidPort()
+
+    subs = [port1]
+
+    eqns = Equation[]
+
+    push!(eqns, port1.P ~ P)
+    push!(eqns, port1.h_outflow ~ h)
+
+    compose(ODESystem(eqns, t, vars, pars; name=name), subs)
+end
+
+
+# N1M1 model and test code.
+function N1M1(;name,
+        P_in=1e6,
+        h_in=400e3,
+        kwargs...,
+    )
+
+    @named port_a = TwoPhaseFluidPort()
+    @named source = SmallBoundary_Ph(P_in=P_in, h_in=h_in)
+
+    subs = [port_a; source]
+
+    eqns = Equation[]
+
+    push!(eqns, connect(source.port1, port_a))
+
+    sys = ODESystem(eqns, t, [], [], name=name)
+    sys = compose(sys, subs)
+end
+
+@named n1m1 = N1M1()
+@named pipe = AdiabaticStraightPipe()
+@named sink = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
+
+streams_a = [n1m1.port_a, pipe.port_a]
+streams_b = [pipe.port_b, sink.port]
+
+eqns = [
+        connect(n1m1.port_a, pipe.port_a)
+        connect(pipe.port_b, sink.port)
+       ]
+
+@named sys = ODESystem(eqns, t)
+@named n1m1Test = compose(sys, n1m1, pipe, sink)
+@test_nowarn structural_simplify(n1m1Test)
+
+
+# N1M2 model and test code.
+function N1M2(;name,
+        P_in=1e6,
+        h_in=400e3,
+        kwargs...,
+    )
+
+    @named port_a = TwoPhaseFluidPort()
+    @named port_b = TwoPhaseFluidPort()
+
+    @named source = SmallBoundary_Ph(P_in=P_in, h_in=h_in)
+
+    subs = [port_a; port_b; source]
+
+    eqns = Equation[]
+
+    push!(eqns, connect(source.port1, port_a))
+    push!(eqns, connect(source.port1, port_b))
+
+    sys = ODESystem(eqns, t, [], [], name=name)
+    sys = compose(sys, subs)
+end
+
+@named n1m2 = N1M2()
+@named sink1 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
+@named sink2 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
+
+eqns = [
+        connect(n1m2.port_a, sink1.port)
+        connect(n1m2.port_b, sink2.port)
+       ]
+
+@named sys = ODESystem(eqns, t)
+@named n1m2Test = compose(sys, n1m2, sink1, sink2)
+@test_nowarn structural_simplify(n1m2Test)
+
+
+@named n1m2 = N1M2()
+@named pipe1 = AdiabaticStraightPipe()
+@named pipe2 = AdiabaticStraightPipe()
+@named sink1 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
+@named sink2 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
+
+eqns = [
+        connect(n1m2.port_a, pipe1.port_a)
+        connect(pipe1.port_b, sink1.port)
+
+        connect(n1m2.port_b, pipe2.port_a)
+        connect(pipe2.port_b, sink2.port)
+       ]
+
+@named sys = ODESystem(eqns, t)
+@named n1m2AltTest = compose(sys, n1m2, pipe1, pipe2, sink1, sink2)
+@test_nowarn structural_simplify(n1m2AltTest)
+
+
+# N2M2 model and test code.
+function N2M2(;name,
+        kwargs...,
+    )
+
+    @named port_a = TwoPhaseFluidPort()
+    @named port_b = TwoPhaseFluidPort()
+    @named pipe = AdiabaticStraightPipe()
+
+    streams_a = [port_a, pipe.port_a]
+    streams_b = [pipe.port_b, port_b]
+
+    subs = [port_a; port_b; pipe]
+
+    eqns = Equation[]
+
+    push!(eqns, connect(port_a, pipe.port_a))
+    push!(eqns, connect(pipe.port_b, port_b))
+
+    sys = ODESystem(eqns, t, [], [], name=name)
+    sys = compose(sys, subs)
+end
+
+@named n2m2 = N2M2()
+@named source = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
+@named sink = SmallBoundary_Ph(P_in=1e6, h_in=400e3)
+
+eqns = [
+        connect(source.port, n2m2.port_a)
+        connect(n2m2.port_b, sink.port1)
+       ]
+
+@named sys = ODESystem(eqns, t)
+@named n2m2Test = compose(sys, n2m2, source, sink)
+@test_nowarn structural_simplify(n2m2Test)