A revised model of fluid transport optimization
in Physarum polycephalum
Vincenzo Bonifaci∗

arXiv:1606.04225v1 [q-bio.TO] 14 Jun 2016

Abstract
Optimization of fluid transport in the slime mold Physarum polycephalum has been the subject of
several modeling efforts in recent literature. Existing models assume that the tube adaptation mechanism
in P. polycephalum’s tubular network is controlled by the sheer amount of fluid flow through the tubes.
We put forward the hypothesis that the controlling variable may instead be the flow’s pressure gradient
along the tube. We carry out the stability analysis of such a revised mathematical model for a paralleledge network, proving that the revised model supports the global flow-optimizing behavior of the slime
mold for a substantially wider class of response functions compared to previous models. Simulations also
suggest that the same conclusion may be valid for arbitrary network topologies.

1

Introduction

Physarum polycephalum is an amoeboid slime mold that exhibits remarkable information processing capabilities. In controlled experiments, the slime mold’s abilities have been leveraged to determine the shortest
path between two locations in a network [Nakagaki et al., 2000, Tero et al., 2006] and, more generally, to
adaptively form efficient transport networks [Tero et al., 2010]. The question remains, however, of detailing
and analyzing the underlying mechanisms and goals of such an optimization process, which have been only
partially explored [Tero et al., 2007, Miyaji and Ohnishi, 2007, Tero et al., 2010, Ito et al., 2011, Bonifaci
et al., 2012, Ma et al., 2013, Bonifaci, 2013].
P. polycephalum, like other Myxomycetes, has a somewhat complicated life cycle consisting of several
stages [Stephenson and Stempen, 2000]. During its mature plasmodium stage it forms a single, giant
multinucleate acellular structure. This acellular structure takes the form of a tubular vein network, through
which protoplasm is periodically transported, driven by the gradient of hydrostatic pressure.
In a remarkable experiment by Nakagaki, Yamada, and Tóth [2000], the plasmodium has been placed
on a preexisting artificial network structure and two food sources (oat flakes) have been laid over two nodes
s0 , s1 of the network. P. polycephalum reacted by adapting its shape dynamically, by controlling the width
of the tubular structures forming its veins, based on feedback from the protoplasmic flow. Gradually, several
branches of the vein network collapsed, eventually leaving the mass of the slime mold only along the shortest path between the two food sources s0 and s1 , thus exhibiting evidence of fluid transport optimization.
Tero, Kobayashi, and Nakagaki [2007] were the first to propose a mathematical model for the transport optimization dynamics of P. polycephalum. For a critical value of the model parameters, the model’s
dynamics indeed provably converge to the shortest path between the two food-source terminals of the underlying network. Such a convergence to the shortest path has been analytically proven by Tero, Kobayashi,
and Nakagaki [2007] for a ring-shaped network, by Miyaji and Ohnishi [2007] for a Wheatstone bridgeshaped network, and finally by Bonifaci, Mehlhorn, and Varma [2012] for arbitrary complex networks.
However, the global convergence of the dynamics to the shortest path fails for other values of the model’s
power-law exponent, called µ by Tero et al. [2007], as well as for nonzero values of the parameter α controlling saturation of the tubular dynamics. In this article we propose a revised model that, while enlarging
the class of admissible response functions, has the property of converging to the shortest path in the network
from any initial condition, independently of the specific details of the response dynamics or of the values of
∗ Istituto di Analisi dei Sistemi ed Informatica, Consiglio Nazionale delle Ricerche, Rome, Italy. Email: vincenzo.bonifaci@
iasi.cnr.it

1

its parameters. We prove this analytically for a parallel-edge network, and test our claim by simulation on
more complex topologies.
Fluid transport in P. polycephalum is based on the positive feedback mechanism between the width of the
veins and the fluid flow. On a longer time scale, it is driven by network peristalsis, which reverses the flow
velocities periodically [Alim et al., 2013, Baumgarten and Hauser, 2013]. It has been found that the time
between reversals is much larger than the time required for the veins to adapt their widths under a steady
flow. In this work we focus on the vein adaptation dynamics, without attempting to model the peristaltic
oscillations. The reader is referred to Kobayashi et al. [2006] for one possible model of the rhythmic
protoplasmic movement. Other details on the physical underpinnings of Physarum’s tubular dynamics are
discussed by Tero et al. [2005]. Finally, the slime mold network dynamics have also been considered from a
combinatorial optimization perspective in the discrete algorithms literature [Becchetti et al., 2013, Straszak
and Vishnoi, 2016].

