Adaptive rewiring in nonuniform coupled oscillators

The anatomical connectivity of the brain network is shaped dynamically through structural plasticity (Butz, Wörgötter, & van Ooyen, 2009). A variety of structural plasticity mechanisms serve the adaptive role of accommodating the functional connectivity, the mutual statistical dependencies between the network components (Avena-Koenigsberger, Misic, & Sporns, 2018; Rubinov, Sporns, van Leeuwen, & Breakspear, 2009). A common principle underlying these mechanisms has become known as adaptive rewiring (Gong & van Leeuwen, 2003, 2004; Papadopoulos, Kim, Kurths, & Bassett, 2017). Adaptive rewiring implements the Hebbian principle of “what fires together, wires together” at the level of network dynamics (Bi & Poo, 2001; Hebb, 1949).

In networks of which the units may represent spiking model neurons (Kwok, Jurica, Raffone, & van Leeuwen, 2007) or neural mass oscillators (Rubinov et al., 2009), adaptive rewiring involves establishing connections between dynamically synchronized units at the expense of currently non-synchronized ones. Over time, adaptively rewiring networks dynamically evolve into complex architectures, showing the characteristics of small-worldness, modularity, and the rich-club effect (Gong & van Leeuwen, 2003, 2004; Hellrigel, Jarman, & van Leeuwen, 2019; Rubinov et al., 2009). Since small-worldness (Sporns & Zwi, 2004), modularity (Meunier, Lambiotte, & Bullmore, 2010), and the rich-club effect (van den Heuvel & Sporns, 2011) are characteristics of large-scale brain connectivity, adaptive rewiring could represent a driving force for the dynamic evolution of brain connectivity structure.

Adaptive rewiring models typically are highly simplified models, of which the units are identical nonlinear oscillators with uniform connection weights and coupling strengths (Rubinov et al., 2009). These simplifications severely reduce the functionality of such models. To implement information processes in these networks, it is vital to enable nonidentical oscillators and nonuniform connection strength. However, the question should be asked whether the properties of adaptive rewiring are robust to such nonuniformities. We therefore aimed to explore adaptive rewiring of coupled oscillators with nonuniform amplitude and coupling strengths and compare their evolution with that of uniform networks. The results will be crucial for the utility of adaptive rewiring in biological and artificial neural networks with information-processing functions such as pattern recognition and learning.

Most adaptive rewiring studies have represented network activity by a model known as coupled logistic maps (Kaneko, 1992). This choice was reached based on a succession of abstractions. Mass activity was described by chaotic attractor dynamics, shown in Figure S1A of the Supporting Information (Breakspear, Terry, & Friston, 2003). Further simplification was obtained via the Poincaré section, yielding the mapping in Figure S1B. This mapping, in turn, can approximately be described, minus the noise, by a logistic map (Figure S1C). Thus, the logistic map is the most straightforward possible abstract representation of neural mass activity (Rubinov et al., 2009). Note that in coupled map networks, the couplings provide noise perturbations to the oscillators, yielding mappings again more similar to Figure S1B.

In this equation, symbols ⊙ and ⊘ denote Hadamard (i.e., element-wise) multiplication and division, respectively. The right-hand side of Equation 3 constitutes the vector form of the logistic map, in which α is the vector of amplitudes and X_t is the vector of node’s activities at time t. The map is Hadamard-multiplied by a coupling term. In the coupling term, ε is the vector of coupling strengths, A_t is the connectivity matrix at time t, and 1_N denotes a vertical unit vector of size N. In the coupling term, A_tε is Hadamard-divided by A_t1_N, normalizing the former by the sum of the weights of the edges connected to each node. For binary networks, the term in the denominator counts the number of connections for each node. We consider only binary graphs for convenience. For adaptive rewiring in weighted networks, see Hellrigel et al. (2019).

