Chaos in Nonlinear Random Walks with Non-Monotonic Transition Probabilities

Random walks serve as important tools for studying complex network structures, yet their dynamics in cases where transition probabilities are not static remain under explored and poorly understood. Here we study nonlinear random walks that occur when transition probabilities depend on the state of the system. We show that when these transition probabilities are non-monotonic, i.e., are not uniformly biased towards the most densely or sparsely populated nodes, but rather direct random walkers with more nuance, chaotic dynamics emerge. Using multiple transition probability functions and a range of networks with different connectivity properties, we demonstrate that this phenomenon is generic. Thus, when such non-monotonic properties are key ingredients in nonlinear transport applications complicated and unpredictable behaviors may result.


I. INTRODUCTION
Random walks have long served as a powerful tool for studying and understanding the structural properties of complex networks due to the wealth of information they provide with remarkably simple dynamics [1,2]. Examples of the widespread utility of using random walks and diffusion dynamics in the context of complex networks include Google's PageRank centrality [3][4][5], network search and exploration [6][7][8], modeling transport, diffusion, and movement processes [9][10][11], detecting communities and other network structures [12,13], and finding geometric properties and topological embeddings for dimensionality reduction [14]. The Markovian, i.e., memory-less, property of a random walk coupled with the static nature of transition probabilities (or transition rates) yields a linear dynamical system where, assuming the relatively mild condition of a network structure being primitive, the dynamics are guaranteed to converge to a unique, globally-attracting fixed point or stationary distribution [15,16]. However, in an effort to generalize these conditions to better model more realistic and complex phenomenon, the static nature of transition probabilities may be relaxed, giving rise to nonlinear random walks [17][18][19]. In particular, by focusing on the case of discrete-time where nonlinearity arises from the the dependence of transition probabilities on the current system state in a manner that biases random walkers towards or away from nodes that are heavily populated, a rich landscape of nonlinear phenomena has been observed, including period-doubling bifurcations, multistability, and quasi-periodic dynamics [20,21].
In this Letter we extend this paradigm to consider nonlinear random walks with non-monotonic transition probabilities. In particular, while the nonlinear random walks with monotonic transition probabilities studied in Refs. [20,21] operate on the modeling premise that random walkers are preferentially biased towards nodes with either proportionally more or fewer random walkers present, allowing for non-monotonic transition probabilities allows this bias to be generalized. For instance, depending on the shape of the function that defines * persebastian.skardal@trincoll.edu the transition probabilities, random walkers might be biased strongly towards both nodes where very many and very few random walkers are present, while being biased away from moderately populated nodes. Alternatively, random walkers might be biased more strongly towards nodes where a moderate number of random walkers are present, while being biased away from nodes that are both heavily and sparsely populated. In general, such non-monotonic transition probabilities may arise in nonlinear transport applications where the decision-making process that informs movement is nuanced. Examples may include personal travel, where an individual or family may prefer to visit either very remote (e.g., camping) or very densely populated (e.g., large cities) vacation destinations, and finances and investments, where individuals may seek to invest in or choose to purchase commodities from institutions that are neither exceedingly large nor small. Here we explore the dynamics that emerge from nonlinear random walks with non-monotonic transition probabilities and show that they give rise to chaotic dynamics that are not present when transition probabilities are monotonic. We show that this phenomenon is widespread by exploring multiple choices for functions that define the transition probabilities of a system as well as a range of different network topologies. Our results suggest that in nonlinear transport processes where nonmonotonicity in transition probabilities or transition rates is a key ingredient, the resulting behavior may be extremely complicated, making prediction and forecasting difficult even with very accurate system information.
The remainder of this paper is organized as follows. In Sec. II we present the governing equations and present our main results, illustrating the emergence of chaotic dynamics in nonlinear random walks with non-monotonic transition probabilities. In Sec. III we explore dynamics on different network topologies, using both a minimal model and a large real-world network. In Sec. IV we present an initial exploration of mean return times. In Sec. V we conclude with a discussion of our results.