1.1

Outline of the paper

We recall the Tero-Kobayashi-Nakagaki network adaptation model in Section 2.1, and present a revised
model in Section 2.2. A stability analysis of the new model’s equilibria is carried out in Section 3 for networks consisting of parallel edges; it is proved that the unique stable fixed point corresponds to the shortest
path in the graph, independently of the power-law exponent and of the saturation parameter, or, indeed, of
the particular shape of the tubular response functions. Moreover, any nontrivial trajectory approaches such
a fixed point. In Section 4 we report on simulations for more complex network topologies; these suggest
that the globally optimal behavior of the dynamics may hold in arbitrary topologies, again independently of
the particular shape of the tubular response functions. We close by summarizing and discussing our findings
in Section 5.

2
2.1

Mathematical model
The Tero-Kobayashi-Nakagaki model

Let G be an undirected multigraph with node set N, edge set E, edge lengths l ∈ RE>0 and two distinguished
nodes s0 , s1 ∈ N. The graph models P. polycephalum’s vein network; the edges represent the tubular channels, and the nodes represent junctions between the tubes. The two distinguished nodes s0 , s1 are two
junctions corresponding to the location of the food sources.
As food is absorbed by the organism, and protoplasmic flow is distributed through the network, the
widths of the tubular channels adapt to the flow. In our discussion, x ∈ RE>0 will be a state vector representing
the fourth powers of the radii of the tubular channels of the slime mold. For an edge e ∈ E, the value xe
is called the capacity of e. The transport optimization process described in the introduction has previously
been modeled [Tero et al., 2007] as a system of coupled, nonlinear ordinary differential equations,
ẋe = f (|qe |) − xe

for all e ∈ E.

(1)

Equation (1) is called the adaptation equation. The function f models the response dynamics of
the tubular channels to the fluid flow qe along the edge. Proposed forms include f (y) = yµ and f (y) =
(1 + α)yµ /(1 + αyµ ), where µ, α > 0 are parameters of the model. The dynamic vector q ∈ RE , called
the (fluid) flow, is determined at any time by the capacity and length of the edges, by solving a network
Poisson equation, as follows. Without loss of generality, assume that N = {1, 2, . . . , n}, E = {1, 2, . . . , m}
and assume an arbitrary orientation of the edges. Let B = (Bve )v∈N,e∈E be the incidence matrix of G under
this orientation, that is,


+1 if v is the tail of e
def
Bve = −1 if v is the head of e


0
otherwise.
Then q is defined as the unit-value flow from s0 to s1 of minimum energy, that is, as the unique optimal

2

solution to the following continuous quadratic optimization problem:
min q> Rq

(2)

s.t. Bq = b.
Here, R ∈ RE×E is the diagonal matrix with value re
b ∈ RN is the vector defined by


0
def
bv = +1


−1

def

= le /xe for the e-th element of the main diagonal, and
if v ∈
/ {s0 , s1 },
if v = s0 ,
if v = s1 .

It is well-known that a vector q is optimal for system (2) if and only if it satisfies Kirchhoff’s circuit
laws [Bollobás, 1998, Chapter IX]. In particular, Kirchhoff’s current law is expressed by the constraint
Bq = b, which, in words, requires that the flow has zero divergence everywhere except at nodes s0 and s1 .
Kirchhoff’s voltage law is implicit in the optimality condition for (2), which implies that there exist values
p1 , . . . , pn ∈ R (the node potentials) satisfying the hydrodynamic analogue of Ohm’s law [Bollobás, 1998,
Section II.1]:
qe = (pu − pv )/re ,
whenever edge e is oriented from u to v.
(3)
The sum of potential differences along any cycle of the network is thus zero.
Node s0 is the source of the flow, node s1 the sink. It is very important to remark that while the flow
has been somewhat arbitrarily directed from s0 to s1 , the opposite choice yields exactly the same dynamics,
because of the absolute value in (1); the only effect would be to replace q with −q. In other words, a flow
reversal has no effect on the veins’ dynamics.
The value re is called the resistance of edge e, while bv is the divergence of the flow q at v. The constant
bs0 (here, bs0 = 1) is the flow’s value, that is, the divergence of the flow at the source.
def

