Physarum Can Compute Shortest Paths∗
Vincenzo Bonifaci†

Kurt Mehlhorn‡

Girish Varma§

arXiv:1106.0423v3 [] 12 Oct 2011

October 13, 2011

Abstract
Physarum Polycephalum is a slime mold that is apparently able to solve shortest path
problems. A mathematical model has been proposed by biologists to describe the feedback
mechanism used by the slime mold to adapt its tubular channels while foraging two food
sources s0 and s1 . We prove that, under this model, the mass of the mold will eventually
converge to the shortest s0 -s1 path of the network that the mold lies on, independently
of the structure of the network or of the initial mass distribution.
This matches the experimental observations by the biologists and can be seen as an
example of a “natural algorithm”, that is, an algorithm developed by evolution over
millions of years.

Contents
1 Introduction

2

2 Related Work

6

3 Discussion and Open Problems

6

4 Parallel Links

7

5 Electrical Networks and Simple Facts

9

6 Convergence

11

7 Rate of Convergence for Stable Flow Directions

18

8 The Wheatstone Graph

25

9 The Uncapacitated Transportation Problem

28

∗
An extended abstract of this paper appears in SODA (ACM-SIAM Symposium on Discrete Algorithms)
2012.
†
Istituto di Analisi dei Sistemi ed Informatica “Antonio Ruberti” – CNR, Rome, Italy. Most of the work
was done at the MPI for Informatics, Saarbrücken, Germany. Email: bonifaci@mpi-inf.mpg.de.
‡
MPI for Informatics, Saarbrücken, Germany. Email: mehlhorn@mpi-inf.mpg.de.
§
Tata Institute of Fundamental Research, Mumbai, India. Most of the work was done at the MPI for
Informatics, Saarbrücken, Germany. Email: girish@tcs.tifr.res.in.

Figure 1: The experiment in [NYT00] (reprinted from there): (a) shows the maze uniformly
covered by Physarum; the yellow color indicates the presence of Physarum. Food (oatmeal)
is provided at the locations labelled AG. After a while, the mold retracts to the shortest path
connecting the food sources as shown in (b) and (c). (d) shows the underlying abstract graph.
The video [You10] shows the experiment.

1

Introduction

Physarum Polycephalum is a slime mold that is apparently able to solve shortest path problems. Nakagaki, Yamada, and Tóth [NYT00] report on the following experiment, see Figure 1:
They built a maze, covered it with pieces of Physarum (the slime can be cut into pieces that
will reunite if brought into vicinity), and then fed the slime with oatmeal at two locations.
After a few hours, the slime retracted to a path that follows the shortest path connecting the
food sources in the maze. The authors report that they repeated the experiment with different mazes; in all experiments, Physarum retracted to the shortest path. There are several
videos available on the web that show the mold in action [You10].
Paper [TKN07] proposes a mathematical model for the behavior of the slime and argues
extensively that the model is adequate. We will not repeat the discussion here but only define
the model. Physarum is modeled as an electrical network with time varying resistors. We
have a simple undirected graph G = (N, E) with distinguished nodes s0 and s1 , which model
the food sources. Each edge e ∈ E has a positive length Le and a positive diameter De (t);
Le is fixed, but De (t) is a function of time. The resistance Re (t) of e is Re (t) = Le /De (t).
We force a current of value 1 from s0 to s1 . Let Qe (t) be the resulting current over edge
e = (u, v), where (u, v) is an arbitrary orientation of e. The diameter of any edge e evolves
according to the equation
Ḋe (t) = |Qe (t)| − De (t),
(1)
where Ḋe is the derivative of De with respect to time. In equilibrium (Ḋe = 0 for all e), the
flow through any edge is equal to its diameter. In non-equilibrium, the diameter grows or
shrinks if the absolute value of the flow is larger or smaller than the diameter, respectively. In
the sequel, we will mostly drop the argument t as is customary in the treatment of dynamical
systems.
The model is readily turned into a computer simulation. In an electrical network, every
vertex v has a potential pv ; pv is a function of time. We may fix ps1 to zero. For an edge
e = (u, v), the flow across e is given by (pu − pv )/Re . We have flow conservation in every

2

vertex except for s0 and s1 ; we inject one unit at s0 and remove one unit at s1 . Thus,
bv =

X pv − pu
,
Ruv

(2)

u∈δ(v)

where δ(v) is the set of neighbors of v and bs0 = 1, bs1 = −1, and bv = 0 otherwise. The
node potentials can be computed by solving a linear system (either directly or iteratively).
Tero et al. [TKN07] were the first to perform simulations of the model. They report that
the network always converges to the shortest s0 -s1 path, i.e., the diameters of the edges on
the shortest path converge to one, and the diameters on the edges outside the shortest path
converge to zero. This holds true for any initial condition and assumes the uniqueness of the
shortest path.
Miyaji and Ohnishi [MO07, MO08] initiated the analytical investigation of the model.
They argued convergence against the shortest path if G is a planar graph and s0 and s1 lie
on the same face in some embedding of G.
Our main result is a convergence proof for all graphs. For a network G = (V, E, s0 , s1 , L),
where (Le )e∈E is a positive length function on the edges of G, we use G0 = (V, E0 ) to denote
the subgraph of all shortest source-sink paths, L∗ to denote the length of a shortest sourcesink path, and E ∗ to P
denote the set of all source-sink flows of value one in G0 . If we define
the cost of flow Q as e Le Qe , then E ∗ is the set of minimum cost source-sink flows of value
one. If the shortest source-sink path is unique, E ∗ is a singleton. The dynamics are attracted
by a set A ⊆ RE if the distance (measured in any Lp -norm) between D(t) and A converges
to zero.
Theorem [Theorem 2 in Section 6] Let G = (V, E, s0 , s1 , L) be an undirected network
with positive length function (Le )e∈E . Let De (0) > 0 be the diameter of edge e at time zero.
The dynamics (1) are attracted to E ∗ . If the shortest source-sink path is unique, the dynamics
converge to the flow of value one along the shortest source-sink path.
We conjecture that the dynamics converge to an element of E ∗ but only show attraction
to E ∗ . A key part of our proof is to show that the function
V =

1
minS∈C CS

X

Le De + (C{s0 } − 1)2

(3)

e∈E

decreases along all trajectories that start in a non-equilibrium configuration. Here,
P C is the
set of all s0 -s1 cuts, i.e., the set of all S ⊆ N with s0 ∈ S and s1 6∈ S; CS = e∈δ(S) De
is the capacity of the cut S when the capacity of edge e is set to De ; and minS∈C CS (also
abbreviated by C) is the capacity of the minimum cut. The first term in the definition of V
is the normalized hardware cost; for any edge, the product of its length and its diameter may
be interpreted as the hardware cost of the edge; the normalization is by the capacity of the
minimum cut. We will show that the first term decreases except when the maximum flow
F in the network with capacities De is unique, and moreover, |Qe | = |Fe |/C for all e. The
second term decreases as long as the capacity of the cut defined by s0 is different from 1. We
show that the capacity
Pof the minimum cut converges to one and that the derivative of V is
upper bounded by − e (Lmin /4)(De /C − |Qe |)2 , where Lmin is the minimum length of any
edge. Since V is non-negative, this will allow us to conclude that |De − |Qe || converges to zero
for all e. In the next step, we show that the potential difference ∆ = ps0 − ps1 between source
and sink converges to the length L∗ of a shortest-source sink path. We use this to conclude
3

s0

s1

e1

u
c

a
e2

e4

e3
u

s0

s1

e

v
e5

b

e6

d
v

w
(a)

(b)

Figure 2: Part (a) illustrates the path decomposition. All edges are assumed to have length
1; P0 = (e1 ), P1 = (e2 , e3 , e4 ), P2 = (e5 , e6 ), p∗s0 = 1, p∗s1 = 0, p∗v = 1/3, p∗u = 2/3, p∗w = 1/2,
f (P1 ) = 1/3, and f (P2 ) = 1/6.
Part (b) shows the Wheatstone graph. The direction of the flow on edge {u, v} may change
over time; the flow on all other edges is always from left to right.
that De and Qe converge to zero for any edge e 6∈ E0 . Finally, we show that the dynamics
are attracted by E ∗ .
We found the function V by analytical investigation of a network of parallel links (see
Section 4), extensive computer simulations, and guessing. Functions decreasing along all
trajectories are called Lyapunov functions in dynamical systems theory [HS74]. The fact that
the right-hand side of system (1) is not continuously differentiable and that the function V
is not differentiable everywhere introduces some technical difficulties.
The direction of the flow across an edge depends on the initial conditions and time. We
do not know whether flow directions can change infinitely often or whether they become
ultimately fixed. Under the assumption that flow directions stabilize, we can characterize
the (late stages of the) convergence process. An edge e = {u, v} becomes horizontal if
limt→∞ |pu − pv | = 0, and it becomes directed from u to v (directed from v to u) if pu > pv
for all large t (pv > pu for all large t). An edge stabilizes if it either becomes horizontal or
directed, and a network stabilizes if all its edges stabilize. If a network stabilizes, we partition
→
−
its edges into a set Eh of horizontal edges and a set E of directed edges. If {u, v} becomes
→
−
directed from u to v, then (u, v) ∈ E .
We introduce the notion of a decay rate. Let r ≤ 0. A quantity D(t) decays with rate at
least r if for every ε > 0 there is a constant A such that ln D(t) ≤ A + (r + ε)t for all t. A
quantity D(t) decays with rate at most r if for every ε > 0 there is a constant a such that
ln D(t) ≥ a + (r − ε)t for all t. A quantity D(t) decays with rate r if it decays with rate at
least and at most r.
Lemma [Lemma 20 in Section 7] For e ∈ Eh , De decays with rate −1 and |Qe | decays
with rate at least −1.
We define a decomposition of G into paths P0 to Pk , an orientation of these paths, a slope
f (Pi ) for each Pi , a vertex labelling p∗ , and an edge labelling r. P0 is a1 shortest s0 -s1 path in
G, f (P0 ) = 1, re = f (P0 )−1 for all e ∈ P0 , and p∗v = dist(v, s1 ) for all v ∈ P0 , where dist(v, s1 )
1

We assume that P0 is unique.

4

is the shortest path distance from v to s1 . For 1 ≤ i ≤ k, we have2 Pi = argmaxP ∈P f (P ),
where P is the set of all paths P in G with the following properties: (1) the startpoint a
and the endpoint b of P lie on P0 ∪ . . . ∪ Pi−1 , p∗a ≥ p∗b , and f (P ) = (p∗a − p∗b )/L(P ); (2) no
interior vertex of P lies on P0 ∪ . . . ∪ Pi−1 ; and (3) no edge of P belongs to P0 ∪ . . . ∪ Pi−1 .
If p∗a > p∗b , we direct Pi from a to b. If p∗a = p∗b , we leave the edges in Pi undirected. We set
re = f (Pi ) − 1 for all edges of Pi , and p∗v = p∗b + f (Pi ) distPi (v, b) for every interior vertex v
of Pi . Figure 2(a) illustrates the path decomposition.
Lemma [Lemma 21 in Section 7] There is an i0 ≤ k such that
f (P0 ) > f (P1 ) > . . . > f (Pi0 ) > 0 = f (Pi0 +1 ) = . . . = f (Pk ).

