Using Pareto optimality to explore the topology and dynamics of the human connectome

Graph theory has provided a key mathematical framework to analyse the architecture of human brain networks. This architecture embodies an inherently complex relationship between connection topology, the spatial arrangement of network elements, and the resulting network cost and functional performance. An exploration of these interacting factors and driving forces may reveal salient network features that are critically important for shaping and constraining the brain's topological organization and its evolvability. Several studies have pointed to an economic balance between network cost and network efficiency with networks organized in an ‘economical’ small-world favouring high communication efficiency at a low wiring cost. In this study, we define and explore a network morphospace in order to characterize different aspects of communication efficiency in human brain networks. Using a multi-objective evolutionary approach that approximates a Pareto-optimal set within the morphospace, we investigate the capacity of anatomical brain networks to evolve towards topologies that exhibit optimal information processing features while preserving network cost. This approach allows us to investigate network topologies that emerge under specific selection pressures, thus providing some insight into the selectional forces that may have shaped the network architecture of existing human brains.


Introduction
The emergence of network science and the increasing availability of brain connectivity data have recently opened up a network-based perspective on brain function. Studies in this area use diverse mathematical and computational tools to study the architecture of brain networks and its role in the dynamics of information processing [1,2]. In the human brain, the use of diffusion imaging techniques to detect white matter pathways connecting anatomical brain regions has enabled the mapping and analysis of structural brain networks. While descriptive studies have identified a number of characteristic topological attributes [3 -5], the fundamental selectional forces and factors that have shaped human brain network topology remain poorly understood. Three candidate factors explored here are network cost [6][7][8], network communication efficiency [9][10][11] and dynamic complexity [12][13][14].
It has been recognized for over a century that one fundamental factor shaping neuronal morphology and connectivity is that brain networks are embedded in space. A major consequence of spatial embedding is that the generation, maintenance and use of connections incur a cost, as connectivity consumes various resources such as wiring length [6][7][8] and metabolic energy [15]. biological system, these resources are limited and the need to conserve such resources places strong constraints on the topology of the system. Underscoring the importance of spatial embedding [16,17], many studies have shown that the topology of connections between individual neurons is strongly influenced by the spatial distance between them [18] and that cortical regions that are spatially close have high probability of being connected to each other [19][20][21].
However, spatial distance alone is insufficient to fully explain the connectivity patterns observed in brain networks. Structural brain networks are characterized by the existence of specific long-range pathways [22], highly connected regions (hubs) [23] and community structure [24,25], features that violate the concept of wiring minimization. In general, these connectional attributes are found in systems that achieve highly efficient global communication and whose components are highly clustered, e.g. small-world networks [26]. Many studies have shown that human brain networks are organized in an economical small-world manner, which tends to minimize wiring costs while supporting a few long-range connections that are thought to ensure high efficiency in global communication [22,27,28].
Another aspect of network organization relates to the brain's capacity to support a great diversity of dynamic patterns which are highly complex and essential to sustain a large number of competing functional demands. This diversity of dynamic patterns has been conceptualized as a 'functional repertoire' of network states that enables flexibility across a broad range of cognitive functions [29]. From a network perspective, this aspect points to the importance of local or specialized information processing in the cortex as well as the integration of information between different specialized regions [12,13]. It has been suggested that both aspects of information processing, integration and segregation, underlie the complex dynamics taking place in the network and that the patterns of structural connectivity found in the brain promote such complex dynamics [1,12,30].
None of these three factors alone is sufficient to account for all aspects of human brain network architecture. Instead, this architecture appears to represent a trade-off between these (and possibly other) competing factors, enabling economic information processing within a small-world topology [31]. Here, we explore the extent to which the structure of the human cortical brain network is optimally organized in order to achieve efficient and economical information processing. By defining a network morphospace and a multi-objective evolutionary algorithm that operates on the morphospace, we investigated the capacity of brain networks to evolve towards distinct biologically feasible topologies. This approach allows us to address important questions about what topological features emerge as a result of applying different types of selection pressure on the evolution of brain networks. For instance, how different are the networks selected for efficient information processing from networks selecting for certain dynamical properties? Furthermore, the analysis of a brain network morphospace aids in the understanding of actual brain network structure and provides a framework to study the structural variations to which brain networks are subject, due to individual differences or neurological degeneration.
Our analysis was performed over three different structural brain networks, obtained independently and from different imaging techniques. The analysis proceeds in three steps. First, we study the behaviour of our measures of information processing and dynamical complexity as the empirical networks are rewired towards two null models, namely, randomized and latticized networks. Both null models preserve crucial features of the corresponding empirical brain networks, such as network size, density, degree sequence and wiring cost. Second, we create a set of proximal networks by minimally perturbing the structure of the empirical networks. These proximal networks allow us to (i) quantify the sensitivity of our measures to small structural perturbations and (ii) depict the distribution of the proximal morphospace. Third, we use a multi-objective evolutionary algorithm to do a local exploration of the efficiencycomplexity morphospace of brain-like networks. The aim of a local exploration is to investigate alternative biologically feasible topologies for structural brain networks. Through this kind of analysis, we are able to portray how and how much a subregion of the morphospace (the region surrounding the empirical brain network) is filled; this in turn, can provide a picture of the underlying rules and constraints pervading the organization of brain networks.

