Network Diffusion

This example explains the use of the basic functions and constructors in NetworkDynamics.jl by modeling a simple diffusion on an undirected network.

This page will:

• walk you through the theoretical background of a simple diffusion propagating in an undirected network
• explain the network dynamics of the system
• explain how to program them

This example can be downloaded as a normal Julia script here.

Theoretical background

Diffusion processes appear in phenomena as diverse as heat conduction, electrical currents, and random walks. Generally speaking they describe the tendency of systems to evolve towards a state of equally distributed entropy (the entropy being in the form of e.g. heat, charge or concentration). If we assume a thermal system, the temperature of a specific spot changes depending on the temperature gradient between the temperature of the spot itself and that of its neighborhood.

We will build a graph $g$ with $N$ nodes and an adjacency matrix $A$, where $v = (v_1, \dots, v_n)$ is the vector of the (abstract) temperatures at each node $i = 1, \dots, N$. The rate of change of state $v_i$ will be described by the difference between the temperature of the node and that of its neighbours. From the above we obtain the following ordinary differential equation:

\[\dot v_i = \sum_{j=1}^N A_{ji} (v_j - v_i).\]

The sum on the right hand side plays the role of a (discrete) gradient. If the temperature at node $i$ is higher than at its neighboring node $j$ it will decrease along that edge.

If a node were to be disconnected from the rest of the network its state would never change, because its adjacency matrix would be $A_{ji} = 0 \; \forall j$ and hence $\dot v_i = 0$. So because its state would remain unchanged the model will be comprised of nodes that have no internal dynamics. This means that the evolution of a node will depend only on its interaction with its neighbors. In NetworkDynamics.jl, interactions between a node and its neighbors are described using edge equations.

In order to bring this equation into the form required by NetworkDynamics.jl we need to split the dynamics into edge dynamics and vertex dynamics and bring them into the correct input-output formulation.

Vertex Dynamics:

All vertices have one internal state $v$, which is also the vertex output. It is the sum over all incoming edge flows connected to the vertex. This directly corresponds to the component model definition outlined in Mathematical Model:

\[\begin{aligned} \dot x^\mathrm{v} &= f^{\mathrm v}(u^{\mathrm v}, \sum_k^{\text{incident}} y^{\mathrm e}_k, p^{\mathrm v}, t) &&= \sum_k^\mathrm{incident} y^\mathrm{e}_k \\ y^\mathrm{v} &= g^{\mathrm v}(u^{\mathrm v}, \sum_k^{\text{incident}} y^{\mathrm e}_k, p^{\mathrm v}, t) &&= x^\mathrm{v} \end{aligned}\]

Edge Dynamics:

The edge dynamics on the other hand do not have any internal states. Thus we can define the edge output as the difference between the source and destination vertices:

\[\begin{aligned} y^{\mathrm e}_{\mathrm{dst}} &= g_\mathrm{dst}^{\mathrm e}(u^{\mathrm e}, y^{\mathrm v}_{\mathrm{src}}, y^{\mathrm v}_{\mathrm{dst}}, p^{\mathrm e}, t) &&= y^{\mathrm v}_{\mathrm{src}} - y^{\mathrm v}_{\mathrm{dst}}\\ y^{\mathrm e}_{\mathrm{src}} &= g_\mathrm{src}^{\mathrm e}(u^{\mathrm e}, y^{\mathrm v}_{\mathrm{src}}, y^{\mathrm v}_{\mathrm{dst}}, p^{\mathrm e}, t) &&= y^{\mathrm v}_{\mathrm{dst}} - y^{\mathrm v}_{\mathrm{src}} \end{aligned}\]

Modelling dynamics in NetworkDynamics.jl

To model the vertex and edge dynamics we need to create a VertexModel and an EdgeModel, respectively.

As a first step we need to install Julia and the necessary packages (Graphs, NetworkDynamics, OrdinaryDiffEqTsit5, StableRNGs and Plots). Then we need to load them:

using Graphs
using NetworkDynamics
using OrdinaryDiffEqTsit5
using StableRNGs
using Plots

Defining the EdgeModel

Then we must define the EdgeModel. To define it we use the function diffusionedge_g! which takes as inputs the current state of the edge e, its source vertex v_src, its destination vertex v_dst, a vector of parameters p and the time t. In order to comply with the syntax of NetworkDynamics.jl we must always define functions for edges with exactly these arguments. (In the case of this example, the values for p and t are not used since there are no parameters and there is no explicit time dependency in the system).

After the function call the edge's output value e equals the difference between its source and its destination vertex (i.e. the discrete gradient along that edge).

function diffusionedge_g!(e_dst, v_src, v_dst, p, t)
    # e_dst, v_src, v_dst are arrays, so we use the broadcasting operator .
    e_dst .= v_src .- v_dst
    nothing
end

Notably, this function only models $g_\mathrm{dst}$. However we can wrap this single-sided output function in an AntiSymmetric output wrapper to construct the EdgeModel:

nd_diffusion_edge = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow])

EdgeModel :StaticEdgeM PureFeedForward()
 ├─   0 states:  []  
 └─ 1/1 outputs: src=[₋flow] dst=[flow]