→
−
Theorem [Theorem 3 in Section 7] If a network stabilizes, E = ∪i≤i0 E(Pi ), the orien→
−
tation of any edge e ∈ E agrees with the orientation induced by the path decomposition, and
Eh = ∪i>i0 E(Pi ). The potential of each node v converges to p∗v . The diameter of each edge
e ∈ E \ P0 decays with rate re .
We cannot prove that flow directions stabilize in general. For series-parallel graphs, flow
directions trivially stabilize. The Wheatstone graph, shown in Figure 2(b), is the simplest
graph, in which flow directions may change over time.
Theorem [Theorem 6 in Section 8] The Wheatstone graph stabilizes.
The uncapacitated transportation problem generalizes the shortest
P path problem. With
each vertex v, a supply/demand bv is associated. It is assumed that v bv = 0. Nodes with
positive bv are called supply nodes, and nodes with negative bv are called demand nodes.
In the shortest path problem, exactly two vertices have non-zero supply/demand. A feasible
solution to the transportation problem is a flow f satisfying the mass balance constraints,
i.e.,
P
for every vertex v, bv is equal to the net flow out of v. The cost of a solution is e Le fe . The
Physarum solver for the transportation problem is as follows: At any fixed time, the potentials
are defined by (2) and the currents (Qe )e∈E are derived from the potentials by Ohm’s law.
The dynamics evolve according to (1). The equilibria, i.e., |Qe | = De for all e, are precisely
the flows with the following equal-length property. Orient the edges in the direction of Q and
drop the edges of flow zero. In the resulting graph, any two distinct directed paths with the
same source and sink have the same length. Let E be the set of equilibria.
Theorem [Theorem 8 in Section 9] The dynamics (1) are attracted to the set of equilibria
E. If any two equilibria have distinct cost, the dynamics converge to an optimum solution of
the transportation problem.
The convergence statement for the transportation problem is weaker than the corresponding statement for the shortest path problem in two respects. There, we show attraction to
the set of equilibria of minimum cost (now only to the set of equilibria) and convergence to
the optimum solution if the optimum solution is unique (now only if no two equilibria have
the same cost).
This paper is organized as follows: In Section 2, we discuss related work, and in Section 3,
we put our results into the context of natural algorithms and state open problems. The
2

We assume that Pi is unique except if f (Pi ) = 0.

5

technical part of the paper starts in Section 4. We first treat a network of parallel links; this
situation is simple enough to allow an analytical treatment. In Section 5, we review basic
facts about electrical networks and prove some simple facts about the dynamics of Physarum.
In Section 6, we prove our main result, the convergence for general graphs. In Section 7,
we prove exponential convergence under the assumption that flow directions stabilize, and
in Section 8, we show that the Wheatstone network stabilizes. Finally, in Section 9, we
generalize the convergence proof to the transportation problem.

2

Related Work

Miyaji and Ohnishi [MO07, MO08] initiated the analytical investigation of the model. They
argued convergence against the shortest path if G is a planar graph and s0 and s1 lie on the
same face in some embedding of G. Ito et al. [IJNT11] study the dynamics (1) in a directed
graph G = (V, E); they do not claim that the model is justified on biological grounds. Each
directed edge e has a diameter De . The node potentials are again defined by the equations
X pv − pu
for all v ∈ V .
bv =
Ruv
u∈δ(v)

The summation on the right-hand side is over all neighbors u of v; edge directions do not
matter in this equation. If there is an edge from u to v and an edge from v to u, u occurs
twice in the summation, once for each edge. The dynamics for the diameter of the directed
edge uv are then Ḋuv = Quv − Duv , where Quv = Duv (pu − pv )/Luv . The dynamics of this
model are very different from the dynamics of the model studied in our paper. For example,
assume that there is an edge vu, no edge uv, and pu > pv always. Then Qvu < 0 always and
hence Dvu will vanish at least with rate −1. The model is simpler to analyze than our model.
Ito et al. prove that the directed model is able to solve transportation problems and that the
De ’s converge exponentially to their limit values.

3

Discussion and Open Problems

Physarum may be seen as an example of a natural computer, i.e., a computer developed by
evolution over millions of years. It can apparently do more than compute shortest paths and
solve transportation problems. In [TTS+ 10], the computational capabilities of Physarum are
applied to network design, and it is shown in lab and computer experiments that Physarum
can compute approximate Steiner trees. No theoretical analysis is available. The book [Ada10]
and the tutorial [NTK+ 09] contain many illustrative examples of the computational power of
this slime mold.
Chazelle [Cha09] advocates the study of natural algorithms; i.e., “algorithms developed
by evolution over millions of years”, using computer science techniques. Traditionally, the
analysis of such algorithms belonged to the domain of biology, systems theory, and physics.
Computer science brings new tools. For example, in our analysis, we crucially use the maxflow min-cut theorem. Natural algorithms can also give inspiration for the development of new
combinatorial algorithms. A good example is [CKM+ 11], where electrical flows are essential
for an approximation algorithm for undirected network flow.
We have only started the theoretical investigation of Physarum computation, and so many
interesting questions are open. We prove convergence for the dynamics Ḋe = f (|Qe |) − De ,
6

where f is the identity function. The biological literature also suggests the use of f (x) =
xγ /(1 + xγ ) for some parameter γ. Can one prove convergence for other functions f ? We
prove that flow directions stabilize in the Wheatstone graph. Do they stabilize in general?
We prove, but only for stabilizing networks, that the diameters of edges not on the shortest
path converge to zero exponentially for large t. What can be said about the initial stages of
the process? The Physarum computation is fully distributed; node potentials depend only
on the potentials of the neighbors, currents are determined by potential differences of edge
endpoints, and the update rule for edge diameters is local. Can the Physarum computation
be used as the basis for an efficient distributed shortest path algorithm? What other problems
can be provably solved with Physarum computations?

4

Parallel Links

We discovered the Lyapunov function used in the proof of our main theorem through experimentation. The experimentation was guided by the analysis of a network of parallel links.
In such a network, there are vertices s0 and s1 connected with m edges
P of lengths L1 < L2 <
. . . < Lm . Let Di be the diameter of the i-th link, and let D = i Di . Let ∆ = ps0 − ps1
be
source and sink. Then, Qi = ∆/Ri = Di ∆/Li . Since
P the potential difference between
P
i Qi = 1, we have ∆ = 1/
i Di /Li .

Lemma 1 The equilibrium points are precisely the single links.

Proof: In an equilibrium point, Qi = Di for all i. Since Qi = Di ∆/Li , this implies ∆ = Li
whenever Qi 6= 0. Thus, in an equilibrium there is exactly one i with Qi 6= 0. Then, Qi = 1.

Lemma 2 Let D =

P

i

Di . Then, D converges to 1.

P
P
P
Proof: We have Ḋ = i Ḋi = i Qi − i Di = 1 − D. The claim follows by directly solving
the differential equation: D(t) = 1 + (D(0) − 1) exp(−t).
For networks of parallel links, there are many Lyapunov functions.
Lemma 3 Let xi = Di /D, and let L be such that 1/L =
X

Di /D,

i≥2

X
i

xi Li , L,

X
i

Qi Li , ∆

X
i

P

Di Li , and

j

xj /Lj . The quantities

X
(Li ln Di − L1 ln D1 )
i≥2

decrease along all trajectories, starting in non-equilibrium points.
Proof: Clearly,

P

j

xj = 1 and ∆ = L/D. The derivative ẋi of xi computes as:





L
(Di ∆/Li − Di )D − Di (1 − D)
1
Ḋi D − Di Ḋ
1 L
=
=
−
−
1
xi .
x
=
i
D2
D2
Li D D
D Li
P
We have L > LP
1 iff
j≥2 xj > 0. Thus, the derivative of x1 is zero if x1 = 1 and positive if
x1 < 1. Thus, i≥2 xi decreases along all trajectories, starting in non-equilibrium points.
ẋi =

7

Let V =

P

i

xi Li . Then,

X 1 L
1 X
V̇ =
(L − Li )xi .
− 1 xi L i =
D Li
D
i

i

P

P

P
P
So, it suffices to show i Li xi ≥ L = 1/ i xi /Li , or equivalently, ( i Li xi )( i xi /Li ) ≥ 1.
This is an immediate consequence of the Cauchy-Schwarz inequality. Namely,
!
!
!2
Xp
Xp
Xp
p
( xi /Li )2 .
1=
xi Li xi /Li
( xi Li )2 ·
≤
i

i

i

Now, let V = 1/L =

P

xj /Lj . We show that V is increasing. We have




X ẋi
xi
1 X L
1 X Lxi 1
xi
=
−1
=
−
.
V̇ =
Li
D
Li
Li
D
Li Li Li
j

i

i

i

Let zi = Lxi /Li . Then, zi ≥ xi if L ≥ Li , and zi ≤ xi if L ≤ Li . Also
D · V̇ =

X zi − xi
i

Li

=

P

i zi

= 1. Thus,

X zi − xi
X zi − xi
X zi − xi
X zi − xi
+
≥
+
= 0.
Li
Li
L
L

i:L≥Li

i:L<Li

i:L≥Li

i:L<Li

Moreover, V̇ = 0 if and only if zi = xi for all i if and only if x is a unit vector.
P
Consider next the function i Qi Li . Then,
X

Qi Li =

i

X
i

P

∆

Di
1
D
Li = ∆D = P Di = P xi = L;
Li
i Li
i L
i

hence, i Qi Li is decreasing.
P
P
The function ∆ i Di Li = L · i xi Li is the product of decreasing functions and hence
decreasing.
P
Finally, let V = i≥2 (Li ln Di − L1 ln D1 ). Then
!

X  Qi − Di
X
Ḋ1
Q1 − D1
Ḋi
− L1
− L1
=
Li
V̇ =
Li
Di
D1
Di
D1
i≥2
i≥2
 X
X  Di ∆/Li − Di
D1 ∆/L1 − D1
− L1
=
Li
=
(L1 − Li ) < 0.
Di
D1
i≥2

The Lyapunov function

i≥2

P

i≥2 (Li ln Di

− L1 ln D1 ) was already considered in [MO07].

Theorem 1 (Miyashi-Ohnishi [MO07]) For a network of parallel links, the dynamics
converge against D1 = 1 and Di = 0 for i ≥ 2.
Proof: x1 = D1 /D is monotonically increasing and bounded by 1. Hence, it converges.
Assume
that the limit x∗1 is less than one. Clearly, x∗1 > 0. For x1 ≤ x∗1 , we have 1/L =
P
∗
∗
∗
i xi /Li ≤ x1 /L1 + (1 − x1 )/L2 . Moreover, for large enough t, x1 ≥ x1 /2 and D ≤ 2 (Lemma
2), and hence, x˙1 ≥ ε for some ε > 0. Thus, x∗1 < 1 is impossible.
8

P
Some of the
Lyapunov
functions
have
natural
interpretations:
i Qi Li is the total cost
P
P
of the flow; ( i Di Li )/ i Di is the total hardware cost normalized
P by the total diameter,
where a link of length L and diameter D has cost DL; and ∆ i Di Li is the potential
difference between source and sink multiplied by total hardware cost. These functions are
readily generalized to general networks by interpreting the summations as summations over
all edges of the network. Our computer simulations showed that none of these functions is a
Lyapunov function
for general networks.
P
However, i Di can also be interpreted as the minimum capacity of a source-sink
cut in
P
P
a network where Di is the capacity of edge i. With this interpretation, ( i Di Li )/ i Di
becomes
P
e De Le
,
minS∈C CS
where C is the set of all s0 -s1 cuts and CS is the capacity of the cut C. Our computer
simulations suggested that this function may serve as a Lyapunov function for general graphs.
We will see below that a slight modification is actually a Lyapunov function.

