# When does the Physarum Solver Distinguish the Shortest Path from other
Paths: the Transition Point and its Applications

Yusheng Huang, Dong Chu, Joel Weijia Lai, Yong Deng, Kang Hao Cheong Y. Huang,
D. Chu and Y. Deng are with the Institute of Fundamental and Frontier Science,
University of Electronic Science and Technology of China, Chengdu, 610054,
China (e-mail<EMAIL_ADDRESS>prof.deng@hotmail.com). J.W. Lai and
K.H. Cheong is with Science and Math Cluster, Singapore University of
Technology and Design (SUTD), S487372, Singapore
(e-mail:kanghao_cheong@sutd.edu.sg).

###### Abstract

Physarum solver, also called the physarum polycephalum inspired algorithm
(PPA), is a newly developed bio-inspired algorithm that has an inherent
ability to find the shortest path in a given graph. Recent research has
proposed methods to develop this algorithm further by accelerating the
original PPA (OPPA)’s path-finding process. However, when does the PPA
ascertain that the shortest path has been found? Is there a point after which
the PPA could distinguish the shortest path from other paths? By innovatively
proposing the concept of the dominant path (D-Path), the exact moment, named
the transition point (T-Point), when the PPA finds the shortest path can be
identified. Based on the D-Path and T-Point, a newly accelerated PPA named
OPPA-D using the proposed termination criterion is developed which is superior
to all other baseline algorithms according to the experiments conducted in
this paper. The validity and the superiority of the proposed termination
criterion is also demonstrated. Furthermore, an evaluation method is proposed
to provide new insights for the comparison of different accelerated OPPAs. The
breakthrough of this paper lies in using D-path and T-point to terminate the
OPPA. The novel termination criterion reveals the actual performance of this
OPPA. This OPPA is the fastest algorithm, outperforming some so-called
accelerated OPPAs. Furthermore, we explain why some existing works
inappropriately claim to be accelerated algorithms is in fact a product of
inappropriate termination criterion, thus giving rise to the illusion that the
method is accelerated.

###### Index Terms:

Physarum polycephalum inspired algorithm, Bio-inspired algorithm, Physarum
Solver, Shortest path problem

##  I Introduction

Bio-inspired algorithm (BA) has attracted attention and intrigued many due to
its wide range of applications for a long time. Numerous kinds of BAs, such as
the differential evolution algorithm[1], grey wolf optimizer [2], artificial
bee colony algorithm [3], genetic algorithm [4], ant colony algorithm
(ACO)[5], and whale optimization algorithm [6] have been developed. BAs have
played important roles in solving optimization problems, such as large-scale
multiobjective detection and optimization [7, 8], parameter estimation [9],
global numerical optimization [10], feature selection [11], numerical function
optimization [3], and boosting support vector machine [12], just to name a
few. In the past decade, a new BA, named the physarum polycephalum inspired
algorithm (PPA) [13, 14], also called the physarum solver, has attracted great
attention.

The physarum polycephalum, a kind of slime mold, is a single-celled amoeboid
organism [13]. In 2000, Nakagaki et al. [15] observed that the amoeboid
organism would gradually change its shape to improve its efficiency of
foraging, which leads to the result that its body would only cover the
shortest path in the mazes after some time. The maze-solving ability of the
slime mold was then considered as a kind of primitive intelligence [16]. It is
believed that the maze-solving mechanism which is also the path-finding
mechanism was related to the contraction waves in the organism [17]. The
contraction waves influence the thickness of the tube, resulting in a positive
mechanism. If the flow through a certain tube persists or increases for a
certain period, the tube will become thicker; otherwise, the tube will
contract [18].

In 2007, a mathematical model of the physarum polycephalum, which is the
original physarum polycephalum inspired algorithm (OPPA) was proposed [14]. In
the OPPA, the tube-liked body of the physarum polycephalum is considered as a
network with the growing point and the food source serving as the starting
point and ending point, respectively; the Poiseuille flow is adopted to model
the flow flowing through the tube; the conservation law of flux (also called
continuity of flux) is obeyed using network Poisson equation [13]. At the core
of the OPPA is the use of a novel adaptive equation to model the dynamics of
tube thickness. By using the techniques mentioned above, the OPPA perfectly
models the path-finding behavior of the amoeboid organism and was then applied
to solve the minimum-risk path-finding problems [14] and adaptive network
design problems [13]. The effectiveness of the OPPA has been mathematically
proven. Bonifaci et al. provided a mathematical proof of the convergence of
the OPPA on the shortest path problems [19] and provided the time complexity
bounds [20]. Karrenbauer et al. developed on the former proof and proposed a
more general proof for the effectiveness of the OPPA [21].

Many relevant research have been proposed after the initial introduction of
the OPPA. Some studies adopt the PPA to improve the performance of other BAs.
For example, PPA is utilized to initialize the ant colony algorithm to solve
some NP-Hard problems [22]; Gao et al. combined some traditional BAs with PPA
in network community detection; a combined PPA framework with the genetic
algorithm is proposed to handle the traveling salesman problem [23], just to
name a few. Other studies aim to expand the PPA’s applications. For example,
Sun et al. [24] proposed a physarum-inspired approach to identify sub-networks
for drug re-positioning; Xu et al. [25] developed a multi-sink-multi-source
PPA to tackle the traffic assignment problem; Tsompanas et al. [26], and Zhang
et al. [27] also worked on finding equilibrium state of traffic assignments
using PPA; Jiang et al. [28] applied the PPA in routing protocol design; Song
et al. [29] designed a PPA-based optimizer for minimum exposure problem, etc.

Besides the research mentioned above, some researchers are also interested in
improving the convergence rate of the OPPA. The computational complexity for
solving the network Poisson equation using PPA is $O(n^{3})$. While it is
polynomial time, the typical number of variables is of a much higher order,
thus resulting in an inefficient allocation of time resource before converging
to the shortest path [18]. Thus, several studies related to accelerated OPPA
have been developed. Zhang et al. [18] proposed an improved OPPA by
introducing the concept of energy. In the improved PPA, it is assumed that the
foraging of the physarum polycephalum consumes energy and the law of the
conservation of energy should be satisfied. In another work, the OPPA is
accelerated by eliminating some near vanished tubes of the OPPA [30]. Wang et
al. [31] believed that the tubes in the shortest path would be more
competitive than others when OPPA converges. Thus they proposed an
anticipation mechanism to terminate the OPPA earlier. Cai et al. [32] combined
the OPPA with the Bayesian rule to achieve a higher convergence rate and also
proposed a method to tackle negative-weighted edges in the OPPA. Gao et al.
[33] developed an accelerated OPPA by excluding the inactive nodes and
collecting the near-optimal paths. The above researchers share the same
motivation—while the OPPA is effective, its efficiency can be greatly
improved; hence accelerating the OPPA is necessary. These work have sped up
the OPPA using different methods. Yet, despite significant achievements being
made by previous studies, there remain some questions that need to be
answered:

  1. 1)

