Adaptive rewiring based on diffusion balances stability and plasticity in weighted networks while evolving ‘brain-like’ structure

Activity dependent plasticity is the brain’s mechanism for reshaping neural connections. Representing activity by graph diffusion, we model plasticity as adaptive rewiring. The rewiring involves adding shortcut connections where diffusion on the graph is intensive while pruning underused ones. This process robustly steers initially random networks to high-levels of structural complexity reflecting the global characteristics of brain anatomy: modular or centralized small world topologies, depending on overall diffusion rate. We extend this result, known from binary networks, to weighted ones in order to evaluate the flexibility of their evolved states. Both with normally- and lognormally-distributed weights, networks evolve modular or centralized topologies depending on a single control parameter, the diffusion rate, representing a global homeostatic or normalizing regulation mechanism. Once settled, normally weighted networks lock into their topologies, whereas lognormal ones allow flexible switching between them, tuned by the diffusion rate. For a small range of diffusion rates networks evolve the largest variety of topologies: modular, centralized or intermediate. Weighted networks in the transition range show topological but not weighted rich-club structure matching empirical data in the human brain. The simulation results allow us to propose adaptive rewiring based on diffusion as a parsimonious model for activity-dependent reshaping of the brain’s connections. Author Summary The brain is adapting continuously to a changing environment by strengthening or adding new connections and weakening or pruning existing ones. This forms the basis of flexible and adaptable behaviors. On the other hand, uncontrolled changes to the wiring can compromise the stability of the brain as an adaptive system. We used an abstract model to investigate how this basic problem could be addressed from a graph-theoretical perspective. The model adaptively rewires an initially randomly connected network into a more structured one with properties akin to the human brain, such as small worldness and rich club structure. The adaptive changes made to the network follow the heat diffusion, an abstract representation of brain functional connectivity. Moreover, depending on a parameter of the model, the heat diffusion rate, either modular or centralized connectivity patterns emerge, both found across different regions of the brain. For a narrow range of intermediate heat diffusion rates, networks develop a full range from modular to centralized connectivity patterns. Once settled into a connectivity pattern networks with normally distributed weights lock into that state, whereas networks with lognormally distributed weights show greater flexibility to adjust, while maintaining their small-world and rich club properties. Networks with lognormally distributed weights, therefore, show the combination of stability and flexibility needed to address the fundamental requirements of adaptive networks.

An increase in a network's average clustering coefficient (C) or a decrease in its average 147 path length (L) induces by definition an increase in its small worldness. We probed the dynamics 148 of C and L that drove up the small worldness of the networks. We found that along τ, S shows a 149 strong positive correlation with C (ρ binary = 0.91, ρ normal = 0.96, ρ lognormal = 0.97) and a moderate 150 one with L (ρ binary = 0.42, ρ normal = 0.52, ρ lognormal = 0.31). In the Watts and Strogatz model (24), 151 the authors start with a structured network and parametrically swap ordered connections with 152 random ones. This induces a significant decrease in L but a more moderate one in C, resulting in 153 an increase in the small worldness of the rewired network. Our model is in a way the reverse 154 process, in which rewiring induces a random network configuration to morph into a more 155 structured one, and thus becoming small world due to an increase in C (See S1 Appendix for a 156 more detailed description). Binary, normally and lognormally weighted networks develop modularity as a function of 191 τ and p random in a similar fashion ( Fig 3A). Generally, Q initially increases as a function of τ, 192 reaching a plateau, but eventually drops off to a value close to that of a random network (τ = 0) 193 ( Fig 3B). Binary networks reach their maximum Q faster (τ=10 -15 ) compared to normal (τ=2), 194 and lognormal networks (τ = 4) ( Fig 3B) We probed the heavy tail degree and strength distributions of representative centralized 245 networks (τ binary = 5, τ normal = 5, τ lognormal = 7, in all cases p random = 0. 2) The power-law and 246 lognormal distributions were a better fit than the exponential one ( Fig S2). Power law 247 distribution functions are of the form P(k) ~ k -α . For the degree distributions the exponent α that 248 fitted best the data was 1.7, for the strength distributions it varied between 2 and 2.6 for the 249 different types of networks (See S2 Appendix for a more detailed description of the analysis).
250 We used publicly available code for this analysis (26).