Defining the VertexModel

Next we need to define the VertexModel. For undirected graphs, the edgefunction! specifies the coupling from a source to a destination vertex. The contributions of the connected edges to a single vertex are "aggregated". By default, the aggregation is the summation of all incident edge states. The aggregated edge state is made available via the esum argument of the vertex function.

function diffusionvertex_f!(dv, v, esum, p, t)
    # dv, v and esum are arrays, so we use the broadcasting operator .
    dv .= esum
    nothing
end

Just like above, the input arguments v, esum, p, t are mandatory for the syntax of vertex functions. The additional input dv corresponding to the derivative of the vertex' state is mandatory for vertices described by ordinary differential equations. The outputs of the vertex function are just a subset of the internal states. Therefore, we can use the StateMask helper function g = StateMask(1:1) to access them.

nd_diffusion_vertex = VertexModel(; f=diffusionvertex_f!, g=StateMask(1:1), dim=1)

VertexModel :VertexM PureStateMap()
 ├─ 1 state:  [s]
 └─ 1 output: [s]

Constructing the network

With the components defined, we can now define the topology and assemble the network dynamics.

The first step is to define the number of nodes the network is comprised of:

N = 20 # number of nodes
k = 4  # average degree distribution

Then we define the network model. We will use the Barabási–Albert model which generates a scale-free random graph.

g = barabasi_albert(N, k)

We then create the network of nodes and edges by using the constructor Network. It combines the component model with the topological information contained in the graph g and returns an Network compatible with the solvers of DifferentialEquations.jl.

nd = Network(g, nd_diffusion_vertex, nd_diffusion_edge)
rng = StableRNG(1)

StableRNGs.LehmerRNG(state=0x00000000000000000000000000000003)

We then generate random initial conditions

x0 = randn(rng, N)

20-element Vector{Float64}:
 -0.5325200748641231
  0.098465514284785
  0.7528865221245234
 -0.8435374038742546
 -2.0138293556311067
 -0.26980764119712713
  1.0322934723039892
  0.6755660991191906
  0.08216082754396425
  0.24692664035915107
 -1.5005291473189497
  0.32863200921212465
 -1.8531987065432076
 -0.33719082654695276
  0.03566251637914036
  0.23214355889044092
  0.556345368618102
 -1.0278210092684674
  1.206403573149015
 -0.6205921994412764

We then solve the diffusion problem on the time interval $[0, 2]$ with the Tsit5() algorithm, which is recommended by the authors of DifferentialEquations.jl for most non-stiff problems. We first create an ODEProblem:

ode_prob = ODEProblem(nd, x0, (0.0, 2.0))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 2.0)
u0: 20-element Vector{Float64}:
 -0.5325200748641231
  0.098465514284785
  0.7528865221245234
 -0.8435374038742546
 -2.0138293556311067
 -0.26980764119712713
  1.0322934723039892
  0.6755660991191906
  0.08216082754396425
  0.24692664035915107
 -1.5005291473189497
  0.32863200921212465
 -1.8531987065432076
 -0.33719082654695276
  0.03566251637914036
  0.23214355889044092
  0.556345368618102
 -1.0278210092684674
  1.206403573149015
 -0.6205921994412764

And then solve the network using solve:

sol = solve(ode_prob, Tsit5());

We plot the results using the plot command, which takes two parameters as input sol and idxs. The first parameter is the solved network sol. The second parameter is idxs, which provides a set of indices. We can use "symbolic" indices, which specify specific components and their symbols directly (see Symbolic Indexing for more details). In order to collect multiple indices we can use the helper function vidxs and eidxs, which help us collect all symbolic indices matching specific criteria. Here, we want to plot all vertex states from the network.

plot(sol; idxs=vidxs(nd, :, :))

In the plot we see, how the individual vertex states evolve over time. We start at a random configuration, and then develop towards an equilibrium where all nodes have the same state, which is to be expected for a simple diffusion process.

Two Dimensional Extension

To illustrate a very simple multi-dimensional case, we simulate two independent diffusions on an identical graph. The first diffusion uses the symbol x and is initiated with initial conditions drawn from the standard normal distribution $N(0,1)$, while the second diffusion uses the symbol ϕ with squared standard normal initial conditions.

The symbols have to be passed with the keyword sym to VertexModel.

Once again we define the number of nodes and the type of network we want to generate:

N = 10 # number of nodes
k = 4  # average degrees distribution
g = barabasi_albert(N, k)

{10, 24} undirected simple Int64 graph

We have two independent diffusions on the network, hence dim = 2. We now define the VertexModel and the EdgeModel

nd_diffusion_vertex_2 = VertexModel(; f=diffusionvertex_f!, g=1:2, dim=2, sym=[:x, :ϕ])
nd_diffusion_edge_2 = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow_x, :flow_ϕ])
nd_2 = Network(g, nd_diffusion_vertex_2, nd_diffusion_edge_2)

Network with 20 states and 0 parameters
 ├─ 10 vertices (1 unique type)
 └─ 24 edges (1 unique type)