In the case of ensuring accuracy, is there an inherent limit so that no matter
how the termination criteria changes, the OPPA cannot increase the convergence
rate? Both [31] and [33] stated the necessity of defining an appropriate
stopping condition for the OPPA and claimed that the proposed solutions could
terminate the OPPA earlier. This leads to a secondary question, how much
earlier? Is there a limit for the early termination, prior to which the OPPA
have not converged to the shortest path?

  2. 2)

Is there a better evaluation method that could more intuitively evaluate how
much improvement the accelerated methods are making? All of the studies
mentioned above used either the running time or the number of iterations to
compare different accelerated OPPAs. However, the two metrics are highly
dependent on the chosen stopping criteria, in other words, if the evaluation
method is dependent on the stopping criteria, the comparison may not be
objective. Hence, is there an evaluation method that is not dependent to the
stopping criteria of the OPPA?

The two questions above ultimately lead to the main question: when does the
OPPA distinguish the shortest path from the other paths?

The contributions of this work are listed below:

  * •

We have defined two concepts, i.e., the dominant path and the transition point
of the OPPA. By combining the use of these two concepts, we identify the exact
moment that the OPPA starts to detect the shortest path for the first time.
These concepts answer the first and main questions.

  * •

Combining the defined concepts, we propose a new stopping criterion for the
OPPA. By combining the proposed criterion and the OPPA, a new algorithm named
the OPPA-D is developed. According to the experimental results in Section
IV-A, the OPPA is the fastest among all of the tested algorithms. The validity
and the superiority of the proposed criterion are also demonstrated in Section
IV-B.

  * •

By adopting the concept of the transition point, we have developed a novel
method to compare different accelerated OPPAs more objectively. By doing so,
the second question can be answered. Using the proposed methods, we compare
our algorithm with existing accelerated OPPAs in Section IV-C where we also
present some key findings.

Broadly, the paper is developed as follows: Section II introduces the
background and development of the PPA; Section III provides the defined
concepts and the proposed algorithm; the proposed algorithm is then
experimentally tested in Section IV; Section IV also provides the proposed
evaluation method and the corresponding analysis. Section V discusses and
concludes the paper.

##  II Brief introduction of the original PPA

This section is mainly based on [13, 14], where more details could be found.
The pseudo-code of the OPPA is provided in Section II-B for the readers’
reference.

###  II-A The mathematical model of the physarum polycephalum

As mentioned, the slime mold’s body would form a tube-like network when
foraging. Consider the network as a graph $G(N,E)$ where $N$ is the set of
nodes and $E$ represents the set of edges. Note that $G$ is assumed to be an
undirected connected graph without cycles and negative-weight edges, except
where specifically mentioned. Other abstract mappings are as follows: the
tubes of the network are treated as edges in $G(N,E)$; the junctions between
the tube are regarded as nodes in $G(N,E)$; the growing point of the slime
mold is mapped to the ending node of the $G(N,E)$ while the food source is the
starting node. For simplicity, in this paper, the scenario with one starting
node and one ending node is considered. For scenarios with multiple ending
nodes please refer to [34] and as for scenarios with multiple starting nodes
please refer to [25].

In this paper, $N_{1}$ and $N_{2}$ denotes starting and the ending nodes,
respectively; the other nodes are labeled as $N_{3},N_{4},N_{5},\cdots$;
$E_{ij}$ represents the edge between nodes $N_{i}$ and $N_{j}$; $Q_{ij}$
represents the flux flowing through the edge $E_{ij}$ from node $N_{i}$ to
node $N_{j}$.

In [14], the flow in the slime mold’s network is regarded as laminar. Thus,
the Hagen-Poiseuille equation is observed, giving us

$Q_{ij}=\frac{D_{ij}}{L_{ij}}(p_{i}-p_{j}),$ (1)

where $p_{i}$ denotes the pressure at node $N_{i}$; $D_{ij}$ represents the
conductivity of the edge $E_{ij}$ and $L_{ij}$ is the length of the edge.

The whole network is driven by the inflow at the food source (the starting
node). Thus, at $N_{1}$, the source term of flux is

$\sum_{i\not=1}{Q_{i1}=IN_{0}}$ (2)

where $IN_{0}$ is the inflow of the network.

The outflow of the network is considered equal to the inflow which is
$IN_{0}$. Then, at the ending node ($N_{2}$) of the graph, we have

$\sum_{i\not=2}{Q_{i2}=-IN_{0}}$ (3)

The inflow and outflow at other nodes should be balanced, that is

$\sum_{j\not=1,2;i\not=j}{Q_{ij}=0}.$ (4)

By combining (1) to (4), the network Poisson equation of the graph $G(N,E)$ is

$\displaystyle\sum_{i\not=j}{\frac{D_{ij}}{L_{ij}}(p_{i}-p_{j})}=\begin{cases}+IN_{0}&\text{for}\
j=1,\\\ -IN_{0}&\text{for}\ j=2,\\\ 0&\text{otherwise}.\end{cases}$ (5)

All $p_{i}$ can be calculated by solving the network Poisson equation, i.e.,
(5), by setting $p_{1}$ to $0$. After obtaining $p_{i}$, $Q_{ij}$ can be
determined using (1). However, we still need $D_{ij}$ to solve (5). We now
discuss the method to update $D_{ij}$ at each iteration which is the main
feature of the OPPA.

In [13, 14], the following adaptation equation is adopted to model the
dynamics of the thickness of tubes:

$\frac{d}{dt}D_{ij}=f(|Q_{ij}|)-\alpha D_{ij},$ (6)

where $f(|Q_{ij}|)=|Q_{ij}|$ and $\alpha=1$ are typically used. Further
discussion on different parameter settings of the adaptation equation can be
found in previous literature [14, 21]. This adaptation equation suggests that
the conductivity of tube tends to decrease exponentially while it increases
linearly with flux along this tube.

In order to implement the adaptation equation, we discretize (6) by performing
linear approximations to get

$\frac{D_{ij}^{n+1}-D_{ij}^{n}}{\Delta t}=|Q_{ij}^{n}|-D_{ij}^{n+1},$ (7)

where $Q_{ij}^{n}$ is the flow flowing through edge $E_{ij}$ at $n$th
iteration; $\Delta t$ is normally considered as $1$ (please refer to [33] for
various other considerations); $D_{ij}^{n}$ and $D_{ij}^{n+1}$ represents the
conductivity of edge $E_{ij}$ at $n$-th and $(n+1)$-th iteration,
respectively. A more concise form of the above formula is