Based on the network activity defined by Equation 3, adaptive rewiring takes the following form. After several updates to the network activity, a rewiring step is made. At each rewiring step, the connections of a random node are updated as follows: The node is disconnected from the neighbor most dissimilar in activity and is connected to the one most similar in activity among those it is presently unconnected to. Dissimilarity of two nodes at a given time is defined as the absolute value of the difference of their activity values. Note that although rewiring steps are local, the most dissimilar unconnected node is obtained through a global search. For algorithms using local or, instead, regional information for this purpose, see Jarman, Steur, Trengove, Tyukin, and van Leeuwen (2017); Jarman, Trengove, Steur, Tyukin, and van Leeuwen (2014).

All adaptive rewiring studies using coupled logistic maps have used fixed, uniform values for the amplitude parameters α and the coupling strengths ε. In logistic map networks with fixed architectures, these parameters have previously been allowed to vary in order to implement information-processing functions. In a perceptual organization model (van Leeuwen, Steyvers, & Nooter, 1997), sensory input function was realized by local modulation of oscillator amplitudes. The presence of sensory input to certain units brought their amplitude parameter values down to impose a more stable regime on the oscillators. As a result, these oscillators showed a bias to become mutually synchronized. The synchronization was understood as representing perceptual grouping. Grouping preferences followed the Gestalt proximity principle and showed spontaneous switching in case of stimulation with bistable patterns.

In a memory model (van Leeuwen & Raffone, 2001), connectivity parameter values were locally incremented to represent the presence of a memory trace. This led to synchronization biases along the incremented connections, resulting in spontaneous rehearsal and relearning of stored pattern information. These examples demonstrate that coupled maps are capable of performing cognitive functions based on ongoing dynamics.

The question is whether the nonidentical parameter settings involved in these functions would dovetail with adaptive rewiring. To answer this question, we partitioned networks into two subsets of units, a majority and a minority, each receiving a different parameter value for the amplitude α or the coupling strength ε. Fixing these parameters at different values involves a hardship test for adaptive rewiring. For instance, with one partition having increased coupling strength (or decreased amplitude), we might expect connections to be established preferentially within this partition and less within the other one. We may expect the opposite when coupling strength within the partition is decreased (or amplitude increased). In addition, we may expect fewer connections to be established between the partitions. This might induce network disintegration or otherwise interfere with the evolution of the network structure. On the other hand, if these perturbations fail to interfere with the emergence of complexity (i.e., small-wordness, modularity, and rich-club effect), adaptive rewiring could be used alongside pattern detection and learning in sparse neural networks. This has implications for the functioning of biological networks, as well as for the sparsification of artificial (i.e., deep) neural networks.

In what follows, the Method section describes details of the composition and initialization of the models, the rewiring algorithm, and the qualitative and quantitative measures of network structures used to describe, characterize, and compare models. In the Results section, we describe our findings, mainly that nonuniformity of parameters is shown not to interfere with the evolution of brain-like structure while giving rise to distinguishable network structures suitable for cognitive functions. A discussion and several concluding remarks end the paper.

An unlabeled binary graph (or network) G = (V, E) with n vertices (or nodes) and m edges (or couplings) is defined mathematically by a set of nodes V = {1, 2, …, n} and a set of edges E = {(i, j) ∈ V² : i is coupled with j}. E is also known as the adjacency list of the network. We use undirected graphs, that is, (i, j) ∈ E ⇔ (j, i) ∈ E, and self couplings are not allowed, that is, (i, i) ∉ E. The size of the set V (i.e., the number of its members) is denoted by |V|. From E, we may construct the adjacency matrix A of G as a square matrix of the size |V|, the elements of which can take values of 0 or 1. The element on its ith row and jth column (i.e., a_ij) is equal to 1 if (i, j) ∈ E and is 0 otherwise. Since G is undirected, A is symmetrical around the main diagonal, and since there are no self connections, its diagonal elements are zeros.

We may partition V into two nonoverlapping subsets, minority and majority, such that |V_minority| = |V| − |V_majority| and |V_minority| ≪ |V_majority|. The edges among members of these subsets form subgraphs within G and are henceforth called minority and majority partitions. A third subgraph comprises all of V but only edges between minority and majority nodes. Such a subgraph is called interpartition.