II. GOVERNING EQUATIONS AND CHAOTIC DYNAMICS
We focus our attention on the case of discrete-time nonlinear random walks on complex networks. In the thermody- namic limit of infinitely-many random walkers, the dynamics on a network with N nodes evolves according to where p i (t) is the fraction of random walkers present at node i at time t for i = 1, . . . , N and π ij (p(t)) denotes the transition probability of a random walker moving from node j to node i at time t. Importantly, the transition probabilities in the transition matrix Π(p(t)) depend on the current system state, denoted by the vector p(t) = [p 1 (t), . . . , p N (t)] T , yielding a nonlinear dynamical system. Since the entries of Π(p(t)) may also be viewed as conditional probabilities, the columns must sum to one, conserving the total probability in the vectors p(t) under the dynamics of Eq. (1), which may also be written in vector form p(t + 1) = Π(p(t))p(t). The most natural formulation for the transition probabilities π ij (p(t)) is that they depend proportionally on the set of states p i (t) at destination nodes i that stem from node j, yielding where a ij is the entry in the adjacency matrix A corresponding to the existence/strength of the link j → i and f is a function that maps non-negative values to positive values, f : [0, ∞) → (0, ∞). Unless otherwise noted, here we focus on the case undirected, binary networks such that a ij , a ji = 1 if nodes i and j are connected, and otherwise a ij , a ji = 0, although these conditions may be easily relaxed to explore dynamics on weighted or directed networks. Previous work investigating the dynamics of Eqs. (1) and (2) focused on the monotonic function f (p) = exp(αp), where α is a biasing parameter that resulted in random walkers being preferentially moved towards heavily or sparsely populated nodes for positive and negative values of α, respectively. However, to relax this condition and allow for more nuanced biasing properties we consider the following functions: For positive α the function f 1 (p) increases from f (0) = 1, reaches a local maximum at p = z/2, then decreases with Similarly, f 2 (p) initially increases, reaches a local maximum (roughly at p = 0.43z), and begins to decrease, but then reaches a local minimum before leveling off with f 2 (p) → 1 − as p → ∞. [Both f 1 (p) and f 2 (p) are plotted in Fig. 1(f) using α = 3 and z = 1.288 and 2 for f 1 and f 2 , respectively.] Noting that α controls the range of both f 1 and f 2 while z controls the location of the local maxima and minima, we treat α as the primary bifurcation parameter that will be varied in this work, and use z as a scaling parameter that primarily is chosen to suit the size of the underlying network topology. Note also that setting α = 0 yields f (p) = 0, which in turn yields π ij = a ij /k j in Eq. (2), i.e., the classical static, unbiased random walk. In this case, provided that the network structure is primitive (i.e., the adjacency matrix is irreducible and aperiodic) the dynamics converge to a unique stationary state given by the normalized degree vector, In contrast to these simple linear dynamics, transition probabilities defined by the function f in Eq. (3) yields richer dynamics.
We begin by illustrating our main result using a small network of size N = 8 with the function f 1 given in Eq. (3). Setting z = 1.288/N , we plot in Fig. 1(a) the bifurcation diagram using the quantity p(t)−p * 1 as α is increased from 0 to 60. For sufficiently small α the dynamics converge to a fixed point near p ≈ p * 1 , but soon undergo a period-doubling bifurcation at α ≈ 0.6. Beyond this point we observe a cascade of further period-doubling bifurcations, eventually giving way to chaotic dynamics with subsequent bands of higher-order peri- odic and chaotic behavior, as is typical in other generic maps displaying chaotic dynamics [22]. To confirm that the dynamics are in fact chaotic, we plot in Fig. 1(b) the largest Lyapunov exponent λ, which agrees nicely with the bifurcation diagram. We note here that the conservation of probability (or more generally, the conservation of the sum of p(t)) in Eqs. (1) and (2) implies that the dynamics lie on the simplex characterized by Moreover, the Jacobian matrix for Eq. (1) always has at least one eigenvalue given precisely by λ = 1, making this simplex a center manifold that is marginally stable in the transverse direction away from the simplex. Since perturbations in this direction neither grow or decay, this guarantees that at least one Lyapunov exponent will always be zero, so the largest Lyapunov exponent will be λ = 0 for non-chaotic (e.g., periodic and quasi-periodic) dynamics and λ > 0 for chaos. Next we probe the chaotic dynamics further, setting α = 52.2 and plotting the chaotic attractor in the ( p(t)− p * 1 , p(t)− p * 2 ) plane, where p * 2 = 1/N is the normalized constant vector, in Fig. 1(c). The stretching and folding that is typically present in strange attractors can be more easily identified by zooming in on, e.g., the upper-right hand side of the attractor, as shown in Fig. 1(d). To confirm the fractal nature of the structure we perform a fractal analysis that reveals a correlation dimension [22] given by the non-integer value of d corr ≈ 1.1015. In Fig. 1(e) we illustrate the network structure used in Figs. 1(a)-(d).
Before continuing on, we pause to discuss the implications of these new findings. While the applicability of classical random walks (i.e., those with static transition probabilities, either biased or not) cannot be overstated, the relaxation of static transition probabilities is a key ingredient in modeling transport dynamics when decision-making is more nuanced. The results presented above show that when non-monotonicity is a feature in such a process, then not only may the the system not converge to a simple fixed point or periodic orbit of low order, but possibly a higher-order periodic orbit or even a chaotic state. Moreover, in the case of a chaotic state we lose long-term predictability of the system behavior, posing a hurdle for forecasting beyond the very short-term, even with very accurate system information.