Material and methods (a) Graph theory
To study brain connectivity, we apply methods from a branch of mathematics called graph theory [32 -34]. In the context of graph theory, an anatomical brain network with N interconnected neural elements is modelled as a graph G ¼ (V, E), where V is the set of vertices (or nodes) representing brain regions and E is the set of edges (links), representing white matter pathways. The size of a network is given by the number of nodes and connections composing the network, whereas the network density is defined as the number of existing connections divided by the maximum possible number of connections that the network can support. Formally, G is described by an N Â N adjacency matrix A G ¼ fa ij g, where a ij ¼ 1 if nodes i and j are connected and a ij ¼ 0 otherwise. In addition, G is a weighted graph if there is a scalar associated with every connection, such that if a ij ¼ 1, then there is a non-zero weight w ij assigned to the connection fi, jg. In the case of anatomical brain networks, we focus on two weighted matrices associated to the connections: a matrix of fibre densities W G ¼ fw ij g and a matrix of fibre lengths L G ¼ fl ij g estimated as in [3]. Thus, the degree of a node is given by the number of edges incident to the node, whereas the weighted degree is given by the sum of the fibre densities of the edges incident to the node. Finally, in this study, all the connections of anatomical brain networks are undirected, that is a ij ¼ a ji for all pairs fi,jg and thus A G , W G and L G are all symmetric matrices.

(b) Brain networks
Our analyses were carried out over three anatomical human brain networks (labelled LAU1, LAU2 and UTR; datasets are provided in the electronic supplementary material) constructed from data acquired independently in different imaging centres, using different acquisition protocols and different subject cohorts.

(i) LAU1
Five healthy right-handed male subjects (mean age 29.4 years, s.d. 3.4) were scanned on a 3-T Philips Achieva scanner. A highresolution T1-weighted gradient echo sequence was acquired in a matrix of 512 Â 512 Â 128 voxels of isotropic 1 mm resolution. Diffusion spectrum imaging (DSI) was performed using a diffusion-weighted single-shot echoplanar imaging sequence (TR ¼ 4200 ms; TE ¼ 89 ms) encoding 129 diffusion directions rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 369: 20130530 over a hemisphere. The maximum diffusion gradient intensity was 80 mT m 21 , the gradient duration was 32.5 ms and the diffusion time was 43.5 ms, yielding a maximal b-value of 9000 s mm 22 . The acquisition matrix was 112 Â 112, with an in-plane resolution of 2 Â 2 mm. Following diffusion spectrum and T1-weighted MRI acquisitions, the segmented grey matter was partitioned into 998 regions of interest (ROIs). Following white matter tractography, connectivity was aggregated across all voxels within each of the 998 ROIs. Further details are available elsewhere [3].

(ii) LAU2
Forty healthy subjects (24 males and 16 females, 25.3 + 4.9 years old) underwent an MRI session on a 3-T Siemens Trio scanner with a 32-channel head-coil. The magnetization-prepared rapid gradient-echo (MPRAGE) sequence was 1 mm in-plane resolution and 1.2 mm slice thickness. The DSI sequence included 128 diffusion-weighted volumes þ 1 reference b_0 volume, maximum b-value of 8000 s mm 22 and 2.2 Â 2.2 Â 3.0 mm voxel size. The echo planar imaging (EPI) sequence was 3.3 mm in-plane resolution and 0.3 mm slice thickness with TR 1920 ms. DSI and MPRAGE data were processed using the Connectome Mapping Toolkit [35]. Segmentation of grey and white matter was based on MPRAGE volumes. Cerebral cortex was parcellated into 1000 equally sized ROIs [36] followed by whole-brain streamline tractography [37].

(iii) UTR
Imaging data were acquired from 25 subjects (17 males and eight females, 29.4 + 7.7 years old). Diffusion-weighted imaging (DWI) was performed at 3 T, with two sets of 30 weighted diffusion scans (b ¼ 1000 s mm 22 ), each set consisting of five unweighted B0 scans (b ¼ 0 s mm 22 ) and 30 weighted scans (SENSE, p-reduction 3; gradient set of 30 weighting directions, TR ¼ 7035 ms, TE ¼ 68 ms, EPI factor 35; FOV 240 Â 240 mm, 2 mm isotropic, 75 slices, second diffusion set acquired with a reversed k-space readout). Preprocessing of the DWI involved the following steps: (i) diffusion images were realigned, corrected for eddy currents and susceptibility distortions; (ii) diffusion profiles were fitted with a single tensor and deterministic streamline tractography was used to reconstruct streamlines; and (iii) streamlines were used to build subject-specific structural brain networks among 1170 equally sized randomly partitioned cortical parcels (nodes). For a detailed description see [38].
For all three datasets, all subsequent analyses and modelling were carried out on group consensus matrices, built by averaging over all existing connections (expressed as fibre densities) that were present in at least 25% of participants in each dataset. For this study, we limit the analysis to networks containing nodes and connections in the right hemisphere of the brain, for two reasons. First, inter-hemispheric connections are less reliably captured by diffusion imaging and more difficult to reconstruct with tractography [39]. Second, the computational cost of running the analyses and simulations proposed here in whole-brain networks was prohibitive.