To each v_i ∈ V, an activation value is assigned according to Equation 3. The corresponding parameter values, that is, coupling strength ε and amplitude α, remain fixed in our model simulations. Models with identical parameter sets are called families. Five families of models are simulated, each with 10 model instantiations, comprising a total of 50 model instantiations. Each network is run for 20 million iterations. All the simulations and analyses are conducted in R programming language version 3.6.0 (R Core Team, 2019) using computational resources provided by VSC (Flemish Supercomputer Center).

In our models, all Gs have |V| = 300 nodes and 5,200 edges, a connectivity density providing robust evolution of small-world structure with uniform parameter setting (van den Berg, Gong, Breakspear, & van Leeuwen, 2012). The model structure is initialized by randomly assigning 5,200 × 2 = 10,400 values “1” symmetrically to nondiagonal entries of A, and zeros to the remaining entries. Each node in the network is randomly and independently assigned an initial value, uniformly distributed between 0 and 1, that is, x_i1 ∼ Unif (0, 1).

Previous adaptive rewiring studies have been using values of α and ε in the ranges of [1.7–1.9] and [0.3–0.5], respectively (Gong & van Leeuwen, 2003; Hellrigel et al., 2019; van den Berg & van Leeuwen, 2004). Here the midpoints of these ranges, that is, α = 1.8 and ε = 0.4, are used for the parameters in the baseline (BL) condition. In the BL condition, all nodes have identical parameter values.

The same applies to the majority (250 nodes) of the other conditions. However, depending on the condition, the minority subset (i.e., the first 50) of nodes could have either lowered or increased values of either the α or ε parameters. As shown in Figure S2, higher values of the amplitude α tend to yield greater divergence in activity values; reducing the coupling parameter ε has a similar effect (Hellrigel et al., 2019). Conditions with lowered values of α are called less chaotic (LC), and those with increased values more chaotic (MC); conditions with lowered ε values are called sub-coupled (SC) and those with increased values hyper-coupled (HC). While keeping the parameters of the majority at the baseline level (α_i∈51:300 = 1.8, ε_i∈51:300 = 0.4), five different combinations of parameters were assigned to the minorities, each combination called a “family”: The baseline family (BL; α_i∈1:50 = 1.8, ε_i∈1:50 = 0.4), and the families with less chaotic minority (LC; α_i∈1:50 = 1.7, ε_i∈1:50 = 0.4), more chaotic family (MC; α_i∈1:50 = 1.9, ε_i∈1:50 = 0.4), sub-coupled minorities (SC; α_i∈1:50 = 1.8, ε_i∈1:50 = 0.3), and hyper-coupled minority (HC; α_i∈1:50 = 1.8, ε_i∈1:50 = 0.5). In the Results section, we identify model instantiations by the two capitals indicating their family, together with a serial number [1–10], for example, BL7, MC10. The 10 model instantiations within each condition are run with different initializations, which are identical across conditions to allow matched comparison between families.

A rewiring takes place after every 20 updates of the logistic maps, meaning that 1 million rewiring attempts are performed over the 20 million updates of the model. At each rewiring attempt, at time t, a node i is selected randomly from V, a vector of its distance from other nodes is calculated as d_i,t = |X_t − x_i,t1_|V||, and another vector of similarities is defined as s_i,t = 1_|V| − d_i,t.

The state of each model at any given time t is described by adjacency matrix A_t (henceforth, “anatomical connectivity”), which is subject to adaptive rewiring, and the vector of activation values, X_t. A model’s “functional connectivity” at t (represented by the |V| × |V| matrix F_t) is defined by the momentary pairwise absolute differences of its node activation values.

Network structure can be qualitatively assessed by means of visual inspection of the wiring diagram or the adjacency matrix. Using the package seriation (Hahsler, Hornik, & Buchta, 2008), the adjacency matrix is serialized by ordering rows and columns according to the projection of the matrix on its first principal component. Seriation maximizes the visual identifiability of modules within the network.

After each rewiring attempt, we calculate six network connectivity measures (Costa, Rodrigues, Travieso, & Villas Boas, 2007): edge density, clustering coefficient, average path length, small-world index, modularity, and assortativity. Furthermore, we calculate an additional measure, namely, the rich-club coefficient, after the final rewiring. All measures are calculated separately for the entire graph and the three subgraphs (viz., the minority, majority, and interpartition subgraphs).