$D_{ij}^{n+1}=\frac{|Q_{ij}^{n}|+D_{ij}^{n}}{2}.$ (8)

###  II-B The pseudo-code of the OPPA

Having introduced the mathematical model for the physarum solver, we will now
provide the pseudo-code of the OPPA. For the rest of the paper, ‘OPPA’ will
refer to the following algorithm which uses the original PPA to solve the
shortest path problem. The pseudo-code of the OPPA is demonstrated in
Alg.II-B.

Algorithm 1 The OPPA [13]

1: //Initialization part

2: Input: the statistics of graph $G(N,E)$ and the corresponding matrix of the
weight of each edges W.

3: $D_{ij}\Leftarrow 0.5\ (\forall E_{ij}\in E,)$ // Initialize the
conductivity of each edge.

4: $Q_{ij}\Leftarrow 0\ (\forall E_{ij}\in E)$ // Initialize the flow of each
edge.

5: $L_{ij}\Leftarrow W_{ij}(\forall E_{ij}\in E)$ //The length of the edges
equal to their weights.

6: $p_{i}\Leftarrow 0\ (\forall i=1,2,\cdots,N)$ // Initialize the pressure at
each node.

7: $count\Leftarrow 1$ // Initialize the counting variable.

8: //Iterative part

9: repeat

10: $p_{1}\Leftarrow 0$ // The pressure of the starting node (Node $1$) is set
to 0.

