We propose a working hypothesis supported by numerical simulations that brain networks evolve based on the principle of the maximization of their internal information flow capacity. We find that synchronous behavior and capacity of information flow of the evolved networks reproduce well the same behaviors observed in the brain dynamical networks of Caenorhabditis elegans and humans, networks of Hindmarsh-Rose neurons with graphs given by these brain networks. We make a strong case to verify our hypothesis by showing that the neural networks with the closest graph distance to the brain networks of Caenorhabditis elegans and humans are the Hindmarsh-Rose neural networks evolved with coupling strengths that maximize information flow capacity. Surprisingly, we find that global neural synchronization levels decrease during brain evolution, reflecting on an underlying global no Hebbian-like evolution process, which is driven by no Hebbian-like learning behaviors for some of the clusters during evolution, and Hebbian-like learning rules for clusters where neurons increase their synchronization.
A plethora of phenomena in nature can be effectively described by networks. Neuroscientists have used tools for the analysis of complex networks that help realize even more deeply the functionality and structure of the brain. It was found that many aspects of brain network structures are typical of a wide range of non-neural or non-biological complex networks [1,2]. One of the main findings in neuroscience is the modular organization of the brain, which in turn implies an inherent parallel nature of brain computations [1]. Modular processors have to be sufficiently isolated and dynamically differentiated to achieve independent computations, but also globally connected to be integrated in coherent functions [1]. It has been revealed that the cortical network is a hierarchical and clustered network with a complex connectivity [3]. A possible network description for this modular organization is that brain networks may be small-world structured [4] with properties similar to many other complex networks [5]. This viewpoint has been driven by the systematic finding of small-world topology in a wide range of human brain networks derived from structural [4], functional [6], and diffusion tensor MRI [7] studies. Small-world topology has also been identified at the cellular-network scale in functional cortical neural circuits in mammals [8] and also in the nervous system of the nematode Caenorhabditis elegans (C. elegans) [9]. Moreover, this topology seems to be relevant for the brain function because it is affected by diseases [10], normal ageing, and by pharmacological blockade of dopamine neurotransmission [11]. Synchronization is ubiquitous in nature. Insightful findings regarding synchronization in complex networks were reviewed recently in Ref. [12]. Recently, synchronization in complex modular or clustered networks has been investigated [13,14]. It appears as the interplay between the intrinsic dynamics associated to the nodes of the network and its graph topology and connecting functions. In this work, synchronization will be considered as a mean to quantify functional behaviors of the brain dynamical networks (BDNs) studied. By BDN we mean a network that represents the connectome equipped with neural dynamics for its nodes to account for their time evolution. Mathematical and computational approaches have a long tradition traced back to the early mathematical theories of perception [15] and of current integration by a neural cell membrane [16]. Hebb’s idea on assembly formation [17] inspired simulations on the largest computers available at that time (1956) [18] to understand the relation between neural connectivity strength and response, i. e. behavior. It is a learning rule which proposes an explanation for the adaptation of neurons during the learning process. The simultaneous activation of pairs of neurons leads to pronounced increases in their synaptic strength. There is also the possibility of many other kinds of learning [19]. In this work we find evidence of Hebbian-like and no Hebbian-like learning rules characterized by the relationship between the addition of synapses and the increase or decay in the synchronization behavior of neurons in the evolved BDNs. Several brain models have been studied so far that do not take into account particular behaviors of isolated neurons, but higher level functions of the external stimulus and the behavioral response attributed to ensembles of neurons of cortical areas. There are two classes: Those based on collective functional dynamics of local groups of neurons, such as the Wilson-Cowan model for the cortical and thalamic nervous tissue [20], and those based on the conditional probabilities (as well as information and mutual information (MI) ) of stimuli and responses, prominent examples of which are the Bayesian brain hypothesis [21] and the infomax theory [22]. Recently, it was shown that most of the probabilistic brain models can be unified under a single free-energy principle [23], the one that interprets the brain as a system that tries to minimize surprises of the sensations from the world. If the brain learns by maximizing the MI between stimuli and response [22], or by updating the internal model of probabilities by using Bayesian techniques [21], or by minimizing the free-energy [23], the surprise of the response provides little insight about the dynamical mechanisms appearing in a brain network when it evolves based on such principles. The drawback however of the probabilistic brain models is that they describe little about the underlying dynamical structure of what is really happening in the neural level and the representation of the stimuli [24]. For deterministic systems with correlations, an appropriate quantity for measuring the transfer of information is the Mutual Information Rate (MIR), the MI per time unit. In Ref. [25], the authors have developed alternative methods to overcome problems that stem from the definition of probabilities from time series and derived an upper bound for the MIR, Ic, between two nodes of a complex dynamical network from time averages that do not rely on probabilities, but instead on the two largest Lyapunov exponents l1, l2 of the subspace of the network formed by the two nodes (see Eq (10) in Materials and Methods). In our study, Ic stands for the upper bound for the information transferred per time unit between any two nodes of the BDN, what represents the information flow capacity of the BDN. We discuss more on these in Materials and Methods, Section Upper Bound for MIR. Inspired by the infomax theory and by theoretical studies that proposes that the maximization of information transmission between subsystems can be used as a principle for understanding the development and evolution of complex brain networks (see Ref. [2] and references therein, and Ref. [26]) and, aiming at elucidating the microscopic dynamic mechanisms associated to the evolution of brain networks, we propose a working hypothesis and, provide evidence that a brain dynamical network may evolve based on the maximization of the information flow capacity it can internally handle at each step of its evolution process, i. e. a new inter-neuron connection is established if it leads to a subsequent increase of the information flow capacity of the new brain circuitry. Our hypothesis is based on the internal neural network dynamics and on a plausible model for brain structure per se without the need to resort to probabilistic models based on the external influence in the brain and its response. We have been able to show that our evolved BDNs present similar synchronization and information flow capacity behaviors with those found for the simulated dynamical networks for the brain structure of the C. elegans and humans. Moreover, we show that BDNs evolved with coupling strengths that maximize the information flow capacity are the ones with the smallest spectral graph distance from the BDNs of the C. elegans and humans, and that, during the growing process, their information flow capacity increase is related to moderately low amounts of global neural synchronization. This work provides ample evidence that brain networks may grow by maximizing the capacity of information flow they can internally handle, driven by global no Hebbian-like evolution processes, according to which the addition of interconnections between clusters during the evolution process leads to a decrease in the global synchronization level of the BDN. This behavior is accompanied by a similar no Hebbian-like learning process for neurons in some of the clusters and by Hebbian-like processes for neurons in the remaining clusters, leading to an increase in the synchronization level between these neurons during brain network evolution. Effectively, the no Hebbian-like mechanism is akin to the unlearning anti-Hebbian mechanism of Crick and Mitchison [27] that refers to the elimination of unnecessary connections to prevent overload and to render the network more efficient. In our evolution we do not delete links. However, both mechanisms lead to a decrease of synchronization and to more efficient networks, being in our case the evolved networks able to maximize their information flow capacity. Global synchronization takes into consideration the synchronous behavior of all neurons in the BDN whereas local synchronization of neurons in a cluster of the BDNs. Finally, our work shows that optimising information flow capacity leads to evolved networks that are heterogeneous, in accordance with the line of research in Ref. [2], where the authors report on a mathematical model for the evolution of heterogeneous modules in the brain based on the maximization of bidirectional information flow transmission.
We present in panels (A), (B) of Fig 1 the global synchronization measure ρ and upper bound of information flow Ic for the BDN of the C. elegans in the parameter space of chemical coupling gn in [0,2] and electrical coupling gl in [0,2]. We refer the reader to Materials and Methods, Subsection C. elegans Data. A direct comparison between the two panels reveals a number of conclusions for the different parameter space regions. At first, for relatively high chemical and electrical couplings almost full global synchronization can be achieved (yellow and red regions in Fig 1 (A) ). For the same region, panel (B) shows an almost absence of capability of information transmission as the upper bound for MIR, Ic ≈ 0 (dark blue region). High levels of global synchronization accompanied by low values of Ic indicate that not only neural activities are similar but also they have very low entropy since both Lyapunov exponents λ1, λ2 are practically zero leading to their difference Ic ≈ 0 (see Materials and Methods, Subsection Upper Bound for MIR). Second, for chemical couplings smaller than 0. 3 (i. e. gn ∈ [0,0. 3]) and electrical couplings gl ∈ [0,2], we observe a multitude of different functional behaviors: There are regions of high synchronization (red region in panel (A) ) and low Ic (blue region in panel (B) ) and others with exactly the opposite behavior (i. e. low global synchronization accompanied by high information flow capacity). These different functional behaviors will become even more evident in Fig 2, where we plot the regions between the left vertical axes and the white dotted lines of Fig 1 in a finer resolution. Our findings suggest that for the C. elegans BDN, for most of the chemical and electrical couplings, global brain synchronization is roughly speaking inversely related to the information flow capacity, Ic. As we shall see in Materials and Methods, Section Similarities between C. elegans and human BDNs for a more detailed analysis of the C. elegans and human BDNs, high synchronization as depicted by ρ implies small information flow capacity. The relation between ρ and Ic can be better understood in the basis that typically ρ ∝ 1/λ2, where λ2 is the second largest Lyapunov exponent of the BDN, for the so-called non-excitatory networks [28,29]. The non-excitatory character is expected to be prominent when the electrical coupling has a dominant contribution to the behavior of the network with respect to the chemical, a situation that promotes global neural synchronization. We have performed a similar study for the global synchronization ρ and upper bound of information flow Ic for the human brain networks in Materials and Methods, Subsection Human Subjects Data, shown in Fig 1 (C) and 1 (D). Particularly, we first prepared parameter spaces for all six human subject BDNs with the same coupling ranges with those of the first two panels of the C. elegans and then computed their average, presented in Fig 1 (C) and 1 (D). The averaged ρ and Ic quantities from the six subjects are indicated by ⟨ρ⟩6 and ⟨Ic⟩6, respectively. A direct comparison between panels (C) and (D) of Fig 1 for the humans reveals a number of parameter space regions associated to different functional behaviors, similar to those observed for the C. elegans. For relatively high chemical and electrical couplings almost full global synchronization is achieved (yellow and red regions in Fig 1 (C) ) since ⟨ρ⟩6 ≈ 1. For the same region, Fig 1 (D) shows an almost absence of information flow capacity as ⟨Ic⟩6 ≈ 0, the reason being the same as for the C. elegans BDN. For chemical couplings smaller than 0. 2 (i. e. gn ∈ [0,0. 2]) and electrical couplings gl ∈ [0,2] we observe regions of high synchronization (red regions in Fig 1 (C) ) and low ⟨Ic⟩6 (blue region in Fig 1 (D) ), as well as others with exactly the opposite property, having low synchronization and large ⟨Ic⟩6. In order to look deeper into the details of the functional behaviors of these BDNs and, to understand the structural and functional similarities between them, we study in Materials and Methods, Section Similarities between C. elegans and human BDNs, zoom-in plots of the previous parameter spaces. There has been enormous research devoted on the C. elegans worm which has revealed its ability to learn about mechano, chemo and thermosensory inputs and stimuli [30,31]. It was also shown that its neural system has the ability to distinguish between tastes, odours or any indication related to the presence or absence of food. It also shows different kinds of learning behavior, such as associative (classical conditioning and differential classical conditioning), and non-associative forms of learning, such as habituation and dishabituation [32]. These properties are reminiscent of the human brain ability to adapt to different stimuli and environments. In Fig 2, we present a finer resolution version of the parameter spaces of Fig 1 for the C. elegans and humans. They allow us to reveal the extraordinary similarity on the functional level, the global synchronization and Ic patterns between the C. elegans and human BDNs. In these plots, wherever we observe high synchronization as evidenced by orange and red regions in panel (A) for the C. elegans and (C) for the human brain network, the upper bound for MIR, Ic, that stands for the internal information flow capacity of the BDN, is small (blue region) and vice-versa. Since the C. elegans brain connectivity network is about four times smaller in size than the human BDNs [33] studied here, we should expect that the global synchronization and upper bound for MIR patterns could occur for different ranges of chemical and electrical couplings as compared to those of Fig 1 (see for example Ref. [34]). Therefore, a rescaling of the coupling strengths was employed to allow for both BDNs to have the possibility to produce equivalent dynamical behaviors. This rescaling is described in Materials and Methods, Subsection Rescaling of Chemical and Electrical Couplings for Parameter Spaces of Networks with Different Eigenvalue Spectra. It is worth noting that there is an optimal coupling range for both BDNs that allows for large information flow capacity in the brain networks (orange and yellow regions in panels (B), (D) of Fig 2), a coupling range that promotes moderately low global synchronization! In Materials and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information Flow Capacity, we propose an artificial brain network evolution model that presents important structural and functional properties of the BDNs of the C. elegans and humans. It is based on the combined effect of chemical and electrical synapses and, on a topology reminiscent of interconnected brain communities found in these BDNs. We use chemical synapses for the communication of neurons of different clusters (inter-cluster connections) and electrical for the communication of neurons within each cluster (intra-cluster connections). Here, we compare the functional properties of this brain network evolution model with those of the C. elegans and humans. By functional we mean the properties of the dynamics of the BDN such as local, global synchronization and information flow capacity. Particularly, we report on the similarities we found for the functional properties of the BDNs of the C. elegans and humans, and for those of the proposed model for brain network evolution. To support further the validity of the presented results for the model for brain network evolution and to show its independence on the particular initial small-world cluster configuration, we computed its parameter spaces averaging over the functional measurements obtained from five different initial small-world cluster configurations (for a discussion about the creation of small-world networks or clusters see Materials and Methods, Subsection Analysis of Networks and Communities). In all cases, we used the same evolution process as described in Materials and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information Flow Capacity. The results of this study are shown in Fig 3, where ⟨ρ⟩5 stands for the average of the global synchronization of the evolved BDNs for the five realizations and ⟨mMIR⟩5 for the average maximal value of the MIR for the same realizations. We first provide evidence in panel (A) and (C) to (H) that the finally evolved BDN of small-world clusters captures similar local and global synchronization properties to those of the C. elegans BDN. The local synchronization ρci of the i-th community (Fig 3 (C) –3 (H) ) also reproduces similarly the global synchronization patterns of Fig 2 (A) and 2 (C) for the C. elegans and human BDNs. Comparing these panels, we conclude that the averaged synchronization measure ⟨ρ⟩5 of the brain network evolution model (Fig 3 (A) ) attains almost similar values in the same coupling regions to those of the C. elegans in Fig 2 (A). The different ranges on the horizontal axes of the chemical coupling strength gn can be explained by the fact that the finally evolved network consists of 60 neurons whereas the C. elegans of 277 neurons (for more details see Materials and Methods, Subsection Rescaling of Chemical and Electrical Couplings for Parameter Spaces of Networks with Different Eigenvalue Spectra). It is remarkable that the finally evolved BDN exhibits almost identical synchronization and information flow capacity features as the BDN of the C. elegans and humans for almost all coupling ranges considered. Particularly, focusing on panels (A), (B) (for the C. elegans) and (C), (D) (for the humans) of Fig 2 and, on (A), (B) of Fig 3 (for the brain network evolution model), we observe an almost identical pattern of global synchronization and information flow capacity as depicted by Ic and its averages. Again, here we have used the rescaling in Materials and Methods, Subsection Rescaling of Chemical and Electrical Couplings for Parameter Spaces of Networks with Different Eigenvalue Spectra, to create networks that can potentially reproduce similar dynamical behaviors for the finally evolved BDNs, in agreement with those identified for the C. elegans and human BDNs earlier. The upper bound for MIR, Ic, depends on various factors, such as network topology and connectivity patterns, coupling strengths and types (chemical, electrical), synchronicity, etc. The values of Ic, the information flow capacity, can be comparable (or different) for networks with different topologies. This is due to the fact that Ic is a function of the chemical and electrical coupling strengths. The surprising fact is that as we evolve a network with an initial small-world network configuration, by maximizing information flow capacity, the final network exhibits not only similar topological (structural network characteristics) but also similar functional or behavioral (synchronization and upper bound for MIR) properties as those found for the brain dynamical networks of the C. elegans and human subjects (see Materials and Methods, Subsection Spectral Similarity of C. elegans and Human Brain Networks with those of the Model for Brain Network Evolution). For completeness, in S2 Fig, we present a similar analysis to the one of Fig 3 based on Erdős-Rényi random networks (panels (A), (B) ), scale-free (Barabási-Albert) (panels (C), (D) ) and star topologies perturbed by 20% for the clusters of the model for brain evolution of 60 neurons and 6 small-world clusters (panels (E), (F) ) and found out that the evolution model fails to capture similar functional (local and global synchronization and, information flow capacity patterns in the parameter spaces) as the same model equipped with small-world topologies for its clusters (see Fig 3). Here, we are interested in studying and proposing a model for brain network evolution that is able to reproduce not only similar functional properties such as information flow capacity and, local and global synchronization properties, but also importantly to reproduce similar structural properties for the finally evolved full brain network and of its clusters. Based on our results so far, we show next that the initial cluster configuration able to fulfil both requirements is the small-world cluster topology. The discussion about the results of Fig 4 in Materials and Methods suggest that the evolution process of a basic small-world clustered network is capable of generating an evolved one with similar structural properties with those for the C. elegans and human BDNs (see Materials and Methods, Subsection Structural Properties of the Model for Brain Network Evolution). The structural similarity to the human brain network is even more remarkable if the couplings of the network to be evolved are within the range that promotes high levels of information flow capacity and low neural synchronization, a prominent example of which is case 𝓑 (for a definition and discussion about cases 𝓐, 𝓑, see Materials and Methods, Subsection Brain Network Evolution Promotes Global no Hebbian-like and, Local Hebbian-like and no Hebbian-like Evolution Learning). Here we study how close the normalized Laplacian spectral plots of the networks considered in Fig 4 (Q) are from those of evolved BDNs. Particularly, we examined the spectral similarity by comparing spectral plots and computing their average Euclidean distance following Ref. [35]. We refer the reader to Materials and Methods, Subsection Normalized Laplacian Spectra, for the details of our computations. We compared the normalized Laplacian spectral plot of the C. elegans and of the six human brain networks with the spectral plots of all finally evolved networks of the five realizations used to compute the averaged parameter spaces of the model for brain network evolution of Fig 3 (A) and 3 (B). We provide additional support by demonstrating in the last four panels of Fig 5 results from a similar comparison between the C. elegans and humans with a double-sized version of 120 neurons of the model for brain network evolution of the Materials and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information Flow Capacity. The normalized Laplacian spectra for the C. elegans and for the averaged over the six human BDNs are shown in Fig 4 (Q). As it can be seen, both spectra share some common features: Both show a left-skewed distribution in which the largest eigenvalue is closer to one in agreement with results in Ref. [35]. Also, the distributions show peaks around one and the eigenvalues are scattered at the beginning of the spectra, suggesting similarities in their community structure, reminiscent of their small-worldness. Although the spectral plots of both cases of the model for brain network evolution shown in Fig 4 (R) do not exhibit all properties of the spectral plots of the C. elegans and humans of Fig 4 (Q), maybe due to the considerably smaller size of the former networks, they do exhibit interesting similarities: They are both left-skewed distributions with a peak around 1. 3, being closer to 1 than to 2. We also observe low relative frequency eigenvalues at the beginning of both spectra. Both spectral properties suggest similarities in their community structure [35] as well (i. e. their small-worldness). This is in accordance with the spectral plots of the averaged human and C. elegans brain networks of Fig 4 (Q). We argue that this similarity comes from the small-worldness of the communities. The close relation between the normalized Laplacian spectra of Fig 4 suggests the existence of common underlying structural properties of the neural networks of the C. elegans, the humans and the model for brain network evolution. We measured the similarity between the spectral plots of the C. elegans and the averaged human brain network with the averaged model for brain network evolution and, plot in panels (A), (B) of Fig 5 their spectral distance D for different chemical and electrical coupling ranges. In this framework, the closer D is to zero, the closer structurally the compared networks are. In particular, panel (A) is the parameter space for the spectral distance between the C. elegans and the averaged model for brain network evolution and, panel (B) is a similar plot for the spectral distance between the averaged human brain network and the same model for brain network evolution. We pinpoint case 𝓐 by a ▲ and case 𝓑 by ●. Such results allow one to draw interesting relations between structure and function of the proposed model for brain network evolution with the structural properties of the brain networks of C. elegans and humans. Panels (A) and (B) of Fig 5 already reveals that the smallest mean spectral distance happens for the pair of chemical and electrical couplings that gives rise to BDNs that present small amount of synchronization and high information flow capacity, in other words to cases such as 𝓑. In contrast, one of the largest mean spectral distances was found for case 𝓐 that promotes high amount of neural synchronization and small information flow capacity in the brain network!
The complexity of the circuitry of the nervous system of the human brain is still a big challenge to be resolved as it contains about 86 billion neurons and thousands times more synapses. A synapse is a junction between two neurons and it is a mean through which neurons communicate with each other. There are electrical and chemical synapses: An electrical synapse is a physical connection between two neurons which allows electrons to pass through neurons by a very small gap between nerve cells. Electrical synapses are bidirectional and of a local character, happening between neurons whose cells are close. They are believed to contribute to the regulation of synchronization in the brain network. In contrast, chemical synapses are special junctions through which the axon of the pre-synaptic neuron comes close to the post-synaptic cell membrane of another neuron or non-neural cell. In Ref. [37], the authors report on the self-consistent gap junctions and chemical synapses in the connectome of the C. elegans. In our work, we use both kinds of synapses. Following Ref. [34], we endow the nodes (i. e. neurons for the c. elegans and neural ensembles for the humans) of the networks with Hindmarsh-Rose brain dynamics [38]: p ˙ = q - a p 3 + b p 2 - n + I ext, q ˙ = c - d p 2 - q, n ˙ = r [ s (p - p 0) - n ], (1) where p is the membrane potential, q is associated with the fast current, Na+ or K+, and n with the slow current, for example Ca2+. The rest of the parameters are defined as a = 1, b = 3, c = 1, d = 5, s = 4, p0 = −1. 6 and Iext = 3. 25 for which the system exhibits a multi-scale chaotic behavior characterized as spike bursting. r modulates the slow dynamics of the system and was set to 0. 005 so that each neuron or neural ensemble be chaotic. For these parameters, the HR model enables the spiking-bursting behavior of the membrane potential observed in experiments made with a single neuron in vitro. It is also a relatively simple model that provides a good qualitative description of the many different patterns empirically observed in neural activity. We couple the HR system to create an undirected BDN of Nn neurons connected simultaneously by electrical (linear diffusive coupling) and chemical (nonlinear coupling) synapses: p ˙ i = q i - a p i 3 + b p i 2 - n i + I ext - g n (p i - V syn) ∑ j = 1 N n B i j S (p j) - g l ∑ j = 1 N n G i j H (p j), q ˙ i = c - d p i 2 - q i, n ˙ i = r [ s (p i - p 0) - n i ], ϕ ˙ i = q ˙ i p i - p ˙ i q i p i 2 + q i 2, i = 1, …, N n. (2) In our study, ϕ. i is the instantaneous angular frequency of the i-th neuron and ϕi is the phase defined by the fast variables (pi, qi) of the i-th neuron. We consider H (pi) = pi and: S (p j) = 1 1 + e - λ (p j - θ syn), (3) with θsyn = −0. 25, λ = 10, and Vsyn = 2 to create excitatory BDNs. Eq (3) is a sigmoidal function that acts as a continuous mechanism for the activation and deactivation of the chemical synapses and, also allows for analytical calculations of the synchronous modes and synchronization manifolds of the coupled system of Eq (2) [34]. In Eq (2), gn is the coupling strength associated to the chemical synapses and gl to the electrical synapses. For the chosen parameters, we have ∣pi∣ < 2 and that (pi − Vsyn) is always negative for excitatory networks. If two neurons are connected under an excitatory synapse then, when the presynaptic neuron spikes, it induces the postsynaptic neuron to spike. We adopt only excitatory chemical synapses in our analysis. We use as initial conditions for each neuron i: pi=−1. 30784489+ηir, qi=−7. 32183132+ηir, ni=3. 35299859+ηir and ϕi = 0, where η i r is a uniformly distributed random number in [0,0. 5] for all i = 1, …, Nn, following Ref. [34]. These initial conditions place the trajectory quickly in the attractor of the dynamics and thus, there is less need to consider longer transients. Gij accounts for the way neurons are electrically (diffusively) coupled and it is a Laplacian matrix (i. e. Gij = Kij − Aij, where A is the binary adjacency matrix of the electrical connections and K is the degree identity matrix based on A), and so ∑ j = 1 N n G i j = 0. By a binary adjacency matrix, we mean an adjacency matrix with entries either 0 (no connection) or 1 (connection). Bij is a binary adjacency matrix and describes how the neurons are chemically connected and therefore its diagonal elements are equal to 0, giving thus ∑ j = 1 N n B i j = k i, where ki is the degree of the i-th neuron, i. e. it represents the number of chemical links that neuron i receives from all other j neurons in the network. A positive (i. e. 1) off-diagonal value in both matrices A, B in row i and column j means that neuron i perturbs neuron j with an intensity given by glGij (electrical diffusive coupling) or gnBij (chemical excitatory coupling), respectively. Therefore, the binary adjacency matrices C of the BDNs considered in this work are given by: C = A + B. We numerically integrated Eq (2) in Fortran 90 using the Euler integration method (order one) with time step δt = 0. 01. We decided to do so to reduce the numerical complexity and CPU time of the required simulations to feasible levels as a preliminary comparison of trajectories computed for the same parameters (i. e. δt, initial conditions, etc.) with integration methods of order 2,3 and 4, revealed similar results. We evolve the dynamics of the brain networks and calculate their two largest Lyapunov exponents λ1, λ2 for the estimation of the upper bound for MIR, Ic. We use the well-known method of Refs. [39,40] to compute the Lyapunov exponents needed for the estimation of the upper bound Ic for MIR (see also Materials and Methods, Subsection Upper Bound for MIR). The numerical integration of the HR system of Eq (2) for the C. elegans and human BDNs was performed for the final integration time tf = 5000 and the computation of the different quantities such as the order parameter ρ of Materials and Methods, Subsection Synchronization Measures in BDNs and the Lyapunov exponents, starts after the transient time tt = 300 to make sure that orbits converged to the attractor of the dynamics. The same parameters for the model of brain network evolution of Materials and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information Flow Capacity were set to tf = 2500 and tt = 300 to reduce the numerical complexity and CPU time to feasible levels, retaining similar results. We have been careful to check that using exactly the same values as for the C. elegans and humans, the conclusions were practically the same. In this work we propose an artificial evolution model for brain network connectivity that captures important structural and functional properties of the BDNs of the C. elegans and humans. Our idea is reminiscent of modular processors that are sufficiently isolated and dynamically differentiated to achieve independent computations, but also globally connected to be integrated in coherent functions [1]. The model is based on the consideration of the combined effect of chemical and electrical synapses between neurons based on a topology reminiscent of interconnected brain network communities found in the BDNs of C. elegans and humans. In this study, we consider chemical synapses solely for the communication of neurons of different clusters of the network (inter-cluster connections) and electrical synapses for the communication of neurons within each cluster (intra-cluster connections). This idea comes from the biological local and non-local nature of these connections. Particularly, we consider a starting network topology for brain network evolution, where Nc clusters of electrically coupled neurons are connected in a closed ring as shown in S1 Fig of the Supporting Information. We endow each cluster with a small-world topology [9] as this is what we found to be more plausible to happen on the brain networks of the C. elegans and humans (see Subsection. Analysis of Networks and Communities in Materials and Methods). We also use in our evolution model, for simplicity but this is not mandatory, clusters of the same number of neurons. We denote the total amount of neurons in the network by Nn. Each small-world cluster is connected by an inter-cluster connection with its two nearest neighbour clusters by only chemical excitatory connections. In S1 Fig of the Supporting Information, one can see such an example of a small-world network topology that comprises Nc = 6 clusters and Nn = 60 neurons, where the red links denote the chemical inter-cluster connections and black the electrical intra-cluster connections. For this model network, we subsequently compute the two largest Lyapunov exponents λ1, λ2 of the BDNs following Refs. [39,40] to estimate the upper bound Ic = λ1 − λ2 for the MIR of the network, i. e. the maximum amount of information per time unit that can be exchanged between the neurons of the basic (not yet evolved) network, aka its information flow capacity (for the details see Subsection. Upper Bound for MIR in Materials and Methods). We then evolve the starting network, such as the one in S1 Fig, by adding new chemical excitatory inter-cluster connections to simulate the creation of new chemical synapses between neurons of different clusters. The electrical connections, topology and the values of the chemical and electrical coupling strengths are not modified during the evolution process. We adopt the following evolutionary rule to imitate brain plasticity [19]: If the newly added inter-cluster chemical connection leads to an increase of Ic prior to the addition, the new synapse is retained. If, instead, it is found not to increase Ic then it is deleted from the network and the random search for another one starts, being this procedure iterative. We choose the nodes of the different small-world clusters so that to simulate the addition of new inter-cluster connections in a random fashion (i. e. the candidate links are randomly chosen from a uniform distribution). The iterative procedure is repeated until the maximum number of possible pairs of neurons from different clusters is exhausted. We denote by mMIR the value of Ic of the finally evolved BDN which is always bigger or equal than the Ic of the starting BDN. For different values of the coupling strengths, mMIR can be achieved for different numbers of added interconnections. In all cases studied, we also compute the global synchronization measure ρ of the finally evolved BDN (for the details, see Materials and Methods, Subsection Synchronization Measures in BDNs) to allow for direct comparisons with mMIR and for the identification of relations between synchronization and information flow capacity in the BDNs. Following Ref. [34], a rough estimation on the range of chemical and electrical couplings based on those used for the parameter space of another network that is capable of reproducing similar synchronous behaviors [28] and similar amounts of Kolmogorov-Sinai entropy [29], can be computed as following: Suppose we have produced a parameter space such as those of Fig 2 showing behaviors of the BDN of the C. elegans as a function of gn and gl. Let us denote the maximum electrical coupling as g l C, the maximum chemical as g n C, the smallest positive eigenvalue of the Laplacian matrix of the electrical connections as ω m C and the average degree of the chemical connections as d ‾ C. Suppose now we want to compute a similar parameter space for another BDN. Let us denote by g l max, g n max, ω m max and d ‾ max the corresponding values of the new parameter space. Then, for the new maximum couplings we have: g n max = (d ¯ C d ¯ max) g n C (chemical coupling), (4) g l max = (ω m C ω m max) g l C (electrical coupling). (5) To arrive at Eqs (4) and (5) using the results of Ref. [34], we have assumed that a network with an average degree d ‾ for its chemical connections behaves similarly to a network with the same degree for its chemical connections. Eqs (4) and (5) provide a rough estimation on the maximum coupling strengths that can be used for the new parameter spaces. S1 Table presents ωm, d ‾, g n max and g l max for the different BDNs considered in our work. Based on these rough predictions, we then identified as best matching ranges, those depicted in the figures of the paper. Based on the maximum values of the parameter space ranges of Fig 2 (A) for the C. elegans (g n C = 0. 3 and g l C = 2), we get for the average humans g n max = 0. 21 and g l max = 1. 71, in good agreement with the maximum values of the ranges in Fig 2 (C) and 2 (D). For the model for brain network evolution we get g n max = 1. 14 to 2. 28 and g l max = 0. 97 to 1. 17 depending on the particular BDN, in accordance with the range of couplings used in Fig 3 and, 2 in S2 Fig. Similarly, for the large version of our model for brain network evolution (Nn = 120, Nc = 6) we estimated g n max = 2. 16 and g l max = 2. 9, consistent with the maximum coupling strengths used in the last four panels of Fig 5. Consequently, our methodology allowed us to identify regions of synchronization ρ and upper bound for MIR, Ic, for the different BDNs of this work with similar functional and structural properties. We initially identified the communities of the networks using the walktrap method [41] of the igraph software with six steps. The algorithm detects communities through a series of short random walks, with the idea that the vertices encountered on any given random walk are more likely to be within a community. The algorithm initially treats all nodes as communities of their own, then merges them into larger communities, and these into still larger, and so on. Essentially, it tries to find densely connected subgraphs (i. e. communities) in a graph via random walks. The idea is that short random walks tend to stay in the same community. Following this procedure we have been able to identify 6 communities in the C. elegans BDN, 10 in human subject A1,5 in A2,9 in B, 6 in C, 10 in D and finally, 7 in E. After this step, we computed various statistical quantities such as the global clustering coefficient, the average of local clustering coefficients, the mean shortest path, the degree pdf of the network and the small-worldness measure. The latter property is characterized by a relatively short minimum path length on average between all pairs of nodes in the network, together with a high clustering coefficient. Even though small-worldness captures important aspects of complex networks at the local and global scale of the structure, it does not provide information about the intermediate scale. Properties of the intermediate scale can be more completely described by the community structure or modularity of the network [42]. The modules of a complex network, also-called communities, are subsets of nodes that are densely connected to other nodes in the same module but sparsely connected to nodes belonging to other communities. Since nodes within the same module are densely intra-connected, the number of triangles in a modular network is larger than in a random graph of the same size and degree distribution, while the existence of a few links between nodes in different modules plays the role of topological shortcuts in the small-world topology. Systems characterized by this property tend to be small-world networks, with high clustering coefficient and short path length with respect to random networks. To infer the small-worldness of a network or community, we first compute the mean local clustering coefficient C of the network or community and the mean of the local clustering coefficients of one hundred randomly created networks ⟨Cr⟩100 of the same degree pdf with the studied network and also, the mean shortest path of the studied network L and the average of the mean shortest paths of the same one hundred random networks ⟨Lr⟩100. Watts and Strogatz [9] measured that many real-world networks have an average shortest path length comparable to those of a random network (L ∼ Lr), and at the same time a clustering coefficient significantly higher than expected by random chance (C ≫ Cr). Then, they proposed a novel graph model, currently named the Watts-Strogatz model, with a small average shortest path length L, and a large clustering coefficient C. We adopt this as a working definition of a small-world network or community. Therefore, small-world networks are in between the limit cases of regular graphs with large L and C and random networks with small L and C. To quantify small-worldness we use the ratios [43]: μ = L ⟨ L r ⟩ 100, γ = C ⟨ C r ⟩ 100, (6) in such a way that, for a small-world network or community, we compute: σ = γ μ > 1, (7) being the small-worldness measure. The higher is σ from unity for a given network or community, the better it displays the small-world property. For completeness, we note that for the human subject D, the walk-trap community analysis with step equal to six detected eleven communities for which the last one comprised only one neuron. Hence, in all computations, we disregarded this community as a trivial case. We present the results of the above analysis in S2 Table. We have performed all structural analyses of this paper using the igraph software. The evolution of the basic network of S1 Fig, under the principle of the maximization of the information flow capacity, is able to capture behaviors of real brain connectivity networks such as those for the communities of the C. elegans BDN, being their small-world structure a prominent reason. We found out that in all networks studied, the small-world measure σ gets values much higher than unity (see S2 Table), clearly indicating that they all display the small-world property, though in different degrees. Here, based on the extraordinary functional similarities found so far, we focus on two characteristic behaviors of the model for brain network evolution defined in Materials and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information Flow Capacity, to reveal important modular behaviors and underlying learning evolution processes: We associate the first one to the case of high global synchronization and low information flow capacity (which we call case 𝓐) and the other one to the opposite situation of low synchronization and high information flow capacity (case 𝓑). For illustration purposes, we select from the proposed model for the brain network evolution of 60 neurons and 6 clusters of Materials and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information Flow Capacity, two finally evolved BDNs with the following coupling strengths: For case 𝓐 the pair gn = 0. 2, gl = 1. 8 from the fifth realization (indicated by ▲ in Figs 3 and 5) and for case 𝓑 the pair gn = 0. 9, gl = 1. 5 from the fourth realization (indicated by ● in Figs 3 and 5). We have checked that the conclusions are valid independently of the realization and other similar pairs of coupling strengths. Also, our results reported here are independent on the initial small-world cluster configuration and initial conditions. In panels (A), (B) of Fig 6 we demonstrate the relation between global synchronization ρ and mMIR for these cases. Case 𝓑 of moderately low global synchronization (dashed black curve in panel (A) ) and high information flow capacity (dashed black curve in panel (B) ) is characterized by a larger number of added interconnections (i. e. 30) present in the finally evolved network with respect to case 𝓐 of only 12 (corresponding to the solid black curves in panels (A), (B) )! In both cases, ρ of Eq (8) for global synchronous behavior, has the tendency to decrease with different slopes (denoted by θ in Fig 6 (A) ) during brain network evolution as is evident by the fitting to the data in dashed blue lines. We attribute this behavior to an evolutionary brain network behavior, where global neural synchronization levels decrease during brain network evolution, reminiscent of a possible underlying global no Hebbian-like evolution process that promotes a decrease in global synchronization levels as new connections are added in the network. We regard this as an important global property of the model for brain network evolution that needs to be further clarified, as we need to account for what happens on the local, cluster level as well. We thus use the following procedure to infer about the underlying learning rules for the clusters: During each step of the evolution process of the model of 60 neurons and 6 small-world clusters (which is a particular BDN), we compute ⟨ρij⟩cl based on Eq (9) in Materials and Methods, to account for the average pair-wise synchronization of cluster cl, l = 1, …, 6. At the end of the time evolution, we have six such values for the clusters, lying in the interval [0,1]. We then record these values if the newly added chemical interconnection leads to an increase of the information flow capacity as depicted by the corresponding Ic and disregard them if not. At the end of the process, we result with a relationship between ⟨ρij⟩cl and the added interconnections that maximizes Ic, for all clusters (see panels (C), (D) of Fig 6). Doing so, we can account for the underlying cluster learning processes. In particular, the terminology “Hebbian-like” is employed to represent a learning rule that not only involves quantities for synchronous events defined between pairs of neurons (see Eq (9) ), but also a direct relation between synapse strength and synchronization increase. The “no Hebbian-like” terminology refers to a learning rule that involves not only a global measure of synchronous behavior (see Eq (9) ), but also refers to a phenomenon where the addition of synapses is accompanied by a decrease in the synchronization levels. We present these results in Fig 6 (C) and 6 (D). Following our considerations, we already know that globally, both finally evolved BDNs are following a no Hebbian-like learning rule. However, panels (C), (D) reveal a substantial difference for the synchronization behavior for pairs of neurons in the clusters in the two cases that allows us to assign different kinds of learning rules to them. For case 𝓑 in panel (D), the slopes θcl of the fitted lines (in blue) to the data ⟨ρij⟩cl, l = 1, …, 6 (in black) show that the third and fifth cluster have the tendency to increase their internal synchronization level as the BDN evolves. In other words, neurons belonging to these two clusters follow a Hebbian-like learning process that comes in contrast to the no Hebbian-like learning behavior of neurons belonging to the other clusters as they show the opposite trend during brain network evolution, exhibiting negative slopes for the synchronization. These findings are in agreement with the conclusions drawn from panel (F) in which neurons belonging to the third and fifth cluster are seen to be more synchronized with respect to the others (in yellow and red), a situation that promotes the modular organization in the brain with different internal levels of synchronization. On the other hand, case 𝓐 demonstrates a completely different situation in which the level of global synchronization decreases and cluster synchronization is very strong between neurons in all clusters during the whole brain network evolution, a behavior that results in a network whose neurons are unable to exchange but very small amounts of information. This is a case where the brain network can not transmit information by increasing its synaptic efficacy and thus, it is unable to learn new information! Panels (E), (F) of Fig 6 show the pair-wise synchronization level of the neurons denoted as ρij (see Eq (9) in Materials and Methods, Subsection Synchronization Measures in BDNs) of the finally evolved BDNs for cases 𝓐 and 𝓑, resulting in networks of 12 and 30 chemical inter-links, respectively. Panel (F) of case 𝓑 for moderately low global synchronization and high information flow capacity reveals a clusterized synchronization behavior. The third and fifth small-world clusters show much higher levels of internal synchrony (yellow and red points), i. e. synchronization of pairs of neurons in the same small-world community, with respect to the blue or dark blue points of low neural pair synchronization in other small-world clusters. The same also happens for a small number of pairs of neurons belonging to different clusters (off-diagonal points). One can notice that given the cluster that has the same internal level of synchronization, there can always be found two subnetworks, each belonging to one cluster, that have an equal amount of synchronization. This means that case 𝓑 corresponds to a BDN that has clustered multilayer synchronization, where different clusters become functionally connected with different common behaviors. A situation reminiscent of findings in neuroscience that demonstrate the modular organization of the brain in which modular processors are sufficiently isolated and dynamically differentiated to achieve independent computations, but also globally connected to be integrated in coherent functions [1]. Our results for case 𝓑 (see Fig 6 (F) ) come in contrast to those of 𝓐 of high global synchronization and low information flow capacity shown in Fig 6 (E) which demonstrates that all clusters and pairs of neurons attain almost the same state of almost complete synchronization, revealing a highly synchronized brain dynamical network that is not able however to exchange but only very small amounts of information between its different parts! The no Hebbian-like mechanism is effectively similar to the unlearning anti-Hebbian mechanism of Crick and Mitchison [27] that proposes the elimination of unnecessary connections to prevent overload and to render the network more efficient. Both mechanisms lead however to more efficient networks, being in our case the evolved networks able to maximize their information flow capacity. Typically, low synchronization implies λ1 ≈ λ2 leading to Ic ≈ 0. In critical points however, Ic tends to be large as long as most of the oscillation modes are stable resulting in a situation where λ1 > 0 and, λ2 ≈ 0 and positive. Thus, maximum Ic at low synchronization corresponds to maximum Ic near the critical point of the dynamics where λ1 > 0 and, λ2 ≈ 0 and positive. Self critical phenomena happen when the network has marginal Lyapunov exponents, in other words at the critical point that results in typically large Ic. In this context, our results are in agreement with the work in Ref. [44], supporting our findings for the existence of Hebbian and no Hebbian-like learning mechanisms in the level of the communities (local synchronization) at self criticality, which is responsible for the maximization of the information flow capacity of the evolved BDNs (low global synchronization), such as in case 𝓑 (for a similar result see Ref. [45]). We do not study our BDNs under different stimuli. Our main hypothesis is that the final topology and behavior of an evolved BDN that maximizes mutual information rate between its neurons is similar to real brain networks, such as those from the C. elegans and human subjects. Since our results for the maximization of the information flow capacity of the evolved BDNs happen when self critical phenomena emerge, they are in agreement with the results reported in Refs. [46,47]. We study here the structural properties of the evolved BDNs and compare with those of the C. elegans and human BDNs. As we have already demonstrated, they share common functional properties. We present the results of this study in Fig 4. Particularly, panel (A) shows the degree probability distribution function (pdf (k ‾ i) ), panel (B) the clustering coefficient CC (k ‾ i) as a function of the normalized degree k ‾ i, panel (C) the average degree k n n (k ‾ i) of the neighbors of nodes with degree k ‾ i and panel (D) the network with its distinct communities depicted by different color-shaded neighborhoods, for case 𝓐. In this context, k ‾ i is the normalized with respect to the maximum, degree. Panels (E) to (H) show similar plots for case 𝓑. Panels (I) to (L) are similar plots for the C. elegans brain network and, panels (M) to (P) for human subject A1. The plots in the second column show the correlation between different normalized degrees of the network whereas those of the third the tendency of the nodes of a certain normalized degree k ‾ to link with other nodes of a given degree. From the plots of the second column we observe that high degree nodes have the tendency to link with low degree nodes following an exponential dependence (with different exponents) implying disassortative mixing by degree. Our results from the second column of Fig 4 suggest that evolving BDNs based on the maximization of the upper bound for MIR (cases 𝓐 and 𝓑) gives rise to disassortative mixing by degree meaning that high degree nodes are preferentially connected to other low degree nodes and low to high degree nodes. It is worth noting that the structural properties of all human subjects are similar. Fig 4 (I) –4 (O) show that human subject A1 has similar structural properties to the C. elegans structure. Cases 𝓐 and 𝓑 seem also to present strong structural similarities, for the quantities considered in Fig 4, despite the profoundly different functional behaviors. We found out that case 𝓑 of moderately low global synchronization and high information flow capacity is characterized by a big number of added interconnections (i. e. 30) present in the final network whereas case 𝓐 by only 12! We have checked that this relationship between a large (small) number of added chemical inter-links with low (high) synchronization and large (low) information flow capacity happens for all similar cases of the parameter spaces. It is also noteworthy that in case 𝓑 (see Fig 4 (H) ), the number of clusters of the finally evolved network is reduced by one (since we start with six and end up with five) and is in contrast with what is happening in case 𝓐! Depending on the couplings which give rise to different functional behaviors, the brain network evolution model is capable of merging different communities, i. e. of restructuring the initial network configuration. Case 𝓑 is a more globally connected network. The previous observation can be quantified in terms of the modularity of the network, i. e. by the strength of the division of the network into modules (groups, clusters or communities). Networks with high modularity have dense connections between the nodes within modules and sparse connections between nodes of different modules, being more clustered. We have used the igraph software for these computations. By applying this idea here, we find that the modularity of the final BDN of case 𝓑 is 0. 596, very close to the average modularity of the six human subjects 0. 588 ± 0. 023. In contrast, the modularity of case 𝓐 is 0. 702. For the sake of completeness, we also report the modularity of the C. elegans topology which is 0. 375, the smaller of all cases we considered. The last result shows a network with sparser connections between the nodes within the modules and denser between nodes of different modules! For the C. elegans, the natural evolution process led to a smaller clusterization of its brain network. The results of Fig 4 suggest that the evolution process of a basic small-world clustered network is capable of generating evolved ones with similar structural properties to those for the C. elegans and human BDNs. The structural similarity to the human brain is even more remarkable if the couplings of the network to be evolved are within the range that promotes high levels of information flow capacity such as in case 𝓑. Synchronous activity has been observed in neural systems and reported to be associated not only with pathological brain states [48] but also with various cognitive functions [49]. It has been found that burst synchronization of neural systems may be strongly influenced by many factors, such as coupling strengths and types [50], noise [51], and the existence of clusters in neural networks [14]. In this paper we use the order parameter ρ to account for the synchronization level of the neural activity of the studied BDNs and of their communities [52]. It is originated from the theory of measures of dynamical coherence of a population of Nn oscillators of the Kuramoto type [53] and, can be computed by a complex number z (t) defined as: z (t) = ρ (t) e i Φ (t) = ∑ j = 1 N n e i ϕ j (t). (8) By taking the modulus ρ (t) of z (t), one can measure the phase coherence of the population of the Nn neurons of the BDN, and by Φ (t) to measure the average phase of the population of oscillators. In this context, ϕi is the phase variable of the i-th neuron of the HR system Eq (2) given by its fourth equation. Actually, one averages over time ρ (t) to obtain the order parameter ρ = ⟨ρ (t) ⟩t, the tendency of ρ in time. A value of ρ = 1 corresponds to complete synchronization of the oscillators, whereas ρ = 0 to complete desynchronization. We use Eq (8), adapted accordingly, wherever in the paper we need to compute the synchronization level of BDNs or clusters. In particular, in the case of BDNs, Nn is the number of neurons of the BDN and j runs through all Nn neurons of that network whereas in the case of clusters, Nn represents the number of neurons of the particular cluster and j refers to the particular neurons which are members of this cluster. We also compute and plot in Fig 6 (C) and 6 (D), for cases 𝓐 and 𝓑 respectively of the model for brain network evolution of 60 neurons and 6 small-world clusters, the pair-wise neural synchronization by looking at the synchronization patterns between all pairs of neurons i, j of the network as: ρ i j = lim Δ t → ∞ | C i j Δ t ∫ τ τ + Δ t e i [ ϕ i (t) - ϕ j (t) ] d t |, (9) where Cij is the adjacency matrix of the brain network and ϕi is the phase variable of the i-th neuron of system Eq (2). ρij are bounded in the interval [0,1], being ρij = 1 when neurons i, j are fully synchronized and 0 when they are dynamically uncorrelated. To correctly compute ρij, we take the averaging time large enough in order to obtain good measurements of the coherence degree of each pair. We similarly compute and plot in Fig 6 (C) and (D) for the same cases, the averaged quantity ⟨ρij⟩cl, l = 1, …, 6 over the 6 clusters, where i, j run through the neurons of each cluster. After Shannon’s pioneering work [54] on information, it became clear that it is a very useful and important concept as it can measure the amount of uncertainty an observer has about a random event and thus provides a measure of how unpredictable it is. Another related concept to the Shannon entropy that can characterize random complex systems is the MI [54] which is a measure of how much uncertainty one has about a state variable after observing another state variable. In Ref. [25], the authors have derived an upper bound for the MIR between two nodes or two groups of nodes of a complex dynamical network that depends on the two largest Lyapunov exponents l1, l2 of the subspace of the network formed by these nodes. In particular, they have shown that: MIR ≤ I c = l 1 - l 2, l 1 ≥ l 2, (10) where l1, l2 are the two finite time and size Lyapunov exponents calculated in the bi-dimensional observation space of the two considered nodes [25,55], which typically should approach the two largest Lyapunov exponents λ1, λ2 of the dynamical network if it is connected and the time considered to calculate l1, l2 is sufficiently small. In our study, the upper bound Ic for the MIR between any two nodes of the BDNs is effectively estimated by Ic = λ1 − λ2 (i. e. l1 = λ1 and l2 = λ2) and will stand for the upper bound for the information transferred per time unit between any two nodes of the BDN (i. e. between the neurons), what represents the information flow capacity of the BDN. The phase spaces of the dynamical systems associated to the neural networks we study here are excessively highly multi-dimensional and thus, estimating an upper bound for the MIR using λ1 and λ2 calculated by the methods of Refs. [39,40] instead of the MIR itself between all pairs of nodes, reduces enormously the computational complexity of the numerical calculations of this work. Besides, parameter changes that causes positive or negative changes in the MIR are reflected in the upper bound with the same proportion [25]. It is well-known that similarities between the structure of networks can be used for their classification [56]. The architecture of brain networks that describe the organization of maps of connections between neurons and brain elements at a systems level can be achieved by examining the eigenvalue spectrum of the normalized Laplacian of the connectome [35,57]. In our study, the connectome is given by the adjacency matrix of the brain network and thus we compute the eigenvalues of the normalized Laplacian based on this matrix. The eigenvalues νi, i = 1, …, Nn of the normalized Laplacian matrix L are in [0,2] which helps to compare networks of different sizes, and is defined as: L i j = { 1 if i = j, - 1 k i if i, j are connected, 0 otherwise, (11) where i, j represent any two nodes of the brain network, Lij the link between nodes i, j and ki the degree of node i. Therefore, the Laplacian spectrum of the network is given by the set of all eigenvalues of L, namely by its eigenspectrum.
In this paper we propose a working hypothesis, and provide evidence, that neural networks that evolve based on the principle of the maximization of their internal information flow capabilities produce networks whose functional behavior and topology are similar to those features observed in dynamical neural networks whose topology is provided by the C. elegans and humans. Our hypothesis goes along the lines of the infomax theory that proposes that the brain evolves by maximizing the mutual information between external stimuli and its response. When maximizing the internal information flow capacity, we are creating a network capable of processing information about external stimuli for which its information content is smaller than the information flow capacity of the evolved network. Notably, the brain evolves by the action of input signals. Our working hypothesis simplifies enormously the complexity of the involved calculations and, allow us to understand function and behavior in the brain, without the need to externally perturb non-autonomous neural networks. We have been able to show that our evolved brain networks present similar synchronization and information flow capacity behaviors with the ones found for the simulated dynamical networks for the structure of the C. elegans and human brain. Moreover, we have shown that BDNs evolved with coupling strengths that maximize the information flow capacity are the ones that have the smallest spectral graph distance from the BDNs of the C. elegans and humans, and that, during the growing process, their MIR increase is related to moderately low amounts of global neural synchronization. Actually, the global neural synchronization levels decrease during brain network evolution, revealing an underlying global no Hebbian-like evolution process (where synapse strength leads to global decay of synchronization) driven by a mix of local no Hebbian-like learning rules for neurons in some clusters and by Hebbian-like learning rules in neurons belonging to other clusters where synapse strength leads to cluster synchronization. In this context, the no Hebbian-like mechanism is effectively similar to the unlearning anti-Hebbian mechanism as both lead to more efficient networks, in the sense that in our case the evolved networks are able to maximize their information flow capacity. We note that if other models than Hindmarsh-Rose will be used, such as the Morris-Lecar, Izhikevich, or spiking map-based neural models, then the parameter regions that maximize information flow capacity and minimize synchronization (or vice versa) will be different, however we expect that this would not change the main results and conclusions of this work in the sense that brain network evolution based on the maximization of information flow capacity will lead to similar topologies, behaviors and relations for the evolved networks. For the human subjects, the graphs represent functionally connected brain regions. The Wilson-Cowan model could be appropriate to model the human brain, whereas it will not be suitable to model the C. elegans brain. As we have done here, using neurons to represent nodes in the human connectome do not reproduce the real dynamics of the brain but gives us a mean to compare results with the C. elegans and the evolved networks. In relation to our work, the maximization of the information flow capacity for low synchronization corresponds to the critical point of the dynamics in which self critical phenomena occur when the second or larger Lyapunov exponents of the BDNs are marginally positive. Finally, our results support further the hypothesis made in Ref. [2] that maximization of the information flow capacity can serve as a principle for the development of heterogeneous structures in brain dynamical networks, such as the neocortex of mammalian brains.
The study of the function of the brain is of primordial importance in neuroscience. Several brain models have been studied so far that take into account higher level functions of the external stimulus and the behavioral response attributed to ensembles of neurons of cortical areas. If the brain learns by maximizing the Mutual Information between stimuli and response, or by updating the internal model of probabilities by using Bayesian techniques, or by minimizing the free-energy, it provides little insight about the dynamical mechanisms appearing in brain networks when evolved based on such principles. In our work we propose a working hypothesis supported by numerical simulations that brain dynamical networks evolve based on the principle of the maximization of their internal information flow capacity, i. e. the upper bound for the information transferred per time unit between any two nodes. We make a strong case to verify our hypothesis by showing that the neural networks with the closest spectral graph distance to the brain networks of Caenorhabditis elegans and humans are the Hindmarsh-Rose neural networks evolved with coupling strengths that maximize information flow capacity. We also find that synchronous behavior and capacity of information flow of the evolved neural networks reproduce well the same behaviors observed in the brain dynamical networks of Caenorhabditis elegans and humans. Finally, we find that global neural synchronization levels decrease during brain network evolution.
Do Brain Networks Evolve by Maximizing Their Information Flow Capacity?