Since the total number of edges remains the same during the adaptive rewiring, this coefficient gives an indication of how strongly each partition has attracted new nodes at every rewiring step.

The normalized index makes the comparison of networks with different sizes and densities possible. Since all networks modeled in this study start off with random networks of equal size and density, for computational reasons, the non-normalized small-world index (Equation 9) is calculated and reported.

This measure requires an a priori labeling of nodes that defines the communities to which the nodes belong. A variety of algorithms have been suggested to discover module, or communities, within a network so that the value of Q is maximized (for a review, cf. Zhang, Ma, Zhang, Sun, & Yan, 2018). The communities discovered by these algorithms can thus be used as labels for calculating modularity of the network. In line with Clauset, Newman, and Moore (2004), we use the fast greedy algorithm to optimally detect communities and thus calculate the modularity based on community membership of the nodes. The igraph package was used to automatically detect communities and calculate modularity.

The assortativity coefficient is a measure of homophily in networks that indicates the preferences of nodes to connect to alike nodes. The likeness can be imposed externally, for example, by assigning categories to the nodes using labels (known as nominal assortativity), or by internal criteria such as node degrees (degree assortativity). Degree assortativity is defined as the Pearson correlation coefficient of degrees of connected nodes, thus taking values in the range [−1, 1].

The absolute value of the rich-club coefficient is hard to interpret and is not comparable among networks of different sizes, densities, and degree distributions. Hence, this coefficient is often normalized by the average rich-club coefficient of random networks of the same size with similar degree sequence. For each model (and its minority and majority partitions), we simulated 200 such networks, for each club size k, and normalized the values of the rich-club coefficients of the network by the average rich-club coefficient among the random networks. Moreover, for each k, we statistically tested whether the non-normalized RC is significantly larger than the rich-club coefficients of randomly generated networks. Since the distribution of RCs of random networks were non-Gaussian, the conventional one-sample t test was not applicable. Instead, we performed one-sample Wilcoxon signed rank test and set α = 0.01 as the significance level for the p values.

For a certain k, a normalized rich-club coefficient larger than 1 indicates that nodes with degree k tend to connect to the “rich” nodes (i.e., those with degrees equal to or greater than k), thus forming “rich clubs.” Conversely, RC_norm(k) < 1 implies that the nodes with degree k show a tendency to connect to nodes with lower degrees. Finally, RC_norm(k) = 1 suggests that nodes with degree k show no preference to connect to nodes with lower or higher degrees.

To compare network families with each other, we assume, in line with Berlingerio, Koutra, Eliassi-Rad, and Faloutsos (2012), that the structural information embedded in networks can be summarized by the distributions of local network measures. Comparison of networks is thus reduced to comparing these distributions. To obtain measures of distributional distances, we use the NetSimile method (Berlingerio et al., 2012) and the Heller–Heller–Gorfine (HHG) algorithm (Heller, Heller, & Gorfine, 2013), both of which are discussed in the Supporting Information. In short, NetSimile indicates the degree of dissimilarity between the distributions attributed to the networks being compared via comparing signature vectors that encapsulate the distributions by their summary statistics.

HHG, on the other hand, provides p values for a test of independence among the distributions; a small p value derived from HHG (e.g., below the threshold of α = 0.05) provides evidence to reject the null hypothesis of distributional independence. Thus, loosely speaking, HHG p value can be regarded as an indicator for dissimilarity; a higher value of this measure entails a smaller “resemblance” (or dependence) between the networks. Nevertheless, interpreting HHG p values as such a measure is somewhat unorthodox and is hardly meaningful unless used along with another dissimilarity measure such as NetSimile.

We first make pairwise comparisons among the 1,225 unique pairs of model instantiations at their final state after 1 million rewiring attempts. Having quantitative measures for dissimilarities among the networks, we quantify the within-family resemblances and between-family contrasts among the models. Finally, we define a score for family distinction in order to compare how families vary with respect to this measure.