(c) Metrics of network performance
Complex network analysis has provided various metrics that aim to characterize different aspects of network topology [32,33]. Here, we selected four measures that jointly capture the performance of a network at combining integrated and segregated information processing in an economical manner.

(i) Wiring cost
This measure quantifies the cost of making and maintaining anatomical connections between neurons [8,15]. By assuming that the wiring cost is proportional to the wiring volume [28,31], we can express the cost of a single connection fi, jg as the product between its fibre density and length. Then, the total wiring cost of a network with N nodes is given by cost ¼ P N i,j w ij l ij : (ii) Efficiency of information processing We approach the measurement of efficiency of information processing from the perspective of two different communication schemes, one based on the routing of information in a network, and the other one based on the diffusion of information within the network [11].

Routing efficiency
In this work, the measure E glob defined in [9] is referred to as E rout [11]. This measure is based on the shortest path length matrix w ¼ [w ij ] where the distance between a pair of nodes is computed in terms of the inverse of the fibre densities of the connections. Then, the routing efficiency is computed as follows:

Diffusion efficiency
We start by defining a transition matrix P as the matrix whose elements p ij represent the probability of a random walker going from node i to node j in one step. If the transition probabilities are proportional to the fibre density of the connections, p ij ¼ w ij /k i , where k i is the weighted degree of node i. Given a transition matrix, the mean first passage time (MFPT) between node i and j is defined as the average number of steps it takes a random walker starting at node i, to arrive at node j for the first time [40]. If the network is connected and has no self-connections, the MFPT between any pair of nodes is finite and can be computed as follows: where the vector v is the left eigenvector associated to the eigenvalue of value unity; Z ¼ [z jj ] is the fundamental matrix, computed as Z ¼ (I 2 P þ W ) 21 , where I is the N Â N identity matrix, P is the transition matrix and W is an N Â N matrix with each column being the vector v such that 8j W ij ¼ v i . Diffusion efficiency is then defined as follows [11]: (iii) Neural complexity Neural complexity (C N ) is a measure that captures the coexistence of functional segregation and functional integration in a neural system [12]. C N is a statistical measure of the dynamics of the system defined in terms of the mutual information between subsystems. C N was originally defined in terms of the integration associated with a system of n neural components and a stationary stochastic process X(t) ; {X i (t)ji ¼ 1, 2, . . . , n}, where X i (t) represents the activity of the ith neural component at time t. Integration is computed as I ¼ P n i¼1 H i À H, where H is the entropy of the entire system and H i is the entropy of the ith individual component. Then, neural complexity is defined as follows: where k . l k denotes an average over all n k subsystems of size k. C N tends to be low for systems whose components are either statistically independent or highly dependent; conversely, C N is high for systems whose components are (on average) independent in small subsets and increasingly dependent in subsets of increasing size.
rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 369: 20130530 Computation of C N Note that in the definition of C N given above, the number of subsystems of size k ¼ 1, 2,. . ., n 2 1 increases exponentially with n as 2 n . The subsequent combinatorial explosion makes it unfeasible to calculate C N for large values of n. To overcome this problem, C N can be approximated by a sampling method, which consists of taking M random samples of subsystems of size k (for k ¼ 1,. . ., n 2 1); as M increases, the approximation approaches the exact value of C N ; hence, M must be sufficiently large to get an accurate approximation and this computation becomes very intensive. Alternatively, a computationally less intensive approximation for C N has been proposed [41], used here to compute C N values. This approximation assumes that the correlation between the activities of neural components is small and that the dynamics of the system are stable; C N can then be approximated as andR is the correlation matrix with elements equal to zero in the diagonal (r ii ¼ 0 8i).
To derive correlation matrices from the structural connectivity matrices, we implemented a linear model of neural activity described in [42]. This model is based on the linearization of a coupled neural system driven around a fixed point by spatially and temporally independent Gaussian noise sources. The correlation of the dynamics, matrix A, can be obtained analytically as The coupling matrix C used here is the structural connectivity matrix, I is the identity matrix and a is the rate of activation leakage per node (here set to a ¼ 2). Finally, the condition for weakly coupled neural elements is met when the spectral radius (the absolute value of the largest eigenvalue) of the covariance matrix is smaller than unity.
The condition of weakly coupled neural elements comes from the assumption that neural dynamics can be approximately characterized by a stationary multivariate Gaussian process. By implementing a Gaussian neural model as a linear process, the computation of C N is remarkably simplified, given that it is possible to express the interactions between neural elements (i.e. entropies and mutual information) in terms of a covariance matrix that can be derived analytically from the network's connectivity matrix. However, the linearization of the neural model is an approximation in the weakly coupled near-linear regime of the nonlinear dynamics of the system. Linear approximations are commonly used in neuroscience to model large-scale neural systems [43], which is what the nodes of the networks used in this work represent.