The quantity E = q> Rq is the (instantaneous) energy of the flow q. By the conservation of energy
principle, the energy of the flow equals the difference between the source and sink potentials, times the
value of the flow [Bollobás, 1998, Corollary IX.4]:
E = (ps0 − ps1 )bs0 = ps0 − ps1 .

(4)

An alternative way to express the fluid flow vector arises from the Laplacian operator of the graph
def
[Strang, 1988, Biggs, 1997]. Let C = R−1 . The Laplacian of G is the symmetric and positive semidefinite
def
matrix L = BCB> . If we represent the potential vector by p ∈ RN , Ohm’s law (3) can be written in matrix
form as
q = CB> p.
(5)
Multiplying both sides by B yields the discrete Poisson equation Lp = b, with solution p = L+ b, where
L+ is the Moore-Penrose pseudoinverse1 of L. Substituting in (5), we get
q = CB> L+ b.

(6)

The fluid flow has been observed to be laminar [Kamiya, 1950]. In this case, Poiseuille’s law expresses
the relation between the flow rate, the tube radius, and the pressure gradient: for an edge e = (u, v),
|qe | =
def 1/4

where Re = xe

πR4e |pu − pv | πxe |pu − pv |
=
,
8η
le
8η
le

(7)

is the radius of the tube, and η is the viscosity constant.

1 The Moore-Penrose pseudoinverse of a matrix A is the unique matrix A+ such that AA+ A = A, A+ AA+ = A+ , and AA+ and
A+ A are both Hermitian.

3

2.2

The revised model

The revised model that we propose resembles closely that discussed in the previous section, following in
particular Equations (2)–(7). However, since in this work we conjecture that pressure gradients, rather than
sheer flow amounts, control the response of the tubular channels, we take the controlling variables to be the
ratios |pu − pv | /le (instead of the edge flows |qe |). After appropriate normalization, Poiseuille’s law implies
that the pressure gradients are equivalent to the ratios |qe | /xe . Therefore, we replace Equation (1) by
 


|qe |
ẋe = xe fe
−1
for all e ∈ E.
(8)
xe
Each edge e ∈ E has its own dimensionless response function fe : R≥0 → R≥0 . These response functions
are assumed to satisfy the following condition.
Definition 2.1. A function f : R≥0 → R≥0 is a standard response function if:
1. f (1) = 1;
2. f is strictly increasing on R>0 ;
3. f is differentiable on R>0 .
Our main result, the stability analysis in Section 3, will not require any other property from the fe , apart
from being standard response functions. However, when contrasting our findings with those of Tero et al.
[2007] (Section 3.4) and in the simulations (Section 4), we consider for concreteness the same types of
response functions that have been considered in earlier literature:
(Type I) Nonsaturating response: f (y) = yµ , for some µ > 0;
(Type II) Saturating response: f (y) = (1 + α)yµ /(1 + αyµ ), for some µ, α > 0.
While Type I functions have a simpler structure, Type II functions have the additional property of saturating
as y → ∞, implying a finite maximum radius for the tubes, and may therefore be considered more realistic.
Note in any case that, by using
(1 + α)yµ
(9)
f (y) =
1 + αyµ
with µ > 0, α ≥ 0, one can capture both Type I (α = 0) and Type II (α > 0) responses.
We observe that, when using response functions of the form (9) with µ = 1 and α = 0, the dynamics (8)
are identical with (1). We will show, however, that in general they have a qualitatively different behavior
(Section 3.4).
For the purpose of analysis, we finally assume that the edge length vector l is such that each s0 -s1 path in
G has a distinct overall length. That is, the configuration of lengths is nondegenerate, which is a physically
realistic assumption. While this assumption is not crucial to our main results, it simplifies their statement:
for example, without this assumption, the fixed points of the system may not be isolated and the shortest
path in the network may not be unique.

2.3

Basic properties of the revised model

We end this section with a couple of simple but useful properties of the revised model.
Proposition 2.1. E = p> Lp.
Proof. By (6) and the definition of L, E = q> Rq = b> L+ BCB> L+ b = p> Lp.
Proposition 2.2. The set [0, 1]E is an attracting set for the dynamics (8).
Proof. Observe that, for any edge e ∈ E, 0 ≤ |qe | ≤ 1, since the flow on any edge cannot exceed the global
flow value, which is bs0 = 1. Therefore, as long as xe > 1, we have |qe | /xe < 1 and ẋe < xe · ( fe (1) − 1) =
0.

4

l1

s0

l1

...

1

2

s0

s1

1

