Identifying Controlling Nodes in Neuronal Networks in Different Scales

Recent studies have detected hubs in neuronal networks using degree, betweenness centrality, motif and synchronization and revealed the importance of hubs in their structural and functional roles. In addition, the analysis of complex networks in different scales are widely used in physics community. This can provide detailed insights into the intrinsic properties of networks. In this study, we focus on the identification of controlling regions in cortical networks of cats’ brain in microscopic, mesoscopic and macroscopic scales, based on single-objective evolutionary computation methods. The problem is investigated by considering two measures of controllability separately. The impact of the number of driver nodes on controllability is revealed and the properties of controlling nodes are shown in a statistical way. Our results show that the statistical properties of the controlling nodes display a concave or convex shape with an increase of the allowed number of controlling nodes, revealing a transition in choosing driver nodes from the areas with a large degree to the areas with a low degree. Interestingly, the community Auditory in cats’ brain, which has sparse connections with other communities, plays an important role in controlling the neuronal networks.


Introduction
Synchronization is widely observed in many fields such as coupled nonlinear systems and complex networks [1][2][3][4][5][6]. Especially, synchronization of distributed brain activity has been found to play an important role in neural information processing [7][8][9][10][11]. The experimentally observed brain activity, characterized by synchronization phenomena over a wide range of spatial and temporal scales, reflects the relevance for cognitive dysfunctions and pathophysiology [8]. Structurally, the analysis of the anatomical connectivity of the mammalian cortex has uncovered that large-scale neuronal networks display both high clustering and short pathlength [12,13]. The cortical network also shows a hierarchy of complex connectivity [12,[14][15][16][17].
Extensive information in mammalian cortex, such as the brains of macaque monkeys and cats, has been collected [18][19][20][21][22]. Recently, hub regions, which are believed to play pivotal roles in the coordination of information flow in brain networks [22][23][24], have been identified using modern tools from complex networks [20,22]. The hub regions of cortical networks are analyzed using node degree, structural motif, path length and clustering coefficient distributions [22]. The results in [20] highlight the influence of the topological connectivity in the formation of synchronization, revealing a few cortical areas forming a Rich-Club connectivity pattern.
Control of complex networks is a hot topic, which is closely related to synchronization of complex networks [25][26][27]. Some vertices in complex networks serve as reference sites, leaders or pacemakers [28] and drive all the other vertices toward desired targets or evolutions and thus synchronization is achieved. It is valuable to study the controllability of complex networks, especially for cortical networks due to the technical [29,30] and neuroscience backgrounds [8,16,20]. By fully utilizing the structure of the networks, Lu et al. [27] found out the minimum number of controllers for the pinning synchronization control of complex network with general topology and derived some efficient criteria to judge the success of the designed pinning controllers, which are illustrated by small-world and scale-free networks to be valid and efficient for large-scale networks.
Recently, controllability of complex networks has been studied using control theory or master stability function (MSF) [25,29,31]. Most recently, in [32], the authors reported on a generic procedure to steer a network's dynamics towards a given desired evolution, where techniques from MSF were used in connection with a greedy algorithm to determine a specific, suboptimal, sequence of nodes to be driven in order to control a network toward a desired dynamics. It is shown that there is a striking correlation between the suboptimal ranking and the inverse of the degree sequence [32]. However, it is still not clear how to determine the locations of optimal driving sequences, which is crucial in to achieve the most efficient controlling performance.
Understanding a complex network's structure is beneficial to understanding its function [33,34]. The past decade has witnessed an increasing of methods developed in this cross-disciplinary of physics community [35]. Structural properties in complex networks exist on both the microscopic level, arising from differences between single node properties, and the mesoscopic level resulting from features shared by groups of nodes. In [34], it is shown by benchmark problems how multiscale generative probabilistic exponential random graph models combined with efficient inference techniques can be used to achieve this separation of scales, resulting in an improved detection accuracy of latent classes. In [20,36], extensive numerical evidences are given to confirm the original claims that the microscopic and mesoscopic dynamics of synchronized patterns indeed follow different routes. In [33,37], mesoscopic analysis of networks is applied to exploratory analysis and data clustering.
In this study, we use the cortico-cortical network of cats' brain, which is a weighted and directed network with community structure [22]. We aim to identify controlling regions (driver nodes) of brain networks of a cat, which is equivalent to enhancing controllability of cortical networks. By converting the problem of identification of controlling nodes into a single-objective optimization problem, a recent wellstudied evolutionary computation method, the self-adaptive differential evolution (JaDE), is utilized to uncover the controlling nodes of the neuronal network. By utilizing JaDE, the controlling nodes are identified in microscopic, mesoscopic and macroscopic ways. In addition, the controlling nodes selected by JaDE are compared with the usual hubs [22], which are identified using node degree, betweenness centrality, closeness and node importance. In contrast to the usual hubs, most of the controlling nodes are selected from the nodes with a small degree. Our results reveal that the number of driver nodes plays a key role in the controllability of neuronal networks.

Results
Firstly, several examples are provided to verify the performance of JaDE [38]. JaDE is used to detect the controlling nodes/areas/ regions of the cortical network of cats' brain in microscopic, mesoscopic and macroscopic ways, respectively.
We will analyze three different scales of controlling nodes/ areas/regions in the cortical network: (1) the microscopic scale refers to the mean degree, the mean betweenness centrality (BC) and the mean closeness of driver nodes that are calculated under different numbers l of driver nodes; (2) the mesoscopic scale corresponds to the controlling communities; (3) the macroscopic scale is the controlling nodes sorted according to their total times of serving as driver nodes.
In the following, the reliability of evolutionary computation methods is shown in terms of the convergence speed, the mean value and the best value of ten runs. In order to show that JaDE is suitable for identification of controlling nodes of the cortical network, we compare it with some well-known efficient evolutionary computation approaches CLPSO [39], jDE [40], SaDE [41] and CoDE [42]. Also, JaDE is compared with some methods in complex networks theory.

Parameter Setting
The population sizes NP of all DEs and Particle Swarm Optimizations (PSOs) are set as 20 and the search range in each dimension is set to (0,N (see Materials and Methods). The maximum fitness evaluation f e, max is set as f e, max~N P Ã T Ã D, where T~250 is a constant and D~2 Ã l is the size of problem dimension. If a large T is given, the accuracy of the solutions might be refined and the computation consumption is increased linearly and vice visa. Evolutionary computation algorithms will be repeated 10 times independently for eliminating random discrepancy. Algorithms will be terminated when they achieve f e, max .

Comparison of JaDE with Evolutionary Computation Methods
The best value B and the mean value M of the solutions in ten runs are listed in Table 1. The number of driver nodes is increased from 6 to 48 with a stepsize 6. B is used to describe the best solution of algorithms found in 10 times and M is used to represent the mean value of solutions in 10 times. Note that both the best value and the mean value of solutions are of great significance for measuring the reliability of algorithms, hence we use [43] Q~ffi ffiffiffiffiffiffiffiffiffi ffi where both B and M are involved. Obviously, Q should be made as small as possible. Therefore, we also sort Q of five algorithms in an ascending way under different l and their orders P i . The mean order of each algorithm is calculated as follows 12,18,24,30,36,42,48g), ð2Þ and is also listed in Table 1. Based on the mean order Q m , the final rank of five algorithms is obtained in Table 1 (See ''Score''). Table 1 and Fig. 1 show that JaDE, CoDE and jDE perform better than the other two algorithms in terms of both search speed and convergence rate. From Table 1, JaDE ranks first and has good reliability of finding potential optimum with a satisfactory convergence speed. It is worth mentioning that JaDE is equipped with an elitism approach. Therefore, JaDE is able to find the global optimum when f e, max ??. In reality, it is unreasonable to run an algorithm with infinite generations. However, the performance of JaDE is confirmed by our simulation results (Table 1 and Fig. 1). Furthermore, a series of scientific experiments in [38] reveal that JaDE is a powerful and efficient algorithm for handling real-world optimization problems. In the following, JaDE is adopted to all the following simulations.

Comparison of JaDE with Network-based Methods
JaDE is compared with some other schemes (See Materials and Methods) from complex networks in terms of enhancement of controllability of the cortical network. The best solutions in 10 runs of JaDE under different l are used to produce the following results. It is worth pointing out that one can run JaDE for one time due to its reliability, as confirmed above.
Figs. 2 and 3 show that JaDE always performs better than the other methods. When l is large, the degree descending strategy, the BC descending strategy and the closeness ascending strategy are getting worse. Conversely, the degree ascending strategy, the BC ascending strategy and the closeness descending strategy are becoming better. The U and S-based strategies are intermediate among all the algorithms.
When only minimizing s and neglecting the effect of R, Fig. 3 shows that s (See Materials and Methods) can easily reach zero when applying JaDE, implying that it is easy to enhance controllability in the cortical network in terms of s. This phenomenon supports the finding in [44,45], in which the imaginary part of the eigenvalues of network connection matrix can be neglected to measuring synchronizability of complex networks. When minimizing R and increasing l, the controllability of the cortical network is becoming better using all the methods. However, when minimizing s and increasing l, the controllability of the cortical network is getting better when only using JaDE, which is strongly different from the case of only minimizing R.

Controllability of the Cortical Network -a Microscopic Way
When only minimizing R, Figs. 4, 5, 6, 7, 8 and 9 depicts the mean values of degree, BC and closeness of driver nodes by various methods. Figs. 4, 6 and 8 show that, the driver nodes selected by JaDE are the nodes with a large degree, a small closeness and a large BC at the very beginning. Then, the driver nodes selected by JaDE abruptly change to the nodes with a small degree, a small BC and a large closeness, when increasing l.
Specially, when p~l N Ã 100% is near 20%, the mean value of degree of the controlling nodes selected by JaDE achieves its minimum value. After the mean value of degree of driver nodes reaches its minimum value, it increases gradually and finally attains the mean value of degree of the cortical network. As a whole, Fig. 4 shows that the mean values of degree of driver nodes display a concave shape as a function of l. The standard deviation also becomes gradually larger when increasing l. The observed phenomenon indicates that, when l is not large, driver nodes are usually selected from the nodes with a small degree and nearly no nodes with a large degree are chosen. Some similar phenomena are observed when the BC and closeness of the driver nodes are shown (Figs. 6 and 8). This finding is consistent with the work in [29], in which the nodes with a large degree should be avoided choosing as driver nodes. It is worth mentioning that there exists a major difference with the finding in [29], i. e., when l is very small, the nodes with a large degree should be considered as driver nodes, as illustrated in Fig. 4. Different from optimizing R, when minimizing s, Figs. 5, 7 and 9 show that the mean values (degree, BC and closeness) of driver nodes selected by JaDE fluctuate around the mean values (degree, BC and closeness) of the network. The standard deviations (degree, BC and closeness) of driver nodes selected by JaDE keep stable when l increases. All the findings indicate that one should select the nodes to make the mean values (degree, closeness and BC) of driver nodes around those of the network.
Finally, the relationship between R, l r 1 , l r N and l (See Materials and Methods) is investigated in terms of minimizing R. Fig. 10 shows that R(l)!l {c , which can help to predict R when knowing The measurements of Q and Q m are provided in (1) and (2). ''Order'' is obtained by sorting Q and ''Score'' is obtained by sorting Q m in an ascending way. doi:10.1371/journal.pone.0041375.t001 l. Moreover, in order to minimize R under a small l, l r N should be suppressed near a constant value and l r 1 should be enlarged as much as possible. As l increases, both l r N and l r 1 grow exponentially and the growth of the amplitude of l r 1 is larger than that of l r N . Fig. 10 illustrates that the shape of R largely depends on l r 1 . The observed phenomena indicate that l r 1 plays a more important role in minimizing R than l r N does. When l~N, it is shown that l r 1 &l r N , which makes R&1. In summary, when minimizing R, enlarging l r 1 is more important than reducing l r N .
This finding is similar to our finding in [43], where only undirected complex networks are studied.

Controlling Nodes of the Cortical Network -a Macroscopic Way
By means of JaDE, we control the cortical network under different l in terms of minimizing R and s, respectively. Denote     Table 2. Fig. 11 shows that when minimizing R, the standard deviation of T R,i is large, which means that some nodes in the cortical network, such as VPc, 2 and AMLS, are of great importance to be controlled. Some areas are negligible to be selected as driver nodes, such as 20a, CGp and 5AI. When minimizing s, the standard deviation of T s,i is small and nearly all the areas in the cortical network are important for minimizing s. Hence, the controlling nodes are different from the usual hubs, which are generally selected from nodes with a large degree [22]. In addition, the controlling nodes in the case of minimizing R are different from those in the case of minimizing s (Table 2). In order to show what factors have impacts on selection of controlling nodes, Dk~k in {k out of each area in the cortical network is depicted in Table 2, where k in and k out can be referred to Materials and Methods. Table 2 shows that, when optimizing R, most of the controlling nodes are selected from the nodes with a large k in and a small k out . Therefore, the areas with Dkw0 should be considered as controlling nodes when

Controlling Communities of the Cortical Network -a Mesoscopic Way
In the following, we show which module/community is significant to be controlled in a mesoscopic way. According to Table 2   community, respectively. Tables 3 and 4 show that most of the areas in the community Auditory serve as CN. Specifically, when minimizing R, most of the areas in the community Visual work as CN and ICN, most of the areas in the community Somato-motor belong to ICN and WCN and most of the areas in the community Fronto-limbic serve as ICN and WCN. When minimizing s, most of the areas in the community Visual work as ICN and WCN, most of the areas in the community Somato-motor belong to CN and ICN and most of the areas in the community Fronto-limbic serve as WCN. From the above observations, when minimizing R, the importance of each community is listed in a descending order: Auditory]Visual]Somato{motor]Fronto{limbic. When minimizing s, the importance of each community is listed in a descending order: Auditory]Somato{motor]Visual] Fronto{limbic. Hence, although the community Auditory is sparsely connected with other communities and is the smallest community, it is the most important one to control the cortical network. The observed phenomenon indicates that community with sparse connection to other communities should be paid special attention to control the network efficiently. Figure 10. The relationship between log 10 R, log 10 l r 1 , log 10 l r N and log 10 l by JaDE. doi:10.1371/journal.pone.0041375.g010

Discussion
The cortical hubs are believed to play pivotal roles in the coordination of information processing in cortical networks. In previous studies, the identification and classification of hub regions have been analyzed in terms of node degree, structural motif, path length, clustering coefficient distributions and synchronization [20,22]. In these works, the intrinsic relationship between structural and functional connectivity is analyzed by using ensembles of neurons coupled by a cortical network of cats' brain. By means of statistical methods, the crucial importance of nodes and clusters are revealed to analyze the separation and integration of sensory information in the cerebral cat cortex [24,46].
Additionally, one of the major challenges for human is to control natural systems or networks efficiently. As a typical natural network, identifying controlling nodes of a realistic anatomical network of cat cortical connectivity is of crucial significance to provide insights into avoiding abnormal synchronization in typical neural diseases [8,9,12]. In the light of previous studies, the problem of identification of controlling nodes of cortical networks remains open.
In this study, we have investigated the identification of controlling nodes in a network representing the connectivity among cortical areas in cats' brain. The issue regarding controllability of the cortical network is converted into a combinatorial optimization problem [43]. A representative evolutionary computation method, JaDE, which is a self-adaptive and efficient algorithm to solve real-world optimization problems [38], is used to identify controlling nodes with an appropriate encoding scheme. The comparison with some well-known network-based methods and evolutionary computation methods is presented, revealing JaDE performs best among all the algorithms.
The controlling nodes of the cortical network are detected in microscopic, mesoscopic and macroscopic ways. Using such various scales will help us to understand the controllability of neuronal networks in depth. We have shown a close relationship of the number of driver nodes and the locations of the driver nodes, indicating a concave shape of the mean degree of driver nodes as an increase of the number of driver nodes. For low values of the number of driver nodes, the areas with a large degree govern the coordination dynamics of the network. As a whole, the nodes with a small degree are important to be selected as controlling regions, which is in contrast to the work in [22] and supports the finding in  [29,32]. More importantly, the most prominent community in the cortical network of cats is the community Auditory, which has sparse connections with other communities. The comparative results of two quantities for measuring controllability of complex networks are also investigated in detail.
The model and methods can be extended and improved in several ways. Firstly, it is meaningful to propose more efficient optimization methods to deal with controllability of cortical networks. Secondly, we have only focused on the highest level of cortical networks and thus large subnetworks [14,47,48] with other biologically realistic features [11,49,50] should be considered. Finally, the results should be applied to other realistic natural systems to illustrate controlling rules. The achievements would require further developments in neuroscience, in the theory of dynamical complex networks, in optimization methods as well as in control science.

Cortico-cortical Network of Cats' Brain
The cortico-cortical network of cats' brain is a biological network that describes the anatomical connectivity of cats' brain [18,19]. Here, we use a version of a dataset in [21]. The cat cerebral cortex can be divided into 53 cortical areas, linked by about 830 fibres of different densities into a weighted and directed complex network. It consists of four topological clusters that broadly agree with four functional cortical sub-divisions: visual cortex (16 areas), auditory (7 areas), somato-motor (16 areas) and fronto-limbic (14 areas). We also refer to the topological clusters as communities or modules. The community Auditory is sparsely connected while the communities Visual, Somato-Motor and Fronto-Limbic are densely connected among each other [16].

Model and Problem Formulation
We consider a reference evolution/state as follows: (a(t)):  This equation is general, since many real-world systems such as social networks, biological systems and other natural systems can be modeled as differential equations [30].
Then, the following model of a diffusively coupled array of identical systems is considered as a general complex network: where is a continuous vector function. k is the coupling gain of the network. In the coupling term, the node is connected through a generic output function h(x i (t)). The matrix G stands for the connectivity about the cortical network topology. The graph G is supposed to be directed, weighted and simple (without self-loops and multiple edges). Let weighted and directed matrix G~½g ij N i,j~1 be the adjacency matrix of graph G, which is defined as follows: for any pair i=j,g ij v0 if e(i,j)[E; otherwise, g ij~0 . g ii~{ P N j~1,j=i g ij (i~1,2, Á Á Á ,N). The adjacency matrix G can be converted into the Laplacian matrix L by neglecting the weights over the networks. For any pair i=j,l ij~{ 1 if e(i,j)[E; otherwise, l ij~0 . l ii~{ P N j~1,j=i l ij , (i~1,2, Á Á Á ,N). The output degree k out (i)~{ P N j~1,i=j l ij of a node i is the number of efferent connections that it projects to other nodes, and its input degree k in (i)~{ P N j~1,i=j l ji , is the number of the afferent connections it receives. Denote by m i~m the set of eigenvalues of G and assume that they are ordered in such a way that m r 1 ƒm r 2 ƒ Á Á Á ƒm r N . To control such a cortical network to the reference evolution a(t), feedback controllers are added to (5): where c i are control gains or coupling strengths. Suppose that 1ƒ P N i~1 d M (i)ƒN. We aim to lead the cortical network (5) toward the desired reference evolution a(t), i. e., By linear manipulations, the stability analysis of (6) can be transformed into the dynamics of N independent blocks in the parameters e i~k l i ,i~1,2, Á Á Á ,N [26,51,52], where Jf (a(t)) and J(h(a(t))) are the Jacobians of the functions f and h calculated around the time varying reference evolution a.
Without loss of generality, we assume that l r i are sorted as l r 1 ƒl r 2 ƒ Á Á Á ƒl r N . As pointed out in [25,26], through above transformation, the problem of controllability of complex networks is converted into synchronizability of networks. Similar to the analysis method of checking synchronizability of networks, the enhancement of controllability can be characterized by reducing the eigenratio.  s~maxfl m i g, as small as possible [25,44], i. e., the smaller the R and s are, the easier the network is controllable. Previous works have shown that s can be neglected, since usually s is very small and has only minor effects on synchronizability/controllability of networks [45]. We also consider s and illustrate the impact of s on controllability, since s is important when one considers some special graphs, e. g., normalized Laplacian graph. It should be noted that the selection of driver nodes is a typical combinatorial optimization problem [43], where the locations of driver nodes are discrete variables, and the design of control gains is a continuous optimization problem. Taking the locations of driver nodes and their control gains into account together, the controllability of networks can be viewed as a multimodal optimization problem.
Here, minimizing R and s by determining locations of driver nodes d M (i) and designing c i (i[M) can be formulated as follows: From the above equations, we study the controllability of cortical networks by minimizing R and s, respectively. Evolutionary computation methods are employed to study the controllability and identify controlling regions.

The Strategies for Determining the Locations of Driver Nodes
Several well-known strategies for determining the locations of driver nodes or controlling nodes are illustrated as follows [43].
(1) Degree-based strategies. Degree-based pinning schemes are the most popular methods to select potential driver nodes, in which the locations of driver nodes are chosen according to degree information of networks in a decreasing or an ascending way [25,30,53]. Here, the two schemes are called ascending and descending degree-based strategies, respectively. The output degree k out is used to provide degree information.
(2) Betweenness centrality (BC)-based strategies. Similar with the degree-based scheme, we consider descending and ascending BC-based strategies. (3) Closeness-based strategies. Two kinds of closeness-based strategies, i. e. descending and ascending closeness-based strategies are taken into account. (4) Node importance-based strategies. Since the controllability of the cortical network is mainly related to its eigenvalues, it is interesting to determine the locations of driver nodes by considering their importance in the network [54]. We analyze two measures of node importance for the cortical network. The first one is to minimize strategy. The other one is to minimize S~maxfm m i g of G upon sequential removal of nodes, which is called S-based strategy. It should be noted that U and S are usually used to measure synchronizability performance of complex networks [44]. (5) Evolutionary algorithm-based strategies. Using an appropriate encoding scheme, differential evolution (DE) is used to select driver nodes and design control gains. Evolutionary algorithms have been successfully used in the synchronization of two coupled systems in [55], the coordination of unmanned aircraft vehicle [43] and networks topology with optimal synchronizability [56]. Here, adaptive differential evolution is adopted to identify the controlling nodes [38].
In the degree-based, the BC-based, the closeness-based and the node importance-based strategies, control gains in all the nodes are considered to be identical and one can tune the control gains of driver nodes in the cortical network gradually with a step size 0.1, like [25,26].

Differential Evolution and its Encoding Scheme
In order to determine the locations of driver nodes in the cortical network and design their control gains, an appropriate encoding scheme is used according to [43]. In addition, equipped with this encoding scheme, JaDE [38] is used to detect the controlling nodes/areas/regions of the cortical network of cats' brain in microscopic, mesoscopic and macroscopic ways, respectively.