5

Electrical Networks and Simple Facts

In this section, we establish some more notation, review basic properties of electrical networks,
and prove some simple facts.
Each node v of the graph G has a potential pv that is a function of time. A potential
difference ∆e between the endpoints of an edge e induces a flow on the edge. For e = (u, v),
Qe = De ∆e /Le = De (pu − pv )/Le = (pu − pv )/Re

(4)

is the flow across e in the direction from u to v. If Qe < 0, the flow is in the reverse direction.
The potentials are such that there is flow conservation in every vertex except for s0 and s1
and such that the net flow from s0 to s1 is one, that is, for every vertex u, we have
X
Qu,v = b(u),
(5)
v:(u,v)∈E

where b(s0 ) = 1 = −b(s1 ) and b(u) = 0 for all other vertices u. After fixing one potential
to an arbitrary value, say ps1 = 0, the other potentials are readily determined by solving a
linear system. This means that each Qe can be expressed as a function of R only.
For the main convergence proof, we will use some fundamental principles from the theory
of electrical networks (for a complete treatment, see for example [Bol98, Chapters II, IX]).
Thomson’s Principle. TheP
flow Q is uniquely determined as a feasible flow that minimizes
2
the total energy dissipation
e Re Qe , with Re = Le /De . In other words, for any flow x
satisfying (5),
X
X
Re Q2e ≤
Re x2e .
(6)
e

e

Kirchhoff ’s Theorem. For a graph G = (N, E) and an oriented edge e = (u, v) ∈ E, let
• Sp be the set of all spanning trees of G, and let
• Sp(u, v) be the set of all spanning trees T of G, for which the oriented edge (u, v) lies
on the unique path from s0 to s1 in T .
9

For a set of trees S, define Γ(S) =
is

P

Quv =

T ∈S

Q

e∈T

De /Le . Then, the current through the edge e

Γ(Sp(u, v)) − Γ(Sp(v, u))
.
Γ(Sp)

(7)

Gronwall’s Lemma. Let α, β ∈ R and let x be a continuous differentiable real function on
[0, ∞). If αx(t) ≤ ẋ(t) ≤ βx(t) for all t ≥ 0, then
x(0) eαt ≤ x(t) ≤ x(0) eβt
Proof:

for all t ≥ 0.

ẋeβt − βxeβt
x(t)
x(0)
d x
=
≤ 0 ⇒ βt ≤ β0 = x(0).
βt
dt e
e2βt
e
e

A similar calculation establishes x(t) ≥ x(0)eαt .
The next lemma gives some properties that are
Peasily derived from (1), (4), and (5).
Recall that C is the set of s0 -s1 cuts and CS =
e∈δ(S) De . Also, let Lmin = mine Le ,
Lmax = maxe Le , n = |N |, and m = |E|.
Lemma 4 The following hold for any edge e ∈ E and any cut S ∈ C:
(i) |Qe | ≤ 1.
P
(ii)
e∈δ({s0 }) |Qe | = 1.

(iii) De (t) ≥ De (0) exp(−t) for all t,
(iv) De (t) ≤ 1 + (De (0) − 1) exp(−t) for all t.
(v) Re ≥ Lmin /2 for all sufficiently large t.
(vi) CS (t) ≥ 1 + (CS (0) − 1) exp(−t) for all t, with equality if S = {s0 }.
(vii) C{s0 } → 1 as t → ∞.
(viii) Orient the edges according to the direction of the flow. For sufficiently large t, there is
a directed source-sink path, in which all edges have diameter at least 1/2m.
(ix) |∆e | ≤ 2nmLmax for all sufficiently large t.
(x) Ḋe /De ∈ [−1, 2nmLmax /Lmin ] for all sufficiently large t.
Proof:
(i) Since Q is a flow, it can be decomposed into s0 -s1 flow paths and cycles. If |Qe | > 1,
since b(s0 ) = 1, there exists a positive cycle in this decomposition, a contradiction to the
existence of potential values at the nodes. The claim is also an immediate consequence
of (7).
(ii) It follows
P from equations
P (4) and (5) that ps0 = maxv pv , so Qs0 ,v ≥ 0 for all {s0 , v} ∈ E,
and e∈δ({s0 }) |Qe | = e∈δ({s0 }) Qe = 1.
10

(iii) From the evolution equation (1), Ḋe ≥ −De . The claim follows by Gronwall’s Lemma.
(iv) |Qe | ≤ 1 for any edge e, so Ḋe ≤ 1 − De from (1), and the claim follows as before.
(v) From (iv), De ≤ 2 for all sufficiently large t, so Re = Le /De ≥ Lmin /2 for the same t’s.
P
P
(vi) ĊS = e∈δ(S) Ḋe = e∈δ(S) (|Qe | − De ) ≥ 1 − CS , with equality if S = {s0 }.

(vii) Follows by noting that the inequality in (vi) becomes tight for the cut {s0 }, due to (ii).
(viii) From (vi), eventually CS ≥ 1/2 for all S ∈ C, so there is an edge of diameter at least
1/2m in every cut. Thus, there is a s0 -s1 path in which every edge has diameter at least
1/2m.
(ix) Consider a source-sink path in which every edge has diameter at least 1/2m. By (4) the
total potential drop ps0 − ps1 is at most 2nmLmax .
(x) Ḋe /De = (|Qe | − De )/De = |∆e |/Le − 1, and the bound follows from (ix).

6

Convergence

We will prove convergence for general graphs. Throughout this section, we will assume that
t is large enough for all the claims of Lemma 4 requiring a sufficiently large t to hold.

6.1

Properties of Equilibrium Points.

Recall that D ∈ RE
+ is an equilibrium point, when Ḋe = 0 for all e ∈ E, which by (1) is
equivalent to De = |Qe | for all e ∈ E.
Lemma 5 At an equilibrium point, minS∈C CS = C{s0 } = 1.
Proof:
1 ≤ min
S∈C

X

e∈δ(S)

|Qe | = min CS ≤ C{s0 } =
S∈C

X

|Qe | = 1.

e∈δ({s0 })

Lemma 6 The equilibria are precisely the flows of value 1, in which all source-sink paths
have the same length. If no two source-sink paths have the same length, the equilibria are
precisely the simple source-sink paths.
Proof: Let Q be a flow of value 1, in which all source-sink paths have the same length.
We orient the edges such that Qe ≥ 0 for all e and show that D = Q is an equilibrium
point. Let E1 be the set of edges carrying positive flow, and let V1 be the set of vertices lying
on a source-sink path consisting of edges in E1 . For v ∈ V1 , set its potential to the length
of the paths from v to s1 in (V1 , E1 ); observe that all such paths have the same length by
assumption. Let Q′ be the electrical flow induced by the potentials and edge diameters. For
any edge e = (u, v) ∈ E1 , we have Q′e = De ∆e /Le = De = Qe . Thus, Q′ = Q. For any edge
e 6∈ E1 , we have Qe = 0 = De . We conclude that D is an equilibrium point.
11

Let D be an equilibrium point and let Qe be the corresponding current along edge e,
where we orient the edges so that Qe ≥ 0 for all e ∈ E. Whenever De > 0, we have
∆e = Qe Le /De = Le because of the equilibrium condition. Since all directed s0 -s1 paths
span the same potential difference, all directed paths from s0 to s1 in {e ∈ E : De > 0} have
the same length. Moreover, by Lemma 5, minS CS = 1. Thus, D is a flow of value 1.
Let E ∗ be the set of flows of value one in the network of shortest source-sink paths. If the
shortest source-sink path is unique, E ∗ is a singleton, namely the flow of value one along the
shortest source-sink path.

6.2

The Convergence Process

Lemma 7 Let W = (C{s0 } − 1)2 . Then, Ẇ = −2W ≤ 0, with equality iff C{s0 } = 1.
Proof: Let C0 = C{s0 } for short. Then, since
Ẇ = 2(C0 − 1)

X

P

e∈δ({s0 }) |Qe |

= 1,

(|Qe | − De ) = 2(C0 − 1)(1 − C0 ) = −2(C0 − 1)2 ≤ 0.

e∈δ({s0 })

The following functions play a crucial role. Let C = minS∈C CS , and
VS =

1 X
Le De for each S ∈ C,
CS
e∈E

V = max VS + W, and
S∈C

h=−

1 X
1 X
Re |Qe |De + 2
Re De2 .
C
C
e∈E

e∈E

Lemma 8 Let S be a minimum capacity cut at time t. Then, V̇S (t) ≤ −h(t).
Proof: Let X be the characteristic vector of δ(S), that is, Xe = 1 if e ∈ δ(S) and 0
otherwise. Observe that CS = C since S is a minimum capacity cut. We have
V̇S =

X ∂VS
Ḋe
∂De
e

X 1
=
C2
e

Le C −

X

Le′ De′ Xe

e′

1
1 X
Le |Qe | − 2
=
C e
C

X

!

Le′ De′

e′

(|Qe | − De )
!

X
e

!

Xe |Qe | +

!
!
X
1 X
1 X
−
Le De + 2
Xe De
Le′ De′
C e
C
e
e′
1 X
1 X
1 X
1 X
Re |Qe |De − 2
Le De +
Le De
≤
Re De2 −
C e
C e
C e
C e

= −h.

12

The only inequality follows from Le = Re De and
one unit current must cross S.

P

e Xe |Qe |

≥ 1, which holds because at least

Lemma 9 Let f (t) = maxS∈C fS (t), where each fS is continuous and differentiable. If f˙(t)
exists, then there is S ∈ C such that f (t) = fS (t) and f˙(t) = f˙S (t).
Proof: Since C is finite, there is at least one S ∈ C such that for each fixed δ > 0, f (t + ε) =
fS (t+ε) for infinitely many ε with |ε| ≤ δ. By continuity of f and fS , this implies f (t) = fS (t).
Moreover, since
maxS ′ fS ′ (t + ε) − maxS ′ fS ′ (t)
lim
ε→0
ε
exists and is equal to f˙(t), any sequence ε1 , ε2 , . . . converging to zero has the property that
maxS ′ fS ′ (t + εi ) − maxS ′ fS ′ (t)
→ f˙(t)
εi

for i → ∞.