2

s1

lm

l2

(a) A parallel-edge network

(b) A ring-shaped network

Figure 1: Examples of parallel-edge network topologies. Labels inside the nodes denote their identifier.
Labels s0 and s1 denote the source and sink node, respectively. The label on each edge describes its length.

3
3.1

Stability analysis for parallel-edge networks
Network structure

To allow analytical tractability, in this section we limit our discussion to simple networks consisting of two
nodes and a set of parallel edges between them (Figure 1(a)); a special case is the ring-shaped network
(Figure 1(b)). While admittedly a simplification, parallel-edge and ring-shaped networks already exhibit
a wide range of dynamical properties and have been the departure point of previous analyses; see, for
example, the discussion by Tero et al. [2007].
The incidence matrix B of a parallel-edge network has the simple structure


+1 +1 · · · +1
B=
.
(10)
−1 −1 · · · −1
Using this structure, the instantaneous energy of the system is easily derived.
Proposition 3.1. In a parallel-edge network G = (N, E),
!−1
−1

E = (tr C)

=

∑ xe /le

.

(11)

e∈E

Proof. The matrix B given by (10) yields the Laplacian

tr C
>
L = BCB =
− tr C

− tr C
tr C.



and combining (4) and Proposition 2.1 yields
E = p> Lp = (p1 − p2 )2 tr C = E 2 tr C.
Solving for E yields the claim.

3.2

Location of fixed points

A fixed point of (8) is a vector x ∈ RE≥0 such that
 


|qe |
xe · fe
− 1 = 0 for all e ∈ E.
xe
Lemma 3.2. The fixed points of (8) in a parallel-edge network G = (N, E) are exactly the standard basis
vectors of RE .

5

Proof. Since the response function fe is assumed to be standard (Definition 2.1), the unique solution to
f (y) = 1 is y = 1. The fixed point condition for (8) is thus equivalent to
(xe = 0) or (|qe | /xe = 1)

for all e ∈ E.

In a parallel-edge network, each term |qe | /xe simplifies to E /le by Ohm’s law (3) and by (4). Moreover,
we assumed that no two source-sink paths have the same length. This implies that, in a fixed point, there
cannot be two distinct xe ’s with xe 6= 0; which in turn implies that in a fixed point, xi = qi = 1 for exactly
one i ∈ E and xe = qe = 0 for all e 6= i. Conversely, it is straightforward to verify that any standard basis
vector χi of RE is a fixed point of (8).

3.3

Nature of fixed points

Recall that our assumption on the response functions fe : R≥0 → R≥0 is that they are increasing, differentiable, and such that fe (1) = 1. After substitution in (8), using |qe | /xe = E /le , the adaptation equation can
be equivalently written as
  

E
ẋe = xe fe
−1
for all e ∈ E,
(12)
le
where, as observed in Proposition 3.1, E = (∑e xe /le )−1 .
Theorem 3.3. System (12) has exactly one stable fixed point, namely, the standard basis vector χi∗ corresponding to the edge i∗ of shortest length in the network. All other fixed points are unstable.
Proof. Let x ∈ RE≥0 and i, j ∈ E (i 6= j). Direct computation of the terms of the Jacobian matrix J(x) of (12)
yields, using the substitution y = E /li ,
 
 
∂ fi E
E
− 1 + xi
Jii (x) = fi
li
∂ xi li
 
 
E
∂ y ∂ fi E
= fi
− 1 + xi
li
∂ xi ∂ y li

 

 
xe −2 0 E
E
1
= fi
− 1 − xi 2 ∑
fi
li
li
li
e le
 
 
2
E
E
E
− 1 − xi 2 fi0
,
= fi
li
li
li
 
∂ fi E
Ji j (x) = xi
∂ x j li
 
∂ y ∂ fi E
= xi
∂ x j ∂ y li
 
E2 0 E
fi
.
= −xi
li l j
li
We evaluate the Jacobian at any standard basis vector χi to obtain (for distinct i, j, k ∈ E)
Jii (χi ) = − fi0 (1)
Ji j (χi ) = −(li /l j ) fi0 (1)
J ji (χi ) = 0
J j j (χi ) = f j (li /l j ) − 1
J jk (χi ) = 0
Jk j (χi ) = 0.

6

After rearranging the rows and columns of J(χi ) so that the ith column and row are swapped with the first
column and row, respectively, we obtain the following matrix J(i) :


