NOLTA, IEICEPaper

A simplistic dynamic circuit analogue ofadaptive transport networks in true slimemold
Hisa-Aki Tanaka 1 ,2 a) , Kazuki Nakada 2 , Yuta Kondo 1 ,Tomoyuki Morikawa 1 , and Isao Nishikawa 2
1

Graduate School of Information Systems,
The University of Electro-Communications
1-5-1 Chofugaoka, Chofu 182-8585, Japan

2

Graduate School of Informatics and Engineering,
The University of Electro-Communications
1-5-1 Chofugaoka, Chofu 182-8585, Japana)htanaka@uec.ac.jp

Received June 6, 2015; Revised October 28, 2015; Published April 1, 2016
Abstract: This paper presents a simplistic dynamic circuit analogue for an adaptive transport network model in true slime mold by Tero et al. This circuit analogue model is derived
from Tero’s model through nontrivial simpliﬁcation under certain assumptions, and it realizes
less computational complexity through a reduction of the number of variables. Despite ofits simplicity, systematic simulations conﬁrm that the shortest path search task is eﬃciently
accomplished with this model; (i) the shortest path is always identiﬁed, for various random
networks; (ii) if there are multiple, competing shortest paths in the network, they are simultaneously identiﬁed; and (iii) for random deletions of a link in the shortest path, a new shortest
path is quickly identiﬁed accordingly. The model proposed here is easily implemented on thecircuit simulator SPICE for instance, and hence the path search time will be further reducedwith certain numerical devices including automatic adaptive numerical integration schemes aswell as an acceleration method proposed in the end of the paper.
Key Words: biology-inspired algorithm, shortest path search methods, nonlinear resistive
circuits

1. Introduction
A pioneering study by Tero et al. [1, 2] has recently uncovered a new shortest path search method,motivated by observation of the amoeba-like organism (plasmodium) of the true slime mold Physarum
polycephalum [3] (hereafter abbreviated to Physarum). They showed that their method (i.e., Tero’s
model) identiﬁes the shortest path and, according to systematic numerical simulations and analysis,
this ability is robust to temporal changes in the network topology [1, 2]. Later on, some of their results

86
Nonlinear Theory and Its Applications, IEICE, vol. 7, no. 2, pp. 86–94c
IEICE2016

DOI: 10.1587/nolta.7.86

are mathematically proved by Miyaji and Ohnishi [4], where a kind of ‘local’ Lyapunov function iscleverly constructed for planar graphs. This result in [4] enabled further analytical investigations. For
instance, convergence to shortest paths are proved for all graphs [5], and for a time-discretization of
the Tero’s model [6].
Incidentally, Aono’s group has explored a more diﬃcult combinatorial optimization; the travelingsalesman problem (TSP), using Physarum [7]. Their systematic experiments showed that Physarum
ﬁnds a high quality solution for the (8-city) TSP with a high probability. Thus, Physarum-inspired
algorithms and their physical implementation are expected to explore new, primitive but sophisticatedcomputing machinery.
In this paper, we propose a nontrivial circuit analogue of Tero’s model [1, 2], which is directly
implemented on circuit simulators such as SPICE. From systematic investigations of this circuit, we
numerically verify the following patterns for all the instances considered here1 : (i) the shortest pathis always identiﬁed, (ii) multiple, competing shortest paths are simultaneously identiﬁed, and (iii) a
new shortest path is identiﬁed immediately after a link is lost in the shortest path. Also, the required
times for this shortest path search method are systematically investigated. As far as we know, suchproperties by a simplistic shortest path search algorithm have never been reported before.
Some results in this paper were presented in NOLTA’07 [8]. The present paper carefully reexamines
all the simulations and further investigates the above item (iii) of the presented method, and hence
it adds new and original contents. The remainder of this paper is organized as follows. In Section 2,
we brieﬂy review the mathematical model by Tero et al. [1, 2] and construct a circuit analogue for it.
In Section 3 we investigate the path search time and the path search process in the proposed circuits.
Then, the discussion and conclusions are presented in Section 4.

2. Mathematical model for transport network in Physarum
polycephalum and its dynamic circuit analogue
Systematic experiments demonstrate that adaptive transport networks in the true slime moldPhysarum have a shortest path search ability [3]. Motivated by this observation, Tero et al. [1,2] constructed a simple mathematical model for adaptive transport networks in Physarum (Tero’s
model, below). They numerically and analytically investigated how their model reaches the shortest path. In this section, we review Tero’s model and construct a nontrivial circuit analogue for it.
Implementation of this circuit on SPICE is also considered.