Taking (εi )∞
i=1 to be a sequence converging to zero such that f (t + εi ) = fS (t + εi ) for all i,
we obtain
fS (t + εi ) − fS (t)
f˙(t) = lim
= f˙S (t).
i→∞
εi
Lemma 10 V̇ exists almost everywhere. If V̇ (t) exists, then V̇ (t) ≤ −h(t) − 2W (t) ≤ 0, and
V̇ (t) = 0 if and only if Ḋe (t) = 0 for all e.
Proof: V is Lipschitz-continuous since it is the maximum of a finite set of continuously
differentiable functions. Since V is Lipschitz-continuous, the set of t’s where V̇ (t) does not
exist has zero Lebesgue measure (see for example [CLSW98, Ch. 3], [MN92, Ch. 3]). When
V̇ (t) exists, we have V̇ (t) = Ẇ (t) + V̇S (t) for some S of minimum capacity (Lemma 9). Then,
V̇ (t) ≤ −h(t) − 2W (t) by Lemmas 7 and 8.
The fact that W ≥ 0 is clear. We now show that h ≥ 0. To this end, let F represent a
maximum s0 -s1 flow in an auxiliary network, having the same structure as G, and where the
capacity on edge e is set equal to De . In other words, F is an s0 -s1 flow satisfying |Fe | ≤ De
for all e ∈ E and having maximum value. By the max-flow min-cut theorem, this maximum
value is equal to C = minS∈C CS . But then,
1 X
1 X
Re De2
Re |Qe |De − 2
−h =
C e
C e
!1/2
!1/2
X
1 X
1 X
2
2
≤
− 2
Re De
Re Q e
Re De2
C
C
e
e
e
!1/2
!1/2
X
1 X Fe2
1 X
Re 2
Re De2
≤
− 2
Re De2
C
C
C
e
e
e
!1/2
!1/2
X
1 X
1 X
≤ 2
Re De2
Re De2
Re De2
− 2
C
C
e
e
e
= 0,

where we used the following inequalities:
13

- the Cauchy-Schwarz inequality

P

1/2
1/2
e (Re |Qe |)(Re De )

P
P
≤ ( e Re Q2e )1/2 ( e Re De2 )1/2 ;

- Thomson’s Principle (6) applied to the unit-value flows Q and F/C; Q is a minimum
energy flow of unit value, while F/C is a feasible flow of unit value;
- the fact that |Fe | ≤ De for all e ∈ E.
Finally, one can have h = 0 if and only if all the above inequalitiesPare equalities, which
implies that |Qe | = |Fe |/C = De /C for all e. And, W = 0 iff
e∈δ({s0 }) De = 1 =
P
e∈δ({s0 }) |Qe |. So, h = W = 0 iff |Qe | = De for all e.
The next lemma is a necessary technicality.

Lemma 11 The function t 7→ h(t) is Lipschitz-continuous.
Proof: Since Ḋe is continuous and bounded (by (1)), De is Lipschitz-continuous. Thus, it
is enough to show that Qe is Lipschitz-continuous for all e.
First, we claim that De (t+ε) ≤ (1+2Kε)De for all ε ≤ 1/4K, where K = 2nmLmax /Lmin .
For if not, take
ε = inf{δ ≤ 1/4K : D(t + δ) > (1 + 2Kδ)D(t)},
then ε > 0 (since Ḋe (t) ≤ KDe (t) by Lemma 4) and, by continuity, De (t+ε) ≥ (1+2Kε)De (t).
There must be t′ ∈ [t, t + ε] such that Ḋe (t′ ) = 2KDe (t). On the other hand,
Ḋe (t′ ) ≤ KDe (t′ ) ≤ K(1 + 2Kε)De (t)
≤ K(1 + 2K/4K)De (t) < 2KDe (t),
which is a contradiction. Thus, De (t + ε) ≤ (1 + 2Kε)De for all ε ≤ 1/4K. Similarly,
De (t + ε) ≥ (1 − 2Kε)De .
Q
Consider now a spanning tree T of G. Let γT = e∈T De /Le . Then γT (t + ε) ≤ (1 +
2Kε)n γT (t) ≤ (1+4nKε)γT (t) for sufficiently small ε. Similarly, γT (t+ε) ≥ (1−4nKε)γT (t).
By Kirchhoff’s Theorem,
P
P
T ∈Sp(v,u) γT
T ∈Sp(u,v) γT −
P
,
Quv =
T ∈Sp γT

and plugging the bounds for γT (t + ε)/γT (t) shows that Qe (t + ε) = Qe (t)(1 + O(ε)), where
the constant implicit in the O(·) notation does not depend on t. Since |Qe | ≤ 1, we obtain
that |Qe (t + ε) − Qe (t)| ≤ O(1) · ε, that is, Qe is Lipschitz-continuous, and this in turn implies
the Lipschitz-continuity of h.

Lemma 12 |De − |Qe || converges to zero for all e ∈ E.
Proof: Consider again the function h. We claim h → 0 as t → ∞. If not, there is ε > 0 and
an infinite unbounded sequence t1 , t2 , . . . such that h(ti ) ≥ ε for all i. Since h is Lipschitzcontinuous (Lemma 11), there is δ such that h(ti + δ′ ) ≥ h(ti ) − ε/2 ≥ ε/2 for all δ′ ∈ [0, δ]
and all i. So by Lemma 10, V̇ (t) ≤ −h(t) ≤ −ε/2 for every t in [ti , ti + δ] (except possibly a
zero measure set), meaning that V decreases by at least εδ/2 infinitely many times. But this
is impossible since V is positive and non-increasing.
14

Thus, for any ε > 0, there is t0 such that h(t) ≤ ε for all t ≥ t0 . Then, recalling that
Re ≥ Lmin /2 for all sufficiently large t (Lemma 4.v), we find
2 X 
2
X Lmin  De
De
Re
− |Qe | ≤
− |Qe |
2
C
C
e
e
X
2 X
1 X
Re De2 +
Re Q2e −
Re |Qe |De
= 2
C e
C e
e
2 X
2 X
≤ 2
Re De2 −
Re |Qe |De
C e
C e

= 2h ≤ 2ε,
P
P
2
2
2
where we used once more the inequality
e Re De /C , which was proved in
e Re Qe ≤
Lemma 10. This implies that for each e, De /C − |Qe | → 0 as t → ∞. Summing across
e ∈ δ({s0 }) and using Lemma 4.ii, we obtain C{s0 } /C − 1 → 0 as t → ∞. From Lemma 4,
C{s0 } → 1 as t → ∞, so C → 1 as well.
To conclude, we show that De /C − |Qe | → 0 and C → 1 together imply De − |Qe | → 0.
Let ε > 0 be arbitrary. For all sufficiently large t, |De /C − |Qe || ≤ ε, |1 − C| ≤ ε, De ≤ 2,
and C ≥ 1/2. Thus,
|De − |Qe || ≤ |De − De /C| + |De /C − |Qe || ≤ De

|C − 1|
+ |De /C − |Qe || ≤ 5ε.
C

Lemma 13 Let ∆ = ps0 −ps1 be the potential difference between source and sink. ∆ converges
to the length L∗ of a shortest source-sink path.
Proof: Let L be the set of lengths of simple source-sink paths. We first show that ∆
converges to a point in L and then show convergence to L∗ .
Orient edges according to the direction of the flow. By Lemma 4.viii, there is a directed
source-sink path P of edges of diameter at least 1/2m. Let ε > 0 be arbitrary. We will show
|∆ − LP | ≤ ε. For this, it suffices to show |∆e − Le | ≤ ε/n for any edge e of P , where ∆e
is the potential drop on e. By Ohm’s law, the potential drop on e is ∆e = (Qe /De )Le , and
hence, |∆e − Le | = |Qe /De − 1|Le = |(Qe − De )/De |Le ≤ 2mLmax |Qe − De |. The claim follows
since |Qe − De | converges to zero.
The set L is finite. Let ε be positive and smaller than half the minimal distance between
two elements in L. By the preceeding paragraph, there is for all sufficiently large t a path Pt
such that |∆ − LPt | ≤ ε. Since ∆ is a continuous function of time, LPt must become constant.
We have now shown that ∆ converges to an element in L.
∗
′
We will next show that ∆ converges to
P L . Assume otherwise, and let P be a shortest
undirected source-sink path. Let WP ′ = e∈P ′ Le ln De . This function was already used by
Miyaji and Ohnishi [MO08]. We have
ẆP ′ =

X
X
X Le
(|Qe | − De ) =
|∆e | −
Le ≥ p s 0 − p s 1 − LP ′ = ∆ − L∗ .
D
e
′
′
′

e∈P

e∈P

e∈P

Let δ > 0 be such that there is no source-sink path with length in the open interval (L∗ , L∗ +
2δ). Then, ∆ − L∗ ≥ δ for all sufficently large t, and hence, ẆP ′ ≥ δ for all sufficiently large
15

t. Thus, WP ′ goes to +∞. However, WP ′ ≤ nLmax for all sufficiently large t since De ≤ 2 for
all e and t large enough. This is a contradiction. Thus, ∆ converges to L∗ .

Lemma 14 Let e be any edge that does not lie on a shortest source-sink path. Then, De and
Qe converge to zero.
Proof: Since |De − |Qe || converges to zero, it suffices to prove that Qe converges to zero.
Assume otherwise. Then, there is a δ > 0 such that |Qe | ≥ δ for arbitrarily large t.
Consider any such t and orient the edges according to the direction of the flow at time
t. Let e = (u, v). Because of flow conservation, there must be an edge into u and an edge
out of v carrying flow at least Qe /n. Continuing in this way, we obtain a source-sink path
P in which every edge carries flow at least Qe /nn ≥ δ/nn ; P depends on time and LP > L∗
always. We will show |∆ − LP | ≤ (LP − L∗ )/4 for sufficiently large t, a contradiction to the
fact that ∆ converges to L∗ . For this, it suffices to show |∆g − Lg | ≤ (LP − L∗ )/(4n) for
any edge g of P , where ∆g is the potential drop on g. By Ohm’s law, the potential drop
on g is ∆g = (Qg /Dg )Lg , and hence, |∆g − Lg | = |Qg /Dg − 1|Lg = |(Qg − Dg )/Dg |Lg ≤
Lmax |Qg − Dg |/Dg . For large enough t, |Qg − Dg | ≤ min(δ/(2nn ), δ(LP − L∗ )/(8nn+1 Lmax )).
Then, Dg ≥ Qg − |Qg − Dg | ≥ δ/(2nn ), and hence, Lmax |Qg − Dg |/Dg ≤ (LP − L∗ )/(4n).
Theorem 2 The dynamics are attracted by E ∗ . If the shortest source-sink path is unique,
the dynamics converge against a flow of value 1 on the shortest source sink path.
Proof: Q is a source-sink flow of value one at all times. We show first that Q is attracted
to E ∗ . Orient the edges in the direction of the flow. We can decompose QP
into flowpaths. For
an oriented path P , let 1P be the unit flow along P . We can write Q = P xp 1P , where xP
is the flow along the path P . This decomposition is not unique. We group the flowpath into
two sets, the paths running inside G0 and the paths using an edge outside G0 , i.e.,
X
Q = Q0 + Q1 , where Q0 =
xP 1P .
P is a path in G0

Q0 is a flow in G0 , and each flowpath in Q1 is a non-shortest source-sink path.3 We show
that the value of Q0 converges to one.
Assume otherwise. Then, there is a δ > 0 such that the value of Q1 is at least δ for
arbitrarily large times t. At any such time, there is an edge e 6∈ E0 carrying flow at least δ/m;
this holds since source-sink cuts contain at most m edges. Since there are only finitely many
edges, there must be an edge e 6∈ E0 for which Qe does not converge to zero, a contradiction
to Lemma 14.
We have now shown that the distance between Q and E ∗ converges to zero. By Lemma 12,
|De − |Qe || converges to zero for all e, and hence, the distance between Q and D converges
to zero. Thus, D is attracted by E ∗ .
3