Using the seriation and igraph packages, we plot the raw (unserialized) and ordered (serialized) adjacency matrices and the wiring diagrams of the networks in their final state. In the plots, the minority and majority nodes are colored sky blue and pink, respectively. In both matrix visualizations and wiring diagrams, the within-minority and within-majority edges are colored blue and red, respectively. The interpartition edges, connecting nodes of minority subset to nodes of majority, are colored green. Although there are variations among models, either within- or between-families, in all models (except for the MC2, MC3, SC1, and SC3; see below), several densely coupled sets of nodes, that is, modules, emerged. These modules can be identified as squares in the serialized adjacency matrices. Figure 1 shows two representative networks per family. As evident in the plots, the networks manifest a range of different structures. Yet, similarities can be observed among models belonging to the same families.

The baseline models (BL2 and BL8 in Figure 1) typically include three densely coupled modules and a few larger, sparser sets of nodes. The modules are not isolated from the rest of the network, as intermodular edges keep them connected to other nodes. The unserialized adjacency matrices show that the density of edges is relatively uniform over subsets of nodes.

In the HC family, wherein the minority nodes have relatively higher coupling strengths, the edge density is higher in the minority subgraph. Moreover, the HC family networks have more distinct modules than other families. HC2 (Figure 1), for instance, has only two modules, both of which are highly connected. Higher edge density in the minority subset and highly distinct modules are also observed for the LC family, which has lower amplitude in the minority nodes (LC5 and LC8, Figure 1). The similarity between the HC and LC models was expected, as the nodes with lower amplitudes and higher coupling strengths can synchronize more easily. However, it is worth noting that the effect is not limited to the minority nodes; highly connected modules also emerge among the majority nodes.

In the MC family, the edge density in the minority partition is reduced. Moreover, the higher level of amplitude for minority nodes resulted in highly connected modules among the majority nodes (see MC4 in Figure 1). In the SC family, the edge density of the minority (which had lower coupling strengths) is lower than that of the majority. Additionally, the lower coupling strength of the minority prevented minority nodes from forming modules, and they were absorbed into modules formed mainly by the majority nodes. See SC5 and SC7 in Figure 1.

The evolution of clustering coefficient, modularity, edge density, small-world index, degree assortativity, and average path length were plotted for the structural network connectivity of all models and their subnetworks (viz., within-minority, within-majority, and interpartition). Figures 2–6 show these plots for the minority, majority, and whole networks, grouped per family. Furthermore, the normalized rich-club coefficient of the final states of the minority, majority, and whole networks are plotted in Figure 7.

Let us first consider the evolution of network statistics for the whole network. As evident in the plots, modularity, clustering coefficients, and small-world index of all models increase noticeably after 60,000 rewiring attempts (roughly by a factor of 4.7, 5.3, and 4.6, respectively) and degree assortativity approaches 0.53. Meanwhile, the average path length increases no more than 15%. The network statistics fluctuate afterward, and the fluctuation is more pronounced for degree assortativity. The majority subnetworks show similar trends. In the minority subnetworks—although with higher fluctuations and instability in HC, LC, and SC—we observe an increase in these statistics after 60,000 rewiring attempts. The network statistics of minority subnetworks of MC and BL families resemble those of the whole networks in the same families. Yet, clustering coefficient and small-world index show higher growth in the minority subnetworks (more than six- and sevenfold, respectively). The exception to these improvements is the modularity, which increases only moderately (roughly 35%) in this subnetwork. Given the finite number of nodes in this partition, this was not unexpected; a larger subnetwork would likely have shown more improvement in modularity as well. The edge densities in the minority subnetworks, although unstable, are often large, roughly 2.4 times that of the whole network. In line with elevated edge density, average path length within the minority subgraph drops by almost 7%. This means that the minority nodes tend to connect more strongly to each other than to other nodes. Tables 1–3 summarize the mean and standard deviation of changes in network statistics after 60,000 rewiring attempts for the whole network and the majority and minority subnetworks.