(d) Network morphospace
The concept of morphospace originated in the context of evolutionary biology [44]. It provides a framework to map all the possible biological forms that can result by varying the parameter values of a geometrical or mathematical model of form. Most importantly, it allows the identification of forms that have been produced in nature and forms that have not. The parameters of a model of form define the axes or dimensions of a morphospace, where different locations within each dimension specify the parameter values and are associated with a particular biological form. Here, we extend the concept of morphospace to encompass network structure; hence, the dimensions of a network morphospace are given by network structural measures and positions within the morphospace correspond to characteristic aspects of network topology (figure 1a). Specifically, in this paper we define a threedimensional morphospace with axes given by E diff , E rout and C N ; therefore, networks are placed within the morphospace according to the previously defined measures.
In the context of this study, a morphospace exploration is the process of simulating brain-like networks and identifying their locations within the efficiency-complexity morphospace by measuring their corresponding values of E diff , E rout and C N . Simulated brain-like networks are generated by incremental rewiring of a population of minimally perturbed empirical brain networks, which is implemented through a multi-objective evolutionary algorithm that approximates a Pareto-optimal set within the morphospace.
Given that the brain networks constructed from each dataset display differences in size and density, we define three distinct morphospaces, one for each dataset [46]. Within each morphospace, we constrained the set of possibly simulated networks to preserve the number of nodes and edges, degree sequence and network cost invariant with respect to the values measured from the corresponding empirical network. Finally, graph metrics are normalized in each morphospace with respect to the values of the corresponding empirical networks. Therefore, each empirical network is located within its morphospace in the coordinates (E diff , E rout , C N ) ¼ (1, 1, 1).

(e) Morphospace analysis
In a theoretical morphospace, there is a distinction between possible and impossible forms (topologies), and a second distinction between functional and non-functional forms (topologies). The former distinction refers to topologies that are impossible because the combination of parameters (structural traits) is meaningless or infeasible and no network can satisfy such combination. The later distinction defines a subset of the formally possible topologies, which consists of the functionally feasible topologies. This means that within the space of possible topologies, there are networks that are not functionally feasible and therefore, such networks will not be found in the real world. For instance, all disconnected networks belong to the set of impossible brain network topologies. In addition, crucial properties of brain networks such as density, length and volume of neuronal connections, are subject to functional constraints that define the subset of biologically feasible networks within the set of possible networks in the brain network morphospace. In our work, we perform a morphospace exploration that implements functional constraints through a rewiring algorithm, whose objective is to preserve the total cost of the network connections. This strategy of morphospace exploration is based on evolving a population of networks by repeatedly carrying out two steps: network selection and network variation.
Network selection occurs according to Pareto optimality, a concept used in economics and engineering to describe a set of solutions that optimize multiple objectives simultaneously [47]. In general, a solution is said to be Pareto-optimal if an improvement of any single objective cannot be achieved without negatively affecting some other objective (figure 1b). In the context of a population or ensemble of networks that are being evaluated by multiple objective functions, a network G belongs to the Pareto-front set if and only if (1) G is not worse than any other network within the population, with respect to all objectives; (2) G is strictly better than any other network in the population, with respect to at least one objective [45] (figure 1c).
Network variation refers to small structural changes that are implemented with a rewiring algorithm. The algorithm is based on a random rewiring algorithm whose elementary moves are the so-called 'edge-swaps' [48]. An edge swap consists of the following steps: (a) Randomly, select four distinct nodes, namely (i, j, k, l ). (b) If (a ik , a jl , a il , a jk ) ¼ (1, 1, 0, 0) then swap the edges, so that the adjacency matrix entries become (a ik , a jl , a il , a jk ) ¼ (0, 0, 1, 1). (c) Otherwise, go back to (a).
This rewiring procedure ensures that the rewired network always remains connected and that the degree of each node is rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 369: 20130530 unchanged [49]. In this study, we require that each edge swap satisfies two additional conditions. The first condition is that the total wiring cost of the network must remain constant. The second condition is that the value of the fibre densities of all connections remains confined to the interval (0, w max ), where w max is the maximum fibre density found among all the connections of the empirical network. Therefore, when an edge swap is performed, the fibre densities w il , w jk corresponding to the added edges fi, lg and f j, kg are randomly chosen, provided that they satisfy the equation w ik l ik þ w jl l jl ¼ w il l il þ w jk l jk subject to 0 , w il , w jk , w max (figure 1d ). We refer to an edge swap that satisfies these conditions as a rewiring step.