11: Calculate the pressure of all the nodes using (5):  $\sum_{i\in
N}{\frac{D_{ij}}{L_{ij}}(p_{i}-p_{j})=\left\\{\begin{gathered}+IN_{0}\ \ for\
j=1,\\\ -IN_{0}\ \ for\ j=2,\\\ 0\ \ \ otherwise.\\\ \end{gathered}\right.}$

12: Calculate the flux using (1):  $Q_{ij}\Leftarrow
D_{ij}\cdot(p_{i}-p_{j})/L_{ij}.$

13: Calculate the conductivity of the iteration using (8):
$D_{ij}^{n+1}=\frac{|Q_{ij}^{n}|+D_{ij}^{n}}{2}.$

14: $count\Leftarrow count+1$

15: until The given termination criterion is met.

16: Output: The flow matrix Q.

Given a graph, after the execution of the OPPA, according to the flow matrix
Q, one would find a path that contains the most amount of inflow in the graph,
this is the shortest path.

##  III Defining concepts and the proposed OPPA-D algorithm

###  III-A The dominant path

Dominant path (D-Path): Given a flow matrix generated by the PPA-based
algorithms, the D-path is the path found by the following iterative
procedures:

  1. i.

Set up a node list $List_{Node}$ initialized as an empty list. The starting
node $N_{1}$ of the graph would be the first element of $List_{Node}$.

  2. ii.

To find the next element of $List_{Node}$ in the graph, find all the nodes and
their corresponding edges that are connect to the current last element of
$List_{Node}$.

  3. iii.

Among the edges found in the second step, find the edge that contains the
maximum flow. Since one of the two nodes that relate to this edge is the last
element of the node set $List_{Node}$, set the other node as the next element
of $List_{Node}$.

  4. iv.

Delete the current last element of $List_{Node}$ and its corresponding edges
from the graph. Push the next element of $List_{Node}$ into the list
$List_{Node}$. This element would be the last element of $List_{Node}$ in the
next iteration.

  5. v.

Repeat Step ii to Step iv until the next element we found in Step iv is the
ending node of the graph. The nodes in $List_{Node}$ form one of the paths in
the graph. This path is defined as the D-Path.

Given a flow matrix generated by the PPA-based algorithms, the D-Path can be
generated following the above procedures, and the length $L_{D-Path}$ of the
path can be calculated.

###  III-B The transition point of the physarum solver

With the definition of D-Path, it is interesting to note that after the
convergence of the OPPA, the algorithm will direct most (or even all) of the
inflow to the shortest path; thus, when the OPPA converges, the D-Path
generated is exactly the shortest path of the graph, and the $L_{D-Path}$ will
be equal to the length of the shortest path. One can also deduce that the
process of the OPPA finding the shortest path is indeed the process of OPPA
gradually converging the D-Path to the shortest path.

However, the means to quantify the convergence of the OPPA is still an open
problem. The convergence of the OPPA is typically regarded as the iteration
when the conductivity of each edge in the graph does not change with
significantly with the subsequent iterations. Thus, the stopping (or
convergence or termination) criterion is usually set to $\sum_{i}\sum_{j\neq
i}D_{ij}\leq\epsilon$. In [33] $\epsilon$ is set to $10^{-5}$ while in [18,
30] it is set to $10^{-2}$. There is still no conclusive method or consensus
for determining the parameter $\epsilon$.

In most of the current studies, we allow the algorithm to run until it
achieves the pre-determined stopping criterion. This would often result in
long computational time; yet it remains a widely implemented method for
deciding whether the OPPA finds the shortest path. The question remains, does
the OPPA need such a long computational time to achieve convergence?

With our definition of D-Path, the answer is no. The OPPA distinguishes the
shortest path before the stopping criterion is met. To demonstrate this, four
complete graphs, i.e., ’Instance 1’ to ’Instance 4’, are randomly generated
with five nodes and the length of edges varying from 1 to 10000. The OPPA is
then implemented in the four instances to find the shortest path. The
termination criterion is set to $\sum_{i}\sum_{j\neq i}D_{ij}\leq 10^{-3}$. At
each iteration, the $L_{D-Path}$ is recorded.

Figure 1: The transition point.

Fig.1 demonstrates the results. According to Fig.1, in all of the instances,
the $L_{D-Path}$ fluctuates at the first few iterations, and then decreases
rapidly, after which it converges to the length of the shortest path. Of
importance is that the convergence of $L_{D-Path}$ happens long before the
OPPA’s termination. Thus, from Fig.1, one could tell that the OPPA
distinguishes the shortest path long before the termination condition is met.

Transition point (T-Point): The T-point of the PPA-based algorithms is defined
as the iteration when the $L_{D-Path}$ starts to converge to the length of the
shortest path. For example, the four points with green arrows pointing at in
Fig.1 are the T-Points of the OPPA in each of the four instances,
respectively.

Given the definition of T-Point, the main question proposed in Section I,
i.e., “When does the OPPA distinguish the shortest path from the other
paths?”, can be answered. We would like to provide our opinion: after the
transition point, the PPA-based algorithms start to distinguish the shortest
path from the others. For the first question in Section I, i.e., “In the case
of ensuring accuracy, is there an inherent limit, so that no matter how the
termination criteria changes, the OPPA cannot increase the convergence rate?”,
again, we would like to propose a possible answer: the inherent limit would be
the T-Point, before which the OPPA does not achieve the shortest path.

###  III-C The proposed algorithm: OPPA-D

Further developing from the previous sub-sections, we now introduce our
application of the proposed concepts, i.e., the OPPA-D algorithm.

The OPPA-D algorithm builds on the OPPA with a newly proposed termination
criterion. This criterion is defined as follows: Execute the OPPA until the
$L_{D-Path}$ does not change for $K$ iteration. $K$ is a pre-defined
parameter. Thus, after the $L_{D-Path}$ remains steady for $K$ iterations, the
occurrence of the T-Point is assumed, which also implies the convergence of
the OPPA.

The pseudo-code of the OPPA-D algorithm is provided in Alg.III-C. For purpose
of illustration, $K$ is set to $10$. The algorithm’s sensitivity to $K$ will
be subsequently tested in Section IV-B.

Algorithm 2 The OPPA-D

1: //Initialization part

2: Input: the statistics of graph $G(N,E)$, the corresponding matrix of the
weight of each edges W, and the pre-defined parameter $K$.

3: $D_{ij}\Leftarrow 0.5\ (\forall E_{ij}\in E,)$ // Initialize the
conductivity of each edge.

4: $Q_{ij}\Leftarrow 0\ (\forall E_{ij}\in E)$// Initialize the flow of each
edge.

5: $L_{ij}\Leftarrow W_{ij}(\forall E_{ij}\in E)$ //The length of the edges
equal to their weights.

6: $p_{i}\Leftarrow 0\ (\forall i=1,2,\cdots,N)$ // Initialize the pressure at
each node.

7: $count\Leftarrow 1$ // Initialize the counting variable.

8: $Length_{D-Path}^{0}\Leftarrow 0$ // Variable used to record the length of
D-Path.

9: $COUNT\Leftarrow 0$ // Variable used to count the number of iterations that
the $Length_{D-Path}$ remains unchanged.

10: //Iterative part

11: while true do

12: $p_{1}\Leftarrow 0$ // The pressure of the starting node (Node $1$) is set
to 0.

13: Calculate the pressure of all the nodes using (5):  $\sum_{i\in
N}{\frac{D_{ij}}{L_{ij}}(p_{i}-p_{j})=\left\\{\begin{gathered}+IN_{0}\ \ for\
j=1,\\\ -IN_{0}\ \ for\ j=2,\\\ 0\ \ \ otherwise.\\\ \end{gathered}\right.}$

14: Calculate the flux using (1):  $Q_{ij}\Leftarrow
D_{ij}\cdot(p_{i}-p_{j})/L_{ij}.$

15: Calculate the conductivity of the iteration using (8):
$D_{ij}^{n+1}=\frac{|Q_{ij}^{n}|+D_{ij}^{n}}{2}.$

16: Find the D-Path and calculate its length $Length_{D-Path}^{count}$
according to Section III-A.

17: if $Length_{D-Path}^{count}$==$Length_{D-Path}^{count-1}$ then

18: $COUNT=COUNT+1$

19: if $COUNT\geq K$ then

20: break

21: end if

22: else

23: $COUNT=0$

24: end if

25: $count\Leftarrow count+1$

26: end while

27: Output: The shortest path length $Length_{D-Path}^{count}$.

##  IV Experiments & Analysis

In this section, we present our experimental methodologies and our findings.
In particular, we compare the proposed OPPA-D with other state-of-the-art
accelerated OPPAs in randomly generated complete networks, randomly generated
small-world networks, and real-world networks, to demonstrate its
effectiveness and efficiency in Section IV-A. Since the OPPA-D is a modified
algorithm of the OPPA with the newly proposed termination criterion, we will
verify the validation of the proposed stopping criterion. Hence, the OPPA-D is
further compared with other commonly adopted criteria in Section IV-B. In
Section IV-C, three accelerated OPPAs will be evaluated using a new
transition-point based evaluation method. All the experiences are done through
computer simulations using Matlab R2018a on an Intel Core i5-8500 CPU (3GHz)
with 8 GB RAM under Windows 10.

###  IV-A Comparing the OPPA-D with other accelerated OPPAs

Data Set 1: There are three data sets. 1) ‘Da-Com-1’ to ‘Da-Com-8’: Eight
complete graphs are randomly generated with size ranging from 100 to 5000, and
length of the edges varying from 1 to 10000. 2) ‘Da-SW-1’ to ‘Da-SW-8’: Eight
small-world networks are randomly generated with size ranging from 100 to
5000, and length of the edges varying from 1 to 10000. All the small-world
networks are created using a Matlab function ’WattsStrogatz’ 111More details
of the function ’WattsStrogatz’ could refer to
https://ww2.mathworks.cn/help/matlab/math/build-watts-strogatz-small-world-
graph-model.html?lang=en.. 3) ‘Da-R-1’ and ‘Da-R-2’: Two real-world networks,
i.e., the Sioux-Falls network and the Anaheim network 222The meta-data of
networks could refer to https://github.com/bstabler/TransportationNetworks.,
are adopted. Some basic information about the data sets is provided in Table
I.

TABLE I: Basic information of Data Set 1

Instance Name | $|N|^{a}$ | $|E|^{a}$ | Instance Name | $|N|^{a}$ | Mean node degree | Instance Name | $|N|^{a}$ | $|E|^{a}$
---|---|---|---|---|---|---|---|---
Da-Com-1 | 50 | 1.22E+03 | Da-SW-1 | 50 | 6 | Da-R-1b | 24 | 76
Da-Com-2 | 100 | 4.95E+03 | Da-SW-2 | 100 | 12 | Da-R-2b | 416 | 914
Da-Com-3 | 250 | 3.11E+04 | Da-SW-3 | 250 | 30 |  |  |
Da-Com-4 | 500 | 1.25E+05 | Da-SW-4 | 500 | 60 |  |  |
Da-Com-5 | 750 | 2.81E+05 | Da-SW-5 | 750 | 90 |  |  |
Da-Com-6 | 1000 | 4.99E+05 | Da-SW-6 | 1000 | 120 |  |  |
Da-Com-7 | 2000 | 2.00E+06 | Da-SW-7 | 2000 | 240 |  |  |
Da-Com-8 | 5000 | 1.25E+07 | Da-SW-8 | 5000 | 600 |  |  |

a$N$ is the number of nodes and $E$ is the number of edges.
bInstances ‘Da-R-1’ and ‘Da-R-2’are the Sioux-Falls network and the Anaheim
network, respectively.