The decomposition into Q0 and Q1 can be constructed as follows: Initialize Q0 to Q and Q1 to the empty
flow. Consider any edge e 6∈ E0 carrying positive flow in Q0 , say ε. Let P be an oriented source-sink path
carrying ε units of flow and using e. Add ε1P to Q1 and subtract it from Q0 . Continue until Q0 is a flow in
G0 .

16

Finally, if the shortest source-sink path is unique, E ∗ is a singleton, and hence, D converges
to the flow of value one along the shortest source-sink path.
Lemma 15 If the shortest source-sink path is unique, pv converges to dist(v, s1 ) for each
node v on the shortest source-sink path, where dist(v, s1 ) is the shortest path distance from v
to s1 .
Proof: Let P0 be the shortest source-sink path. For any e ∈ P , De converges to one and
|De − Qe | converges to zero. Thus, ∆e converges to Le .

6.3

More on the Lyapunov Function V

P
In this section, we study V = e Le P
De /C + (C{s0 } − 1)2 as a function of D. Recall that
C = C(D) = minS∈C CS , where CS = e∈δ(S) De .
Lemma 16 Let D 0 and D 1 be two equilibrium points. Define
D λ = (1 − λ)D 0 + λD 1 ,

λ ∈ [0, 1].

If V (D 0 ) < V (D 1 ), then V (D λ ) is a linear, increasing function of λ.
Proof: By Lemma 5, C(D 0 ) = C(D 1 ) = 1, and C{s0 } (D 0 ) = C{s0 } (D 1 ) = 1. Since CS (D)
is linear in D for any fixed cut S, one has CS (D 0 ) ≥ 1 and CS (D 1 ) ≥ 1, so CS (D λ ) ≥ 1
for all S. Thus, C(D λ ) ≥ 1. On the other hand, C{s0 } (D λ ) = 1. Thus, C(D λ ) = 1, and
P
V (D λ ) = e Le Deλ , that is, V (D λ ) is a linear function of D λ .
Lemma 17 The problem of minimizing V (D) for D ∈ RE
+ is equivalent to the shortest path
problem.
Proof: By introducing an additional variable C = minS CS > 0, the problem of minimizing
V (D) is equivalently formulated as
2

X
X
1
min
De − 1
Le De + 
C e
e∈δ({s0 })

s.t. CS ≥ C

∀S ∈ C

C>0
D ≥ 0.
Substituting xe = De /C, we obtain
min

X
e

s.t.



Le xe + C 1/2 

X

X

e∈δ({s0 })

xe ≥ 1

∀S ∈ C

e∈δ(S)

x ≥ 0, C > 0,
17

xe −

2

1
C

which is easily seen to be equivalent to the (fractional) shortest path problem.
Lemma 17 was the basis for the generalization of our results to the transportation problem
(Section 9). We first generalized the above Lemma to Lemma 33 and then used the Lyapunov
function suggested by the generalization.

7

Rate of Convergence for Stable Flow Directions

The direction of the flow across an edge depends on the initial conditions and time. We
do not know whether flow directions can change infinitely often or whether they become
ultimately fixed. In this section, we assume that flow directions stabilize and explore the
consequences of this assumption. We will be able to make quite precise statements about the
convergence of the system. We assume uniqueness of the shortest source-sink path and add
more non-degeneracy assumptions as we go along.
An edge e = {u, v} becomes horizontal if limt→∞ |pu − pv | = 0, and it becomes directed
from u to v (directed from v to u) if pu > pv for all large t (pv > pu for all large t). An edge
stabilizes if it either becomes horizontal or directed, and a network stabilizes if all its edges
stabilize. If a network stabilizes, we partition its edges into a set Eh of horizonal edges and
→
−
→
−
a set E of directed edges. If {u, v} becomes directed from u to v, then (u, v) ∈ E .
We already know that the diameters of the edges on the shortest source-sink path (we
assume uniqueness in this section) converge to one. The diameters of the edges outside G0
converge to zero. The potential of a vertex v ∈ G0 converges to dist(v, s1 ). For stabilizing
networks, we can prove a lot more. In particular, we can predict the decay rates of edges, the
limit potentials of the vertices, and for each edge the direction in which the flow will stabilize.
Definition 1 (Decay Rate) Let r ≤ 0.
A quantity D(t) decays with rate at least r if for every ε > 0 there is a constant A > 0
such that for all t
D(t) ≤ Ae(r+ε)t ,

or equivalently,

ln D(t) ≤ (ln A) + (r + ε)t.

A quantity D(t) decays with rate at most r if for every ε > 0 there is a constant a > 0
such that for all t
D(t) ≥ ae(r−ε)t ,

or equivalently,

ln D(t) ≥ (ln a) + (r − ε)t.

A quantity D(t) decays with rate r if it decays with rate at least and at most r.
We first establish a simple Lemma that, for any edge, connects the decay rate of the flow
across the edge and the diameter of the edge.
Lemma 18 Let −1 ≤ a < 0 and let e, g ∈ E. If Qe decays with rate at least a, then so does
De . De decays with rate at most −1. If ||Qe |− |Qg || decays with rate at least a, then |De − Dg |
decays with rate at least a.
Proof: Assume first that Qe decays with rate at least a, where −1 ≤ a < 0. Then, for any
ε > 0, there is an A > 0 such that Qe ≤ Ae(a+ε)t for all t. Consider f with f˙ = Ae(a+ε)t − f .

18

This has solution f = f0 e−t + αe(a+ε)t , where α = A/(1 + a + ε) and f0 is determined by the
value of f at zero, namely, f (0) = f0 + α. Consider De − f . Then,
d
(De − f ) = |Qe | − De − (Ae(a+ε)t − f ) ≤ −(De − f ).
dt
Thus, De − f ≤ C ′ e−t for some constant C ′ by Gronwall’s Lemma, and hence,
De ≤ (f0 + C ′ )e−t + αe(a+ε)t ≤ C ′′ e(a+ε)t
for some constant C ′′ . Thus, De decays with rate at least a.
Ḋe = |Qe | − De ≥ −De . Thus, De decays with rate at most −1 by Gronwall’s Lemma.
Finally, assume that ||Qe | − |Qf || decays with rate at least a. Then,
d
(De − Dg ) = |Qe | − |Qf | − (De − Dg ) ≤ ||Qe | − |Qf || − (De − Dg ),
dt
and therefore, De − Dg decays with rate at least −a. The same argument applies to Dg − De .

P
For a path P , let W (P ) := e∈P Le ln De be its weighted sum of log diameters, and let
∆(P ) = pa − pb be the potential difference between its endpoints. The function W (P ) was
introduced by Miyaji and Ohnishi [MO07, MO08].
Lemma P
19 Let P be an arbitrary path, let ∆(P ) be the potential drop along P , and let
W (P ) = e∈P Le ln De . Then,
Ẇ (P ) = ∆(P ) − L(P ) + 2

X

|∆(e)|.

e∈P : ∆(e)<0

If ∆(P ) ≤ ∆ and ∆(e) ≥ −δ for some δ ≥ 0, all e ∈ P and for all sufficiently large t, then
W (P )(t) ≤ C + (∆ − L(P ) + 2nδ)t
for some constant C and all t. If ∆(P ) ≥ ∆ for all sufficiently large t, then
W (P )(t) ≥ C + (∆ − L(P ))t
for some constant C and all t.
Proof: The first claim follows immediately from the dynamics of the system.
X
X
Ẇ (P ) =
|∆(e)| − L(P ) = ∆(P ) − L(P ) + 2
|∆(e)|.
e∈P

e∈P : ∆(e)<0

Let t0 be such that ∆(P ) ≤ ∆ and ∆(e) ≥ −δ for all t ≥ t0 . We integrate the equality
from t0 to t and obtain
Z t
Ẇ (P )dt ≤ (∆ − L(P ) + 2nδ)(t − t0 ).
W (P )(t) − W (P )(t0 ) =
t0

19

This establishes the claim for t ≥ t0 . Choosing C sufficiently large extends the claim to all t.
Let t0 be such that ∆(P ) ≥ ∆. We integrate the equality from t0 to t and obtain
W (P )(t) − W (P )(t0 ) =

Z

t

Ẇ (P )dt ≥ (∆ − L(P ))(t − t0 ).

t0

This establishes the claim for t ≥ t0 . Choosing C sufficiently large extends the claim to all t.

Edges that do not lie on a source-sink path never carry any flow, and hence, their diameter
evolves as De (0) exp(−t). From now on, we may therefore assume that every edge of G lies
on a source-sink path.
Lemma 20 For e ∈ Eh , De decays with rate −1, and |Qe | decays with rate at least −1.
Proof: We certainly have De ≤ 2 for all large t. Let e = {u, v}, and let ε > 0 be arbitrary.
Then, |pu − pv | ≤ εLe for all large t, and hence, |Qe | = (De /Le )|pu − pv | ≤ εDe for all large
t. Thus, Ḋe ≤ (ε − 1)De for all large t, and hence, (d/dt) ln De ≤ −1 + ε. Thus, De decays
with rate at least −1. Since Ḋe ≥ −De , De decays with rate at most −1.
|Qe | = (De /Le )|pu − pv | ≤ ADe for some constant A. Thus, |Qe | decays with rate at least
−1.
We define a decomposition of G into paths P0 to Pk , an orientation of these paths, a slope
f (Pi ) for each Pi , a vertex labelling p∗ , and an edge labelling r. P0 is a4 shortest s0 -s1 path
in G, f (P0 ) = 1, re = f (P0 ) − 1 for all e ∈ P0 , and p∗v = dist(v, s1 ) for all v ∈ P0 , where
dist(v, s1 ) is the shortest path distance from v to s1 . For 1 ≤ i ≤ k, we have5
Pi = argmax f (P ),
P ∈P

where P is the set of all paths P in G with the following properties:
- the startpoint a and the endpoint b of P lie on P0 ∪ . . . ∪ Pi−1 , p∗a ≥ p∗b , and f (P ) =
(p∗a − p∗b )/L(P );
- no interior vertex of P lies on P0 ∪ . . . ∪ Pi−1 ; and
- no edge of P belongs to P0 ∪ . . . ∪ Pi−1 .
If p∗a > p∗b , we direct Pi from a to b. If p∗a = p∗b , we leave the edges in Pi undirected. We set
re = f (Pi ) − 1 for all edges of Pi , and p∗v = p∗b + f (Pi )distPi (v, b) for every interior vertex v of
Pi . Here, distPi (v, b) is the distance from v to b along path Pi . Figure 3 illustrates the path
decomposition.
Lemma 21 There is an i0 ≤ k such that
f (P0 ) > f (P1 ) > . . . > f (Pi0 ) > 0 = f (Pi0 +1 ) = . . . = f (Pk ).
4
5

We assume that P0 is unique.
We assume that Pi is unique except if f (Pi ) = 0.

20

s0

s1

e1

e2

e4

e3
u

v
e5

e6
w