(i) Fibre-length interpolation
Note that many times, when a pair of edges are swapped, the values of l il and l jk are not defined in the original connectivity matrices extracted from the neuroimaging data, simply because there is no actual connection between the pairs of nodes fi, lg and f j, kg. To assign fibre length values to edges created during the rewiring process, for each dataset a fully connected matrix LI was constructed combining existing and interpolated fibre length values between all pairs of nodes. This was done by making the assumption that two fibres whose starting and ending points are close in space should follow similar trajectories in the brain and thus have similar lengths. Under this assumption, an estimated fibre length based on similar existing fibres can be assigned to pairs of nodes that are not connected in the empirical network. We consider two fibres to be similar if we can define two neighbourhoods-one containing the starting nodes of the fibres, and another containing the ending points of the fibres-such that the radius of each neighbourhood is smaller than d, where d is defined as 20% of the Euclidean distance that separates the centres of both neighbourhoods. For every pair of unconnected nodes fi, jg in the empirical network, we defined such neighbourhoods and looked for fibres whose endpoints are in each of the neighbourhoods. If such fibres were found, we assigned their average length to the length of a fibre between nodes i, j. Finally, for all pairs of nodes for which we could not find similar fibres connecting their neighbourhoods, we used a polynomial interpolation of degree 2 to fit the fibre lengths as a function of Euclidean distance. For the right hemisphere sub-network of each dataset used in this study, the fraction of fibre lengths estimated by polynomial interpolation was 88%, 87% and 78% for LAU1, LAU2 and UTR respectively; the correlation between the interpolated fibre length matrix and the Euclidean distance matrix is 0.980, 0.982 and 0.992 for LAU1, LAU2 and UTR, respectively. Figure 1e shows the matrix of Euclidean distances between all pairs of nodes of the LAU1 network, together with the full fibre length matrix LI after the interpolation process.
(ii) Evolutionary process All morphospace explorations start with an initial population of M 0 ¼ 500 networks, derived from the empirical brain network by performing three rewiring steps. During the evolutionary process, all networks preserve the number of nodes and edges, degree sequence and network cost of the corresponding empirical network. For every epoch of the evolution (defined by a single iteration of network selection, followed by network variation), the objective functions are evaluated on all of the population members. Selection according to Pareto optimality is applied to define the set of networks that pass unchanged to the next epoch. The set of networks that do not belong to the Pareto front are eliminated and substituted with a random sample of the Pareto-front members, which is subjected to minimum variation by carrying out one rewiring step on each network. We avoid falling in to local maxima by introducing noise in the evolutionary process as follows. At a given epoch, if more than 90% of the population belongs to the Pareto-front set, then half of the population is randomly selected and subjected to one rewiring step. Then the simulation carries on.
In order to explore a greater extent of the sub-region of the morphospace surrounding the empirical brain network, we carried out eight independent runs of the evolutionary process. All runs start with the same initial population but implement distinct objective functions that aimed to drive the population of simulated brain networks towards the eight quadrants of the three-dimensional morphospace. The eight objective functions are defined as follows: where max() and min() stand for the maximization and minimization function. The maximum number of iterations (epochs) for each objective function was set to 2000; however, all eight runs of the evolutionary process (one per objective function) required different CPU times to compute, and CPU times varied depending on the objective functions and the datasets. Therefore, the stopping condition for each evolutionary process was either 2000 iterations completed or 7 days of computation.
It is worth mentioning that in this work, the terms evolution and selection pressure are used within the context of a computational algorithm, and their usage does not necessarily reflect a direct mapping onto processes studied in the field of evolutionary biology.