Baseline algorithms: The OPPAs with termination criteria $\epsilon\leq
10^{-2}$ and $\epsilon\leq 10^{-5}$, i.e., the OPPA($10^{-2}$) and
OPPA($10^{-5}$), are utilized as the baseline algorithms. Two other state-of-
the-art accelerated OPPAs, i.e., the EHPA[30] and the APS[33] are also
compared against 333The codes of the EHPA and the APS are available at
https://github.com/caigaoub/PhysarumOptimization. We thank the authors for
giving access to their source codes.. Each algorithm will be evaluated 15
times for each graph.

(a) ‘Da-Com-1’ instance . (b) ‘Da-Com-2’ instance . (c) ‘Da-Com-3’ instance .

(d) ‘Da-Com-4’ instance . (e) ‘Da-Com-5’ instance . (f) ‘Da-Com-6’ instance .

(g) ‘Da-Com-7’ instance . (h) ‘Da-Com-8’ instance .

Figure 2: Experimental results in complete graphs.

(a) ‘Da-Sw-1’ instance . (b) ‘Da-Sw-2’ instance . (c) ‘Da-Sw-3’ instance .

(d) ‘Da-Com-4’ instance . (e) ‘Da-Sw-5’ instance . (f) ‘Da-Sw-6’ instance .

(g) ‘Da-Sw-7’ instance . (h) ‘Da-Sw-8’ instance .

Figure 3: Experimental results in small-world networks.

(a) ‘Da-R-1’ instance . (b) ‘Da-R-2’ instance .

Figure 4: Experimental results in real-world networks. TABLE II: Average
running time of tested algorithms in Data Set 1

Instacnes | Running time (Second)
---|---
OPPA($10^{-2}$) | OPPA($10^{-5}$) | EHPA | APS | OPPA-D(Ours)
Da-Com-1 | 0.0020 | 0.0049 | 0.0073 | 0.0023 | 0.0009
Da-Com-2 | 0.0223 | 0.0557 | 0.0131 |  0.0059 | 0.0066
Da-Com-3 | 0.0400 | 0.0425 | 0.0474 | 0.0478 | 0.0398
Da-Com-4 | 0.4839 | 1.0926 | 0.3771 | 0.2293 | 0.2006
Da-Com-5 | 2.7072 | 6.1812 | 1.0617 | 0.6322 | 0.4698
Da-Com-6 | 3.5256 | 8.2404 | 2.8248 | 1.3938 | 0.7721
Da-Com-7 | 5.7222 | 9.4358 | 20.0010 | 16.0612 | 3.4003
Da-Com-8 | 50.1900 | 89.3060 | 444.2410 | 277.3507 | 27.6500
Da-SW-1 | 0.0015 | 0.0035 | 0.0010 | 0.0021 | 0.0008
Da-SW-2 | 0.0086 | 0.0195 | 0.0062 | 0.0066 | 0.0125
Da-SW-3 | 0.0357 | 0.0777 | 0.0466 | 0.0362 | 0.0217
Da-SW-4 | 0.2700 | 0.5818 | 0.2176 | 0.2157 | 0.0883
Da-SW-5 | 1.1505 | 2.7874 | 0.8128 | 0.5841 | 0.3672
Da-SW-6 | 0.9320 | 1.4951 | 1.1180 | 1.1162 | 0.5987
Da-SW-7 | 20.1512 | 53.4956 | 27.4835 | 15.3748 | 3.8334
Da-SW-8 | 121.5098 | 277.3977 | 620.7715 | 279.5441 | 34.8306
Da-R-1 | 0.0012 | 0.0021 | 0.0028 | 0.0039 | 0.0004
Da-R-2 | 0.8085 | 1.8604 | 0.3038 | 0.2507 | 0.1533

aValues presented in the table are the average value of 15 trials.

TABLE III: Average number of iterations of tested algorithms in Data Set 1

Instacnes | Number of iterations
---|---
OPPA($10^{-2}$) | OPPA($10^{-5}$) | EHPA | APS | OPPA-D(Ours)
Da-Com-1 | 26 | 58 | 12 | 16 | 12
Da-Com-2 | 57 | 134 | 21 | 16 | 19
Da-Com-3 | 27 | 49 | 16 | 16 | 18
Da-Com-4 | 58 | 139 | 23 | 17 | 25
Da-Com-5 | 134 | 337 | 26 | 26 | 25
Da-Com-6 | 96 | 238 | 33 | 20 | 22
Da-Com-7 | 31 | 55 | 19 | 16 | 19
Da-Com-8 | 41 | 75 | 26 | 17 | 23
Da-SW-1 | 24 | 47 | 8 | 16 | 11
Da-SW-2 | 27 | 50 | 14 | 16 | 11
Da-SW-3 | 41 | 89 | 20 | 16 | 17
Da-SW-4 | 34 | 73 | 16 | 16 | 11
Da-SW-5 | 62 | 149 | 25 | 18 | 20
Da-SW-6 | 27 | 43 | 16 | 16 | 17
Da-SW-7 | 116 | 305 | 29 | 16 | 22
Da-SW-8 | 105 | 233 | 38 | 36 | 30
Da-R-1 | 27 | 53 | 9 | 16 | 11
Da-R-2 | 133 | 330 | 37 | 44 | 24

aValues presented in the table are the average value of 15 trials.

All the evaluated algorithms successfully solved for the shortest path in all
15 instances. The running time of them are presented in Fig.2, 3, 4. Tables II
and III summarize the running time and the number of iterations, respectively.

In general, the OPPA-D is superior to the baseline algorithms in terms of
efficiency. In the complete graphs, as recorded in Fig.2 and Table II, the
OPPA-D outperforms all other baseline algorithms except for ‘Da-Com-2’ where
the APS is slightly faster than the OPPA-D. In the other instances, the OPPA-D
is 56.40%, 0.37%, 12.48%, 25.69%, 44.61%, 40.58%, and 44.91% faster than the
next fastest algorithm, respectively. In the small-world graphs, according to
Fig.3 and Table II, the OPPA-D is more efficient than all the baseline
algorithms except for ‘Da-Sw-2’ where it ranks fourth. In other instances, the
OPPA-D is 26.21%, 39.19%, 59.06%, 37.14%, 35.76%, 75.07%, and 71.34% faster
than the second-ranked algorithm, respectively. In the real-world network,
according to Fig.4 and Table II, the OPPA-D is the fastest algorithm. It is
63.09% and 38.84% more efficient than the algorithm in second place,
respectively.