2.1 Mathematical model for transport networks in Physarum
We start by reviewing the mathematical model for adaptive transport networks in Physarum per [1,2]. This living network initially contains numerous thin ‘tubes’ ﬁlled with nutrient-transporting liquidthat sense and respond to the environment. If we set two ‘foods’ on N1 and N2 in the network, asshown in Fig. 1(a), Physarum tends to connect N1 and N2 with a shortest path of tubes by adaptively
growing (or degrading) the thickness of these tubes [3]. The dynamics of network adaptation in
Physarum [1, 2] is modeled by experimental results and physiological insights at the molecular level

Fig. 1. Network architecture in (a) adaptive transport networks [1, 2], and(b) time-dependent, nonlinear resistive circuits in the proposed method.
1
We have used 60 instances of randomly generated network with various sizes for shortest path search. See Section 3.1
for more details.

87

(for detail, see [2]), which are summarized as follows.
First, the tube connecting nodes Ni and Nj are denoted as Mij . If there are multiple tubes
1
2
connecting the same pair of nodes, we distinguish these tubes as Mij
and Mij, as shown in Fig. 1(a).
The total numbers of tubes and nodes are M and N , respectively. Then we consider the ﬂux Qij ,which passes through Mij from Ni to Nj . As Poiseuille ﬂow is assumed in the tube, the ﬂux Qij (t)is given by
Dij (t)Qij (t) =[pi (t) − pj (t)],(1)
Lijwhere pi (t) represents the pressure at Ni , Lij represents the length of Mij , and Dij (t) is the conductivity of Mij . As the total amount of liquid is conserved in the network, Kirchhoﬀ’s law holds at eachnode except at nodes with ‘foods’:N
Qij = 0, (j = 1, 2).
(2)
i=1

Assuming N1 and N2 act respectively as the source and the sink, the following two equations hold:N
i=1

Qi1 + I0 = 0,N

Qi2 − I0 = 0,

(3)

i=1

where I0 is the ﬂux from the source node and Q11 and Q22 are assumed to be 0. In accordance with
Refs. [1, 2], we set this I0 to be a positive constant.
To model the adaptation of tubular thickness to the ﬂux Qij , the conductivity Dij is assumed tochange in time asd
Dij = f (|Qij |) − rDij ,(4)dt
where f (Q) is a certain increasing function with f (0) = 0, and r is a positive constant. From thisequation, the conductivity Dij decreases to 0 by itself, but it increases when a certain amount of ﬂux
Qij is retained in the tube Mij . For simplicity, we assume the function f as f (|Qij |) = α|Qij |, where
α is a positive constant. Then, the equation for Dij is simpliﬁed asd
Dij = α|Qij | − rDij .
dt(5)Thus, Dij (t), Qij (t), pi (t), and pj (t) at time t are determined by Eqs. (1), (2), (3), and (5). Moreprecisely, Eqs. (1), (2), and (3) contain M + N linear equations with M + N variables (Qij and pi )and Eq. (5) contains M nonlinear equations with M variables (Dij ), which are well-deﬁned. Then,
Dij (t) of l.h.s. in Eq. (5), a function of Dij (t), Qij (t), and pi (t), is obtained by Eq. (5).
Note that as pi appears in the form of pi − pj in Eq. (1), there remains an arbitrariness in the valueof pi . To remove this arbitrariness, we assume p2 (t) ≡ 0 naturally. This assumption is reasonablebecause the pressure p2 always becomes 0 at the sink node N2 , as shown in the liquid transport
network of Fig. 1(a). In addition, α and r in Eq. (5) are set to 1 without loss of generality after
rescaling pi and t.

2.2 Dynamic circuit analogue of the network
A circuit analogue for Tero’s model is constructed as in Fig. 1(b). First, we regard the ﬂux Qij as−1
as the unit resistance Rij . The change ofthe current Iij , the pressure pi as the voltage Vi , and Dij−1variables is: pi ≡ Vi , Dij ≡ Rij = Vij /Iij , and Qij = (Dij /Lij )Vij = Iij /Lij , in which Vij ≡ Vi − Vj .
Then, the equivalent equations to (1), (2), (3), and (5) are obtained for Iij (t), Vij (t), and Lij . The
circuit analogue of Eqs. (1), (2), and (3) is clear and trivial, since these equations represent Ohm’s law
and Kirchhoﬀ’s law. On the contrary, the circuit analogue of Eq. (5) is nontrivial, because it is notclear how to implement an Eq. (5) analogue with simple circuit elements. However, if we assume thatthe time evolution of Vij is relatively slow, compared with that of Iij , a reduction becomes possible
as follows.

