Heterogeneous Kuramoto Network

Three vertex types on a ring: standard, static, and inertia

Replicates the full NetworkDynamics.jl Heterogeneous System tutorial — an 8-node ring network with three different vertex types:

• Kuramoto (6 nodes): standard phase oscillator \(\dot{\theta} = \omega_0 + \text{esum}\)
• StaticNode (1 node): no dynamics, outputs a fixed phase parameter
• KuramotoInertia (1 node): second-order oscillator with angular velocity \(\dot{\omega} = \omega_0 - \omega + \text{esum}\)

Per-node parameter overrides set heterogeneous natural frequencies.

YAML Specification

# =============================================================================
# Heterogeneous Kuramoto Network
# =============================================================================
# Source: https://juliadynamics.github.io/NetworkDynamics.jl/stable/generated/heterogeneous_system/
#
# Three different vertex types on the same graph:
#   1. Kuramoto oscillator (nodes 2-4, 6-8): dtheta/dt = omega0 + esum
#   2. Static node (node 1): output = theta_fix (no state variables)
#   3. Kuramoto with inertia (node 5): 2D system with damping
#
# All vertices are coupled via K sin(theta_src - theta_dst).
# =============================================================================

id: 104
label: "Heterogeneous Kuramoto Network"
description: >
    Kuramoto oscillators with three vertex types: standard Kuramoto oscillators
    (6 nodes), one static node with fixed phase output, and one oscillator
    with inertia (second-order dynamics).  Demonstrates heterogeneous vertex
    assignment on a Watts-Strogatz ring network.
    Recreates the NetworkDynamics.jl Heterogeneous System tutorial.
references:
    - "https://juliadynamics.github.io/NetworkDynamics.jl/stable/generated/heterogeneous_system/"
    - "Kuramoto, Y. (1975). Self-entrainment of a population of coupled non-linear oscillators."

# -- Default (most common) vertex dynamics ------------------------------------
dynamics:
    name: Kuramoto
    label: "Kuramoto phase oscillator"
    description: >
        Standard Kuramoto oscillator: dtheta/dt = omega0 + esum.
        Output (coupling variable) is theta.
    system_type: continuous
    autonomous: true
    parameters:
        omega0:
            label: "Natural frequency"
            symbol: "omega_0"
            value: 0.0
            unit: "rad/s"
            description: "Intrinsic angular frequency of the oscillator"
    coupling_terms:
        c_coupling:
            description: "Aggregated sinusoidal coupling from incident edges"
    state_variables:
        theta:
            label: "Phase angle"
            symbol: "theta"
            equation:
                rhs: "omega0 + c_coupling"
            initial_value: 0.0
            unit: "rad"
            coupling_variable: true
            variable_of_interest: true

