Transportation Network with Fluctuating Input/Output Designed by the Bio-Inspired Physarum Algorithm

In this paper, we propose designing transportation network topology and traffic distribution under fluctuating conditions using a bio-inspired algorithm. The algorithm is inspired by the adaptive behavior observed in an amoeba-like organism, plasmodial slime mold, more formally known as plasmodium of Physarum plycephalum. This organism forms a transportation network to distribute its protoplasm, the fluidic contents of its cell, throughout its large cell body. In this process, the diameter of the transportation tubes adapts to the flux of the protoplasm. The Physarum algorithm, which mimics this adaptive behavior, has been widely applied to complex problems, such as maze solving and designing the topology of railroad grids, under static conditions. However, in most situations, environmental conditions fluctuate; for example, in power grids, the consumption of electric power shows daily, weekly, and annual periodicity depending on the lifestyles or the business needs of the individual consumers. This paper studies the design of network topology and traffic distribution with oscillatory input and output traffic flows. The network topology proposed by the Physarum algorithm is controlled by a parameter of the adaptation process of the tubes. We observe various rich topologies such as complete mesh, partial mesh, Y-shaped, and V-shaped networks depending on this adaptation parameter and evaluate them on the basis of three performance functions: loss, cost, and vulnerability. Our results indicate that consideration of the oscillatory conditions and the phase-lags in the multiple outputs of the network is important: The building and/or maintenance cost of the network can be reduced by introducing the oscillating condition, and when the phase-lag among the outputs is large, the transportation loss can also be reduced. We use stability analysis to reveal how the system exhibits various topologies depending on the parameter.


Effect of output frequency
The numerical calculations of this system were performed over the frequency range 2π × 10 −2 ≤ ω ≤ 2π × 10 7 . When 2π × 10 −1 ≤ ω ≤ 2π × 10 7 , the conclusions in the main text are valid, except for the oscillatory behavior of the converged D ij (t) detailed below. Interestingly, the network topology still depends on phase-lag ϕ in two oscillating outputs, even when the frequency is very high (ω = 2π × 10 7 ) against time constant of D-degeneration (estimated as 1 from Eq. (4)). This could be caused by nonlinearity of the growth function, Eq. (5), where the dependence of the time-average of the flux Q on ϕ survives in this system at least in this frequency range (see also Section 4 in File S1). Conversely, when the frequency is extremely low, such as ω = 2π ×10 −2 , the system can follow the slow change of flux. Then, the converged network topology differs from that for higher frequency, as shown in Fig. A in File S1. The system exhibits a V-shaped network when ϕ = 0 ( Fig. A(a)) and a Y-shaped network when ϕ = π (Fig. A(b)) in the higher frequency, which is consistent with the results in the main text. However, the system always exhibits a V-shaped network in the lower frequency ω = 2π × 10 −1 To be more specific about the oscillatory behavior of converged D ij (t), the mean value (equivalent toD ij ) does not depend on frequency but the amplitude does, as illustrated in Fig Left panels show time-variable D ij . Right panels represent the converged network topologies. The colors of the plots correspond to those in the network diagrams. Bar represents the amplitude of oscillation in D ij (t). The data are obtained for the links colored red in the inset diagram. The phase lag was set at ϕ = 0.

The converged valueD ij versus time-variable
D ij (t) on estimation of performance functions.
In the numerical calculations, converged valueD ij is used instead of timevariable D ij (t) for application to man-made systems such as practical power grids because wires with time-variable conductance are not realistic. However, for application to natural systems such as slime mold and ant trails, the time-variable conductance is possible. We tested the case of D ij (t) instead ofD ij in the calculation of Eqs. (9) and (10) using the following definition, and obtained virtually the same result. The relative residuals between the two calculation methods were a maximum of |P t − P |/P = 5.6 × 10 −5 for loss P , and a maximum of |B t − B|/B = 2.1 × 10 −3 for cost B. These are sufficiently small that the results in the main text are still valid for the case of the time-variable conductance.

Comparison of performance functions starting from homogeneous and non-homogeneous initial conditions.
Performance functions are calculated for the systems starting from the homogeneous and non-homogeneous initial conditions of D ij . The results in the case of the non-homogeneous condition hold almost the same feature as in the case of the homogeneous condition, as evidenced in Figs. C-E. A small difference is observed only in P . The non-homogeneity causes multiple values of P at a certain combination of µ and ϕ, for example, multiple crosses at µ = 3.5, 4.0, 4.5, 5.0 in Fig. C(a). The multiplicity is derived from the difference in the network topology. Color of plots represents network topology: Black, red, blue, and green denote, respectively, mesh, partial mesh, V-shaped, and Y-shaped networks.
Circles are results starting from the homogeneous initial conditions of D ij . Crosses are results starting from the non-homogeneous initial conditions, where D ij are distributed according to a normal distribution with mean 1.0 and standard deviation 0.1. One hundred samples for each parameter set of µ and ϕ were tested for the non-homogeneous conditions.

Dependence of integration of f (|Q|) on ϕ
In the calculation of Eq. (13) in the main text, the time-averaged value of f (|Q|) over a period depends on ϕ because the function is nonlinear. The flux Q i at each link l i is the sum of those originating from out 1 and out 2 , which are respectively denoted as q out 1 i = a i (1 + sin ωt) and q out 2 The parameters a i and b i are certain constants calculated for each parameter condition of µ and ϕ.
First, assume the linear case: f (|Q|) = |Q|, the integration over a period is as follows: which is not dependent on ϕ. Second, assume the simplest example for the nonlinear case: f (|Q|) = |Q| 2 , the integration over a period is as follows: (S1.4) Now the integration depends on ϕ.