88First, from Eqs. (1) and (5), the following equation is directly obtained after setting α and r tounity:


d Iij
Iij
= |Iij | − Lij
,
(6)
Lij
dt VijVijwhere Lij is a positive constant. If Vij evolves much slower than Iij (i.e., Assumption 1; the validityof this assumption is numerically veriﬁed later in Section 3.2), then the following is valid for a certain(inﬁnitesimally) short time span,
Vijd
Iij (t) =|Iij (t)| − Iij (t),
dt
Lij

(7)

where Vij is safely approximated as a constant parameter and its values is determined with Iij (t) from
Ohm’s law and Kirchhoﬀ’s law (i.e., Eqs. (1), (2), and (3)) at each time span. Then, the solution of
Iij is given by
 
 V
ij
Iij (t) = Iij (0) exp− 1 t > 0,
(8)
Lij
or
 
 −V
ij
Iij (t) = Iij (0) exp− 1 t < 0,(9)
Lijrespectively for a positive or negative initial value of Iij (0). Now the variable Dij (t) in Eq. (5) iseliminated, and Eq. (5) is reduced to Eqs. (8) and (9) via Eq. (6). This reduction makes our pathsearch process computationally more eﬀective, since the number of variables (Iij and Vi ) is now M +N .
This M +N is reduced from the original 2M +N variables (Qij , Dij , and pi ), and the reduced number
M (= 2M + N − M − N ) is often huge when the links are dense.
Here we note, in identifying the unknown shortest path in the network, the ﬁnal direction of thecurrent on each link is not given in advance. This requires some modiﬁcations for Eqs. (8) and (9).
A simple but nontrivial construction is obtained by superimposing both currents Eqs. (8) and (9):
 
 
 V
 −V
ij
ij
− 1 t − I0 exp−1 t ,(10)
Iij = I0 exp
Lij
Lij
in which I0 and −I0 respectively correspond to Iij (0) in Eqs. (8) and (9). In this modiﬁcation, both
directions of currents are implicitly assumed in each link, which tolerates the uncertainty of the ﬁnal
currents directions. In conjunction with this uncertainty, Iij in Eq. (10) becomes 0 at t = 0 since
we have set Iij (0) in Eqs. (8) and (9) as I0 and −I0 respectively. Therefore, this time-dependent,
nonlinear voltage-dependent current source Iij is well-deﬁned for t ≥ 0, and can be stably simulated.
The identiﬁed shortest path is characterized as follows. In Eq. (10), Iij converges to I0 or −I0 if andonly if Vij converges to Lij or −Lij respectively2 , and otherwise Iij dies out to 0. This implies all|Iij | → I0 on the shortest path, which is consistent to Eqs. (2), (3) (Kirchhoﬀ’s law) and it is veriﬁed
in all simulations in this study.
Even though the above construction is somewhat heuristic, the shortest paths (and the second
shortest paths) for all 60 cases given in Section 3.1 are successfully obtained in this circuit, withSPICE3 . This result and the underlying mechanism will be investigated in detail in the next section.

3. Path search process and elapsed times
Thus far, we have derived a dynamic circuit analogue from Tero’s model [1, 2]. In simulating thiscircuit we set a constant current source (∼ 10 A) between the nodes N1 and N2 , and we connect alarge resistance R (∼ 100 kΩ) in parallel with the nonlinear current sources of Eq. (10), as shownin Fig. 1(b). The reason why this resistance R is introduced is understood as follows. Firstly, at
t = 0 there is no current on every link because Iij becomes 0 at t = 0. However, as we assumed a2
In Eq. (10) the speed of convergence in Vij is fast enough, compared with the increasing speed of t, which is not yet
proved but verified in all simulations in this study.
3
An illustrative example on SPICE (netlist) is available from the following site: http://synchro3.ee.uec.ac.jp/
netlist2015.pdf, which generates the data shown in Fig. 5(a).