Results (a) Randomization and latticization of the brain networks
We explored the behaviour of the efficiency and complexity measures when the empirical networks are rewired towards two canonical models, namely a spatial lattice-like network and a random network. Both processes involving the randomization and latticization of the empirical networks use the edge swapping algorithm (see §2e) iteratively to rewire the networks; this guarantees that the latticized and randomized networks preserve the number of nodes and edges, degree sequence and wiring cost. The latticization process takes into account the spatial positions (Euclidean distances) of the network nodes in order to create a lattice-like network where nodes tend to be connected to their spatially nearest neighbours. Figure 2a shows E diff , E rout and C N as a function of the number of rewiring steps carried out on all three empirical networks during their randomization (blue dots) and latticization (red dots), respectively. All values are averages over 40 repetitions of the randomization and the latticization processes.
In all three datasets, the three measures reach a stable regime after 2 15 rewiring steps, suggesting that additional rewiring steps do not further change the topology of the networks. Qualitatively, the efficiency measures behave similarly across the three datasets during the randomization process: E diff increases while E rout decreases as a function of the number of rewiring steps towards randomization. C N decreases in the LAU1 and UTR datasets, in agreement with earlier studies showing that randomly rewiring cortical networks tends to decrease their complexity [14]. However, we found inconsistent behaviour of C N for the LAU2 network, as well as larger variance of C N across repetitions of the randomization and latticization. In order to test whether the observed variability of C N is intrinsic to the LAU2 network or if it is produced by the approximation method used to calculate C N (see §2c(iii)), we computed C N values of randomized and latticized networks with the sampling method (see §2c(iii)); the correlation between approximated and sampled C N values is 0.989 ( p , 0.01) for the randomized networks and 0.995 ( p , 0.01) for the latticized networks, thus suggesting that variability in C N is due to differences in network topology across the datasets.
An interesting finding is that randomization decreases E rout . By contrast, high E rout is a typical characteristic of random networks [9] when these networks are not constrained to preserve network cost. The study of the fibre length and fibre density distributions of the empirical, randomized and latticized networks (figure 2b) reveals topological changes that have an effect on communication efficiency. While randomization tends to increase the amount of long-range connections in the networks, it also has the effect of decreasing the fibre density of the majority of the connections; overall, such thinning of connections diminishes E rout . Regarding the observed increase of E rout during the latticization, we note that latticized networks tend to have higher fibre density on short-distance connections, which promotes shorter path length and thus favours E rout . The second aim of our work is to characterize the effects of small perturbations on the structure of brain networks. To do so, we generated three populations of 10 000 minimally rewired network variants, each created by carrying out three rewiring steps on each of the three empirical networks, LAU1, LAU2 and UTR. We used three rewiring steps because it is the minimum number of rewiring steps that allows us to distinguish numerical differences in all three dimensions of the morphospace. In this way, we explore the proximal morphospace, that is, the space that contains the closest neighbouring network elements of each empirical network. Interestingly, in all three datasets we found that 99% of the neighbouring networks are in a region of the morphospace defined by E diff .  Figure 3 shows the distributions of proximal networks embedded in the three morphospaces corresponding to each dataset. The shape of the region occupied by these networks shows that the accessibility of the sub-region of morphospace surrounding the coordinates (E diff , E rout , C N ) ¼ (1, 1, 1) is not uniform and that there are 'preferred' directions along each axis in which networks are located.
(c) Exploring the efficiency-complexity morphospace of brain networks We implemented an evolutionary algorithm to explore a subregion of the efficiency-complexity morphospace in search for alternative biologically feasible brain-like networks (see §2e(ii)). In order to characterize the structural traits of the networks simulated through distinct selection pressures, we used eight objective functions to drive a population of networks towards the eight quadrants of the three-dimensional morphospace. For each dataset, we carried out 10 repetitions of the exploratory process. Each repetition of the process includes: (i) creating an initial population of 500 networks by minimally rewiring the empirical network and (ii) applying the evolutionary process eight times independently, once for each objective function (see §2e(ii)). Therefore, the completion of an exploratory process yields eight distinct populations of 500 networks that have been subject to different selection pressures. We will refer to each one of these populations as a front: front 1 is the population evolved by optimizing the objective function f 1 ; front 2 is the population evolved by optimizing f 2 , and so on. The stopping condition for the evolutionary process was either 2000 iterations completed or 7 days of computation. During the 10 repetitions of the morphospace exploration, for the datasets LAU1 and LAU2, the evolutionary process of all eight fronts was stopped after 2000 iterations. For the UTR dataset, the evolutionary process of fronts 3, 5, 6, 7 and 8 completed 2000 iterations, while the evolution of fronts 1, 2 and 4 was stopped after 7 days of computation, during which the processes had simulated on average 1460, 1315 and 1670 iterations, respectively. Note that all networks belonging to any front are embedded in a sub-region of the morphospace that is still fairly close to the empirical brain network; structural changes in the simulated networks account for, on average, 20% and 8% of the total number of edges present in the networks in fronts 1 through 4 and fronts 5 through 8, respectively. Figure 4 shows the regions of the morphospace explored during the evolution of the eight fronts, starting with a population of networks derived from the LAU1 network (see the electronic supplementary material, figures S1 and S2, for LAU2 and UTR datasets, respectively). Although the shape and extent of the regions explored by each evolving front vary across datasets, we find three important aspects that are consistent, regardless of the dataset used to derive the initial population. First, the evolutionary algorithm is unable to find solutions within the region defined by fE diff , 1g. Second, none of the fronts follows the trajectory of a randomization or latticization process. This demonstrates that the evolution of the network populations towards different regions of the morphospace is driven by distinct selection pressures, and not by the random nature of the rewiring algorithm. Third, the evolutionary process is able to generate brain-like networks within the region fE diff . 1, E rout . 1, C N . 1g; that is, all three topological aspects of brain networks can simultaneously increase, while preserving wiring cost.

(d) Characterization of Pareto-optimal brain networks
To allow comparisons across datasets, brain networks were down-sampled into a commonly used low-resolution partition of the human cortex, composed of 66 anatomical areas [50], with 33 areas representing the right cortical hemisphere of the brain (see figure 5a). For each dataset, networks evolved towards eight fronts (10 repetitions per front) and the final populations of 500 evolved networks were downsampled to the low-resolution partition and then aggregated according to front membership. Thus, we obtained eight populations (one for each front) of low-resolution brain networks, each population containing 5000 evolved brain networks. As the four fronts driving networks towards fE diff , 1g failed to advance, their evolved networks were not investigated further. For the remaining four fronts, we identified anatomical pathways whose fibre density and/or cost has significantly changed during the evolutionary process to favour particular topological traits. For each front, all final populations of evolved networks were aggregated into a single average network, representative of the corresponding front. The differences between the average networks of each front and the corresponding empirical network are shown in figure 5c, together with the corresponding plots recording the consistency with which connections increased or decreased in strength (figure 5d). Each of the fronts is associated with a characteristic pattern of changes in connection weights, and visual inspection suggests greater similarity in the patterns for fronts 1 and 3, and for patterns for fronts 2 and 4, respectively. Analysis of the pairwise  cosine angles between average networks of each front confirms this observation, with fronts 1 and 3 (both maximizing C N ) and fronts 2 and 4 (both minimizing C N ) exhibiting the greatest similarity across all three datasets.
Other aspects of changes in connection patterns were consistently observed across all three datasets. First, the density of evolved networks in all four fronts increases significantly (figure 6a), indicating that areas originally unconnected   have a strong tendency to become weakly connected. These new projections appear as a result of rewiring of edges away from denser pathways, thus sculpting their overall pattern into a new topology and steering the population towards one of the four Pareto fronts. Second, most of these newly formed projections extend over long spatial distances, while most of the projections that become weakened involve nodes that are spatially close, including node pairs that belong to the same anatomical region (figure 6b). Finally, a cost analysis suggests that high-cost connections, i.e. connections that contribute strongly to the overall cost of the network (which is conserved in our simulations) are principal targets for rewiring (figure 6c). Their rewiring results in a dispersal of their contribution to network cost to a larger set of connections spanning a greater number of anatomical regions.