Figure 7 shows the normalized rich-club coefficient RC_norm(k) of the minority, majority, and whole networks, grouped by families, as a function of club size k. RC_norm(k) above 1 (dashed line) indicates rich clubs. Values significantly larger than 1 (p < 0.01, based on one-sample Wilcoxon signed rank test) are marked by solid circles. Despite differences among families, they all show significant rich-club structures for larger club sizes, both in the whole network and in the majority partition. Let us first consider the plots for the whole networks, depicted in the left column of this figure. The HC family has clubs of sizes 30 < k < 70 (involving 10–23% of the nodes), most of them significant, with consistent RC_norm values among members. In the LC family, all models form rich clubs with 25 < k < 55 (roughly 8–18% of the nodes), and some models form larger rich clubs up to the sizes of 65–85 (roughly 22–28% of the nodes). Other families, that is, MC, BL, and SC, form clubs with sizes ranging from 30 to roughly 65 (10–22% of the nodes), with relatively higher (but less consistent) RC_norm values compared with HC.

In the majority partition (middle column of Figure 7), we consistently observe rich clubs larger than 25, involving 10% of the nodes within that partition. The HC family has rich clubs as large as 40 (16% of nodes), and other families form even larger rich clubs, with sizes spanning 45 to 60 (making up 18–24% of the nodes). In general, we observe more diverse values among models for the rich-club coefficients in the majority partition. Finally, in the minority partition (the right column of the same figure), we observe a relatively remarkable—and consistent—rich-club effect in the HC family (12 < k < 30; 24–60% of nodes), while other families have relatively larger rich-club coefficients for smaller club sizes (5 < k < 20; 10–40% of nodes). Some models of the LC family form much larger rich clubs with sizes up to 40 (i.e., 80% of the nodes).

As can be observed in the evolution plots, four models (viz., MC2, MC3, SC1, and SC3) stop evolving before 10,000 rewiring attempts. Visual inspection reveals that, at some point in their evolution, one node reaches maximum degree and is connected to all other nodes (cf. Figure 8). This indicates that these models reached the end point of a pathological development, something we chose not to prevent in our algorithm. Consequently, the rewiring algorithm naturally terminates, as division by 0 looms in Equation 3. Terminated models failed to form any modules and were omitted from family-wise comparisons. We consider these incidental losses as a sign that robustness against perturbation is not always guaranteed for our models.

The matrices of Contrast^N are plotted as upper triangular matrices in Figure 10 using the corrplot R package (Wei & Simko, 2021). The cell colors, coded similarly to the heat maps, denote average contrast measures derived from the NetSimile algorithm while average HHG p values (i.e., $ContrastHHGN$⁠) are indicated in each cell. The HC family manifests the least within-family contrast. Based on the HHG test of multivariate independence, except for HC-HC and HC-LC family pairs, no conclusive evidence exists for distributional dependence among families. Finally, as can be seen in this figure, the within- and between-family NetSimile contrasts of both structural and functional connectivity networks show similar patterns. More specifically, both the structural and the functional connectivity of HC-BL, HC-LC, BL-BL, BL-LC, and LC-LC all share close contrast values compared with other family pairs. This is also the case for the BL-SC and LC-SC pairs.

The differentiation scores for structural and functional connectivity between families are plotted in Figure 11. Differentiation values above 1 (dashed line) imply that the within-family resemblance of network structures of family f_i is higher than the average resemblance of its members to the members of other families. We observe elevated differentiation in both structural and functional networks of HC, BL, and LC. This measure is remarkably higher for the structural networks of the HC family while barely exceeding the threshold for the functional network of the same family.

Finally, in order to have both family resemblance and family differentiation in a single frame, we summarized their values in the graphs shown in Figure 12. In these graphs, individual nodes represent families of models. Edge color and size code between-family contrast and node color captures within-family contrast. The size of each node is proportional to the value of the differentiation score of its corresponding family. The families with Differentiation^N(f_i) > 1 are marked with asterisks. It can be noticed that in both structural and functional networks, the families with differentiation scores larger than 1 (i.e., HC, BL, and LC) have lower within-family contrast values. Moreover, the pairwise contrasts among these families (i.e., HC-BL, HC-LC, and BL-LC) are lower than any other family pair. This suggests that models with increased amplitude and decreased coupling strengths lead to structures less resembling the baseline family.