89constant current source between N1 and N2 , there must be a bypass (a constant resistance) in parallel.
Secondly, it is desirable that all the current should pass through only the shortest path eventually.
So this constant resistance should be as large as possible. The most widely used circuit simulator
SPICE is used for all simulations in this study, as its numerical results are reliable and they can beeasily traced by many researchers.

3.1 Elapsed times for shortest path search
As we shall see in Section 3.2, typical path search processes consist of two stages: an early fast-evolving
transient process and a ﬁnal process slowly converging to an shortest path (which is shown later in
Fig. 3(a) for a typical case). Thus, the shortest path is obtained immediately after this fast evolving
stage in practice, which is judged by changing rates of the currents being less than a threshold (∼ 1.0
×10−3 A/s for instance).
To investigate the averaged elapsed times required for the shortest path search, we have systematically generated random networks with N nodes and M (= 4N ) links for N = 25 , 26 , 27 , 28 , 29 , and
210 , respectively4 . For each N , 10 networks are randomly generated and our shortest path searchmethod is applied to each network. Also, Dijkstra’s algorithm for ﬁnding the shortest path [10] is
applied to the same networks to verify the advantage and the correctness of our path search results.
Figure 2 shows the elapsed times for our proposed method, where a comparison is made amongelapsed times for simulations based on Dijkstra’s algorithm (plotted with ) on a PC (Dell Dimension
8300, Pentium 4 3.2 GHz CPU), elapsed times on (virtual) circuits () on the SPICE simulator, andtheir averages (×). By comparing these two elapsed times, we can clarify the essential physics behind
this path search process, which is not inﬂuenced by simulation details. Namely, we observe thefollowing patterns;

Fig. 2. Elapsed times for the proposed method.  (blue) and  (red) respectively represent the elapsed time for numerical simulations on PC (Simulation
Time) and the time for virtual circuits on the SPICE simulator (Circuit Time).
× represents the average of both times for a given number of nodes N . Data
sets ‘a,’ ‘b,’ ‘c,’ ‘d,’ and ‘e’ are exceptional slowly converging cases, as explainedin the main text.

(i) for all cases except for the ﬁve exceptional slowly converging cases ‘a,’ ‘b,’ ‘c,’ ‘d,’ and ‘e,’ thecircuit ﬁnds the shortest path within around 100 s, irrespective of the number of nodes N , and(ii) although the simulation times clearly show an increasing tendency (partly reﬂecting the timecomplexity of the numerical integration in the simulator), the average elapsed times (×) for the circuit
do not apparently show such a tendency with respect to N . Nevertheless, the exceptional cases ‘a,’‘b,’ ‘c,’ ‘d,’ and ‘e’ deviate from the averaged circuit time signiﬁcantly.

3.2 Shortest path search process
To investigate the mechanism behind the above observations (i) and (ii), we have numerically analyzed
the time courses of path search processes in detail, by setting I0 in Eq. (10) to 10 A. Figures 3(a)4Those random networks are downloaded from the DIMACS challenge site [9], which is widely used by researchers ofcombinatorial algorithms and its validity has been tested by them.