Discussion
In this work, we applied a multi-objective evolutionary approach to rewire brain networks and place them within an efficiency-complexity morphospace. Using various Paretooptimal selection criteria, we were able to explore how brain networks evolve when subject to distinct selection pressures. Furthermore, the approach allowed us to investigate relationships and trade-offs between distinct topological traits associated with the principal axes of the morphospace, E diff , E rout and C N . Our results demonstrate that the empirical networks we used as seed points for evolution are surrounded by a large space of variant network topologies, even when holding wiring cost constant, including networks that combine a higher capacity to support efficient communication with higher neural complexity.
Our work attempts to make several methodological contributions. First, building on fundamental work in evolutionary theory [44,51,52], we extend the morphospace analysis framework into the realm of human brain networks. In the past, morphospace analysis has been applied independently in the field of complex networks [11,53,54], and in the field of neuroscience [55]; here, we combine methodology from the three fields, evolutionary biology, neuroscience and complex networks to study the topological features available for biologically feasible brain-like networks. Our work proposes that differences among variants of human brain networks can be investigated by placing these variants into a space formed by several principal axes representing fundamental measures of network organization. Within this space, gradual rewiring of network nodes and edges, for example by applying multi-objective optimization, 'moves' networks towards new topological patterns.
Second, as previously introduced in [11], here we explored two separate measures of network efficiency, one based on communication along shortest paths (routing-based communication) [9] and the other based on diffusion processes. In addition, we considered a dynamic measure of neural complexity that expressed the coexistence of segregation and integration in the network [12].
Third, as we were interested in the trade-off of these efficiency and complexity measures within the cost constraints imposed by human brain size and geometry, we employed a rewiring rule that conserved not only node degree [48,49] but also overall network cost, by adjusting the fibre density of rewired connections according to their wiring length (figure 1d).
Our first interesting finding was that full randomization resulted in networks that were less efficient in routing communication. Notably, this result diverged from the sharp increase in routing efficiency observed when non-cost-conserving randomization models are applied. When conserving cost, randomizing brain connectivity is necessarily accompanied by a thinning out of the fibre densities (shown in the distribution of the fibre density of randomized networks, figure 2), because most randomized connections span greater distances; such low-density connections do not contribute towards efficient routing communication. Full latticization of brain networks does not produce perfect lattices because of the constraints imposed by the rewiring rule; in fact, a comparison of the fibre length and density distributions between the latticized networks and the respective empirical networks reveals great similarity between these networks. Latticized networks differ slightly from the empirical networks in that they tend to have higher fibre density on short-distance connections, which promotes shorter path length and thus favours E rout .
Our next results were derived from studying the effects that minimal structural perturbations have over the measures of communication efficiency and dynamical complexity. The accessibility of the proximal morphospace (the region immediately surrounding the empirical networks) was found to be non-uniform, e.g. many more network variants exhibited greater E rout and E diff . This result held for all three datasets.
A local exploration of the efficiency-complexity morphospace confirmed that, given the network invariants of cost, density and degree sequence, the morphospace region E diff , 1 is restricted, and thus it is all but impossible to further decrease diffusion efficiency in all three datasets. Four separate fronts failed to advance in the direction of E diff , 1, suggesting that the topology of empirical brain networks tends to minimize E diff . Most important, we note that three different concepts play into the discussion of this result. First, what is the relative magnitude of a network's E diff (with respect to a null model) and to what extent can E diff be decreased (increased) through a rewiring process that preserves certain network features? Second, to what extent can the dynamics occurring on top of a network structure be explained or predicted by a diffusion-based model, regardless of whether diffusion-like dynamics are an efficient communication scheme for the system? Third, to what extent have diffusion-based dynamics been a critical evolutionary pressure shaping the structure of brain networks?
The relative value of E diff captures to what degree the structure of the network facilitates the integration of information when information spreads through diffusion-based dynamics. Therefore, in this paper, we can answer the first question by providing evidence that brain networks are close to a minimum of E diff (of course, within the space of networks with a fixed number of nodes and connections, a predetermined degree sequence and an invariant global connection cost). Furthermore, we can conclude that the topology of structural brain networks does not facilitate an efficient communication between all pairs of nodes, provided that information within the network spreads solely as a diffusion process. However, the relative value of a network's E diff does not provide any information about the underlying process through which information actually spreads within the network. Other studies have used different approaches to address this question; for example, Betzel et al. [56] and Goñ i et al. [57] have provided evidence that brain dynamics can be modelled and/or explained by diffusion-like processes. It is important to bear in mind that in these studies diffusion-like processes were used as models for the spreading of perturbations (as generative models for resting-brain functional connectivity), and it is well known from studies in other systems that modular topologies tend to limit the spread of perturbations across module boundaries. Finally, in this paper we can only speculate about the third question, regarding the role of diffusion-based dynamics as an evolutionary pressure. Nonetheless, in our opinion, it seems very unlikely that E diff is not a critical evolutionary factor shaping brain networks given that the value of E diff that brain networks exhibit is not arbitrary, but is actually at a minimum. One possible interpretation is that minimizing E diff has been due to critical evolutionary pressure shaping brain network topology in order to limit passive diffusion (e.g. of noisy perturbations) on global scales while promoting efficient diffusive communication on local scales, such as within network communities [56].
Conversely to the severe constraints found on the E diff axis, we did not find such strict constraints on the other two axes of the efficiency-complexity morphospace. Greater E rout as well as greater C N , singly or in combination, could be achieved through rewiring of specific pathways in all three datasets. Fronts advancing towards greater or lesser C N exhibited greater consistency in rewiring of specific pathways, suggesting that neural complexity depends more strongly on specific network topologies. Instead, there appear to be more structural configurations available for networks belonging to the fronts advancing towards higher or lower E rout .
In addition to these trends that were specific to the multiobjective function employed, we also observed some aspects of the rewiring process that were shared among all fronts. These aspects included a strong tendency to create new (albeit weak) pathways linking previously unconnected anatomical regions by redistributing connections away from node pairs that were spatially close and/or linked by highcost connections. This study does not allow us to determine whether these general tendencies mainly reflect constraints imposed by the cost-conserving rewiring rule or if they point to greater accessibility of parts of the morphospace by structural network variants that are more diffusely or densely connected. We acknowledge that there are several methods to rewire weighted networks and that the constraints imposed by any rewiring method will have certain effects on the topologies of the rewired networks. For instance, a rewiring algorithm that is constrained to preserve the distribution of connection weights will necessarily have the effect of increasing the total connection cost. This is because the rewiring process tends to create long-range connections and the cost measure used in this study is proportional to the connection lengths. Alternatively, to preserve crucial aspects of brain networks throughout the rewiring process, one could enforce preserving both the distribution of connection weights and the distribution of connection lengths. However, such a set of restrictions would have the effect of drastically reducing the space of solutions, imposing severe limitations on the exploration of the morphospace. Hence, for this study, we have opted for an approach that is not as restrictive as the latter but is still conservative, that is to preserve the total cost of the connections of the networks. This approach provides sufficient degrees of freedom to explore the morphospace, while imposing a strong functional constraint that allows us to study biologically feasible brain networks. Furthermore, the rewiring algorithm we present in this study provides an alternative null model to perform statistical tests of graph measures of brain networks. The use of random networks as null models is very common [46]; however, as we have pointed out previously, there are several ways to randomize a network, and hypothesis testing will yield different results depending on the selected null model. Here, we suggest that the appropriate null model is one that preserves most of a brain network's basic features that make it biologically feasible, such as grey/white matter volume, connection density and degree sequence, among others.
Several aspects of the present work require future extensions. First, networks derived from the three datasets employed here shared numerous topological features but were also somewhat variable due to differences in subject cohort, data acquisition and tractography (see §2b). These differences did not allow firm inferences about which changes in specific anatomical pathways were associated with specific Pareto fronts (i.e. due to specific selection pressures). This inference awaits the arrival of more uniformly acquired, normative datasets, for example those collected as part of the Human Connectome Project [58]. Second, while the current work examined the trade-off between different measures of network communication efficiency and complexity, the trade-off of these measures with network cost (held constant in this study) remained unexplored. This could be explored, for instance, by relaxing the invariant cost assumption used here and introducing network cost as a term in the objective function instead (either conforming or being part of one of the morphospace axes). Third, the measures forming the principal axes of the morphospace were chosen on the basis of previous work [11, 12,14] which suggested that they are relevant for various aspects of brain function and dynamics; however, alternative formulations of network morphospace that target other features of network structure and topology may be explored in future work. Furthermore, as opposed to using explicit network measures, one could favour orthogonality or statistical independence (e.g. independent component analysis) when defining the morphospace axes and the topological invariants in the evolutionary algorithm. Finally, further extensions of this work may also include an analysis of local or within-community communication efficiency, because global measures do not capture how communication efficiency is distributed among network communities.
These and other extensions of the current work may become useful for characterizing regions of network morphospace that are occupied by existing topological variants of the human brain. As is the case for biological forms [44], we expect that the majority of the morphospace is empty, i.e. most possible network configurations are either physically or economically infeasible or have been selected against in evolution. Among those variants that do occur, we expect that embedding of individual human brains in network morphospace will highlight important patterns in individual differences of network organization, including those associated with disease-related network disturbances.