The role of noise in self-organized decision making by the true slime mold Physarum polycephalum

In principle, Markov theory offers us powerful means to further analyze this process. However, it is impossible (or at least extremely difficult) to apply such an analysis to the two-dimensional system (Eqs (9) and (10)) with a discontinuous forcing function. We thus aim to reduce it to a one-dimensional system that approximates the main properties of the full system well and is amenable to a formal analysis. A one-dimensional continuous-time continuous-space Markov process for the c value can be specified as an Itô-Diffusion dcdt=μ(c,t)+σ(c,t)ξ(t),(13) where μ describes the deterministic development (so-called drift), and ξ is a Gaussian noise |ξ(t)| = 1, with mean 〈ξ(t)〉 = 0, and uncorrelated in time 〈ξ(t)ξ(t′)〉 = δ(t − t′). σ captures the fluctuation of the noise amplitude.

For a temporally homogeneous process we can attempt to infer a one-dimensional Itô-Diffusion from experimental data or simulation data by a technique known as equation-free analysis (EFA [25]). The idea is the following: We assume the existence of μ(⋅) and σ(⋅) and measure them from simulation data of the full system Eq (6). We then compare the evolution of c(t) measured from this simulation data with the evolution of c(t) obtained by forward integration of Eq (13) for a large number of sample paths. If these agree statistically, we are justified in our choice of μ(⋅) and σ(⋅) and can proceed by analysis of Eq (13) to understand the properties of the full system. We use a variant of EFA [15] that estimates drift and variance parameters for a process X(t) as: μ(x)=⟨X(t+δt)-X(t)|X(t)=x⟩δt(14) σ2(x)=⟨[X(t+δt)-X(t)-μ(x)δt]2|X(t)=x⟩δt(15) where 〈⋅〉 denotes the sample average.