90Fig. 3. Two typical path search processes in the proposed method. Horizontal axes (Circuit Time) represent the elapsed time in the circuits. I0 is set to
10 A. (a) Fast converging case. ◦ and  represent Iij (t) for the shortest pathand the other paths respectively. (b) Slowly converging case. ◦, , and  respectively represent Iij (t) for the shortest path, the second shortest path, andthe other paths. (c) Time evolution of the shortest path solution Vij (red) andIij (black) for the fast converging case (a). (d) Time evolution of the shortest
path solution Vij (red) and Iij (black) for the slowly converging case (b). [Notethat right after t = 0 some Vij become too large to be included in Figs. 3(c),(d), since all current (10 A) passes through the bypass resistance R (100 kΩ).]and (b) respectively shows two typical examples from a fast converging case and a slowly converging
case (i.e., the data point ‘a’ in Fig. 2). In both ﬁgures, we plotted all the currents in the network,
but they quickly coalesced into some groups.
As we observe in Fig. 3(a), in a typical instance of N = 26 = 64, the shortest path (denoted by ◦)can be already identiﬁed at around t = 20 s, while other paths () quickly disappear, namely theirassociated currents go to 0.
In contrast to such fast converging cases, we have investigated slowly converging cases ‘a,’ ‘b,’ ‘c,’‘d,’ and ‘e’ in Fig. 2, as follows. For all these ﬁve cases, we veriﬁed that two groups of solutions are
initially competing in the path search process. As shown in Fig. 3(b), for the case ‘a’ in Fig. 2, twogroups of solutions (◦ and ) coexist for a certain period, and it is conjectured that the appearanceof two such competing groups results in the slow convergence to the shortest path since the second
shortest path, once formed at an early stage of the path search process, requires certain time to be
removed. Actually, it is veriﬁed that these two groups of solutions correspond to the shortest andthe second shortest paths in the network5 . Figures 4(a), (b), (c), (d), and (e) shows a summary ofthese network structures showing the shortest and the second shortest paths for the cases ‘a,’ ‘b,’‘c,’ ‘d,’ and ‘e,’ respectively, where l denotes the distance between nodes. Note these two competing
paths have path lengths that are quite close to each other. From the perspective of application, theabove observation suggests some advantage of the proposed method (and possibly of Tero’s model [1,2] as well) over the conventional combinatorial algorithms; competing multiple paths are identiﬁed
simultaneously at an early stage of the path search process.
Now, we numerically verify the validity of Assumption 1 in Section 2.2 in all simulations for shortest
paths here. In Assumption 1, we have expected Vij evolves slower than Iij , and it is partly true asshown in the typical cases of Figs. 3(c), (d). As observed in Figs. 3(c), (d), this assumption is violated5
By using Dijkstra’s algorithm [10], we have verified that these paths are the shortest and the second shortest paths,
respectively.

91

Fig. 4. Network structures for the cases (a), (b), (c), (d), and (e) with twocompeting shortest paths.

at a very early short span of the fast evolving process (0 ≤ t < 1 s). However, we note this assumption
is satisﬁed for all time after the above very early stage. This implies that the system always satisﬁes
Assumption 1 if we regard the initial conditions of Vij and Iij as their values at t = 1 s for instance.
And, this observation explains the situation how the derivation of Eqs. (8), (9) (and hence Eq. (10))
is validated, although the global convergence to the shortest path from any initial conditions is notclear from this argument.

3.3 Recovering process after random deletions of a link in the network
In Refs. [1, 2], Tero et al. showed that a new shortest path is quickly identiﬁed after random deletionsof a link in the network. To verify that this ability is retained in our path search method, we haveinvestigated the recovering process after random deletions of a link in the shortest path. Figures 5(a)and 6(a) respectively show two typical examples of the path search process for a random deletionof such a link at t = 100 s for instance, after the fast converging process shown in Fig. 3. In theseexamples, the shortest path (◦) is already identiﬁed until t = 100 s, and this path disappears rightafter the deletion of a link on the shortest path. On the other hand, the new shortest path (;veriﬁed with Dijkstra’s algorithm) emerges immediately after the deletion of a link in the shortest
path6 . Although here we only show the examples of N = 25 and N = 26 in Figs. 5 and 6 due tospace limitation, the same pattern is observed in a total of all 60 instances for networks mentioned
in Section 3.1. Thus, we expect that the proposed method can adapt quickly to a new shortest path,suggesting that the proposed method has certain robustness to temporal changes in the network, i.e.,
resilient path ﬁnding ability, at least numerically.

Fig. 5. A typical example with network size N = 25 = 32. START and
GOAL represent the source node and the sink node, respectively. (a) Emergence of the new shortest path (i.e., the second shortest path) after the deletionof a link at t = 100 s. [Note that all currents are 0 at t = 0, although it isdiﬃcult to recognize this from the graph due to their fast movement at the initial moment.] (b) The shortest path (blue) and the new shortest path (green)
obtained in the network.
6
In most cases where random deletions of a link do not destroy the second shortest link, this new shortest path is thesecond shortest path in the original network.

92

Fig. 6. A typical example with network size N = 26 = 64. (a) Emergenceof the new shortest path after the deletion of a link at t = 100 s. [Note that
all currents are 0 at t = 0, although it is diﬃcult to recognize this from thegraph due to their fast movement at the initial moment.] (b) The shortest path(blue) and the new shortest path (green) in the network.