252
Network structure at the transition range 253 Our previous analyses showed that, for a given random rewiring probability, and 254 depending on the value of the control parameter τ, the network could be rewired to be either in a 255 modular or centralized state. At the boundary between those two states the structure of the 256 network is ambiguous. Moving the control parameter τ across the boundary causes a phase 257 transition from modular to centralized connectivity. Our aim in this section is to probe the 258 properties of the network at this transition. Regarding the distribution of modularity values at τ transition , one hypothesis is that the 270 distribution is bimodal with one peak centered at the modular region and the other at the 271 distributed one. A contrasting hypothesis is that the distribution is unimodal with the peak in-272 between the modular and centralized regions. For both binary and normal networks, we found a 273 compromise between those two hypotheses: the modularity distribution at τ transition is uniform, 274 with the modularity values on the left and right boundaries giving rise to modular and centralized 275 networks respectively (Fig 5A and B, center plots). Furthermore, depending on the sign of small 276 perturbations (δτ = 0.1) to τ transition , the resulting modularity distributions show a distinct peak 277 either at the modular region (for τ = τ transition -δτ) or the centralized one (for τ = τ transition +δτ) ( Fig  278 5A and B, the leftmost and rightmost plots). On the other hand, the modularity distribution for 279 lognormal networks is more in line with the second hypothesis: it is unimodal, and its peak is in-280 between the modular and centralized regions. Unlike the other two networks, its peak is not 281 shifted for the small perturbations of its τ value (Fig 5C). Starting from a random configuration, the rewiring process at τ transition can result in vastly 293 different connectivity patterns ranging from highly modular to highly centralized networks (Fig   294 5). On the other hand, at the fringes of its variability range, the initial random network is biased 295 towards a more modular or centralized connectivity pattern. We asked whether small biases in 296 the initial random network are amplified during rewiring, essentially predicting the connectivity 297 pattern of the final network. We tested this hypothesis by comparing the modularity indices of 298 1000 networks before and after rewiring. We found a very weak correlation between the 299 modulation values of the random and rewired networks (ρ binary = 0.10, ρ normal = 0.14, ρ lognormal = 300 0.08; Fig 6A). This indicates that the connectivity pattern of the rewired network at τ transition is 301 nearly independent of the connectivity bias of the initial random network. Our next analyses probed the flexibility of network developed at τ transition , when τ is 319 subsequently set to a different value . An initially random network was rewired 4000 times at 320 τ transition ; next another 4000 rewirings were subsequently performed at τ test . As long as τ test 321 deviates from τ transition only slightly (τ test = τ transition  0.1), the network's final connectivity pattern 322 remains close to the one obtained after the first 4000 rewirings (more so for normal and binary 323 networks than for lognormal ones; Fig S3). Fig 6C shows that with τ test << τ transition , all networks 324 are driven to modular connectivity patterns irrespective of their initial topologies. On the other 325 hand, for τ test >> τ transition , network evolution is strongly affected by the state originally reached at 326 τ transition , at least for the binary and normal networks (Fig 6D). For lognormal networks, however, 327 the correlation coefficient is low, indicating that the explanatory strength of the linear fit is weak.   (Fig 7B and 7C). Thus for τ test > τ transition , binary and normal 365 networks lock in their initial connectivity pattern similarly to when τ test = τ transition , whereas the 366 lognormal network shows a more flexible behavior with a bias nevertheless to switch into a more 367 centralized state. Assortativity and rich club structure 378 We probed the possibility that the rewiring process favors homophily by measuring the 379 topological assortativity coefficient, r. Networks with positive r are assortative, meaning that 380 nodes with similar degrees tend to connect. Networks with negative r are dissasortative; in this 381 case nodes with dissimilar degrees are more prone to connect. We found that for modular 382 networks, r shows weak or zero assortativity; centralized networks dip into the dissasortative 383 realm (Fig S4) 384 385 We further measured the rich club coefficient, Φ(k), of the networks. Topological rich 386 club refers to the tendency of typically high degree nodes to connect with each other; when its 387 normalized counterpart, Φ norm (k), is above the baseline (greater than 1), then the subset of nodes 388 with degree greater than k is more densely interconnected compared to a random control with the 389 same degree distribution. The coefficient is not trivially connected to assortativity, since a 390 dissasortative network could still be rich club and vice versa (27). We found that in the modular 391 and transition range the weighted networks exhibited topological rich club behavior for a range 392 of intermediate degrees (binary networks were not rich club for τ transition ), whilst centralized 393 networks were not rich club for any degree (Fig S5).

395
For normal and lognormal networks, we evaluated the normalized weighted rich club 396 coefficient, Φ w,norm (k). Φ w,norm (k) above the baseline indicates that the edges of the rich club 397 nodes have larger weights compared to a random control network (a network with the same 398 topology but randomly reshuffled weights). We found that for the lognormal networks Φ w,norm (k) 399 did not deviate from 1; for the normal networks, Φ w,norm (k) similarly hovered around 1 in the 400 modular and centralized regimes. Hence for these states the rewiring process does not distribute 401 the larger weights preferentially to the high degree nodes. In general the data show that there is a 402 distinction between the topological and weighted measures; the former being above baseline for 403 a range of degrees, but not the latter.