A complication is that the process under consideration is not temporally homogeneous, because of the forcing function, i.e. μ = μ(c, τ(t)) and σ2 = σ2(c, τ(t)) with τ(t) = t mod 2π. We thus divide the process into three different regimes corresponding to the three different forms of forcing: both paths dark; both paths lit; one path dark and one path lit. Note that the combination lit/dark only occurs in a single form with always the same path lit, depending on br1 and br2. Based on this we estimate the coefficients μi(⋅), σi(⋅) for each of the regimes i separately. Thus, assuming br1 < br2, each of the regimes by itself can be treated as a time-homogenous process. dcdt={μ1(c)+σ1(c)ξ(t),τ≤br1μ2(c)+σ2(c)ξ(t),br1<τ≤br2μ3(c)+σ3(c)ξ(t),τ>br2(16) Fig 5 shows the results of EFA carried out with a bin size of 0.04.

In the first regime both tubes are unlit so that there is no forcing (0 ≤ τ < π). As one would expect based on the underlying deterministic process (Fig 2), the shape of the estimated drift function μ1(c) implies two locally stable equilibria, one globally stable equilibrium at the downgoing zero c ≈ 0, and two further globally stable equilibria at the absorbing interval boundaries (c ≈ ±1).

In the second regime only the first tube is lit (π≤τ<32π), which makes the second one more attractive. Thus the drift is shifted negative. The magnitude of the drift is increased due to the larger difference in risks between tubes.

In the third regime the first tube is more attractive (32π≤τ<2π; both tubes lit, but the second one more strongly than the first one). The drift is shifted to the positive region with the exception of the range c < −0.7 where the absorption of the boundary at c = −1 dominates. The dominating stable equilibrium is, however, the one at the interval boundary c = 1. Diffusion depends only weakly on the forcing regime.

We simulate the one-dimensional Markov process Eq (13) with μ, σ as obtained by EFA. Averages and variance for 500 sample paths are given in Fig 4B and compared to the full two-dimensional system. The figure shows that there is very good agreement between the expected values of the two processes. The variance follows the same pattern for both processes and is lower for the reduced system. The good agreement indicates that the essentially features of the process are captured in the reduced system and we may proceed with further analysis based on the reduced system.

To this end, we construct a temporally-homogeneous Markov process by time-averaging the forcing regimes according to br1, br2. This is equivalent to a linear approximation of the influence of forcing within each period. A similar time-averaging could be achieved by performing an EFA assuming only a single (averaged) forcing regime for the whole process. Fig 6 gives the potential function Φ of this process [26]: Φ(c):=-∫-∞x2μ^(s)σ^2(s)ds(17) where μ^ and σ^2 are the time-averaged drift and variance. Local minima in Φ correspond to meta-stable points of the stochastic process and locally stable equilibria of the underlying deterministic process. It is clearly visible that the meta-stable point at c ≈ 0 has almost become a saddle due to the influence of noise. We are mainly interested in how likely it is that the system will evaluate the time-variant risk correctly and make a decision for Path 1.

This can be computed as the so-called splitting probability [24]. For a given interval [a, b], the splitting probability pa(x) gives the probability that the system initialized at c = x will reach the state c = a before the state c = b, i.e. in the case of a bi-stable process that it will make a decision for a. The corresponding probability pb(x) is defined symmetrically.

Let f(t, y) be the probability density for c(t) to take the value y at time t. The time-development of f(⋅, ⋅) is described by the Kolmogorov-forward or Fokker–Planck Equation (FPE [27]). ∂tf(t,y)=-∂y[μ^(y)f(t,y)]+∂yy[12σ^2(y)f(t,y)](18) Its steady state π(c) = f(c, t) is time-independent, so that Eq (18) reduces to the ODE 0=-ddy[μ^(y)π(y)]+12d2dy2[σ^2(y)π(y)](19) We can thus calculate the steady state probability density function π(x) as the solution of Eq (19) ψ(x)=e∫0x(2μ^(y)/σ^2(y))dyπ(x)=Cψ(x)σ^2(x)(20) where C is a suitable normalization constant [24]. Splitting probabilities for the process to leave the interval through the end at c = 1, i.e. with a correct decisions, can now be calculated for the PDF π(x) as [24, Eq (5.2.190)] pb(x)=∫axψ(s)ds(∫abψ(s)ds)-1(21) The splitting probabilities are shown in Fig 6. The decision point lies at c ≈ −0.5 where p+1 = p−1. Thus, the system will decide correctly with c → 1 for most initialization points in the range. The transition from a wrong decision to a correct one is remarkably sharp considering that outside of the interval [−0.75, −0.25] the residual probability for the system to revise its decision is less than 0.001. From the potential Φ we know that the c-space projection of the attractor for the meta-stable point of no decision D1 ≈ D2 is −0.5 < c < 0.5. The expected probability of the process to decide correctly when initialized at a random position in this range is ∫-0.50.5p+1(c)dc=0.969. In conclusion, the stochastic process will almost always decide correctly unless it is initialized with a very strong bias towards the wrong decision (c < −0.5). This is unlike the underlying deterministic process in which the decision depends almost entirely on the initialization and over a wide range of initializations no decision at all will be achieved.

Based on this analysis, we expect that the model can be tested with the proposed experiment in a straightforward fashion. We expect the success rate of the organism in the biological experiment to be significantly different from the one predicted by the noise-free model (ξ = 0), but we would expect these differences to disappear if ξ ≠ 0 is fitted to the data in a cross-validation procedure.

The mathematical analysis clearly shows that a well-attuned level of noise can enable the organism to correctly assess a time-variant risk, while the corresponding noise-free system fails to do so. This corroborates that noise plays a crucial functional role in self-organized systems. In biological terms, this is of evolutionary significance. Biological systems need to strike a delicate balance between flexible and stable behavior: Stability allows an organism or a group to concentrate its resources and to ignore irrelevant short-term fluctuations in the environment. Yet, adaptation in the case of stronger or more long lasting changes is required. A multi-stable behavior selection mechanism, such as the one analyzed here, can achieve exactly this by exploiting noise. It enables a self-organized system to reliably react to short-term changes in the environment while maintaining a generally stable behavior. The alternative of a control mechanism that follows every change in the environment would potentially be disadvantageous because it leads to unstable behavior.

Related findings have earlier been reported for ant colonies [8, 9]. There are two important differences between these studies and the present one. Firstly, the two biological systems investigated, ants and slime molds, are fundamentally different. Secondly, in the case of ants it was shown that noise enables them to react to changes in the environment by switching between multiple behaviors, whereas our study shows that noise can enable the slime mold to select the correct behavior in an environment that changes too frequently for it to follow individual changes. Instead of attempting to track each change, the organism can adopt a single behavior that maximizes the average long-term benefit. Yet, despite these differences both studies show that the decision making in ants and slime molds can be understood as instances of the same phenomenon: behavior selection as stochastic attractor switching [28].

In the case of ants, the effects have already been experimentally verified for the real biological system [9]. We have proposed a simple and concrete experiment to do this for P. polycephalum. Our mathematical analysis shows that this experiment will yield interesting outcomes whether it verifies our theoretical predictions or reveals that the model does not capture the full spectrum of decision making mechanisms in P. polycephalum.

Generally speaking, the fact that noise facilitates adaptive decision making is not tied to specific physical details of any particular biological system. Instead, it arises from very general mathematical properties of the underlying self-organized processes [10]. Mass recruiting ants and slime molds have very little in common biologically and physically. Yet, despite this, the phenomenological mathematical models that describe their behavior selection are very similar when constructed on the right level of observation. The same holds for a variety of other types of self-organized collective decision-making mechanisms in social organisms and human social systems. For example, food source selection [1] and clustering behavior [15] of honey bees, foraging patterns of bacteria [12], the emergence of fashion trends [13] and the dispersion of innovations [14] all can be and have been described with very similar mathematical models. This explains why some fundamental principles that govern self-organized collective behavior appear to be universal across the range [29, 30]. We may thus expect to also find similar beneficial effects of noise in other instances of self-organized decision making.

Noise has many origins. In any biological system two major influences are fluctuations in the environment and the stochastic nature of the underlying bio-chemical processes themselves. In the behavior of social groups, variations between individuals’ characteristics and the stochastic nature of interactions between group members provide additional sources of noise [9]. Pseudo-randomness in the form of deterministic chaos may also enter into the equation [15]. In fact, no real physical system is noise-free. Usually, however, we expect noise to be a disturbance and to degrade system performance or at best to be irrelevant. It is thus fascinating that evolution seems to have enabled some organisms to make constructive use of noise.

Consider the system ∂Di∂t=-δDi+f(DiD1+D2)(22a) with f(Qi) = (1 + ϵ)Q2/(ϵ + Q2) and Qi = Di/(D1 + D2) for i = {1, 2}. In an equilibrium, it is ∂tD1 = ∂tD2 = 0 and hence 1+ϵδ=ϵ(D1+D2)2+D12D1(22b) 1+ϵδ=ϵ(D1+D2)2+D22D2(22c) Substituting Eq (22b) into Eq (22c) delivers a sixth-order polynomial in D2, i.e. there are a maximum of six equilibria. The first three are easily determined 1)the trivial equlibrium with (δD1, δD2) = (0, 0)2)one equilibrium on the D2 axis: (δD1, δD2) = (0, 1)3)one equilibrium on the D1 axis: (δD1, δD2) = (1, 0)