4. Discussion and conclusions
We have proposed a dynamic circuit analogue for the shortest path search method [1, 2] which shows
unique dynamical characteristics compared with the original Tero’s model [1, 2] as well as the classicalDijkstra’s algorithm [10]. One of the interesting characteristics is that competing multiple shortest
paths (i.e., the shortest path and the second shortest path) are simultaneously identiﬁed during
the path search process. This task is known to be diﬃcult to achieve with Dijkstra’s algorithm or
its modiﬁcations [11]. Although the circuit analogue here is somewhat heuristic and all the results
are obtained numerically, the systematic simulations thus conﬁrm practical utility of the proposed
method, i.e., reduction of huge number of variables (i.e., number of links in the network) as well asa nontrivial nonlinear dynamics behind the system.
The elapsed time for the numerical shortest path search process can be further shorten by making
the time step larger in the numerical integration scheme during the slow converging process in Fig. 3(a),for instance. Also, it is worth considering to replace t in Eq. (10) with t2 by regarding this t here
as a time-dependent parameter, after the fast evolving process shown in Fig. 3(a). Although this
modiﬁcation is non-rigorous, our preliminary results show that it successfully leads to the shortest
path and reduces the simulation time, suggesting one possible practical acceleration method. Suchaccelerations of the proposed method is now ongoing and will be reported elsewhere.
In addition, a proof concerning why the shortest path is obtained with our method is now underconsideration. For instance, in relation to rigorous results in [4–6] and their references within, it isworthwhile to consider if our system or Eqs. (1), (2), (3) and Eq. (10) has a Lyapunov function.
On the other hand, it is also worth investigating if the slow dynamics of the proposed system with
Eq. (10) is reduced to the original Tero’s model with Eq. (5) after the fast evolving process, by using
mathematical techniques such as the method of multiple scales.

Acknowledgments
We thank Drs. A. Tero (Kyushu Univ.), T. Saegusa, T. Nakagaki (Hokkaido Univ.), and R. Kobayashi(Hiroshima Univ.) for stimulating discussions. We also thank Drs. T. Saito (Hosei Univ.), K.
Yamamura (Chuo Univ.), H. Asai (Shizuoka Univ.), K. Okumura (Hiroshima Inst. of Tech.), andS. Umetani (Osaka Univ.) for helpful comments.

93References[1] A. Tero, R. Kobayashi, and T. Nakagaki, “Physarum solver: A biologically inspired method ofroad-network navigation,” Physica A, vol. 363, pp. 115–119, 2006.
[2] A. Tero, R. Kobayashi, and T. Nakagaki, “A mathematical model for adaptive transport networkin path ﬁnding by true slime mold,” J. Theor. Biol., vol. 224, pp. 553–564, 2007.
[3] T. Nakagaki, H. Yamada, and Á. Tóth, “Intelligence: Maze-solving by an amoeboid organism,”
Nature, vol. 407, p. 470, 2000.
[4] T. Miyaji and I. Ohnishi, “Physarum can solve the shortest path problem on Riemannian surface
mathematically rigorously,” Int. J. Pure and Applied Mathematics, vol. 47, no. 3, pp. 353–369,
2008.
[5] V. Bonifaci, K. Mehlhorn, and G. Varma, “Physarum can compute shortest paths,” J. Theor.
Biol., vol. 309, pp. 121–133, 2012.
[6] L. Becchetti, V. Bonifaci, M. Dirnberger, A. Karrenbauer, and K. Mehlhorn, “Physarum cancompute shortest paths: convergence proofs and complexity bounds,” Automata, Languages,
and Programming, pp. 472–483, Springer, Berlin Heidelberg, 2013.
[7] L. Zhu, M. Aono, S.-J. Kim, and M. Hara, “Amoeba-based computing for traveling salesman
problem: Long-term correlations between spatially separated individual cells of Physarum polycephalum,” BioSystems, vol. 112, no. 1, pp. 1–10, 2013.
[8] Y. Kondo and H.-A. Tanaka, “An electric circuit analogue of a mathematical model for adaptive
transport network in true slime mold,” Proc. NOLTA’07, pp. 521–524, Vancouver, Canada, 16–
19 September 2007.
[9] DIMACS, Ninth DIMACS Implementation Challenge – Shortest Paths, 2006. http://www.dis.
uniroma1.it/˜challenge9/.
[10] T. Ibaraki, Algorithms and Data Structures in C, Shokohdo, 1999 (in Japanese).
[11] S. Umetani, private communications, 2015.

94