Though the accelerated OPPAs, i.e., the EHPA and the APS, are faster than the
OPPA($10^{-2}$) and OPPA($10^{-5}$) in most instances, they do not out-perform
the OPPA($\epsilon=10^{-5}$) in the 2000-nodes and 5000-node complete graphs
and 5000-node small-world network. One interesting finding in this experiment
is that, according to Tables II and III, while the number of iterations is
small it does not necessarily translate to a shorter running time. For
example, the OPPA($10^{-2}$) takes 41 iterations to converge while the APS
only take 17 iterations for ‘Da-Com-8’. However, the OPPA($10^{-2}$) is
approximately five times faster than the APS. Therefore, when comparing
different accelerated OPPAs, the running time, instead of the number of
iterations, would be a more appropriate evaluation metric.

###  IV-B Comparing different termination criteria

In the above experiments, we observe that the running time of the OPPA is
highly sensitive to the current stopping criteria, i.e., $\sum_{i}\sum_{j\neq
i}D_{ij}\leq\epsilon$. According to Table II, the OPPA($10^{-2}$) is more
efficient than the OPPA($10^{-5}$); both algorithms solve for the shortest-
path in all the instances successfully. Recall that the OPPA-D is
fundamentally the OPPA with a newly proposed termination criterion. Thus, to
validate the superiority of the proposed criterion, i.e. executing the OPPA
until the $L_{D-Path}$ reamins steady for $K$ iteration, we conduct further
experiments.

The evaluated criteria: 1) The traditional termination criterion with
$\epsilon=10^{-1},10^{-2},\cdots,10^{-5}$ are tested. 2) The proposed
termination criterion with $K=5,10,15,20,30$ are also tested.

Data Set 2: We test the termination criteria in randomly generated complete
graphs. The size of the graphs are 10, 100, and 500. For each size, 50 graphs
are randomly generated with the length of the edge varying from 1 to 10000.

TABLE IV: Experimental results in Data Set 2–Success rate.

Size | The traditional criteria ($\epsilon$) | The proposed criteria ($K$)
---|---|---
($10^{-1}$) | ($10^{-2}$) | ($10^{-3}$) | ($10^{-4}$) | ($10^{-5}$) | 5 | 10 | 15 | 20 | 30
10 | 96%(2$|$0) | 100%(0$|$0) | 100%(0$|$0) | 100%(0$|$0) | 100%(0$|$0) | 96%(2$|$0) | 96%(2$|$0) | 100%(0$|$0) | 100%(0$|$0) | 100%(0$|$0)
100 | 98%(1$|$0) | 100%(0$|$0) | 96%(0$|$2) | 92%(0$|$4) | 88%(0$|$6) | 60%(20$|$0) | 88%(6$|$0) | 90%(5$|$0) | 96%(2$|$0) | 100%(0$|$0)
500 | 100%(0$|$0) | 100%(0$|$0) | 98%(0$|$1) | 96%(0$|$2) | 96%(0$|$2) | 64%(18$|$0) | 92%(4$|$0) | 98%(1$|$0) | 100%(0$|$0) | 100%(0$|$0)

a The values in the table have a format ‘Value1%(Value2$|$Value3)’, where
Value1 is the success rate out of 50 cases, Value2 represents the number of
the cases that the algorithm does not successfully solve, Value3 is the number
of cases where the algorithm could not reach the termination criteria though
sufficient computational time.

TABLE V: Experimental results in Data Set 2–Running time.

Size | The traditional criteria ($\epsilon$) | The proposed criteria ($K$)
---|---|---
($10^{-1}$) | ($10^{-2}$) | ($10^{-3}$) | ($10^{-4}$) | ($10^{-5}$) | 5 | 10 | 15 | 20 | 30
10 | 0.0003 | 0.0013 | 0.0045 | 0.0078 | 0.0111 | 0.0002 | 0.0003 | 0.0003 | 0.0004 | 0.0006
100 | 0.0130 | 0.0410 | 0.0913 | 0.0978 | 0.1061 | 0.0043 | 0.0065 | 0.0084 | 0.0105 | 0.0141
500 | 0.3203 | 0.9352 | 1.4508 | 1.7348 | 2.4032 | 0.1049 | 0.1441 | 0.1806 | 0.2179 | 0.2877

a The values presented in the table are the average values of running time
(seconds) of only the success cases from the 50 cases.

Tables IV and V report the experimental results. From Table IV, different
settings towards parameters $\epsilon$ and $K$ may lead to different success
rates. One interesting thing of note is that when the graph size is 100 or
500, as $\epsilon$ decreases, despite the decrease in the number of failed
attempts, the algorithm could not reach the traditional termination criterion
with a low $\epsilon$ in certain cases. However, the proposed stopping
criterion does not have the above issue.

According to Table V, one could conclude that, from the perpective of running
time, the traditional termination criterion is more sensitive to the parameter
$\epsilon$, compared with the proposed stopping criterion’s sensitivity to
$K$. For example, when the graph size is 500, when $\epsilon$ decreases from
$10^{-1}$ to $10^{-5}$, the running time increases from 0.32 seconds to 2.40
seconds (by 650%); however, when $K$ increases from $5$ to $30$, the running
time increases from 0.10 seconds to 0.29 seconds (by 190%).

In general, in terms of of success rate, the traditional termination criterion
with $\epsilon=10^{-2}$ and the proposed termination criterion with $K=30$
would be good choices, both of which achieve 100% success rate for three
different graph sizes (each of the graph size has 50 cases). After comparing
our experimental results of the running time for both criteria, one would find
that the proposed criterion with $K=30$ would be the better choice.

###  IV-C Analyzing the improvements made by accelerated OPPAs

In the previous sub-section, we demonstrate that the proposed criterion is
better than the traditional criterion. However, the running time of the OPPA
for both stopping criteria is still inevitably dependent on the initial
parameters. It is logical to deduce that the above findings still holds when
these criteria are applied to other accelerated OPPAs. Thus, in this section,
we want to compare the accelerated OPPAs with the OPPA by excluding the
termination criteria.

####  IV-C1 The proposed evaluation method

Using the defined concepts of the D-Path and T-Point, we propose our
evaluation method:

  1. i.

When generating the data set, calculate the shortest path of each graph.

  2. ii.

Evaluate the algorithm in the data set. In each iteration of the algorithm,
calculate the $L_{D-Path}$. If the $L_{D-Path}$ of this iteration is equal to
the shortest path length, record the corresponding running time and proceed to
Step iii; otherwise, repeat Step ii.

  3. iii.

Execute the algorithm for 50 or more iterations to ensure the convergence of
the algorithm. If the $L_{D-Path}$ maintains the value of the shortest path
length for the whole evaluation period, one could deduce that the T-Point has
already occurred. Thus, the recorded running time in Step ii is the running
time used for comparison; otherwise, return to Step ii.