Edge-Aggregation using SequentialAggregator(+)

We define the first diffusion as x ~ $N(0,1)^\mathrm{2}$ and the second diffusion as ϕ ~ $N(0,1)$. So the propagation of the diffusion from node 0 to node 2 is given by creating a vector:

x0_2 = vec(transpose([randn(rng, N) .^ 2 randn(rng, N)]))

20-element reshape(transpose(::Matrix{Float64}), 20) with eltype Float64:
  0.374921338987723
  0.3213483657325198
  2.224257343893195
  0.6882831873751686
  0.45670519258445863
 -0.3899398874160884
  0.5643370627576133
  0.05438833466615296
  0.24399248696097944
  0.7283022866605751
  0.17112343782893982
  1.331489923091551
  0.309857463575159
  0.4986114118380903
  1.5826182355628717
 -0.8584115280138382
  0.2922646350323511
  0.9910360582065267
  0.5301523637072203
  0.03813525217366319

We then define the ODEProblem:

ode_prob_2 = ODEProblem(nd_2, x0_2, (0.0, 3.0))

ODEProblem with uType Base.ReshapedArray{Float64, 1, LinearAlgebra.Transpose{Float64, Matrix{Float64}}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 3.0)
u0: 20-element reshape(transpose(::Matrix{Float64}), 20) with eltype Float64:
  0.374921338987723
  0.3213483657325198
  2.224257343893195
  0.6882831873751686
  0.45670519258445863
 -0.3899398874160884
  0.5643370627576133
  0.05438833466615296
  0.24399248696097944
  0.7283022866605751
  0.17112343782893982
  1.331489923091551
  0.309857463575159
  0.4986114118380903
  1.5826182355628717
 -0.8584115280138382
  0.2922646350323511
  0.9910360582065267
  0.5301523637072203
  0.03813525217366319

Then solve the network using:

sol_2 = solve(ode_prob_2, Tsit5());

To plot the evolution of variables ϕ_i over time we use the command:

plot(sol_2; idxs=vidxs(nd_2, :, :x))

[To write ϕ in the terminal type \phi and press TAB]

Using the eidxs helper function we can also plot the evolution of the flow variables over time:

plot(sol_2; idxs=eidxs(nd_2, :, :flow_x))

As expected, the flows are nonzero first and go towards zero as we reach the equilibrium point.

Putting it all together

using Graphs
using NetworkDynamics
using OrdinaryDiffEqTsit5
using StableRNGs
using Plots

function diffusionedge_g!(e_dst, v_src, v_dst, p, t)
    # e_dst, v_src, v_dst are arrays, so we use the broadcasting operator .
    e_dst .= v_src .- v_dst
    nothing
end

nd_diffusion_edge = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow])

function diffusionvertex_f!(dv, v, esum, p, t)
    # dv, v and esum are arrays, so we use the broadcasting operator .
    dv .= esum
    nothing
end

nd_diffusion_vertex = VertexModel(; f=diffusionvertex_f!, g=StateMask(1:1), dim=1)

N = 20
k = 4
g = barabasi_albert(N, k)

nd = Network(g, nd_diffusion_vertex, nd_diffusion_edge)
rng = StableRNG(1)
x0 = randn(rng, N)
ode_prob = ODEProblem(nd, x0, (0.0, 2.0))
sol = solve(ode_prob, Tsit5());
plot(sol; idxs=vidxs(nd, :, :), fmt=:png)
N = 10 #
k = 4  #
g = barabasi_albert(N, k)
nd_diffusion_vertex_2 = VertexModel(; f=diffusionvertex_f!, g=1:2, dim=2, sym=[:x, :ϕ])
nd_diffusion_edge_2 = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow_x, :flow_ϕ])
nd_2 = Network(g, nd_diffusion_vertex_2, nd_diffusion_edge_2)
x0_2 = vec(transpose([randn(rng, N) .^ 2 randn(rng, N)]))
ode_prob_2 = ODEProblem(nd_2, x0_2, (0.0, 3.0))
sol_2 = solve(ode_prob_2, Tsit5());
plot(sol_2; idxs=vidxs(nd_2, :, :x), xlabel = "x", ylabel = "nd_2",)
plot(sol_2; idxs=eidxs(nd_2, :, :flow_x))

Appendix: The network Laplacian $L$

The diffusion equation on a network can be rewritten as:

\[\dot v_i = \sum_{j=1}^N A_{ji} v_j - d_i v_i = e_i^T A v - d_i v_i\]

where $d_i$ is the degree distribution of node $i$ and $e_i^T$ is the $i$-th standard basis vector.

By introducing the diagonal matrix $D$, that has the degree of node $i$ in its $i$-th row and the Laplacian matrix $L = D - A$ we arrive at:

\[\dot v = e_i^T(A - D) v\]

and finally

\[\dot v = - L v\]

This is a linear system of ordinary differential equations (ODEs) and its solution is a matrix exponential. To study the asymptotic behaviour of the system it suffices to analyze the eigenspectrum of $L$. For this reason $L$ is an important construction in network science.

This page was generated using Literate.jl.