Figure 3: All edges are assumed to have length 1; P0 = (e1 ), P1 = (e2 , e3 , e4 ), P2 = (e5 , e6 ),
p∗s0 = 1, p∗s1 = 0, p∗v = 1/3, p∗u = 2/3, p∗w = 1/2, f (P1 ) = 1/3, and f (P2 ) = 1/6.
The path (e2 , e5 , e6 , e4 ) has f -value 1/4.
Proof: It suffices to show: if there is an i such that f (Pi+1 ) ≥ f (Pi ), then f (Pi ) = f (Pi+1 ) =
0. If no endpoint of Pi+1 is an internal vertex of Pi , then f (Pi+1 ) = f (Pi ); otherwise Pi+1
would have been chosen instead of Pi . By assumption, equality is only possible if the f -values
are zero. So we may assume that at least one endpoint of Pi+1 is an internal node of Pi ; call
it c and assume w.l.o.g. that it is the startpoint of Pi+1 . Split Pi at c into Pi1 and Pi2 , and
let d be the other endpoint of Pi+1 ; d may lay on Pi .
Assume first that d does not lie on Pi and consider the path Pi1 Pi+1 . The f -value of this
path is
p∗a − p∗c + p∗c − p∗d
p∗a − p∗d
=
.
L(Pi1 ) + L(Pi+1 )
L(Pi1 ) + L(Pi+1 )
Next, observe that (p∗a − p∗c )/L(Pi1 ) = f (Pi ) since p∗c is defined by linear interpolation and
(p∗c − p∗d )/L(Pi+1 ) = f (Pi+1 ) ≥ f (Pi ). In case of inequality, Pi1 Pi+1 is chosen instead of Pi .
In case of equality, there are two paths with the same f -value. By assumption, this is only
possible if the f -values are zero.
Assume next that d also lies on Pi . We then split Pi into three paths Pi1 , Pi2 , and Pi3 and
consider the path Pi1 Pi+1 Pi3 . We then argue as in the preceding paragraph.
→
−
→
−
Theorem 3 If a network stabilizes, then E = ∪i≤i0 E(Pi ), the orientation of any edge e ∈ E
agrees with the orientation induced by the path decomposition, and Eh = ∪i>i0 E(Pi ). The
potential of each node v converges to p∗v . The diameter of each edge e ∈ E \ P0 decays with
rate re .
Proof: We use induction on i to prove:
- for every vertex v ∈ P0 ∪ . . . ∪ Pi , the node potential pv converges to p∗v ;
- for every edge e ∈ P0 ∪ . . . ∪ Pmin(i,i0 ) , the flow stabilizes in the direction of the path Pj
containing e;
- for every edge e ∈ P1 ∪ . . . ∪ Pi , the diameter converges to zero with rate re , and the
flow converges to zero with rate at least6 re . If e ∈ Pi and i ≤ i0 , the flow converges to
zero with rate re .
6
If for an edge e = {u, v}, pu − pv = 0 always, then Qe = 0 always. Thus, for horizontal edges, Qe may
converge to zero faster than with rate −1.

21

Lemma 15 establishes the base of the induction, the case i = 0. Assume now that the
induction hypothesis holds for i − 1; we establish it for i. Let P≤i−1 = P0 ∪ . . . ∪ Pi−1 .
For e ∈ E \ P≤i−1 , let
 ∗

pa − p∗b
′
fe = max
; P ∈ Pe ,
L(P ′ )

where Pe is the set of paths P ′ in G \ P≤i−1 from some a ∈ P≤i−1 to some b ∈ P≤i−1
with p∗a ≥ p∗b and containing e. Then, maxe6∈P≤i−1 fe = f (Pi ). For i ≤ i0 , we have further
f (Pi ) > maxe6∈P≤i fe ≥ f (Pi+1 ). In general, the last inequality may be strict; see Figure 3.
Lemma 22 For e ∈ E \ P≤i−1 , |Qe | and De decay with rate at least fe − 1.
Proof: According to Lemma 18, it suffices to prove the decay of |Qe |. Let e ∈ E \ P≤i−1
and let ε > 0 be arbitrary. We need to show
ln |Qe (t)| ≤ C + (fe + ε − 1)t
for some constant C and all sufficiently large t.
If Qe (t) = 0, the inequality holds for any value of C. So assume Qe (t) 6= 0 and also
assume that the flow across e = {u, v} is in the direction from u to v. We construct a path
R(t) containing uv. For every vertex, except for source and sink, we have flow conservation.
Hence there is an edge (v, w) carrying a flow of at least Qe /n in the direction from v to w.
Similarly, there is an edge (x, v) carrying a flow of at least Qe /n in the direction from x to v.
Continuing in this way, we reach vertices in P≤i−1 . Any edge on the path R(t) carries a flow
of at least Qe /nn .
Since potential differences are bounded by B := 2nmLmax (Lemma 4.ix), any edge e′ on
R(t) must have a diameter of at least Qe Le /(nn B) ≥ (Lmin /(nn B))Qe . Let c = Lmin /(nn B).
Then,
X
W (R(t)) =
Le′ ln De′ ≥ L(R(t))(ln c + ln |Qe (t)|).
e′ ∈R(t)

The path R(t) depends on time. Let a(t) and b(t) be the endpoints of R(t). Since e does not
belong to P≤i−1 ,
p∗a(t) − p∗b(t)
≤ fe .
f (R(t)) =
L(R(t))
For large enough t, we have ∆(R(t)) ≤ ∆∗ (R(t)) + εL(R)/2. Every edge e ∈ R(t) either
→
−
belongs to E or to Eh due to the assumption that the network stabilizes. In the former case,
→
−
R must use e in the direction fixed in E , in the latter case, the potential difference across
e converges to zero. We now invoke Lemma 19 with δ = εL(R)/(4n). It guarantees the
existence of a constant C1 such that
W (R(t))(t) ≤ C1 + (∆∗ (R(t)) + εL(R)/2 − L(R) + εL(R)/2)t
for all t. The constant C1 depends on the path R(t). Since there are only finitely many
different paths R(t), we may use the same constant C1 for all paths R(t).
Combining the estimates, we obtain, for all sufficiently large t,
L(R(t))(ln c + ln |Qe (t)|) ≤ C1 + (∆∗ (R(t)) + εL(R(t)) − L(R(t)))t,
22

and hence,
ln |Qe (t)| ≤ C1 /L(R(t)) − ln c + (fe + ε − 1)t.

Corollary 4 For e ∈ E \ P≤i−1 , |Qe | and De decay with rate at least f (Pi ) − 1. If i ≤ i0 ,
then for any e ∈ E \ P≤i , |Qe | and De decay with rate at least f (Pi ) − δ − 1 for some δ > 0.
Proof: If i ≤ i0 , and hence, f (Pi ) > 0, fe < f (Pi ) for any edge e ∈ E \ P≤i . The claim
follows.

Lemma 23 Let e ∈ Pi . Then, De decays with rate f (Pi ) − 1. If i ≤ i0 , then |Qe | decays with
rate f (Pi ) − 1.
Proof: We distinguish the cases f (Pi ) = 0 and f (Pi ) > 0. If f (Pi ) = 0, the diameter of all
edges e ∈ Pi decays with rate at least −1 (Lemma 19). No diameter decays with a rate faster
than −1.
We turn to the case f := f (Pi ) > 0. The flows across the edges in E \ P<i decay with rate
at least f − 1, and the flows across the edges edges in E \ P≤i decay faster, say with rate at
least f − δ − 1 for some positive δ (Corollary 4). We first show
W (Pi ) ≤ C + L(Pi ) · max(ln De , (f − δ − 1)t)

(8)

for sufficiently large t and some constant C. If Pi consists of a single edge e, W (Pi ) =
Le ln De (t) and (8) holds. Assume next that Pi = e1 . . . ek with k > 1. Consider any interior
node u of the path. The flow into u is equal to the flow out of u, and u has two incident
edges7 in Pi . The flow on the other edges incident to u decays with rate at least f − δ − 1.
Thus for any two consecutive edges on Pi , | |Qej | − |Qej+1 | | decays with rate at least f − δ − 1.
By Lemma 18, this implies that |Dej − Dej+1 | decays with rate at least f − δ − 1. Thus, we
have Dej = De + gej , where |gej | ≤ C1 e(f −δ−1)t for some constant C1 and all j. Plugging into
the definition of W (Pi ) yields
X

W (Pi ) ≤
Lej ln 2 max(De , gej )
ej ∈Pi

≤ L(Pi ) ln 2 + L(Pi ) max(ln De , ln C1 e(f −δ−1)t ),

and we have established (8).
Let t0 be large enough such that |∆(Pi ) − ∆∗ (Pi )| ≤ δL(Pi )/2 for all t ≥ t0 . Then, by
Lemma 19,
W (Pi ) ≥ A + L(Pi )(f − δ/2 − 1)t
(9)
for some constant A and all t.
Combining (8) and (9) yields
A + L(Pi )(f − δ/2 − 1)t ≤ C + L(Pi ) · max(ln De , (f − δ − 1)t).
7

Here, we need uniqueness of Pi . Otherwise we would have a network of paths with the same slope.

23

Thus, for every t we have either
A + L(Pi )(f − δ/2 − 1)t ≤ C + L(Pi ) · ln De
or
A + L(Pi )(f − δ/2 − 1)t ≤ C + L(Pi ) · (f − δ − 1)t.
The latter inequality does not hold for any sufficiently large t. Thus, the former inequality
holds for all sufficiently large t, and hence, De decays with rate at most f (Pi ) − 1. By
Lemma 18, |Qe | cannot decay at a faster rate if f (Pi ) > 0.
Lemma 24 For v ∈ Pi , the potentials converge to p∗v . For e ∈ Pi and i ≤ i0 , the flow
direction stabilizes in the direction of Pi .
Proof: Assume i ≤ i0 first. Let Pi = e1 . . . ek . The flows and the diameters of the edges in
Pi decay with rate f (Pi ) − 1 (Lemma 23). The flows and diameters of the edges incident to
the interior vertices of Pi and not on Pi decay faster, say with rate at least f (Pi )−δ −1, where
δ > 0. For large t and any interior vertex of Pi , one edge of Pi must, therefore, carry flow into
the vertex, and the other edge incident to the vertex must carry it out of the vertex. Thus,
the edges in Pi must either all be directed in the direction of Pi or in the opposite direction.
As current flows from higher to lower potential, they must be directed in the direction of Pi .
Because the flow and the diameters of the edges not on Pi and incident to interior vertices
decay faster, we have for any ε > 0 and sufficiently large t
Qej = Qe1 (1 + εj ) and

Dej = De1 (1 + ε′j ),

where |εj |, |ε′j | ≤ ε. The potential drop ∆ej on edge ej is equal to
∆ej =

Qe1 (1 + ε′j )
Qe j Le j
Le ,
=
Dej
De1 (1 + εj ) j

and hence, the potential drop along the path is
pa − pb =

X

∆ej =

j

Qe 1
L(Pi )(1 + ε′′ ),
De1

where ε′′ goes to zero with ε. The potential drop along the path converges to p∗a − p∗b .
Thus, Qe1 /De1 converges to f (Pi ), and therefore, the potential of any interior vertex v of Pi
converges to p∗v .
We turn to the case i > i0 . The potentials of the endpoints of Pi converge to the same
value. Thus, the potentials of all interior vertices of Pi converge to the common potential of
the endpoints.
We have now completed the induction step.

24

s0
b

a
e

L

R
c

d
s1

Figure 4: The Wheatstone graph.

8

The Wheatstone Graph