# -- Network -------------------------------------------------------------------
network:
    label: "8-node ring (Watts-Strogatz k=2, p=0)"
    description: >
        Regular ring lattice: each node connected to 2 nearest neighbours.
        Equivalent to watts_strogatz(8, 2, 0) in Graphs.jl.
    number_of_nodes: 8
    graph_generator:
        name: watts_strogatz
        type: WattsStrogatz
        parameters:
            k: { name: k, value: 2 }
            p: { name: p, value: 0 }

    # -- Additional vertex dynamics models ----------------------------------------
    dynamics:
        StaticNode:
            name: StaticNode
            label: "Static vertex"
            description: >
                No internal dynamics.  Outputs a fixed phase value theta_fix as
                its coupling output.  Used for boundary conditions or fixed
                reference oscillators.  Has no state variables — it is a pure
                algebraic output node (ND.jl: ff=NoFeedForward()).
            system_type: continuous
            autonomous: true
            parameters:
                theta_fix:
                    label: "Fixed phase output"
                    symbol: "theta_fix"
                    value: -0.4375
                    unit: "rad"
                    description: "Constant phase output of this vertex"
            coupling_terms: {}
            state_variables: {}

        KuramotoInertia:
            name: KuramotoInertia
            label: "Kuramoto oscillator with inertia"
            description: >
                Second-order Kuramoto oscillator with damping:
                dtheta/dt = omega, domega/dt = omega0 - 1.0*omega + esum.
                Only theta is the coupling output (g=1:1 in ND.jl).
            system_type: continuous
            autonomous: true
            parameters:
                omega0:
                    label: "Natural frequency"
                    symbol: "omega_0"
                    value: 0.0
                    unit: "rad/s"
                    description: "Intrinsic angular frequency"
            coupling_terms:
                c_coupling:
                    description: "Aggregated sinusoidal coupling from incident edges"
            state_variables:
                theta:
                    label: "Phase angle"
                    symbol: "theta"
                    equation:
                        rhs: "omega"
                    initial_value: 0.0
                    unit: "rad"
                    coupling_variable: true
                    variable_of_interest: true
                omega:
                    label: "Angular velocity"
                    symbol: "omega"
                    equation:
                        rhs: "omega0 - 1.0*omega + c_coupling"
                    initial_value: 5.0
                    unit: "rad/s"
                    variable_of_interest: true

    nodes:
        # Node 1: static vertex (fixed phase output)
        - id: 0
          label: "Static node"
          dynamics: StaticNode
          parameters:
              - name: theta_fix
                value: -0.4375

        # Nodes 2-4: standard Kuramoto
        - id: 1
          dynamics: Kuramoto
          parameters:
              - name: omega0
                value: -0.3125
          state:
              theta: -0.3125
        - id: 2
          dynamics: Kuramoto
          parameters:
              - name: omega0
                value: -0.1875
          state:
              theta: -0.1875
        - id: 3
          dynamics: Kuramoto
          parameters:
              - name: omega0
                value: -0.0625
          state:
              theta: -0.0625

        # Node 5: Kuramoto with inertia (2D)
        - id: 4
          dynamics: KuramotoInertia
          parameters:
              - name: omega0
                value: 0.0625
          state:
              theta: 0.0625
              omega: 5.0

        # Nodes 6-8: standard Kuramoto
        - id: 5
          dynamics: Kuramoto
          parameters:
              - name: omega0
                value: 0.1875
          state:
              theta: 0.1875
        - id: 6
          dynamics: Kuramoto
          parameters:
              - name: omega0
                value: 0.3125
          state:
              theta: 0.3125
        - id: 7
          dynamics: Kuramoto
          parameters:
              - name: omega0
                value: 0.4375
          state:
              theta: 0.4375

    coupling:
        KuramotoCoupling:
            name: KuramotoCoupling
            label: "Kuramoto sinusoidal coupling"
            description: >
                Antisymmetric coupling: K sin(theta_src - theta_dst).
                The AntiSymmetric wrapper ensures g_src = -g_dst.
            delayed: false
            parameters:
                K:
                    label: "Coupling strength"
                    value: 3.0
                    unit: "1"
                    description: "Global coupling strength"
            pre_expression:
                rhs: "K * sin(x_j - x_i)"

# -- Integration ---------------------------------------------------------------
integration:
    description: "Explicit Runge-Kutta (Tsit5) with dt=0.05 over t in [0, 10]"
    method: Tsit5
    step_size: 0.05
    duration: 10.0

Model Report

Heterogeneous Kuramoto Network

Kuramoto oscillators with three vertex types: standard Kuramoto oscillators (6 nodes), one static node with fixed phase output, and one oscillator with inertia (second-order dynamics). Demonstrates heterogeneous vertex assignment on a Watts-Strogatz ring network. Recreates the NetworkDynamics.jl Heterogeneous System tutorial.

---

1. Brain Network: 8-node ring (Watts-Strogatz k=2, p=0)

Regular ring lattice: each node connected to 2 nearest neighbours. Equivalent to watts_strogatz(8, 2, 0) in Graphs.jl.

• Regions: 8

Coupling: Kuramoto sinusoidal coupling

Antisymmetric coupling: K sin(theta_src - theta_dst). The AntiSymmetric wrapper ensures g_src = -g_dst.

Receives \(\theta\) from connected regions.

Pre-synaptic: \(c_{\text{pre}} = - K \cdot \sin{\left(x_{i} - x_{j} \right)}\)

Parameter	Value	Description
\(K\)	3.0	Global coupling strength

2. Local Dynamics: Kuramoto phase oscillator

Standard Kuramoto oscillator: dtheta/dt = omega0 + esum. The model comprises 1 state variables.

2.1 State Equations

\[\dot{\theta} = c_{coupling} + \omega_{0}\]

2.2 Parameters

Parameter	Value	Unit	Description
\(\omega_{0}\)	0.0	rad_per_s	Intrinsic angular frequency of the oscillator

3. Numerical Integration

• Method: Tsit5
• Time step: \(\Delta t = 0.05\) ms
• Duration: 10.0 ms

---

References

Cabral, J., Hugues, E., Sporns, O., & Deco, G. (2011). Role of local network oscillations in resting-state functional connectivity. NeuroImage, 57(1), 130-139.

Kuramoto, Y. (1975). Self-entrainment of a population of coupled non-linear oscillators. Lecture Notes in Physics, 420-422.