Assuming that for the latter three D1 ≠ 0 and D2 ≠ 0, we may let D2 = aD1 where physical realizability implies a∈R+. With D2 = aD1, Eqs (22b) and (22c) read as δD1=f(D1D1(1+a))=f(11+a)=1+ϵϵ(1+a)2+1(23a) δaD1=f(aD1D1(1+a))=f(a1+a)=(1+ϵ)a2ϵ(1+a)2+a2(23b) eliminating D1 from these equations delivers the following relation for a which is independent of δ: (1+ϵ)aϵ(1+a)2+1=(1+ϵ)a2ϵ(1+a)2+a2(23c) ⇒0=ϵa3+(ϵ-1)a2-(ϵ-1)a-ϵ(23d) withtherootsa1=1anda2/3=1-2ϵ2ϵ±2ϵ-12ϵ-1.(23e) The criterion for a2/3 to be real is ϵ<14. Such that we find the three remaining equilibria 4)a = a1 = 1: (δD1,δD2)=(1+ϵ1+4ϵ,1+ϵ1+4ϵ)5)a = a2: (δD1,δD2)=(1+ϵϵ(1+a2)2+1,a2(1+ϵ)ϵ(1+a2)2+1)6)a = a3: (δD1,δD2)=(1+ϵϵ(1+a3)2+1,a3(1+ϵ)ϵ(1+a3)2+1)