405
This distinction becomes more accentuated for normal networks in the transition zone, 406 where Φ w,norm (k) tilts below baseline for nodes with larger degrees (k>20) (Fig 8). Taken 407 together, these topological and weighted coefficient profiles bode well with anatomical data on 408 the human brain. Specifically, van den Heuvel and colleagues (23)  Graph diffusion and other rewiring models 441 We showed how in graph-theoretical terms, the stability-plasticity dilemma can be 467 a wide range of diffusion rates and random rewiring probabilities (Fig 1). This is by and large 468 because diffusion shapes the initial random network into a more structured connectivity pattern 469 with high clustering coefficient while at the same time it maintains small path lengths ( 487 For smaller τ values the emerging network is modular (Fig 3), with dense connectivity within 488 clusters but sparse between them. For larger τ values, we obtain centralized topologies where a 489 subset of nodes are acting as hubs (Fig 4C). Both of these patterns are present in the brain. Furthermore, the centralized structures in the transition zone showed rich club behavior 515 (Fig 9), but the centralized networks for larger τ values did not (Fig S5). The combination of 516 centrality and rich club behavior is a signature of connectivity patterns in the brain (5). Rich club 517 connections in the brain could constitute its backbone, communicating information from diverse 518 functional modules. Empirical evidence has indicated that rich club connections extend over long 519 distances forming a high capacity structure that facilitates global communication (23).

520
Flexibility and stability of topology in and out of the transition zone 522 In the transition zone between modular and centralized structures, we observe the highest 523 variability in network topology. For example, in contrast to the τ values producing robustly 524 either modular or centralized topologies-when starting from a randomly connected one-, at the 525 center of the boundary between them, τ transition , a network's connectivity pattern cannot be 526 predefined before rewiring since it is equally likely that it will have either topology and anything 527 in-between in a continuum (Fig 5). Whereas this transition zone is similar in range for networks 528 with binary and normally distributed connection weights, it covers a considerably broader 529 parameter range for networks with lognormally distributed weights (Fig 4A).
530 531 Very small connectivity biases in the initial random network do not influence the final 532 topology of the rewired network (Fig 6A), however, larger biases lock the network at the 533 connectivity state it already is (Fig 6B). This 'point of no return' is not observed for networks 534 with lognormal weights. Moreover, lognormal networks show a more controllable behavior as 535 the diffusion rate moves away from the transition range in either direction: for diffusion values 536 smaller than the transition ones networks will most likely evolve to become modular and for 537 greater values they develop into a centralized topology irrespective of their initial connectivity 538 pattern (Fig 7). The graph Laplacian matrix 569 The graph Laplacian is defined as , where D is a diagonal matrix having the = -

585
; consequently, its eigenvalues are nonnegative and real and its eigenvectors ℒ ≥ 0, ∀ ∈ ℝ 586 form an orthonormal set. The spectral decomposition of the Laplacian matrix is , where is a diagonal matrix, having in its diagonal entries the eigenvalues Λ = ( 0 , 1 … ) 588 of ( ) and is an nXn matrix having as its columns the where is a function of n spatial variables and time variable , is a positive ℎ 1 ,…, 608 constant and is the Laplacian operator, or else the divergence of the gradient. With the adaptive rewiring algorithm and its variations employed in this paper, we seek 633 to probe the effects that a simple self-organizing rule will have on the properties of an initially 634 randomly connected network. The core algorithm is essentially the same as the one in Jarman 635 and colleagues (9), extended for application to both binary and weighted networks.

637
Starting with a random Erdös-Rényi type with vertices and connections-the latter guaranteeing the random network is fully connected (59)-the ( -1)) 639 rewiring algorithm proceeds as follows: 640 641 Step 1: Select, with uniform random probability, a vertex, from the set of vertices in 642 the graph that are of nonzero degree but also not connected to all other vertices ( ∈ | 0 < 643 . < -1) 644 645 Step 2: With probability p, go to step 2.1 (random rewiring), otherwise go to step 2.2 646 (heat diffusion rewiring). In both cases, select vertex from the set of vertices that are not Step 2.1: For both and the selection is random and uniform among the elements of 1 2 653 each set. 654 655 Step 2.2: Calculate the heat kernel, , of the adjacency matrix. is selected such that ℎ( ) 1 656 from all the vertices not connected to k, it is the one with the highest heat transfer with k. is 2 657 selected such that from all the vertices connected to k, it is the one with the lowest heat transfer 658 with k. In mathematical terms, we can express the statements above as follows: Step 3. Go back to Step 1 until k rewiring iterations are made.