We investigated the effect of nonuniformity of oscillator amplitude and connection strength on the evolution of network structure under adaptive rewiring. The adaptive rewiring was driven by synchronization in coupled logistic maps. The amplitude and coupling parameters of the oscillators govern their synchronization behavior. It has amply been shown that networks evolve to complex, brain-like structures when these parameters were fixed to uniform values. We replicated this behavior for our baseline parameter setting. Additionally, we studied how nonuniform parametrization gets reflected in the evolving network structure and whether adaptive rewiring is robust to these perturbations.

To this aim, a minority subset of network units had either the amplitude reduced (LC) or increased (MC), or the connectivity reduced (SC) or increased (HC). These perturbations may be compared with perceptual input (modeled by increased or decreased divergence in node activities) or memory processes (modeled by varying degrees of coupling). In all these conditions, adaptive rewiring establishes the evolution of random networks into complex structures. These structures are modular and small-world networks and show a fivefold increase in the clustering coefficient, at the cost of an increase of less than 15% in average path length compared with the original random networks. Moreover, the models form rich clubs as large as 8–28% of network size. This evolution is in accord with previous adaptive rewiring studies. Adaptive rewiring, therefore, is generally robust against the symmetry-breaking perturbations of system parameters and thus can accommodate—and integrate—function-specific modules. Note, however, that robustness is not always guaranteed. In a small proportion of our models (four cases in total), we observed pathological development in the perturbed models, resulting in network evolution breakdown.

On top of these findings about the whole network, the subsets of our models also developed brain-like structures. More specifically, across families, the majority subnetwork shows improved properties akin to the whole network. The minority subnetwork of perturbed models—with an exception for modularity—manifests these properties to a much greater degree, especially in the relative sizes of its rich clubs. The latter observation is indicative of more specialized structures in the minority subnetwork as a direct result of heterogeneity in model parameters.

Despite the overall consistency in structural improvements of networks, we observed considerable variability among models, both in their evolution and in the final states of structural and functional connectivity. These variations are partly due to the family-specific parameterization of models and partly due to the random initialization of the models. We quantified the within- and between-family contrasts and defined a measure of family differentiation to score, on average, how well models belonging to one family distinguish from others.

We found that the baseline network family and families with increased coupling strengths or decreased amplitude (respectively, BL, HC, and LC) differentiate themselves most from other families, in both structural and functional connectivity. The differentiation score is the highest for the structural connectivity of the HC family. Moreover, we observe a relatively high resemblance between HC and LC for both structural and functional connectivity. Based on the HHG test, there is relatively strong evidence (p = 0.02) for distributional dependence between structural networks of HC and LC. This means that elevated coupling strength and reduced amplitude have similar effects on both the structural and the functional connectivity of adaptively rewiring oscillator models.

Altogether, perturbation to the parameters of coupled oscillators yields structural and functional differences that are fundamental for implementing cognitive functions in evolving networks. Moreover, for both the structural and the functional connectivity, perturbation leads to differentiation from the baseline. Different perturbations show specific differentiations. From a cognitive neuroscience perspective, this implies that functional and structural differentiation can be used to identify functional components in a network, in accord with the use of structural and functional connectivity measures in neuroimaging.

Further research may extend our findings in several ways. First, the incidental occurrence of network evolution breakdown was unexpected. Technically speaking, the matrix algebraic implementation of coupled logistic maps is sensitive to minor computational errors such as undefined division for one node (which may stem from pathological development of models), as was the case for MC2, MC3, SC1, and SC3. A different implementation could have prevented breakdown by isolating the problematic node and proceeding with the adaptive rewiring algorithm for the remaining network. However, rather than searching for an even more robust version of our algorithm, we choose to keep the current implementation for the sake of consistency with previous work, and as a warning signal that pathologies may arise as a by-product of neural network evolution. Future work may focus more specifically on such pathological development, in relation to developmental disorders in the human brain. This may require a systematic study of the probability distribution of network breakdown under a range of parametrization conditions, network sizes, and connectivity densities to find factors contributing to the pathological development of adaptively rewiring networks.