III. MINIMAL MODEL AND VARYING NETWORK TOPOLOGY
We now demonstrate that the chaotic dynamics illustrated in the example above are generic in the sense that they occur for multiple transition probability functions as well as over a wide range of network topologies. We begin with the fundamental task of finding a minimal system that exhibits chaos. Note that, since we do not allow for self-links, i.e., random walkers may not remain at the same node from one time step to the next, N = 3 is the smallest network with non-trivial dynamics. In fact, chaos can be observed by breaking the symmetry of a complete 3-node network by increasing the strength of one directed link in the network, as illustrated in Fig. 2(a), where the thick arrow indicated a directed link with strength two while the rest of the directed links have strength one. Using f 1 and z = 1.288/N , we plot in Fig. 2(b) the bifurcation diagram for the probability p i (t) describing the upper-most node in Fig. 2(c) and in Fig. 2(c) we plot the corresponding largest Lyapunov exponent. Thus, with only three nodes we observe hallmarks of chaos such as period-doubling bifurcations, bands of periodicity and chaos, and positive Lyapunov exponents. We note here that we found the symmetry breaking described above to be necessary to destabilize the stationary state fixed point, although this does not rule out the possibility for multistability with a chaotic orbit, although we found no such structures in our exploration.
Next we demonstrate that the chaotic dynamics exhibited in the previous examples persists for different probability transition functions and for different network structures. First, using f 2 (p) given in Eq. (4) with z = 2/N on the same network structure depicted in Fig. 1(e), we plot the bifurcation diagram and largest Lyapunov exponent as α is increased from 0 to 20 in Figs. 3(a) and (b), respectively. Again, we observe a cascade of period-doubling bifurcations that eventually give way to chaotic dynamics just above α = 10 that agree with a positive Lyapunov exponent. We also consider the same dynamics on a much larger network structure generated from real data describing the interactions between the 500 busiest airports in the United States in 2002 [23], with nodes representing airports and links existing between airports if they have a direct flight going from one to the other. (Here we consider the undirected, unweighted version of this network.) In addition to being a larger, real-world network, this structure is much more heterogeneous, with several of the N = 500 nodes having a nodal degree larger than 100 while many others have only one or a handful of links. Again using f 2 with z = 2/N we plot the bifurcation diagram and largest Lyapunov exponent as α is increased from 0 to 20 in Figs. 3(c) and (d), respectively. (Here we have plotted the bifurcation diagram using the probability p i (t) corresponding to the node of largest degree, given by k i = 145.) With the increased size, heterogeneity, and overall complexity of the network topology, some signature effects such as period-doubling bifurcations and intermittent bands of periodic and chaotic dynamics are lost, however the largest Lyapunov exponent indicates chaos shortly after α = 5.
In addition to the network structures used so far in this Letter, namely the first toy network (N = 8), the minimal network (N = 3), and the US airport network (N = 500), we have simulated the random walk dynamics of Eqs. (1) and (2) for a number of number of random networks with different parameters. In the Appendix we present these results, showing several examples for each combination of network structures (i) of size N = 8 with mean degree k = 2, (ii) N = 24 with k = 4, and (iii) N = 100 with k = 10 with both f 1 and f 2 , given by Eqs. (3) and (4). We also explore the effect of different choices of non-monotonic functions, specifically presenting results using two additional functions that differ from f 1 and f 2 in that they begin by decreasing before hitting a pair of local extrema. This larger collection of simulations show qualitatively similar results as those presented above and demonstrate that the chaos we observe is generic in the sense that it occurs over a wide range of network structures and for multiple transition functions.