Do edge directions stabilize? We do not know. We know one graph class for which edge
directions are unique, namely series-parallel graphs. The simplest graph which is not seriesparallel is the Wheatstone graph shown in Figure 4. We use the following notation: We have
edges a to e as shown in the figure. For an edge x, Rx = Lx /Dx denotes its resistance and
Cx = Dx /Lx denotes its conductance.8 For edges a, b, c, and d, the direction of the flow is
always downwards. For the edge e, the direction of the flow depends on the conductances.
We have an example where the direction of the flow across e changes twice.
A shortest path from source to sink may have two essentially different shapes. It either
uses e, or it does not. If e lies on a shortest path, Lemma 19 suffices to prove convergence as
observed by [MO08]. If (a, e, d) is a shortest path9 , let P = (a, e) and P ′ = (b). Then,
d
(W (P ) − W (P ′ )) ≥ ∆(P ) − L(P ) − (∆(P ′ ) − L(P ′ )) = L(P ′ ) − L(P ) > 0.
dt
Since W (P ) is bounded, this implies W (P ′ ) → −∞. Thus, Db converges to zero. Similarly,
Dd must converge to zero. More precisely, W (P ′ ) goes to −∞ linearly, and hence, Db and
similarly Dd decay exponentially.
The non-trivial case is that the shortest path does not use e. We may assume w.l.o.g. that
the shortest path uses the edges a and c. The ratio
xa =

Ra
1
1
Cc
=
=
=
Ra + Rc
1 + Rc /Ra
1 + Ca /Cc
Ca + Cc

is the ratio of the resistance of a to the total resistance of the right path; define xb , xc , and
xd analogously. Observe xa + xc = 1 and xb + xd = 1. Let
x∗a =

La
;
La + Lc

define x∗b , x∗c , and x∗d analogously. Without edge e, the potential drop on the edge a is xa
times the potential difference between source and sink. If Da = Dc , which we expect in the
limit, xa = x∗a .
8
9

Observe that we use the letter C with a different meaning than in preceding sections.
For simplicity, we assume uniqueness of the shortest path in this section.

25

Lemma 25 Let S = Ca Cb (Cc + Cd ) + (Ca + Cb )Cc Cd + (Ca + Cb )(Cc + Cd )Ce . Then,


 ∗
Ca Cc
xb
xa
∗
x˙a =
−
(Cb + Cd + Ce )(La + Lc )(Ca + Cc )(xa − xa ) + Ce Cb Lc
SLa Lc (Ca + Cc )2
x∗
xd
 ∗c


x
x
Cb Cd
a
∗
b
(Ca + Cc + Ce )(Lb + Ld )(Cb + Cd )(xb − xb ) + Ce Ca Ld
.
−
x˙b =
SLb Ld (Cb + Cd )2
x∗d
xc
Proof: The derivatives of Ca to Ce were computed by Miyaji and Ohnishi [MO07]:
Ca
C˙a =
(Cb Cc + Cc Cd + Cc Ce + Cd Ce ) − Ca
SLa
Cc
Ċc =
(Ca Cd + Ca Cb + Ca Ce + Cb Ce ) − Cc .
SLc
The derivatives of Cb and Cd can be obtained from the above by symmetry (exchange a with
b and c with d). We now compute x˙a :
Cc
−(C˙a Cc − Ca Ċc )
d
=
dt Ca + Cc
(Ca + Cc )2


Ca
− SL
Cc
(C
C
+
C
C
+
C
C
+
C
C
)
−
C
c
c
c
e
e
a
b
d
d
a
=
+
(Ca + Cc )2


Cc
(C
C
+
C
C
+
C
C
+
C
C
)
−
C
Ca SL
a d
a b
a e
c
b e
c
+
2
(Ca + Cc )


Ca Cd + Ca Cb + Ca Ce + Cb Ce Cb Cc + Cc Cd + Cc Ce + Cd Ce
Ca Cc
−
=
S(Ca + Cc )2
Lc
La





Ca Cc
Ca
Cc
Cb Cd
=
−
−
(Cb + Cd + Ce )
+ Ce
2
S(Ca + Cc )
Lc
La
Lc
La
Ca Cc
((Cb + Cd + Ce )(Da − Dc ) + Ce (Cb La − Cd Lc ))
=
SLa Lc (Ca + Cc )2



Ca Cc
Lb /Db
La
=
−
(Cb + Cd + Ce )(Da − Dc ) + Ce Cb Lc
.
SLa Lc (Ca + Cc )2
Lc
Ld /Dd
Finally, observe
x∗a − xa =

Da − Dc
Cc
La (Ca + Cc ) − Cc (La + Lc )
La
=
.
−
=
L a + L c Ca + Cc
(La + Lc )(Ca + Cc )
(La + Lc )(Ca + Cc )

We draw the following conclusions:
- if Ce = 0, then sign(x˙a ) = sign(Da − Dc ) = sign(x∗a − xa ). Thus, xa converges monotonically against x∗a .
- From xb + xd = 1 and x∗a + x∗c = 1, we conclude
 ∗

xa
xb
sign
−
= sign(x∗a − xb ).
x∗c
xd
26

xb

S
S
xa

M

L

RL

RL

RL

M

LR

RL
LR

L

LR

LR

Figure 5: The transition diagram under the assumption x∗a < x∗b .
- if s = sign(x∗a − xb ) = sign(x∗a − xa ), then sign(x˙a ) = s.
- if xa , xb > x∗a , then xa decreases.
- if xa , xb < x∗a , then xa increases.
- if xd , xc > x∗d , then xd decreases (equivalent to: if xa , xb < x∗b , then xb increases).
- if xd , xc < x∗d , then xd increases (equivalent to: if xa , xb > x∗b , then xb decreases).
Theorem 5 Assume x∗a < x∗b , that is, La /Lc < Lb /Ld . Then,
1. The regime xa , xb > x∗b cannot be entered. By symmetry, the regime xa , xb < x∗a cannot
be entered.
2. In the regime xa , xb ∈ [x∗a , x∗b ], xa decreases and xb increases. Hence, in this regime, the
direction of the middle edge e can change at most once.
3. If the dynamics stay in the regime xa , xb ≥ x∗b forever, xa and xb converge.
4. If the dynamics stay in the regime xa , xb ≤ x∗a forever, xa and xb converge.
Proof: At (1): In the regime xa , xb > x∗b , xa and xb both decrease, and hence, the dynamics
cannot enter the regime from the outside. More precisely, we consider two cases: xb ≥ x∗b and
xa = x∗b , or xa > x∗b and xb = x∗b .
If xb ≥ x∗b and xa = x∗b , xa is non-increasing, and hence, we cannot enter the regime.
If xa > x∗b and xb = x∗b , xb is non-increasing, and hence, we cannot enter the regime.
At (2): Obvious from the equations.
At (3): Then, xa and xb are monotonically decreasing and hence converging. The derivative of xb clearly goes to zero if xb and xa converge to x∗b .
At (4): Symmetrically to (3).

27

In Figure 5, we use S, M , and L to denote the three ranges: S = [0, x∗a ], M = [x∗a , x∗b ],
and L = [x∗b , 1]. The box M × M is divided into the triangles xa < xb and xa > xb . The figure
also shows that the boxes S × S and L × L cannot be entered and that the latter triangle
cannot be entered from the former.
We conclude the following dynamics: Either the process stays in S × S or L × L forever
or it does not do so. If it leaves these sets of states, it cannot return. Moreover, there is no
transition from the set of states RL to the set of states LR. Thus, if the process does not stay
in S × S or L × L forever, the direction of the middle edge stabilizes.
Assume now that the dynamics stay forever in S × S, or in L × L. Then, xa and xb
∞
converge. Let x∞
a and xb be the limit values. If the limit values are distinct, the direction of
the middle edge stabilizes. If the limit values are the same, the edge is horizontal and hence
stabilizes. We summarize the discussion.
Theorem 6 The dynamics of the Wheatstone graph stabilize.

9

The Uncapacitated Transportation Problem

The uncapacitated transportation problem generalizes the shortest path
P problem. With each
vertex v, a supply/demand bv is associated. It is assumed that
v bv = 0. Nodes with
positive bv are called supply nodes and nodes with negative bv are called demand nodes. In
the shortest path problem, exactly two vertices have non-zero supply/demand. A feasible
solution to the transportation problem is a flow F satisfying the mass balance constraints,
P
i.e., for every vertex v, bv is equal to the net flow out of v. The cost of a solution is e Fe Le .
The Physarum solver for the transportation problem is as follows: At any fixed time, the
current Q is a feasible solution to the transportation problem satisfying Ohm’s law (4). The
dynamics evolve according to (1).
For technical reasons, we extend G by a vertex s0 with bs0 = 1, connect s0 to an arbitrary
vertex v, and decrease bv by one. The flow on the edge (s0 , v) is equal to one at all times.
Our convergence proof for the shortest path problem extends to the transportation problem. A cut S is a set of vertices. The edge set δ(S) of the cut is the set of edges having
exactly one endpoint in S, and the capacity C
S of the cut is the sum of the D-values in the
P
cut. The demand/supply of the cut is bS = v∈S bv . A cut S is non-trivial if bS 6= 0. We
use C to denote the family of non-trivial cuts. For a non-trivial cut S, let FS = CS /bS , and
let F = min {FS ; S ∈ C}. One may view F as a scale factor; our transportation problem has
a solution in a network with edge capacities De /F . A cut S with FS = F is called a most
constraining cut.
Properties of Equilibrium Points. Recall that D ∈ RE
+ is an equilibrium point when
Ḋe = 0 for all e ∈ E, which is equivalent to De = |Qe | for all e ∈ E.
Lemma 26 At an equilibrium point, minS∈C CS /|bS | = C{s0 } /b{s0 } = 1.
Proof:
1 ≤ min
S∈C

X |Qe |
C{s0 }
CS
= min
≤
= 1.
S∈C |bS |
|bS |
b{s0 }

e∈δ(S)

28

Lemma 27 The equilibria are precisely the solutions to the transportation problem with the
following equal-length property: Orient the edges such that Qe ≥ 0 for all e, and let N be the
subnetwork of edges carrying positive flow. Then, for any two vertices u and v, all directed
paths from u to v have the same length.
Proof: Let Q be a solution to the transportation problem satisying the equal-length property. We show that D = Q is an equilibrium point. In any connected component of N , fix
the potential of an arbitrary vertex to zero and then extend the potential function to the
other vertices by the rule ∆e = Le . By the equal-length property, the potential function is
well defined. Let Q′ be the electrical flow induced by the potentials and edge diameters. For
any edge e = (u, v) ∈ N , we have Q′e = De ∆e /Le = De = Qe . For any edge e 6∈ N , we have
Qe = 0 = De . Thus, D is an equilibrium point.
Let D be an equilibrium point and let Qe be the corresponding current along edge e.
Whenever De > 0, we have ∆e = Qe Le /De = Le because of the equilibrium condition. Since
all directed paths between any two vertices span the same potential difference, N satisfies the
equal-length property. Moreover, by Lemma 26, minS CS /bs = 1, and hence, Q is a solution
to the transportation problem with the equal-length property.
Let E be the set of equilibria and let E ∗ be the set of equilibria of minimum cost.
Lemma 28 Let W = (C{s0 } − 1)2 . Then, Ẇ = −2W ≤ 0 with equality iff C{s0 } = 1.
Proof: Let C0 = C{s0 } for short. Then, since
Ẇ = 2(C0 − 1)