Second, a broader range of parameter variation conditions than currently imposed should be attempted. For instance, amplitude and coupling strength parameters could both deviate from baseline values simultaneously, either for the same subset of nodes or for two (overlapping or nonoverlapping) subsets. Random and patterned deviations of parameters can be studied in large-scale systems to implement perceptual and memory functions. Ultimately, the aim is to have these functions implemented in a network that simultaneously maintains its optimal structure.

This study was limited to binary, unweighted, and undirected networks. The effect of nonuniform parameters of logistic maps can be studied, via systematic search, on weighted networks with various edge weight distributions, akin to Hellrigel et al. (2019). For directed networks, see Rentzeperis, Laquitaine, and van Leeuwen (2021).

We followed the strategy to provide the simplest possible model of brain structure and function in order to avoid stacking arbitrary assumptions. Simplification is inevitable in modeling; “all models are wrong—but some are useful.” (Attributed to George Box.) For those who consider the model to pose severe limitations on the generalization of our findings to neurobiological systems, more realistic neural mass models instead of coupled maps could be a viable solution, in particular, ones that have facilities for studying the effects of traveling waves and spatial embedding of the network topology (Calvo Tapia, Makarov, & van Leeuwen, 2020) on network evolution, or apply phase–amplitude–frequency coupling (Chehelcheraghi, Nakatani, Steur, & van Leeuwen, 2016; Chehelcheraghi, van Leeuwen, Steur, & Nakatani, 2017) to networks of neural mass activity (Deschle, Ignacio Gossn, Tewarie, Schelter, & Daffertshofer, 2021).

Finally, despite the current limitations, our study has set an essential first step in the development of adaptively rewiring networks capable of pattern recognition and learning. Showing the principled possibility of such networks offers scope for the study of the developing, functional brain as well as for applying adaptive rewiring to sparsify artificial neural networks (Gale, Elsen, & Hooker, 2019).

Part of the research leading to this review was supported by an Odysseus grant (G.0003.12) from the Flemish Organization for Science (FWO) to Cees van Leeuwen. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by FWO and the Flemish Government. The authors would like to thank Ilias Rentzeperis for his feedback and comments on parts of the results.

The reproducible scripts used in this study are available online on the study’s repository on the Open Science Framework (https://osf.io/625d8; Haqiqatkhah & van Leeuwen, 2021). This repository also includes all model files of this study (500 files amounting to 101 GB of data, generated using VSC computational resources) and additional plots of various qualitative and quantitative network measures.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00211.

MohammadHossein Manuel Haqiqatkhah: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Software; Visualization; Writing – original draft. Cees van Leeuwen: Conceptualization; Funding acquisition; Investigation; Project administration; Resources; Supervision; Writing – review & editing.

Cees van Leeuwen, Fonds Wetenschappelijk Onderzoek (https://dx.doi.org/10.13039/501100003130), Award ID: G.0003.12.

Structural plasticity:

The brain’s ability to modify its structure by means of forming and pruning neuron-to-neuron connections or changes in synaptic strength based on neural activities.

Neural mass oscillators:

Simplified models of coarse-grained mesoscopic or macroscopic oscillatory activities of often large coupled populations of excitatory and inhibitory neurons and their synapses.

Chaotic attractor:

An attractor (a set of values or states that a dynamic system tends to evolve toward) that is highly sensitive to the initial conditions of the system.

Poincaré section:

A lower dimension in a continuous dynamic system’s state space with which the periodic orbit of the system’s flow intersects.

Limit-cycle attractors:

An isolated periodic orbit or closed trajectory of a continuous dynamical system that has a two-dimensional state space.

Perceptual organization:

A process that, by grouping visual elements (such as corners and edges), facilitates determining the meaning of the visual input as a whole.

Degree sequence:

A sorted list of the degrees of nodes (i.e., the number of connections each node has to others) in a network.

Traveling waves:

Spatially coherent neural activities that propagate across brain regions at different scales and that are hypothesized to coordinate brain circuits at longer ranges.

https://doi.org/osf.io/625d8