####  IV-C2 Comparing different accelerated PPAs using the proposed method

Data Set 3: The instances ‘Da-Com-1’ to ‘Da-Com-8’ mentioned in Data Set 1 are
adopted.

Tested algorithms: The OPPA, the OPPA-D, the EHPA and the APS are evaluated.
When tested, each algorithm would be run 30 times in each graph.

TABLE VI: Average running time of tested algorithms in Data Set 3

Instacnes | Running time (seconds)
---|---
OPPA | EHPA | APS | OPPA-D (Ours)
Da-Com-1 | 0.0006 | 0.0025 | 0.0029 | 0.0006
Da-Com-2 | 0.0068 | 0.0038 | 0.0065 | 0.0046
Da-Com-3 | 0.0017 | 0.0040 | 0.0346 | 0.0020
Da-Com-4 | 0.1194 | 0.2675 | 0.2061 | 0.1213
Da-Com-5 | 0.2778 | 0.6585 | 0.5086 | 0.2833
Da-Com-6 | 0.4198 | 1.1887 | 1.0730 | 0.4267
Da-Com-7 | 1.5913 | 11.5127 | 15.3770 | 1.5929
Da-Com-8 | 15.2210 | 216.5146 | 262.6555 | 15.3558

aValues presented in the table are the average value of 30 trials.

TABLE VII: Average number of iterations of tested algorithms in Data Set 3

Instacnes | Number of iterations
---|---
OPPA | EHPA | APS | OPPA-D (Ours)
Da-Com-1 | 2 | 2 | 16 | 2
Da-Com-2 | 9 | 9 | 16 | 9
Da-Com-3 | 1 | 1 | 16 | 1
Da-Com-4 | 15 | 15 | 16 | 15
Da-Com-5 | 15 | 15 | 16 | 15
Da-Com-6 | 12 | 12 | 16 | 12
Da-Com-7 | 9 | 9 | 16 | 9
Da-Com-8 | 13 | 13 | 16 | 13

aValues presented in the table are the average value of 30 trials.

Tables VI and VII presents the experimental results. It is not a surprise that
the running time and the number of iterations of the OPPA-D and the OPPA are
approximately the same, as the difference between the OPPA and the OPPA-D lies
only in the termination criterion. The similar performance of the OPPA and the
OPPA-D under the proposed evaluation method proves that the method developed
by us could exclude the effects brought by different termination criteria.

Though the EHPA share the same number of iterations with the OPPA and the
OPPA-D, the running time of EHPA is approximately 14 times longer than that of
the OPPA and the OPPA-D in instance ‘Da-Com-8’. The above findings indicate
that when the effects of the termination criterion is excluded, the EHPA fails
to accelerate the PPA as expected; instead, the running time per iteration is
significantly longer.

According to Table VI, on average, the running time of the APS is even longer
than the EHPA. Furthermore, it is surprising that, in all instances, the
number of the iterations of the APS is the same, which is 16. We carefully
examined the codes provided by the authors [33] and found out that at the
beginning of the iterative procedure of the APS, the algorithm has a
15-iteration initialization run. One could easily deduce that the APS only
uses one iteration to reach the T-Point after the initialization. Hence, it
suggests that the initialization of the APS costs the majority of the running
time of the algorithm, and further leads to the inefficiency of the APS.

In general, the proposed methods do have the ability to avoid the effects
brought by different stopping criteria, as the performance of the OPPA and the
OPPA-D is approximately the same in the experiment. The most important finding
in this experiment is that the accelerated OPPAs, i.e., the EHPA and the APS
fail to accelerate the PPA as expected. A potential reason might be the
failure of the authors to not consider the possible effects brought about by
the termination criteria.

##  V Discussions & Conclusions

Though many studies have proposed various solutions to accelerate the PPA, few
have considered the core question: when does the physarum solver distinguish
the shortest path from other paths? By defining the concepts of the T-Point
and the D-Path in Section III, we showed that the physarum solver
distinguishes the shortest path at the T-Point. We also propose the use of the
T-Point as the inherent limit, beyond which the algorithm can terminate, hence
improving efficiency while maintaining accuracy. The mathematical proof
remains in the pipelines for further work. We hope that in the future, the
phenomenon of T-Point could also be experimentally confirmed in biological
experiments.

Based on the T-Point and the D-Path, a new termination criterion for the OPPA
and a corresponding algorithm, OPPA-D, are developed. The OPPA-D is actually
the OPPA with the newly proposed termination criterion. From our experiments
in Section IV-A, the OPPA-D outperforms two state-of-the-art and widely used
accelerated OPPAs. In order words, the OPPA significantly outperforms the
accelerated OPPAs despite a slight modification to the stopping criterion.
This begs another question: does the accelerated OPPAs beat the OPPA as
claimed?

In section IV-B, we found that although the proposed termination criterion is
superior to the traditional one, both of them are sensitive to the predefined
parameters, which would further affect the validity of comparing different
accelerated OPPAs. Hence, in order to evaluate the efficiency of accelerated
OPPAs, the effects of the termination criteria need to be excluded. Thus, we
proposed a new evaluation method to compare the accelerated OPPAs more
objectively in Section IV-C1. From our results in section IV-C2 one finds that
the accelerated OPPAs is not as efficient as the OPPA, as previously reported.
Previous research did not evaluate the real performance of the accelerated
OPPA by neglecting the possible effects brought by the termination criteria.
However, with the proposed T-Point and the D-Path, and the developed
evaluation method, experiments in our studies demonstrate that the OPPA is the
fastest algorithm among the tested algorithms. We hope the findings in this
paper would change the standing critique that OPPAs are not efficient.

## Acknowledgment

The work is partially supported by National Natural Science Foundation of
China (Grant No. 61973332).