IV. MEAN RETURN TIMES
Before we conclude, we explore the effect that the complicated dynamics that emerge in nonlinear random walks with non-monotonic transition functions have on characteristic return times, i.e., the typical time it takes for a random walker to return to a particular node. In particular, returning to the toy network illustrated in Fig. 1(e) f 1 with z = 1.288, we consider the four choices α = 3, 10, 15, and 20, yielding period-2, period-four, period-8, and chaotic dynamics. In each scenario, we simulate the trajectory throughout the network of a particular random walk (using 10 6 timesteps) to calculate the mean and standard deviation of return times for each node, which we plot in Figs. 4(a) and (b), using blue circles, green squares, orange triangles, and red diamonds for the four different kinds of dynamics. In particular, as the underlying dynamics become more complicated, the return times for different nodes become more heterogeneous, with mean return times ranging between 4.73 and 2.03 × 10 1 when dynamics are period-2, but 2.04 and 9.65 × 10 3 when dynamics are chaotic. The overall spread of mean return times across the network is described in Fig. 4(c), where we plot the standard deviation of mean return times across the network for each of the four cases. While this results constitute only four specific parameter choices, they tend to agree with other exploratory simulations that suggest that more complicated dynamics result in larger heterogeneity in mean return times throughout a network. The topic of return times (as well as hitting times etc.) in nonlinear random walks on networks deserves further attention in future research.

V. DISCUSSION
In this work we have studied the dynamics that take place in a class of discrete-time nonlinear random walks where transition probabilities depend on the current state of the system. Specifically, we showed that when transition probabilities are non-monotonic, that is, they are defined by a non-monotonic function of the state of the random walk dynamics, complicated behaviors that include chaos emerge. In the case of small networks, the classical signatures of low dimensional chaotic systems such as cascades of period-doubling bifurcations and intermittent bands of chaos and periodicity can be observed. For larger systems, however, these fingerprints become difficult to identify with added complexity, but chaos remains a prominent feature. We also illustrated that these dynamics are generic in the sense that they occur over a wide range of network topologies and with several transition functions with different properties. (In addition to the examples given here, see the Appendix). These results suggest that when the rules governing nonlinear transport throughout a system are made non-static and non-monotonic, for instance, by nuanced decision-making that directs individuals' trajectories, the resulting dynamics may become very rich and complicated, making long-term prediction and forecasting difficult even with very accurate system information.