− fi0 (1) −(li /l2 ) fi0 (1) . . . −(li /lm ) fi0 (1)


0
f2 (li /l2 ) − 1 . . .
0


J(i) = 

..


.
0
0
0
0

...

0

fm (li /lm ) − 1

where, for j ∈
/ {1, i}, the jth element on the main diagonal is f j (li /l j ) − 1, and for j = i it is f1 (li /l1 ) − 1.
By construction, J(i) and J(χi ) have the same eigenvalues; the advantage of J(i) is that it is upper triangular,
and so its eigenvalues can be read off its main diagonal. These eigenvalues are
− fi0 (1), f2 (li /l2 ) − 1, . . . , f1 (li /l1 ) − 1, . . . , fm (li /lm ) − 1.
By our assumptions on the response functions, it holds that − fi0 (1) < 0 and the sign of f j (li /l j ) − 1 is
the same as the sign of li − l j . We conclude that the eigenvalues associated to an equilibrium point χi
are all negative if and only if li < l j for all j 6= i. Otherwise, at least one eigenvalue is positive, and the
fixed point is a source or a saddle. Consequently, there is exactly one stable fixed point of the dynamics
(12), corresponding to the edge with shortest length; all other fixed points are unstable. We remark that
this conclusion holds independently of the concrete form of the response functions, as long as they satisfy
Definition 2.1.
Note that Theorem 3.3 does not rule out the existence of periodic orbits in phase space. To exclude this
possibility, we show convergence to equilibrium for all trajectories.
Lemma 3.4. Every trajectory of (12) converges to an equilibrium as t → ∞.
Proof. We claim that the function
!
def >

V (x) = 1 x + ln E =

∑ xe − ln ∑ xe /le

e∈E

e∈E

is a Lyapunov function for (12). First observe that V is bounded from below, as all xe ≥ 0 and E is also
bounded from below (by Proposition 2.2, we can assume that xe ≤ 1 + ε for each e ∈ E, after some finite
time). Then we compute
V̇ (x) =

∂V

∑ ∂ xe ẋe

e∈E


   

E
E
1
−
x
f
−
1
e
e
∑
le
le
e∈E

  

E
E
= − ∑ xe
−1
fe
−1
le
le
e∈E

=

≤ 0,
where the last inequality follows from the facts that xe ≥ 0 for each e ∈ E, and that E /le − 1 has the same
sign as fe (E /le ) − 1 (since fe is increasing and fe (1) = 1). Moreover, from the same derivation it follows
that V̇ (x) = 0 if and only if, for all e ∈ E, either xe = 0 or E = le , that is, if and only if x is a fixed point of
(12).
Given Theorem 3.3 and Lemma 3.4, it is to be expected that all trajectories starting in the interior of the
positive orthant converge to the stable fixed point. This is indeed formalized in our final result.
Theorem 3.5. Let x(0) > 0. As t → ∞, E (t) converges to li∗ , where i∗ is the shortest edge in E, and x(t)
converges to χi∗ .

7

Figure 2: A typical phase portrait of (12) for a ring-shaped topology with a Type II response. Parameter
values are l1 = 1, l2 = 1.2, µ = 0.8, α = 1. Several solutions are shown, all converging to the shortest-path
equilibrium χ1 = (1, 0).
Proof. By Lemma 3.4, it suffices to prove the first part of the claim. Suppose by contradiction that E
def
converges to the length of some other edge e 6= i∗ . Let δ = (le − li∗ )/2li∗ > 0 and define
def

W (t) = ln xi∗ (t).
Then for all sufficiently large t, we may assume E (t) ≥ le − δ li∗ = (1 + δ )li∗ , and therefore
ẋi∗
xi∗
  

xi∗
E
=
fi∗
−1
xi∗
li∗
 
E
= fi∗
−1
li∗

Ẇ =

≥ fi∗ (1 + δ ) − 1,
so that Ẇ (t) is larger than some positive constant for all sufficiently large t. This implies W (t) → ∞. On
the other hand, Proposition 2.2 implies that for any fixed ε > 0 and all large t, xi∗ (t) ≤ 1 + ε, so that
W (t) ≤ ln(1 + ε), yielding a contradiction.
Figure 2 shows a typical phase portrait of (12) for the ring-shaped topology.

3.4

Comparison with existing models