The Jacobian of the System (22a) reads as J=(-δ+∂D1[f(Q1)]∂D2[f(Q1)]∂D1[f(Q2)]-δ+∂D2[f(Q2)])(24a) where we have (without summation over double-occurring indices) ∂Di[f(Qi)]=2ϵ(1+ϵ)D1D2(D1+D2)[Di2+ϵ(D1+D2)2]2(24b) ∂D1[f(Q2)]=2ϵ(1+ϵ)D22(D1+D2)[D22+ϵ(D1+D2)2]2(24c) ∂D2[f(Q1)]=2ϵ(1+ϵ)D12(D1+D2)[D12+ϵ(D1+D2)2]2(24d) such that 1δJ=(-1+2ϵδ(1+ϵ)D1D2(D1+D2)[D12+ϵ(D1+D2)2]22ϵδ(1+ϵ)D12(D1+D2)[D12+ϵ(D1+D2)2]22ϵδ(1+ϵ)D22(D1+D2)[D22+ϵ(D1+D2)2]2-1+2ϵδ(1+ϵ)D1D2(D1+D2)[D22+ϵ(D1+D2)2]2)(25)

For small perturbations, the equilibrium is unstable along the axes. Consider a small perturbation η > 0 along the direction of D1 and no perturbation along D2 such that Q1 = 1: ∂D1∂t=-δη+1>0(26) and we see that for small perturbations η < 1/δ, these perturbations continue to grow. The same holds for symmetry reasons along the D2 axis.

The normalized Jacobian of the system is 1δJ=(-12ϵδ(1+ϵ)1D1[1+ϵ]20-1)=(-1ϵ1+ϵ0-1)(27) with eigenvalues −1 and −1. The system is hence stable with respect to small perturbations in on of the two variables. For symmetry reasons, the same must hold for the equilibrium (δD1, δD2) = (0, 1).

For D1=D2=1+ϵ1+4ϵ, it is ∂Di[f(Qj)]=2ϵ(1+ϵ)2D3[(1+4ϵ)D2]2=δ4ϵ(1+ϵ)(1+4ϵ)21+4ϵ1+ϵ=δ4ϵ1+4ϵ≡δC(28) 1δJ=(C-1CCC-1)(29) with eigenvalues −1 and 2C − 1. The criterion for linear stability of the equilibrium is hence 2C-1<0⇒8ϵ<1+4ϵ⇒ϵ<0.25

In summary, we find that the bifurcation at ϵ=14 is a sub-critical pitchfork bifurcation (cf. Fig 7): For ϵ < 1/4 (Fig 2a), the system has three basins of attraction. The margins of these basins of attraction in the D1–D2 phase space are given by the lines D2 = a2D1 and D2 = a3D1. These lines cannot be crossed by trajectories since for a2 and a3, it is also ∂tD2 = a∂tD1. The values of a2,3 depend on the parameter ϵ only (they are independent of δ). In between these two lines, the equilibrium D1 = D2 = (1 + ϵ)/(1 + 4ϵ) is attracted; trajectories originating elsewhere in the parameter space converge towards an equilibrium on one of the two axes.For ϵ > 1/4 (Fig 2b), the system has two basins of attraction. The margin between these two basins is given by the line D1 = D2 (which, is not crossed by any trajectory) and initial values on either side of this line converge to the respective equilibrium on that side. The equilibrium with utilization of both tubes continues to exist as a saddle point but looses its stability.