## References

  * [1] R. Storn and K. Price, “Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces,” _Journal of global optimization_ , vol. 11, no. 4, pp. 341–359, 1997.
  * [2] S. Mirjalili, S. M. Mirjalili, and A. Lewis, “Grey wolf optimizer,” _Advances in engineering software_ , vol. 69, pp. 46–61, 2014.
  * [3] D. Karaboga and B. Basturk, “A powerful and efficient algorithm for numerical function optimization: artificial bee colony (abc) algorithm,” _Journal of global optimization_ , vol. 39, no. 3, pp. 459–471, 2007.
  * [4] J. H. Holland, “Genetic algorithms,” _Scientific american_ , vol. 267, no. 1, pp. 66–73, 1992.
  * [5] M. Dorigo, M. Birattari, and T. Stutzle, “Ant colony optimization,” _IEEE computational intelligence magazine_ , vol. 1, no. 4, pp. 28–39, 2006\.
  * [6] S. Mirjalili and A. Lewis, “The whale optimization algorithm,” _Advances in engineering software_ , vol. 95, pp. 51–67, 2016.
  * [7] L. Zhang, H. Pan, Y. Su, X. Zhang, and Y. Niu, “A mixed representation-based multiobjective evolutionary algorithm for overlapping community detection,” _IEEE Transactions on Cybernetics_ , vol. 47, no. 9, pp. 2703–2716, 2017\.
  * [8] Y. Tian, X. Zheng, X. Zhang, and Y. Jin, “Efficient large-scale multiobjective optimization based on a competitive swarm optimizer,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 8, pp. 3696–3708, 2020.
  * [9] Y. Huang, Y. Gao, Y. Gan, and M. Ye, “A new financial data forecasting model using genetic algorithm and long short-term memory network,” _Neurocomputing_ , 2020.
  * [10] A. W. Mohamed, A. A. Hadi, and K. M. Jambi, “Novel mutation strategy for enhancing shade and lshade algorithms for global numerical optimization,” _Swarm and Evolutionary Computation_ , vol. 50, p. 100455, 2019.
  * [11] P. Hu, J.-S. Pan, and S.-C. Chu, “Improved binary grey wolf optimizer and its application for feature selection,” _Knowledge-Based Systems_ , p. 105746, 2020.
  * [12] M. Wang and H. Chen, “Chaotic multi-swarm whale optimizer boosted support vector machine for medical diagnosis,” _Applied Soft Computing_ , vol. 88, p. 105946, 2020.
  * [13] 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_ , vol. 327, no. 5964, pp. 439–442, 2010.
  * [14] T. Nakagaki, M. Iima, T. Ueda, Y. Nishiura, T. Saigusa, A. Tero, R. Kobayashi, and K. Showalter, “Minimum-risk path finding by an adaptive amoebal network,” _Physical review letters_ , vol. 99, no. 6, p. 068104, 2007.
  * [15] T. Nakagaki, H. Yamada, and Á. Tóth, “Maze-solving by an amoeboid organism,” _Nature_ , vol. 407, no. 6803, pp. 470–470, 2000.
  * [16] T. Nakagaki, “Smart behavior of true slime mold in a labyrinth,” _Research in Microbiology_ , vol. 152, no. 9, pp. 767–770, 2001.
  * [17] T. Nakagaki, H. Yamada, and A. Toth, “Path finding by tube morphogenesis in an amoeboid organism,” _Biophysical chemistry_ , vol. 92, no. 1-2, pp. 47–52, 2001.
  * [18] X. Zhang, Q. Wang, A. Adamatzky, F. T. Chan, S. Mahadevan, and Y. Deng, “An improved physarum polycephalum algorithm for the shortest path problem,” _The Scientific World Journal_ , vol. 2014, 2014.
  * [19] V. Bonifaci, “Physarum can compute shortest paths: A short proof,” _Information Processing Letters_ , vol. 113, no. 1-2, pp. 4–7, 2013.
  * [20] L. Becchetti, V. Bonifaci, M. Dirnberger, A. Karrenbauer, and K. Mehlhorn, “Physarum can compute shortest paths: Convergence proofs and complexity bounds,” in _International Colloquium on Automata, Languages, and Programming_. Springer, 2013, pp. 472–483.
  * [21] A. Karrenbauer, P. Kolev, and K. Mehlhorn, “Convergence of the non-uniform physarum dynamics,” _Theoretical Computer Science_ , 2020.
  * [22] Y. Liu, C. Gao, Z. Zhang, Y. Lu, S. Chen, M. Liang, and L. Tao, “Solving np-hard problems with physarum-based ant colony system,” _IEEE/ACM transactions on computational biology and bioinformatics_ , vol. 14, no. 1, pp. 108–120, 2015.
  * [23] X. Chen, Y. Liu, X. Li, Z. Wang, S. Wang, and C. Gao, “A new evolutionary multiobjective model for traveling salesman problem,” _IEEE Access_ , vol. 7, pp. 66 964–66 979, 2019.
  * [24] Y. Sun, P. N. Hameed, K. Verspoor, and S. Halgamuge, “A physarum-inspired prize-collecting steiner tree approach to identify subnetworks for drug repositioning,” _BMC systems biology_ , vol. 10, no. 5, p. 128, 2016.
  * [25] S. Xu, W. Jiang, X. Deng, and Y. Shou, “A modified physarum-inspired model for the user equilibrium traffic assignment problem,” _Applied Mathematical Modelling_ , vol. 55, pp. 340–353, 2018.
  * [26] M.-A. I. Tsompanas, G. C. Sirakoulis, and A. I. Adamatzky, “Evolving transport networks with cellular automata models inspired by slime mould,” _IEEE Transactions on Cybernetics_ , vol. 45, no. 9, pp. 1887–1899, 2015\.
  * [27] X. Zhang and S. Mahadevan, “A bio-inspired approach to traffic network equilibrium assignment problem,” _IEEE Transactions on Cybernetics_ , vol. 48, no. 4, pp. 1304–1315, 2018.
  * [28] H. Jiang, X. Liu, S. Xiao, C. Tang, and W. Chen, “Physarum-inspired autonomous optimized routing protocol for coal mine manet,” _Wireless Communications and Mobile Computing_ , vol. 2020, 2020.
  * [29] Y. Song, L. Liu, H. Ma, and A. V. Vasilakos, “Physarum optimization: a new heuristic algorithm to minimal exposure problem,” in _Proceedings of the 18th annual international conference on Mobile computing and networking_ , 2012, pp. 419–422.
  * [30] X. Zhang, Y. Zhang, Z. Zhang, S. Mahadevan, A. Adamatzky, and Y. Deng, “Rapid physarum algorithm for shortest path problem,” _Applied Soft Computing_ , vol. 23, pp. 19–26, 2014.
  * [31] Q. Wang, X. Lu, X. Zhang, Y. Deng, and C. Xiao, “An anticipation mechanism for the shortest path problem based on physarum polycephalum,” _International Journal of General Systems_ , vol. 44, no. 3, pp. 326–340, 2015.
  * [32] Q. Cai and Y. Deng, “A fast bayesian iterative rule in amoeba algorithm.” _International Journal of Unconventional Computing_ , vol. 14, 2019.
  * [33] C. Gao, X. Zhang, Z. Yue, and D. Wei, “An accelerated physarum solver for network optimization,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 2, pp. 765–776, 2018.
  * [34] X. Zhang, A. Adamatzky, X.-S. Yang, H. Yang, S. Mahadevan, and Y. Deng, “A physarum-inspired approach to supply chain network design,” _Science China Information Sciences_ , vol. 59, no. 5, p. 052203, 2016.