It is useful to contrast our findings with those obtained in the original flow-based model by Tero et al. [2007]
in the case of a ring-shaped network. The flow-based model of Tero et al. undergoes a bifurcation for both
Type I and Type II response functions around the value µ = 1. When µ > 1, there are two stable equilibria
and one unstable equilibrium; when µ < 1, there is a single stable equilibrium (and two unstable ones), but
it lies in the interior of the positive orthant and therefore it does not correspond to any edge of the graph.
Only when µ = 1 the equilibria are two, one for each edge, and the stability of each equilibrium depends
on the length of the corresponding edge.

8

4

3

1.7
0.5
1

5

19

51

1

s0
1

26

s1

s0

4

1

s1
20
28

2

6

9

32
2

2

(a) Wheatstone network

34

3

(b) Another arbitrary topology

Figure 3: Examples of non-parallel-edge network topologies. Labels inside the nodes denote their identifier.
Labels s0 and s1 denote the source and sink node, respectively. The label on each edge describes its length.
The bold edges identify the shortest source-sink path.
In contrast, in the pressure-gradient based model proposed in this article, in a ring-shaped network there
are always two equilibria, one of which is stable and the other one of which is not, independently of the
details of the response function. The stable equilibrium is always the one corresponding to the shortest path
in the network.

4

Simulation of other network topologies

In this section we simulate the dynamics (8) on two more general network topologies, to which the analysis
of Section 3 does not apply.
We consider the bridge-shaped Wheatstone network of Figure 3(a), as well as the more complex network
shown in Figure 3(b). The choice of these topologies is motivated by the fact that they are inherently nonseries-parallel networks. In both figures, the shortest path is highlighted in bold; the length of the shortest
path in the two networks is 2.5 and 90, respectively.
We simulate the dynamics (8) under an Euler discretization scheme with a stepsize h = 0.1:


 
qe [t]
− 1 for all e ∈ E,
(13)
xe [t + 1] − xe [t] = h · xe [t] fe
xe [t]
q[t] = C[t]B> p[t],
+

p[t] = (BC[t]B> ) b,
C[t] = diag(x1 [t]/l1 , . . . , xm [t]/lm ),
where, as in Section 2, (BC[t]B> )+ is the pseudoinverse of the Laplacian L[t] = BC[t]B> .
As the fe , we assume identical response functions of the form (9). As the initial condition, we select the
symmetric state x[0] = 1, to ensure absence of bias towards any specific fixed point. We note that, in analogy
with the case of a parallel-edge network, each s0 -s1 path in the network is in one-to-one correspondence
with a fixed point of (13).
Since we expect each xe to approach either 0 or 1 as t → ∞, we expect l> x[t] to approach the length of
some source-sink path in the network. Therefore, we define the quantity l> x[t] to be the transport cost at
time t. The dynamics have an optimal behavior if, in the limit of large t, the transport cost approaches the
length of the shortest path.
Figures 4(a) and 4(b) plot the value of the transport cost for t ∈ {1, 10, 1000} for the two networks, as a
function of the power-law exponent µ and the saturation parameter α, with ranges µ ∈ [0.5, 1.5], α ∈ [0, 2].
The data confirm that indeed, as time passes, the transport cost approaches the optimal costs given by the
shortest path lengths, independently of the values of the parameters µ and α. However, the convergence
speed is affected by the parameters; namely, the dynamics appear to be faster when the power-law exponent
µ is large and the saturation parameter α is small.

9

1 iter.
10 iter.
1000 iter.

6.5
6
5.5
5
4.5
4
3.5

2

3

1.5

2.5
0.6

1
alpha

0.8

1
mu

0.5
1.2

1.4

0

(a) Transport cost values for the network of Figure 3(a) after 1 (top surface), 10 (middle),
1000 (bottom) Euler steps, as a function of µ and α.

1 iter.
10 iter.
1000 iter.

220
210
200
190
180
170
160
150
140
130
120
110
100
90

2
1.5
0.6

1
alpha

0.8

1
mu

0.5
1.2

1.4

0

(b) Transport cost values for the network of Figure 3(b) after 1 (top surface), 10 (middle),
1000 (bottom) Euler steps, as a function of µ and α.

Figure 4: Transport cost values

5

Discussion and concluding remarks