Strogatz, S. (2000). From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D: Nonlinear Phenomena, 143(1–4), 1-20. - https://juliadynamics.github.io/NetworkDynamics.jl/stable/generated/heterogeneous_system/ - Kuramoto, Y. (1975). Self-entrainment of a population of coupled non-linear oscillators.

Generated Julia Code

from tvbo import SimulationExperiment
exp = SimulationExperiment.from_file("yaml/heterogeneous_kuramoto.yaml")
Code(exp.render_code("networkdynamics"), language='julia')

using Graphs
using NetworkDynamics

using OrdinaryDiffEqTsit5







function Kuramoto_f!(dx, x, esum, (omega0,), t)
    theta = x[1]

    c_coupling = esum[1]
    dx[1] = c_coupling + omega0
    nothing
end

vertex_Kuramoto = VertexModel(;
    f = Kuramoto_f!,
    g = StateMask(1:1),
    sym = [:theta],
    psym = [:omega0 => 0.0],
    name = :Kuramoto,
)








function KuramotoInertia_f!(dx, x, esum, (omega0,), t)

    omega, theta = x

    c_coupling = esum[1]
    dx[1] = c_coupling + omega0 - 1.0 * omega
    dx[2] = omega
    nothing
end

vertex_KuramotoInertia = VertexModel(;
    f = KuramotoInertia_f!,
    g = StateMask(2:2),
    sym = [:omega, :theta],
    psym = [:omega0 => 0.0],
    name = :KuramotoInertia,
)



function StaticNode_g!(out, x, p, t)
    out .= p
    nothing
end
vertex_StaticNode = VertexModel(;
    g = StaticNode_g!,
    outsym = [:theta],
    psym = [:theta_fix],
    ff = NoFeedForward(),
    name = :StaticNode,
)







function KuramotoCoupling_edge_g!(e_dst, v_src, v_dst, (K,), t)
    e_dst[1] = -K .* sin(v_dst[1] - v_src[1])
    nothing
end

edge_KuramotoCoupling = EdgeModel(;
    g = AntiSymmetric(KuramotoCoupling_edge_g!),
    outsym = [:coupling],
    psym = [:K => 3.0],
    name = :KuramotoCoupling,
)




g = watts_strogatz(8, 2, 0)


vertex_array = VertexModel[vertex_Kuramoto for _ in 1:8]
vertex_array[1] = vertex_StaticNode
vertex_array[5] = vertex_KuramotoInertia

nw = Network(g, vertex_array, edge_KuramotoCoupling; dealias=true)




s = NWState(nw)


s.p.v[1, :theta_fix] = -0.4375


s.v[2, :theta] = -0.3125

s.p.v[2, :omega0] = -0.3125


s.v[3, :theta] = -0.1875

s.p.v[3, :omega0] = -0.1875


s.v[4, :theta] = -0.0625

s.p.v[4, :omega0] = -0.0625


s.v[5, :omega] = 5.0

s.v[5, :theta] = 0.0625

s.p.v[5, :omega0] = 0.0625


s.v[6, :theta] = 0.1875

s.p.v[6, :omega0] = 0.1875


s.v[7, :theta] = 0.3125

s.p.v[7, :omega0] = 0.3125


s.v[8, :theta] = 0.4375

s.p.v[8, :omega0] = 0.4375



tspan = (0.0, 10.0)
prob = ODEProblem(nw, uflat(s), tspan, pflat(s))
sol = solve(prob, Tsit5(); saveat=0.05)

adj_matrix = Float64.(adjacency_matrix(g))

using Random: MersenneTwister
function spring_layout(g; seed=42, iterations=50, k=1.0)
    rng = MersenneTwister(seed)
    n = nv(g)
    pos = randn(rng, n, 2)
    for _ in 1:iterations
        disp = zeros(n, 2)
        for i in 1:n, j in (i+1):n
            d = pos[i, :] - pos[j, :]
            dist = max(norm(d), 0.01)
            rep = k^2 / dist
            disp[i, :] .+= d / dist * rep
            disp[j, :] .-= d / dist * rep
        end
        for e in edges(g)
            i, j = src(e), dst(e)
            d = pos[j, :] - pos[i, :]
            dist = max(norm(d), 0.01)
            att = dist^2 / k
            disp[i, :] .+= d / dist * att
            disp[j, :] .-= d / dist * att
        end
        for i in 1:n
            dl = max(norm(disp[i, :]), 0.01)
            pos[i, :] .+= disp[i, :] / dl * min(dl, 0.1)
        end
    end
    return pos
end
using LinearAlgebra: norm
node_positions = spring_layout(g)

using Plots
plot(sol; ylabel="state", xlabel="time", title="Kuramoto on 8-node network")

Run & Plot

ts = exp.run(format="networkdynamics")
ts.plot()

Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython