Random walks on networks with preferential cumulative damage: Generation of bias and aging

For the initial configuration, we consider undirected connected networks with $N$ nodes $i=1,\ldots,N$ described by an adjacency matrix $\mathbf{A}$ with elements $A_{ij}=A_{ji}=1$ if there is an edge between the nodes $i$ and $j$ and $A_{ij}=0$ otherwise; in particular, $A_{ii}=0$ to avoid edges connecting a node with itself. In this structure, we denote the set of nodes as $\mathcal{V}$ and the set of edges (links) as $\mathcal{E}$ with elements $(i,j)$. For each pair in $\mathcal{E}$, the corresponding element of the adjacency matrix is non-null. In the following the link from $i$ to $j$ denoted as $(i,j)$ is independent from $(j,i)$, and we denote as $|\mathcal{E}|$ the total number of different edges in the network.
Additionally to the network structure, the global state of the system at time $T=0,1,2,\ldots$ is characterized by a $N\times N$ matrix of weights $\mathbf{\Omega}(T)$ with elements $\Omega_{ij}(T)\geq 0$ and $\Omega_{ii}(T)=0$ which describe weighted connections between the nodes. The matrix $\mathbf{\Omega}(T)$ contains information of the state of the edges and in general is not symmetric. In order to capture in the model the damage impact affecting the complex system, we introduce the variable $T$ as the measure of the total number of damage hits in the links of the network. The variable $T$ can also be conceived as a time measure by assuming a constant damage impact rate, i.e. successive damage events occur with a constant difference of times $\Delta T=1$. We introduce for each link $(i,j)\in\mathcal{E}$ a stochastic integer variable $h_{ij}(T)$ where $h_{ij}(T)-1$ counts the number of random faults that exist in this link at time $T$. The values $h_{ij}(T)$ for all the edges are numbers that evolve randomly, and a new fault in the link $(i,j)$ appears at time $T$ with a probability $\pi_{ij}(T)$ which is given by

for $\,T=1,2,\ldots$ with the initial condition $h_{ij}(0)=1$, i.e. no faults exist for all the edges at $T=0$. Equation (1) indicates the probability for the event that at time $T$ the number of faults $h_{ij}(T)=h_{ij}(T-1)+1$ are increased by one. In our analysis, the damage is distributed without maintaining the symmetry of the initial undirected network. In this manner, in the general case $h_{ij}(T)$ is independent of the value $h_{ji}(T)$ and also $\pi_{ij}(T)$ from $\pi_{ji}(T)$ which is generating a biased network. With Eq. (1) at $T=1$ the first hit (fault) is randomly generated for any selected link $(i,j)$ with equal probability $\pi_{ij}(1)=\frac{1}{|\mathcal{E}|}$. The occurrence of the second fault at $T=2$ depends on the previous configuration and so on. An asymptotic analysis of the time-evolution of the fault number distribution resulting from Eq. (1) shows that a power-law scaling with features of a stochastic fractal emerge (See Eq. (33) in A and consult also Ref. [22]). An essential feature of the probabilities in Eq. (1) is that they produce preferential damage if a link has already suffered damage in the past. A link has a higher probability to get a fault with respect to a link never being damaged. Such preferential random processes have been explored in different contexts in science (see Ref. [32]), being a key element in our model that generates complexity in the distribution of damage reflected by asymptotically emerging power-law and fractal features (See A).
The choice of generating law of faults in Eq. (1) is based on the observation of aging changes in living beings. The development of aging changes such as age spots is self-amplifying and inhomogeneous. Age spots develop in this way because misrepaired structures have increased damage-sensitivity and reduced reparation-efficiency [19, 31].
Now, we aim to describe how the structure reacts to the damage hits occurring stochastically to the edges. We describe the effects of the damage by using the information in the matrix of weights $\mathbf{\Omega}(T)$. In terms of the values $h_{ij}(T)$, the matrix $\mathbf{\Omega}(T)$ defines the global state of the network containing the complete information on the network topology at time $T$. Its matrix elements

contain the local information on the damaged state of edge $(i,j)$ and $\alpha\geq 0$ is a real-valued parameter that quantifies the effect of the damage in each link. We call $\alpha$ the misrepair parameter since it describes the capacity of the system to repair damage in the links: In the limit $\alpha\to 0$ the system responds with perfect reparation with $\Omega_{ij}(T)\to A_{ij}$ as in a perfect undamaged structure, and the effect of the stochastically generated faults is null. In contrast, in the limit $\alpha\to\infty$, a hit in a link is equivalent to its removal from the network. This limit corresponds to a complex system without repair capacity. The fault accumulation dynamics described by Eq. (1) together with the ‘misrepair equation’ (2) is therefore able to mimic the phenomena related to aging processes observed in living beings [18, 33].