The optimization of transport networks is a commonly occurring feature of several natural (as well as
artificial) systems: blood vasculature and leaf venation are two examples. The fluid transport optimizing
behavior of P. polycephalum may not be surprising in light of the idea that a more efficient use of the
available resources (the size of the tubular structures) enables the organism to achieve a higher fitness.
However, the accuracy achieved by the positive feedback mechanism between the pressure gradients along
the veins and the widths of the tubular channels is somewhat remarkable: the steady state solution is not
only approximately or locally optimal; at least in the case of a parallel-edge network, it is the globally
optimal solution from the point of view of the total length of the tubes.
In previous models, based on sheer amounts of flow, this global optimization behavior was known to
rely on very specific values of the power-law exponent and of the saturation parameter (µ = 1 and α = 0).
We have shown that a model where the controlling variables are the pressure gradients is, instead, able
to support the global optimization behavior of the dynamics for a much wider class of response functions
or range of parameters. It is a natural, though perhaps formidable, open problem to prove this conjecture
analytically for network topologies that go beyond the simple parallel-edge topology we considered in this
article.

10

Acknowledgements The author would like to thank Kurt Mehlhorn and two anonymous reviewers for
suggesting several improvements, as well as Alberto Gandolfi and Carmela Sinisgalli for fruitful discussions
on P. polycephalum’s dynamics.

References
K. Alim, G. Amselem, F. Peaudecerf, M. P. Brenner, and A. Pringle. Random network peristalsis in
Physarum polycephalum organizes fluid flows across an individual. Proceedings of the National Academy
of Sciences, 110(33):13306–13311, 2013.
W. Baumgarten and M. J. Hauser. Functional organization of the vascular network of Physarum polycephalum. Physical Biology, 10:026003, 2013.
L. Becchetti, V. Bonifaci, M. Dirnberger, A. Karrenbauer, and K. Mehlhorn. Physarum can compute shortest
paths: Convergence proofs and complexity bounds. In F. V. Fomin, R. Freivalds, M. Kwiatkowska, and
D. Peleg, editors, Automata, Languages and Programming, volume 7966 of Lecture Notes in Computer
Science, pages 472–483. Springer, 2013.
N. Biggs. Algebraic potential theory on graphs. Bull. London. Math. Soc., 29:641–682, 1997.
B. Bollobás. Modern Graph Theory. Springer, New York, 1998.
V. Bonifaci. Physarum can compute shortest paths: A short proof. Inf. Process. Lett., 113(1-2):4–7, 2013.
V. Bonifaci, K. Mehlhorn, and G. Varma. Physarum can compute shortest paths. Journal of Theoretical
Biology, 309:121–133, 2012.
K. Ito, A. Johansson, T. Nakagaki, and A. Tero.
arXiv:1101.5249v1, Jan 2011.

Convergence properties for the Physarum solver.

N. Kamiya. The rate of the protoplasmic flow in the Myxomycete plasmodium. I. Cytologia, 15(3-4):
183–193, 1950.
R. Kobayashi, A. Tero, and T. Nakagaki. Mathematical model for rhythmic protoplasmic movement in the
true slime mold. Journal of Mathematical Biology, 53(2):273–286, 2006.
Q. Ma, A. Johansson, A. Tero, T. Nakagaki, and D. J. T. Sumpter. Current-reinforced random walks for
constructing transport networks. J. R. Soc. Interface, 10:20120864, 2013.
T. Miyaji and I. Ohnishi. Mathematical analysis to an adaptive network of the Plasmodium system.
Hokkaido Mathematical Journal, 36(2):445–465, 2007.
T. Nakagaki, H. Yamada, and Á. Tóth. Maze-solving by an amoeboid organism. Nature, 407:470, 2000.
S. L. Stephenson and H. Stempen. Myxomycetes: A Handbook of Slime Molds. Timber Press, 2000.
G. Strang. A framework for equilibrium equations. SIAM Review, 30(2):283–297, 1988.
D. Straszak and N. K. Vishnoi. Natural algorithms for flow problems. In R. Krauthgamer, editor, Proceedings of the 27th ACM-SIAM Symposium on Discrete Algorithms, pages 1868–1883. SIAM, 2016.
A. Tero, R. Kobayashi, and T. Nakagaki. A coupled-oscillator model with a conservation law for the
rhythmic amoeboid movements of plasmodial slime molds. Physica D, 205:125–135, 2005.
A. Tero, R. Kobayashi, and T. Nakagaki. Physarum solver: A biologically inspired method of road-network
navigation. Physica A, 363:115–119, 2006.
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, 244:553–564, 2007.
A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Nakagaki. Rules for biologically inspired adaptive network design. Science, 327:439–442, 2010.
11