X

P

e∈δ({s0 }) |Qe |

= 1,

(|Qe | − De ) = 2(C0 − 1)(1 − C0 ) = −2(C0 − 1)2 ≤ 0.

e∈δ({s0 })

The following functions play a crucial role: Let F = minS∈C FS , and
VS =

1 X
Le De for each S ∈ C,
FS
e∈E

V = max VS + W, and
S∈C

h=−

1 X
1 X
Re |Qe |De + 2
Re De2 .
F
F
e∈E

e∈E

Lemma 29 Let S be a most constraining cut at time t. Then, V̇S (t) ≤ −h(t).

29

Proof: Let X be the characteristic vector of δ(S), that is, Xe = 1 if e ∈ δ(S), and Xe = 0
otherwise. Observe that FS = F since S is a most constraining cut. Let C = CS . We have
V̇S =
=

X ∂VS
Ḋe
∂De
e

X |bS |
e

C2

Le C −

X

Le′ De′ Xe

e′

|bS | X
|bS |
=
Le |Qe | − 2
C e
C

X

!

Le′ De′

e′

|bS |
|bS | X
Le De + 2
−
C e
C

≤

(|Qe | − De )
!

X

X

Xe |Qe | +

e

Le′ De′

e′

!

!

X
e

Xe De

!

b2 X
|bS | X
|bS | X
|bS | X
Re |Qe |De − S2
Le De +
Le De
Re De2 −
C e
C e
C e
C e

= −h.

The only inequality follows from Le = Re De and
least bS units of current must cross S.

P

e Xe |Qe |

≥ |bS |, which holds because at

Lemma 30 V̇ exists almost everywhere. If V̇ (t) exists, then V̇ (t) ≤ −h(t) − 2W (t) ≤ 0, and
V̇ (t) = 0 iff ∀e, Ḋe (t) = 0.
Proof: The almost everywhere existence of V̇ is shown as in Lemma 10.
The fact that W ≥ 0 is clear. We now show that h ≥ 0. To this end, let f represent a
solution to the (capacitated) transportation problem in an auxiliary network having the same
structure as G and where the capacity of edge e is set equal to De /F ; f exists by Hoffman’s
circulation theorem [Sch03, Corollary 11.2g]: observe that for any cut T , FT ≥ F , and hence,
|bT | ≤ CT /F . Then,
1 X
1 X
Re De2
Re |Qe |De − 2
F e
F e
!1/2
!1/2
X
1 X
1 X
− 2
Re De2
Re Q2e
≤
Re De2
F
F
e
e
e
!1/2
!1/2
X
1 X
1 X
− 2
Re De2
Re fe2
≤
Re De2
F
F
e
e
e
!1/2
!1/2
X
1 X
1 X
− 2
Re De2
Re De2
Re De2
≤ 2
F
F
e
e
e

−h =

= 0,

where we used the following inequalities:
P
P
P
1/2
1/2
- the Cauchy-Schwarz inequality e (Re |Qe |)(Re De ) ≤ ( e Re Q2e )1/2 ( e Re De2 )1/2 ;
30

- Thomson’s Principle (6) applied to the flows Q and f ; Q is a minimum energy flow
solving the transportation problem, while f is a feasible solution; and
- the fact that |fe | ≤ De /F for all e ∈ E.
Finally, one can have h = 0 if and only if all the above
which
P
P inequalities are equalities,
implies that |Qe | = |fe | = De /F for all e. And, W = 0 iff e∈δ({s0 }) De = 1 = e∈δ({s0 }) |Qe |.
So, h = W = 0 iff |Qe | = De for all e.
Lemma 31 The function t 7→ h(t) is Lipschitz-continuous.
Proof: The proof of Lemma 11 carries over.

Lemma 32 |De − |Qe || converges to zero for all e ∈ E.
Proof: The first and last paragraph of the proof of Lemma 12 carry over. We redo the
second paragraph.
The first paragraph establishes that for any ε > 0, there is t0 such that h(t) ≤ ε for all
t ≥ t0 . Then, recalling that Re ≥ Lmin /2 for all sufficiently large t (by Lemma 4), we find
2 X
2

X Lmin  De
De
Re
− |Qe | ≤
− |Qe |
2
F
F
e
e
X
2 X
1 X
Re |Qe |De
Re De2 +
Re Q2e −
= 2
F e
F e
e
2 X
2 X
≤ 2
Re |Qe |De
Re De2 −
F e
F e

= 2h ≤ 2ε,
P
P
2
2
2
where we used once more the inequality
e Re De /F , which was proved in
e Re Qe ≤
Lemma 30. This implies that for each e, De /F − |Qe | → 0 as t → ∞. Summing across
e ∈ δ({s0 }) and using Lemma 4(ii), we obtain C{s0 } /C − 1 → 0 as t → ∞. From Lemma 4,
C{s0 } → 1 as t → ∞, so C → 1 as well.
We are now ready to prove that the set of equilibria is an attractor.
Theorem 7 The dynamics are attracted by the set E of equilibria.

Proof: Assume otherwise. Then, there is a network and initial conditions for which the
dynamics has an accumulation point D that is not an equilibrium; such an accumulation
point exists because the dynamics are eventually confined to a compact set. Let Q be the
flow corresponding to D. Since D is not an equilibrium, there is an edge e with De 6= |Qe |.
This contradicts the fact that |De − |Qe || converges to zero for all e.
Theorem 8 If no two equilibria have the same cost, the dynamics converge to a minimum
cost solution.
31

Proof: Consider any equilibrium D ∗ , and let Q∗ = D ∗ be the corresponding flow. Let T ∗
be the edges carrying non-zero flow. T ∗ must be a forest, as otherwise, there would be two
equilibria with the same cost. Consider any edge e = (u, v) of T ∗ , and let S be the connected
component of T ∗ \e containing u. Then Q∗e = b(S), and hence, distinct equilibria have distinct
associated forests. We conclude that the setP
of equilibria is finite.
The V -value of D ∗ is equal to the cost e Le Q∗e of the corresponding flow since W = 0
and F = 1 in an equilibrium. If no two equilibria have the same cost, the V -values of distinct
equilibria are distinct.
V is a decreasing function and hence converges. Since the dynamics are attracted to the
set of equilibria, V must converge to the cost of an equilibrium. Since the equilibria are a
discrete set, the dynamics must converge to some equilibrium. Call it D ∗ .
We next show that D ∗ is a minimum cost solution to the transportation problem. Orient
the edges in the direction of the flow Q∗ . If Q∗ is not a minimum cost flow, there is an
oriented path P from a supply node u to a demand node v such that Qe > 0 for all edges of
P , and P is not a shortest path from u to v. The potential difference ∆uv converges to LP .
We now derive a contradiction as in the proof of Lemma 13.
′
∗
′
′
P Let P be a shortest path from u to v in G, let L = LP be its length, and let WP =
L
ln
D
.
We
have
′
e
e
e∈P
X
X
X Le
(|Qe | − De ) =
|∆e | −
Le ≥ pu − pv − LP ′ = ∆uv − L∗ .
ẆP ′ =
D
e
′
′
′
e∈P

e∈P

e∈P

Let δ > 0 be such that there is no path from u to v with length in the open interval (L∗ , L∗ +
2δ). Then, ∆ − L∗ ≥ δ for all sufficently large t, and hence, ẆP ′ ≥ δ for all sufficiently large
t. Thus, WP ′ goes to +∞. However, WP ′ ≤ nLmax for all sufficiently large t since De ≤ 2 for
all e and t large enough. This is a contradiction.
Lemma 33 The problem of minimizing V (D) for D ∈ RE
+ is equivalent to the transportation
problem.

Proof: By introducing an additional variable F = minS CS /|b(S)| > 0, the problem of
minimizing V (D) is equivalently formulated as
2

X
X
1
De − 1
min
Le De + 
F e
e∈δ({s0 })

s.t. CS /|b(S)| ≥ F

∀S ∈ C

F >0
D ≥ 0.
Substituting xe = De /F , we obtain
min

X
e

s.t.



Le xe + F 1/2 

X

X

xe −

e∈δ({s0 })

xe ≥ |b(S)|

e∈δ(S)

x ≥ 0, F > 0
32

∀S ∈ C

2

1
F

which is easily seen to be equivalent to the (fractional) transportation problem.

References
[Ada10]

Andrew Adamatzky. Physarum Machines: Computers from Slime Mold. World
Scientific Publishing, 2010.

[Bol98]

Béla Bollobás. Modern Graph Theory. Springer, New York, 1998.

[Cha09]

Bernard Chazelle. Natural algorithms. In Proc. 20th SODA, pages 422–431, 2009.

[CKM+ 11] Paul Christiano, Jonathan A. Kelner, Aleksander Madry, Daniel A. Spielman, and
Shang-Hua Teng. Electrical flows, Laplacian systems, and faster approximation of
maximum flow in undirected graphs. In Proc. 43rd STOC, pages 273–282, 2011.
[CLSW98] F.H. Clarke, Yu.S. Ledyaev, R.J. Stern, and P.R. Wolenski. Nonsmooth Analysis
and Control Theory. Springer, New York, 1998.
[HS74]

M.W. Hirsch and S. Smale. Differential Equations, Dynamical Systems, and Linear Algebra. Academic Press, 1974.

[IJNT11]

Kentaro Ito, Anders Johansson, Toshiyuki Nakagaki, and Atsushi Tero. Convergence properties for the Physarum solver. arXiv:1101.5249v1, January 2011.

[MN92]

Marko M. Mäkelä and Pekka Neittanmäki. Nonsmooth Optimization. World
Scientific, Singapore, 1992.

[MO07]

T. Miyaji and Isamu Ohnishi. Mathematical analysis to an adaptive network of
the plasmodium system. Hokkaido Mathematical Journal, 36:445–465, 2007.

[MO08]

T. Miyaji and Isamu Ohnishi. Physarum can solve the shortest path problem on
Riemannian surface mathematically rigourously. International Journal of Pure
and Applied Mathematics, 47:353–369, 2008.

[NTK+ 09] T. Nakagaki, A. Tero, R. Kobayashi, I. Ohnishi, and T. Miyaji. Computational
ability of cells based on cell dynamics and adaptability. New Generation Computing, 27:57–81, 2009.
[NYT00]

T. Nakagaki, H. Yamada, and Á. Tóth. Maze-solving by an amoeboid organism.
Nature, 407:470, 2000.

[Sch03]

Alexander Schrijver. Combinatorial Optimization – Polyhedra and Efficiency.
Springer, 2003.

[TKN07]

A. Tero, R. Kobayashi, and T. Nakagaki. A mathematical model for adaptive
transport network in path finding by true slime mold. Journal of Theoretical
Biology, pages 553–564, 2007.

[TTS+ 10] A. Tero, S. Takagi, T. Saigusa, K. Ito, D. Bebber, M. Fricker, K. Yumiki,
R. Kobayashi, and T. Nakagaki. Rules for biologically inspired adaptive network
design. Science, 327:439–442, 2010.
33

[You10]

http://www.youtube.com/watch?v=tLO2n3YMcXw&t=4m43s, 2010.

34

s0
a

b

e
L

R

c

d

s1