In Fig. 1, we illustrate the algorithm for the generation of cumulative damage in a network with $N=10$ nodes formed by a comb-like structure with a ring. For this directed graph, we generate random hits in the links at times $T=1,2,\ldots,100$. The probability to generate a hit at time $T$ in the directed edge $(i,j)$ is proportional to the previous configuration given by $h_{ij}(T-1)$ in Eq. (1). In this way, using Monte Carlo simulations, we depict the damage in two realizations. The values of $h_{ij}$ in the edges are represented with colors codified in the colorbar whereas the respective widths reflect their capacity of transport given by $\Omega_{ij}(T)$ in Eq. (2) with the parameter $\alpha=2$. The results in the final configuration at $T=100$ give us an idea of the variety of structures how the damage can be distributed in the network affecting the nodes in the periphery or in the central ring.

Let us now come back to the general algorithm for cumulative damage described before, at a completely different scale of times (significantly less than the characteristic times for the damage) takes place the movement of a random walker in the network with discrete steps $\Delta t$ at times $t=0,\Delta t,2\Delta t,\ldots$. For each state of the network’s weights, the random walker can move visiting nodes and is capable to visit all the nodes of the structure that changes very slowly with time $T$. In a determined configuration at time $T$, the transition probability matrix $\mathbf{W}(T)$ describing the random walker is defined by the elements $w_{i\to j}(T)$ with the probability to pass from node $i$ to node $j$

We assume a Markovian time-discrete random walker that performs at any time increment $\Delta t$ a random step from one node to another. This process is defined by the master equation [3, 26, 34]

valid for $t\ll\Delta T=1$. In this master equation $P_{ij}(t,T)$ indicates the probability that the walker that starts its walk at node $i$ at $t=0$ occupies node $j$ at the $n$-th time step $t=n\Delta t$.
The one-step transition matrix in Eq. (3) is constructed such that the walker has to change the node at any step (i.e $w_{i\to i}(T)=0$). The canonical representation of the ($t=n\Delta t$) of the $n$-step transition matrix is

Where we use Dirac’s (bra-ket) notation. In Eq. (5), $|\phi_{m}(T)\rangle,\langle{\bar{\phi}}_{m}(T)|$ denote, respectively, the right and left eigenvectors of the transition matrix with the respective eigenvalues $0\leq|\lambda_{m}(T)|\leq 1$. The walk which we assume to take place on a (strongly connected) directed weighted and finite network is ergodic with the unique eigenvalue $\lambda_{1}(T)=1\forall T$ reflecting row stochasticity (with corresponding right-eigenvector having identical components) of the transition matrix $\sum_{j=1}^{N}w_{i\to j}(T)=1$ and $|\lambda_{m}(T)|\leq 1$ maintained $\forall T$ for $m=1,\ldots,N$ (see Ref. [30] for a detailed analysis of undirected networks, and Refs. [35, 36] for directed networks). The stationary distribution $P_{j}^{(\infty)}(T)\equiv\lim_{t\to\infty}\frac{1}{t}\sum_{t^{\prime}=0}^{t}P_{ij}(t^{\prime},T)$, that gives the probability to find the random walker at the node $j$ in the limit $t\to\infty$, is given by [27, 30, 37]

where, since $\left\langle i|\phi_{1}(T)\right\rangle=\mathrm{constant}$ (independent of $i$), the stationary distribution $P_{j}^{(\infty)}(T)$ does not depend on the initial condition.
Additionally, we have the mean first passage time $\left\langle{\cal T}_{ij}\right\rangle$ that gives the average number of time steps (in units of $\Delta t$) the walker needs to travel from node $i$ to node $j$ in the form (see Refs. [30, 35, 37, 38] for a complete derivation)

In this relation for $i=j$ the second term and for $j\neq i$ the first term vanishes. For $i=j$ this relation gives the mean first return time or mean recurrence time [30, 39]

This is the Kac-formula relating the mean recurrence time with the inverse of the stationary distribution $P_{j}^{(\infty)}(T)$.

The results in Eqs. (5) and (7) show that the eigenvalues of the transition matrix $\mathbf{W}(T)$ at different times $T$ contain important information of the diffusive transport of each configuration at times $T=0,1,\ldots$. Also, due to the asymmetry in some edges, in the general case, the eigenvalues $\lambda_{j}(T)$ can take complex values. However, the process of cumulative damage for $\alpha$ finite maintains the ergodicity of the random walker described by $\mathbf{W}(T)$ since the initial structure is connected and none of the edges is completely removed with the hits at different times (for a detailed discussion about ergodicity on directed networks see Ref. [35] and references therein) thus the mean first passage times in Eqs. (7)-(8) are well defined. The relation between the eigenvalues of the transition matrix and node or edge removal strategies in various dynamical processes on networks has been explored numerically by Restrepo et. al. in Refs. [40, 41].
In the following part, we explore the effects of cumulative damage in the network topologies illustrated in Fig. 2. In Fig. 3 we depict the eigenvalues of the transition matrix $\mathbf{W}(T)$ for different initial networks with $N=50$ nodes and two realizations of the random cumulative damage algorithm. Figure 3(a) displays the results for a wheel graph formed by a single node connected to all nodes of a cycle (with $49$ nodes). In Fig. 3(b) we show the eigenvalues for a barbell graph, i.e., a network with two well-defined communities composed of two fully connected networks (with $20$ nodes each) connected by a linear graph (with $10$ nodes) [42]. Conversely, in Figs. 3(c)-(d) the initial structures without damage are random networks obtained with the Watts-Strogatz algorithm [43] generated from a ring with nearest-neighbor and next-nearest-neighbor links and a rewiring probability of $p=0.05$ and the Barabási-Albert algorithm, generated with the preferential attachment rule [44].
The results in Fig. 3 show the evolution of the eigenvalues with the damage at times $T=1,2,\ldots,100000$, the initial eigenvalues are real and are represented with circles in each case. The damage in each directed edge induce asymmetry and the bias generated produces complex eigenvalues depicted with different colors. Here it is worth mentioning that we apply the same approach to several types of trees where we found that in this type of graph the eigenvalues are real. This is an important case for which we do not have an analytical result but we hope this leads to future research in the spectral properties of directed weighted trees. These findings motivate the introduction of a measure to quantify the evolution of bias generated with the accumulation of damage in different networks. We discuss this point in the following part.

Once defined four measures to quantify the asymmetry on transport on networks, we explore the evolution of these quantities with time $T$ for structures with $N=50$ nodes. We analyze three deterministic networks: a ring, a wheel and a barbell graph, and three random structures generated with the Watts-Strogatz (WS), Erdös-Rényi (ER), Barabási-Albert (BA) algorithms. In Fig. 4 we depict the evolution of $\mathcal{S}(T)$ that denotes $\mathcal{S}_{W}(T)$, $\mathcal{S}_{\mathrm{Entropy}}(T)$, $\mathcal{S}_{\mathrm{MFPT}}(T)$, and $\mathcal{S}_{\gamma}(T)$ considering a fractional non-locality with $\gamma=0.5$. The results were obtained from Monte Carlo simulations of the accumulation of damage algorithm with a misrepair parameter $\alpha=1$ in Eq. (2). For each time $T$ we calculate the values of $\mathcal{S}(T)$, we consider 100 realizations to obtain the ensemble average of each $\mathcal{S}(T)$.
In Fig. 4(a) we analyze a ring, for this regular network at $T=0$ all the local and non-local measures register $\mathcal{S}(0)=1$, this value occurs when the network has not suffered any damage, therefore it is the largest value that this parameter can reach. With the action of damage at $T=1,2,\ldots,1000$ the ensemble average of all the $\mathcal{S}(T)$ analyzed decrease with $T$. However, in the results for all the realizations, we see that the ensemble average of each measure presents high dispersion, revealing the susceptibility of the ring to damage, the dispersion is huge in the $\mathcal{S}_{\mathrm{MFPT}}(T)$, and $\mathcal{S}_{\gamma}(T)$ showing how the damage in a single link generates asymmetry that affects the transport on the whole ring. In Fig. 4(b) we depict the results for a wheel graph. In this case, at $T=0$, the local measure $\mathcal{S}_{W}(0)=0.5577$ due to an initial bias in the transition matrix $\mathbf{W}$ due to its normalization, an effect that is reduced in non-local measures $\mathcal{S}_{\mathrm{MFPT}}(0)=0.9627$, and $\mathcal{S}_{\gamma}(0)=0.9755$. We see in the ensemble average the reduction of $\mathcal{S}(T)$ with the damage; but, the values in the realizations are less dispersed in comparison to the results obtained for the ring. This reduced dispersion shows the resistance of the structure at local and global scales, a consequence of the existence of multiple paths that can connect two nodes on the network. On the other hand, the barbell graph analyzed in Fig. 4(c) combines two fully connected networks with 20 nodes with a linear graph with 10 nodes, each two fully connected subgraphs are resistant to damage due to the diverse alternative routes to defective links. We see that the local measures decaying slowly with $T$ and with low dispersion. However, at a global scale, the nodes that form the linear graph make the structure fragile, this condition of the barbell graph is observed in the non-local measures $\mathcal{S}_{\mathrm{MFPT}}(T)$ and $\mathcal{S}_{\gamma}(T)$ with fast decay under damage and with high dispersion in the same way as the results observed for the ring in Fig. 4(a).
Regarding the random networks, in Fig. 4(d) we analyze a Watts-Strogatz network obtained from a regular cyclic structure with four nodes (a ring with two additional links in each node) and a random rewiring with probability $p=0.05$, see Ref. [43]. In this case, the rewiring creates a network with more complexity than the initial structure, the alternative paths connecting two nodes make this structure more resistant to damage. In Fig. 4(d) it is worth to notice that the average of $\mathcal{S}_{W}(T)$ and $\mathcal{S}_{\mathrm{MFPT}}(T)$ have a similar behavior showing that the rewiring creates global connectivity of the network, for this reason, $\mathcal{S}_{W}(T)$ behaves like a non-local measure. In Fig. 4(e), we explore the evolution of the asymmetry in an Erdös-Rényi network at the connectivity threshold $p=\log(N)/N$. Due to the reduced number of links the connectivity of the network can be affected with the removal of a link. We see that $\mathcal{S}_{\mathrm{Entropy}}(T)$ decays faster with $T$ in comparison to the results in other networks. For the scale-free network in Fig. 4(f), $\mathcal{S}_{W}(0)=0.5166$ showing that the heterogeneity of the network already includes a bias that increases the probability to pass to high connected nodes, the slow variation of $\mathcal{S}_{W}(T)$ with the damage shows that the distribution of weights in the links has a similar “complexity” since the damage is generated with a preferential aggregation. Also, we see a low dispersion of the values $\mathcal{S}_{\mathrm{MFPT}}(T)$, showing that the high number of paths connecting two nodes make this structure capable to resist the damage without creating strong non-local asymmetry. Our findings in Fig. 4 for deterministic and random networks show that the variations of the local asymmetry measures are in a good approximation rescaled versions; in this sense, the two local measures are similar in quantifying variations of the asymmetry due to damage. Something similar occurs with the two non-local measures; however, in this case, the variation of the parameter $\gamma$ can be used to consider other non-local effects.
In addition to the temporal evolution, it is important to understand the modifications introduced by the parameter $\alpha$ in Eq. (2) that modifies the effect of damage in the transport. In our ‘misrepair picture’ small $\alpha$ corresponds to high reparation capacity of a damaged link and large $\alpha$ can be related to ‘bad’ reparation capacity in the complex system. In the limit $\alpha=0$, the damage does not alter the transport and $\alpha\to\infty$ is equivalent to the complete removal of the link. In Fig. 5 we analyze the ensemble average of local and non-local symmetry measures $\mathcal{S}$ with $\alpha=0.5,1.0,1.5,2.0$ and $T=1000$ for the networks explored in Fig. 4. The results are generated with 1000 Monte Carlo simulations of the process, the error bars represent the standard deviation of the data. All the asymmetry measures in the six networks analyzed decrease with the increase of $\alpha$. Regarding the dispersion of the data in the error bars, we see that the values $\mathcal{S}_{W}(T)$ are less disperse for all the networks; in contrast, the values $\mathcal{S}_{\mathrm{MFPT}}(T)$ have the largest dispersion revealing that this measure is more susceptible to how damage is distributed. Finally, we see that the non-local measures behave in a similar way with $\alpha$ showing that $\mathcal{S}_{\gamma}(T)$ with $\gamma=0.5$ is a good measure to quantify the average asymmetry in the transport.