This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

1

Robust Multiobjective Controllability of Complex Neuronal Networks Yang Tang, Member, IEEE, Huijun Gao, Fellow, IEEE, Wei Du, Jianquan Lu Athanasios V. Vasilakos, Senior Member, IEEE, and Jurgen Kurths ¨ Abstract—This paper addresses robust multiobjective identification of driver nodes in the neuronal network of a cat’s brain, in which uncertainties in determination of driver nodes and control gains are considered. A framework for robust multiobjective controllability is proposed by introducing interval uncertainties and optimization algorithms. By appropriate definitions of robust multiobjective controllability, a robust nondominated sorting adaptive differential evolution (NSJaDE) is presented by means of the nondominated sorting mechanism and the adaptive differential evolution (JaDE). The simulation experimental results illustrate the satisfactory performance of NSJaDE for robust multiobjective controllability, in comparison with six statistical methods and two multiobjective evolutionary algorithms (MOEAs): nondominated sorting genetic algorithms II (NSGA-II) and nondominated sorting composite differential evolution (NSCDE). It is revealed that the existence of uncertainties in choosing driver nodes and designing control gains heavily affect the controllability of neuronal networks. We also unveil that driver nodes play a more drastic role than control gains in robust controllability. The developed NSJaDE and obtained results will shed light on the understanding of robustness in controlling realistic complex networks such as transportation networks, power grid networks and biological networks, etc. Index Terms—Synchronization, Neuronal networks, Controllability, Robustness, Multiobjective optimization.

F

1

I NTRODUCTION

Neuronal networks, as a special type of directed and weighted complex networks [1]–[9], have been put forward for the investigation of high-level information processing in cortical networks [1], [10]–[12]. By means of data extracted from structural and functional MRI, diffusion tensor imaging, magnetoencephalography and electroencephalography in humans and nonhuman brains, structural and functional topologies have been constructed, which have shown different properties such as small-world topology, scale-free property, highly connected hubs and communities, etc [1], [10], [12]–[14]. The complicated structural network demonstrates the dynamic emergence of coherent physiological dynamics such as phase-locked high-frequency electromagnetic oscillations and synchronization [12], [15]–[21]. Synchronization of neural elements plays a vital role in neural information processing [16]–[18], [20], [22]. In particular, in [17], it was • Yang Tang is with the Key Laboratory of Advanced Control and Optimization for Chemical Processes, Ministry of Education, East China University of Science and Technology, Shanghai 200237, China. E-mail: [email protected]; [email protected]. • Huijun Gao is with the Research Institute of Intelligent Control and Systems, Harbin Institute of Technology, Harbin 150080, China. He is also with the King Abdulaziz University, Jeddah 21589, Saudi Arabia. E-mail: [email protected]. • Wei Du is with the Institute of Textiles and Clothing, The Hong Kong Polytechnic University, Hong Kong, China. E-mail: [email protected]. • Jianquan Lu is with the Department of Mathematics, Southeast University, Nanjing 210096, China. E-mail: [email protected]. • Athanasios V. Vasilakos is with the Deptment of Computer Science, Electrical and Space Engineering, Lulea University of Technology, Lulea 97187, Sweden. E-mail: [email protected]. • Jurgen ¨ Kurths is with the Potsdam Institute for Climate Impact Research, Potsdam, Germany, Institute of Physics, Humboldt University of Berlin, Berlin, Germany and Department of Control Theory, Nizhny Novgorod State University, Gagarin Avenue 23, Nizhny Novgorod 606950, Russia. E-mail: [email protected].

found by a set of experiments that several diseases, such as epilepsy, autism, schizophrenia, Alzheimer’s disease, and Parkinson’s disease, have a close relationship with abnormal neural synchronization. Subsequently, it was uncovered in [23] that phase synchrony has been observed for distinct frequency bands among frontoparietal and visual areas of a brain, which elucidates a coordination behavior of neuronal object representations in visual working memory. As a close topic related to synchronization (leaderless consensus) [5], [7], controllability of complex networks (leader-following consensus) has recently gained increasing attention in natural science and engineering [24]. It has been widely appreciated that controlling the dynamics of a complex network is a common demand in many real networks, such as biological transportation, social and neuronal networks [5]. Up to now, numerous definitions of controllability and control strategies have been developed to assess the controllability of complex networks. For recent advances in this area, one can refer to the recent survey [5] and the references therein. As demonstrated in [23]–[25], it is profound to examine the controllability of a neuronal network, which will provide useful guidelines to control other kinds of realistic weighted and directed networks and grasp the potential mechanisms to keep away from anomalous synchronization so that neural diseases such as schizophrenia, epilepsy, autism, Alzheimer’s disease, and Parkinson’s disease can be inhibited accordingly [17], [26]. In this paper, the problem of robust multiobjective controllability of a neuronal network is examined. It is worth mentioning that robust controllability of complex networks has strong medical meaning [27]. For example, in a therapy, the patient’s recovery largely depends on the proper diagnosis and the suitable dosage of antibiotics, where the proper diagnosis can be viewed as the correct identification of the pathogenic areas (organs) and the input of suitable dosage can be viewed as the right control gains. Apparently, errant diagnosis will result in a delay of the recovery and

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

2

even make the symptoms severely. Besides, the excessive injection of dosage of drugs will bring about the creation of multidrug-resistant bacteria and finally no efficient antibiotics can be found in some severe circumstances. On the other hand, an insufficient injection of dosage will not be favorable to patients’ recovery and drag out the recovery time of patients. Based on these discussions, the keerect diagnosis of pathogenic factors (the correct identification of controlling areas) should be made and a suitable dosage (the appropriate control gain) should be injected, which will make the rehabilitation of patients more efficient. Another example for cancer therapy can be referred to [27]. In searching for a candidate for solving robust multiobjective controllability of neuronal networks, multiobjective optimization evolutionary algorithms (MOEAs) emerge as a promising one [28]–[30]. Among them, the nondominated sorting genetic algorithm-II (NSGA-II) [29] is a paradigm due to its efficiency in solving numerous multiobjective problems. Based on NSGA-II, in [30], robust multiobjective optimization was addressed, whose emphasis is to find a robust Pareto frontier, instead of the global Pareto-optimal frontier. Robust solutions to controllability of a neuronal network are illustrated as robust Pareto fronts (RPFs) by means of a nondominated sorting adaptive differential evolution (NSJaDE). The organization of this paper is as follows. Sec. II reviews the prior work related to robust control and controllability of complex networks. Sec. III presents the models and methods of the robust multiobjective controllability of a neuronal network. In Sec. IV, NSJaDE is presented in detail and a set of computer experiments are carried out to manifest its effectiveness of NSJaDE. The main findings are also reported in this Section. Conclusions are drawn in Sec. V.

2

R ELATED

WORK

To the best of authors’ knowledge, robust controllability of complex networks has been overlooked primarily due to the lack of a suitable definition of uncertainty controllability and the development of an appropriate strategy in optimizing controllability. Here, we first review related work in classical robust control theory. Then, we will review recent progress in controllability of dynamical networks and provide comparisons in detail. As an important component in modern control theory, robust control recently has attracted much attention in the last three decades, see a representative book [31] and a recent survey [32]. Nowadays, it is widely recognized that the inclusion of uncertainty in the mathematical model of a practical system is deemed as robustness. Uncertainty serves as a central part of feedback and cannot be ignored. In control systems, uncertainty usually arises from measurement errors, quantization errors and modelling errors [31], [32]. Controllers which ensure an adequate level of control performance in the context of bounded uncertainty are robust controllers [32]. In addition to the broad applications to control systems, uncertainty also widely exists in physics, statistics, economics, finance, insurance, psychology, sociology, and information science, etc. In [24], the controllability of complex networks with linear dynamics was investigated by using Kalman’s controllability. In [27], the controllability of complex networks with nonlinear dynamics was handled by using the concept of

region of attraction. In [33], the controllability of complex networks was tackled by means of linearization. Following this framework, the controllability of neuronal networks was investigated by a constraint optimization evolutionary algorithm [25], [34] and a multiobjective optimization evolutionary algorithm [19], respectively. In [35], although robust distributed pinning synchronization of agent systems has been investigated, the uncertainty only locates in the system dynamics, in which the robustness of pinning controllers is not considered. A detailed survey of controllability of complex networks has been provided in a recent survey [5]. Up to now, the robust controllability of complex networks or robust leader-selection of agent systems has not been investigated [5], [32]. In detail, there are two momentous factors in controllability of networks, i. e., the selection of driver nodes and the design of control gains. It is natural to ask how the uncertainties in these two ingredients affect control performance. In fact, the uncertainty in choosing driver nodes can be understood as the level of accurate implementation of actuators and the uncertainty in designing control gains can be comprehended as implementation errors in control gains. The misplacement of controllers and overlarge/oversmall design of control gains will inevitably lead to the deterioration of the control performance of the entire network. Consequently, it is of paramount importance to investigate the robust controllability of networks, which is also related to fault-tolerant control [36]. Based on the above literature review, the main purpose of this paper is to tackle robust controllability of neuronal networks, which improves our recent works [19], [25], [34] and provides in-depth insights in robust controllability of complex networks. The contributions of this paper can be summarized as follows: (1) compared with the works in controllability of networks [19], [24], [25], [27], [34], [35], [37]– [39], the problem of robust multiobjective identification of a neuronal network is investigated for the first time, in which the uncertainties in selection of controlling regions and design of control gains are considered; (2) several interesting findings are provided regarding the robust multiobjective identification of a neuronal network, which illuminates the importance of introducing robustness in controllability; (3) a nondominated sorting adaptive differential evolution is proposed to improve the search performance of NSGAII [29] and NSCDE [19].

3

M ODELS

AND METHODS

Notations: In this paper, consider a graph by G = [V, E], where V = {1, · · · , N } and E = {e(i, j)} are the node set and the edge set, respectively. The graph G is supposed to be directed, weighted and without self-loops. Denote the weighted and directed matrix G = [gij ]N i,j=1 to be the coupling matrix of graph G: for any pair i ̸=∑j, gij < 0 if e(i, j) ∈ E; otherwise, gij = 0; gii = − N j=1,j̸=i gij (i = 1, 2, · · · , N ). The Laplacian matrix L is defined as follows: for any pair ∑ i ̸= j, lij = −1 if e(i, j) ∈ E; otherwise, lij = 0. lii = − N j=1,j̸=i lij , (i ∈ V). l ∈ [1, N ] is the number of driver nodes, where N is the size of the network. ID (·) stands for the characteristic function of the set D, i.e., ID (i) = 1 if i ∈ D; otherwise, ID (i) = 0. The output-degree ∑ of node i is kout (i) = − N lij and the input-degree of ∑ j=1,i̸=j node i is kin (i) = − N l . ji j=1,i̸=j

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

3.1

3

Controllability of neuronal networks

3.2 Robust Multiobjective Controllability

Consider the following reference state: (1)

s(t) ˙ = f (s(t)).

Our aim here is to force the state of the neuronal network toward the reference state s(t) described by (1) by using an output feedback controller. Then, the neuronal network with an output feedback controller is listed as follows: x˙ i (t)

f (xi , t) − c

=

N ∑

gij h(xj (t))

j=1

− cID (i)mi (h(s(t)) − h(xi (t))), i ∈ V,

(2)

where xi (t) = [xi1 (t), xi2 (t), · · · , xin (t)]T ∈ Rn (i ∈ V) is the state of the ith area and f (xi , t) = [f1 (xi , t), · · · , fn (xi , t)]T is a continuous vector function; c is the coupling strength; h(xi (t)) is the output function; mi (i ∈ V) is the control gain between the area and the reference state s(t) and G is the coupling matrix of the neuronal network of a cat’s brain. The neuronal network is weighted and directed, which is extracted by cortical parcellation, thalamic parcellation, collation of connection data and translation from database. For more details regarding the statistical and biological properties of the neuronal network, one can refer to [25], [40]–[43] and the references therein. There are numerous models to represent f in the literature [2], [20]. For example, in [20], the FitzHughNagumo model is used for representing f in (2). This way, the dynamics in the cortical brain network of a cat can be then modeled. √ Let ξi = ξir + jξim (j = −1, i ∈ V) be the eigenvalues r of G and ξi are ordered by ξ1r ≤ ξ2r ≤ · · · ≤ ξN . By following the way in [25], [33], an extended network of N +1 dynamical systems yi can be formulated to investigate the controllability of networks, where yi = xi for i ∈ V and yN +1 = s: y˙ i (t)

=

f (yi , t) − c

N +1 ∑

Tij h(yj (t)),

j=1

i ∈ V ∪ {N + 1},

(3)

where T = [Tij ] ∈ R(N +1)×(N +1) having the following form     T =  

T1 g21 .. . gN 1 0

g12 T2 .. . gN 2 0

... ... .. . ... ...

g1N g2N .. . TN 0

−ID (1)v1 −ID (2)v2 .. . −ID (N )vN 0

    ,  

(4)

in which Ti = gii + ID (i)mi . Denote λi = λri + jλm i as the ith eigenvalue of T . Similar to ξi , λi is sorted by λr1 ≤ λr2 ≤ · · · ≤ λrN +1 , where λr1 = 0. Hence, based on the analysis [33], the local controllability can be measured by minimizing the following two objectives at the same time: λr f˜1 = min Nr+1 , λ2

(5)

f˜2 = min max{λm i }.

(6)

and i

By using the framework in Sec. 3.1, the controllability of the neuronal network was investigated by means of a constraint optimization evolutionary algorithm (COEA) [25] and a nondominated sorting composite differential evolution (NSCDE) [19], respectively. In addition, by considering the constraint on control gains, the controllability of the neuronal network has been further investigated in [34]. However, as mentioned in [5], the existing results regarding controllability have not taken into account the robustness despite the fact that the importance of robustness has been witnessed in various areas such as physics, statistics, economics, finance, insurance, psychology, sociology, and information science, etc. The purpose here is to make a first attempt to find robust solutions for controllability, which are less sensitive to uncertainties or perturbations in decision variables, i. e., selection of driver nodes and injection of control gains. Actually, the existence of uncertainty in locating driver nodes may be due to a possible misplacement of actuators and the uncertainty in control gains may arise from the inaccuracy of injection of control inputs. The uncertainties in driver nodes and control gains will deteriorate the control performance. Optimal control may be sensitive to disturbances. Degradation in optimal control may enhance robust controllability in applications. In addition, the robust multiobjective controllability has a relationship with fault-tolerant control, i. e., when the controller is inputted at one node nearby the true driver node due to several physical constraints or fault operations, is the network still controllable? For an example of therapy, the patient’s recovery heavily hinges on the accurate diagnosis and the precise dosage of antibiotics. The proper diagnosis can be regarded as the correct detection of the pathogenic areas (organs) and the dosage is called control gains. As well known, faulty diagnosis would delay the recovery process and even lead to more severe symptoms. In addition, one area may not be operated because of some potential physical constraints. The physical constraints and faults would result in the case that one area or several areas are not suitable for operation and a nearby area is operated on purpose or incorrect manipulation. Besides, an excessive or insufficient injection of dosage of drugs will lead to the creation of multidrug-resistant bacteria or prolong the recovery of patients, respectively. Hence, it is momentous to investigate robust multiobjective controllability of neuronal networks, which would provide deep insights in other practical applications. Another example is from the controllability of a power grid network. Usually, the state of a power grid system can be assumed to be described by the swing equation. In [24], the controllability of a power grid network was investigated. When a single-line fault caused by short circuits to the ground occurs and the fault is cleared by disconnecting the line. However, by the time the fault is cleared, the generators lost synchrony and they continue accelerating away from one another. Fortunately, the perturbed network still admits a stable steady-state synchronous solution with frequency wd . In order to achieve synchronization of the power grid network, some small perturbations are added to generator frequencies to drive frequencies of generators to the target frequency state wd . In a practical setting, the perturbations can be implemented by tuning the damping

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

4

of the generators. In light of [27], [44], one can also apply our techniques to identify the key nodes in a power grid network when uncertainties are taken into account. For example, our results can be employed to show how the locations of perturbations in a power grid network (like the locations of driver nodes in this paper) and the values of perturbations (like the values of control gains here) affect tracking performance of the desired frequency wd . In the following, we will give some preliminaries and methods regarding robust multiobjective controllability of neuronal networks. Refs. [31], [32] systematically discuss the results and advances of uncertainties. Usually, uncertainty can be classified into unstructured and structured uncertainties [32]. Generally, for a matrix A, it suffers from uncertainty q = [q1 , ..., ql ]T , which could be time-varying, is assumed to be unknown but contained in a given set Bq . In control theory and optimization [32], [45], there are several popular ways to characterize uncertainties like norm bounded uncertainties and interval uncertainties. In this paper, according to the encoding scheme used in the following, we consider interval uncertainties to account for robustness in solutions. It should be mentioned that we have considered norm bounded uncertainties in the dynamics of each node of networks [35]. It is promising to take into account norm bounded uncertainties when considering controllability in selection of nodes in the near future. For a matrix A ∈ Rn1×n2 , define AI

=

¯ [A, A]

=

{A = [aij ]n1 ×n2 : aij ≤ aij ≤ a ¯ij , 1 ≤ i ≤ n1 , 1 ≤ j ≤ n2 },

(7)

where A

=

[aij ]n1 ×n2 , A¯ = [¯ aij ]n1 ×n2 , ¯ij , aij ≤ a 1 ≤ i ≤ n1 , 1 ≤ j ≤ n2 .

(8)

Define 1 ¯ A˜ = 1 (A¯ − A) = [˜ (A + A), aij ], 2 2 n1 n2 ∑∑ ei aϵij hTj , |aϵij | ≤ a ˜ij , (9) A = Ao + Aϵ = Ao +

3.3.1 Concepts of Robust Multiobjective Optimization As stated in [30], the authors adopt the mean effective objective function for the robust multiobjective optimization problem instead of the objective function itself. Hence the robust multiobjective controllability problem can be formulated as follows: Problem: A solution x∗ is called a multiobjective robust solution, if it is the global feasible Pareto-optimal solution to the following multi-objective minimization problem of a solution x: min F (x), x ∈ RD , x∈Ω

(10)

where F = {f1 , f2 , · · · , fm }, ∫ 1 fi (x) = f˜i (y)dy, |Bδ (x)| y∈Bδ (x)

(11)

Bδ (x) is a δ-neighborhood of the solution x, |Bδ (x)| is the hypervolume of the neighborhood, and D is the dimension size of an optimization problem, Ω is the feasible decision space, x = {x1 , x2 , · · · , xD } is a vector with a set of decision variables and F = {f1 , f2 , · · · , fm } is the objective vector with m objectives to be minimized. For the robust multiobjective controllability problem, m is 2 and f˜i (i = 1, 2) is given in (5) and (6), respectively. fi is the mean effective objective function for i = 1, ..., m. An example for illustration of this problem and robust solutions can be seen from Fig. 1. For calculation of fi , a practical way is to yield a set of H solutions in a random or structure-constraint manner, which are chosen around a δ−neighborhood Bδ (x) of a solution x in the decision space and the mean function value fi is employed as the fitness value for an EA. However, the introduction of robustness in an objective evaluation would inevitably lead to H times more function evaluations than the usual approach. In the following, we will discuss how to obtain the mean effective objective fi .

Ao =

i=1 j=1

where ei or hj denotes the column vector with the ith element being 1 and the others being 0. aϵij is the actual disturbance in each element. It should be mentioned that the interval uncertainties described above are frequently confronted in control systems, which may result from the fluctuation of operating points, impairment or deterioration of the devices and identification errors, etc. In addition, we can also include statistical information to characterize interval uncertainties. For example, Gaussian distribution or Cauchy distribution can be utilized to model interval uncertainties under a truncated interval. It is promising to study stability analysis and control of dynamical systems under such a characterization by means of the Lyapunov stability theory and/or evolutionary computation. 3.3 Robust Multiobjective Optimization In this subsection, the concept of robust multiobjective optimization and an improved MOEA are presented.

3.3.2 NSJaDE In [29], a multiobjective optimization algorithm called NSGA-II was presented, which is made up of a fast nondominated sorting approach and a diversity preservation mechanism. NSGA-II has been widely acknowledged as a paradigm due to its reduction of computational complexity and its powerful adaptive mechanism. In [30], NSGA-II was further extended to robust multiobjective problems. Although NSGA-II has been manifested its efficient performance in various theoretical and practical problems, its search capability could be improved. In [19], a nondominated sorting composite differential evolution (NSCDE) was developed to handle multiobjective controllability of neuronal networks, which shows its superiority over NSGAII. Nonetheless, the composite differential evolution is not adaptive. Hence NSCDE can be further improved by adding several potential adaptive operations and thus possesses exploration and exploitation search capability. Among the existing well-performed search engines, adaptive differential evolution [46] has been validated its accuracy in dealing with a number of benchmarks and engineering problems, while maintaining satisfactory convergence

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

5

(a)

(b)

Fig. 1. An illustration of robust solutions from the decision space to the objective space is shown. To qualify as a robust solution, each Pareto-optimal solution has to show its insensitivity towards small perturbations in its decision variable values. As a small perturbation of point B causes a large variation in objective values, solution B may not serve as a robust solution. Solution A is less sensitive to variable perturbation than Solution B in a multiobjective optimization context, i. e., solution A is more robust to perturbations in variable than solution B when optimizing the two objectives f1 and f2 .

speed. At the beginning, a population of SP individuals is initialized in a D-dimensional decision space. Each individual can be regarded as a chromosome in genetic algorithms (GAs), standing for a potential solution. JaDE is composed of three procedures: mutation, crossover, and selection operators. In order to detect the global optimum, each procedure is performed at each generation to search the space. The population with its individuals is represented as Pt = (x1,t ; x2,t ; ...; xi,t ), i = 1, 2, ..., SP, t = 0, 1, 2, ..., tmax , and xi,t = (x1i,t , xji,t , ..., xD i,t ), j = 1, 2, ..., D, where t is the generation number. In this paper, JaDE is employed to replace the search engine in NSGA-II for enhancing the ability of exploration and exploitation. A mutation strategy and an external archive are taken into account to supply real-time information of the progress direction. It should be mentioned that, in JaDE, the DE/current-to-pbest scheme is utilized to balance the convergence speed and the diversity of the population by considering multiple best solutions. The nondorminated sorting adaptive differential evolution (NSJaDE) is developed by conserving the nondorminated sorting mechanism, the diversity preservation approach in NSGA-II and the replacement of simulated binary crossover (SBX) operator and polynomial mutation by JaDE. The procedure of the proposed constrained multiobjective optimization algorithm NSJaDE is provided in Algorithm 1, which

Algorithm 1 Robust constrained NSJaDE Begin Create a random population Pt (t = 0) in a feasible solution space with SP individuals. while t ≤ tmax do Spatt = neighborhood-pattern-generation(H, δ); /*generating a neighborhood pattern for the current population, where H is the number of neighboring points, δ is the extent of neighborhood*/; Wt = tournament-selection(Pt ); /*using a binary tournament selection to generate a parent population with ⌊SP/2⌋ individuals*/; fP,t = evaluate-objective(Pt , Spatt ); /*evaluate the objective values of the parent population according to the generated neighborhood pattern*/; cP,t = evaluate-constraint(Pt , Spatt ); /*evaluate the constraint values of the parent population according to the generated neighborhood pattern*/; Qt = adaptive-differential-evolution(Wt ); /*offspring Qt is generated according to JaDE and the population size of Qt is SP */; fQ,t = evaluate-objective(Qt , Spatt ); /*evaluate the objective values of the offspring population according to the generated neighborhood pattern*/; cQ,t = evaluate-constraint(Qt , Spatt ); /*evaluate the constraint values of the offspring population according to the generated neighborhood pattern*/; St = Wt ∪ Qt ; /*combine parent and offspring population*/; F =constrained-nondominatedsort(St , fP,t , cP,t ,fQ,t , cQ,t ); /*F = (F1 , F2 , ...) contains all nondominated fronts of St */; Pt+1 = ∅ and i = 1 while |Pt+1 | + |Fi | ≤ SP do /*|.| is the cardinality of a set*/; crowding-distance-assignment(Fi ); /*calculating crowding distance of Fi */; Pt+1 = Pt+1 ∪ Fi ; /*include ith nondominated front in the population*/; i = i + 1; sort(Fi , ≺); /*sort the individuals in Fi using ≺*/; Pt+1 = Pt+1 ∪ Fi [1 : (SP − |Pt+1 |)]; /*fill Pt+1 with the best (SP − |Pt+1 |) individuals in Fi */; t = t + 1. end while end while End

considers a more general case by including the constraint into NSJaDE. In Algorithm 1, neighborhood-pattern-generation(,.,) is to calculate f˜i , (i = 1, 2) in (5) and (6) to produce the mean function value of fi , (i = 1, 2) in (11). From the expression of neighborhood-pattern-generation(H, δ), it is clear that fi depend on the extent of neighborhood size δ and the number of neighboring points H. There are several ways to generate H neighboring points in the proximity of a solution. A simple and basic way is to generate H points according to Monte Carlo approach. However, this inevitably introduces additional randomness in evaluating

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

6

TABLE 1 Statistical strategies in complex network theory and their abbreviations. Methods

Abbreviation

Ascending degree-based method Descending degree-based method Ascending betweenness centrality-based method Descending betweenness centrality-based method Ascending closeness-based method Descending closeness-based method

ADM DDM ABM DBM ACM DCM

TABLE 2 The effect of γ on the mean (standard deviation) and minimum values of f1 (f2 ) and their corresponding f2 (f1 ) in 20 runs are shown here, where H = 5 and δ = 1. Fig. 2. An example of creating H = 5 neighboring points around a solution x = (a, b).

γ = 150 minimizing f1

l = 5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

the same solution more than once [30]. Here, we follow the way in [30]. A Latin hypercube strategy is utilized to propagate a pattern systematically. For the multiobjective controllability of the neuronal network, the perturbed population Z k , (k = 1, 2, ..., H) of the unperturbed population Pt can be represented by the following matrix: Zk

=

=

(z1 ; z2 ; ...; zSP )  z1,1 z1,2  z2,1 z2,2   .. ..  . . zSP,1 zSP,2

... ... .. . ...

z1,D z2,D .. . zSP,D

l = 5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

    

(12)

We assume that Z k ∈ ZI is subjected to interval uncertainty in Sec. 3.2. That is, for a perturbed matrix Z k ∈ RSP ×D , define ZI

=

¯ [Z, Z]

=

{Z = [zij ]SP ×D : z ij ≤ zij ≤ z¯ij , 1 ≤ i ≤ SP, 1 ≤ j ≤ D}.

(13)

where Z and Z¯ are defined like Sec. 3.2. Hence, Z k = Zo + Z ϵ,k = Zo +

SP ∑ D ∑

ϵ,k T ϵ,k ei zij hj , |zij | ≤ z˜ij ,

l = 5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

minimizing f2 f2 f1 f2 0.5275±0.1186 44.0668±4.2237 0.1722±0.0354 0.5722±0.1509 28.4634±4.3187 0.0757±0.0130 0.5730±0.1183 21.4352±2.2152 0.0640±0.0290 0.5092±0.1069 18.3221±2.1936 0.0370±0.0228 0.5833±0.1125 17.5031±2.7199 0.0620±0.0269 0.6367±0.1450 14.2458±2.6811 0.0738±0.0263 0.6369±0.1108 15.3583±4.2474 0.0917±0.0375 0.6452±0.0713 14.0388±1.1768 0.1076±0.0260 0.6968±0.1500 13.1660±1.4754 0.1066±0.0333 0.6730±0.1599 14.1598±2.4728 0.0933±0.0310 γ = 200 minimizing f1 minimizing f2 f1 f2 f1 f2 35.7263±0.5095 0.5161±0.1274 45.4990±4.9795 0.1657±0.0390 19.5515±0.5312 0.5337±0.1631 28.5864±4.5944 0.0715±0.0196 14.6406±0.4917 0.5717±0.0959 21.4522±2.1600 0.0552±0.0319 12.0156±0.4366 0.5732±0.0974 18.3880±2.2422 0.0357±0.0210 10.2899±0.1988 0.6228±0.1063 17.5520±2.7333 0.0519±0.0196 9.0319±0.2712 0.6176±0.1521 14.4369±2.8500 0.0694±0.0286 8.1508±0.2659 0.6474±0.1330 13.9041±1.8878 0.0801±0.0294 7.5523±0.3178 0.6591±0.0766 13.3921±1.3020 0.0958±0.0247 6.9948±0.1352 0.7342±0.1490 12.1063±1.3223 0.0913±0.0276 6.5065±0.1634 0.6758±0.1846 13.8334±2.4948 0.0894±0.0353 γ = 250 minimizing f1 minimizing f2 f1 f2 f1 f2 36.0045±0.8465 0.5201±0.0979 48.0784±4.4632 0.1285±0.0240 19.3815±0.4993 0.4130±0.0464 25.9031±4.2918 0.0805±0.0250 14.6498±0.4244 0.4419±0.1070 20.6629±2.0280 0.0507±0.0266 11.7706±0.5951 0.5367±0.1625 19.1023±2.8639 0.0422±0.0269 10.2228±0.3251 0.6041±0.0798 17.7630±3.6113 0.0683±0.0321 9.0426±0.4608 0.5979±0.1292 13.3545±1.3263 0.0748±0.0529 8.0764±0.2250 0.6666±0.1943 14.8522±2.8404 0.0469±0.0367 7.6234±0.2857 0.7142±0.1459 15.3182±2.2064 0.0606±0.0320 6.8891±0.2640 0.6560±0.1157 12.7325±0.8642 0.0878±0.0306 6.4454±0.2890 0.6486±0.0849 14.8103±1.0164 0.0516±0.0463 f1 35.8037±0.4741 19.6084±0.5281 14.6924±0.5436 12.0839±0.4826 10.3662±0.2434 9.0646±0.2936 8.1995±0.2762 7.5681±0.3273 7.0398±0.1717 6.5804±0.1694

(14)

i=1 j=1

where z˜ij = δj , which means that in the considered problem, the allowed perturbation region for each individual is the same and may be different from each dimension. According to this formulation, the perturbation domain of each variable is divided into H equal grids. Hence, the δ-neighborhood is separated into H D hyperboxes. H hyperboxes are randomly picked from H D hyperboxes so that all H distinct grids can be represented in some degree. An example of the pattern for a two-variable problem is shown in Fig. 2. Once the hyperboxes are determined, a random point within each hyperbox is chosen and is used for the calculation of the mean effective objective values fi (i = 1, 2). A random pattern is created at the beginning of each generation and all solutions are evaluated by using this pattern during the generation. For each generation, another random pattern is generated. The nondominated-sort(.) function is used to generate different levels of robust Pareto fronts according to Pareto

dominance, while the diversity preservation approach is to make the solutions widely spanned in the obtained set of solutions. It should be mentioned that in this paper, we do not use the number of fitness evaluation to control the termination of NSJaDE or other multiobjective optimization evolutionary algorithms. The reason is that as the inclusion of H is used to calculate the mean effective objectives, the number of fitness evaluation fe relies on the value of H. For simply control the termination of MOEAs, we adopt the generation number t and set the maximum number of generation tmax .

4

C OMPUTER E XPERIMENTS

In this section, in order to illustrate the superiority of NSJaDE, we compare the proposed NSJaDE with some representative statistical methods in complex network theory (shown in Table 1) [25], and two representative MOEAs: NSGA-II [30] and NSCDE [19]. In addition, we illustrate the importance of δ and H on the robust Pareto front. The controlling regions are identified under different levels of δ.

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

7

TABLE 3 Comparison of NSJaDE with statistical methods in Table 1 for l = 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, when H = 5 and δ = 1. The comparison results are shown in minimizing f1 and f2 , respectively. The results of NSJaDE are shown by calculating the mean values (standard deviations) of the minimum f1 (f2 ) in 20 runs. Minimizing f1 f1 76.8747 45.0113 35.2532 32.1015 28.4064 25.2918 24.2035 16.9005 15.4636 14.3258

l = 5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

f1 94.4355 46.1946 36.5891 36.2013 30.4026 27.2327 28.2020 24.5626 20.6922 22.7219

f2 0.759 0.6207 0.6740 0.7287 0.3781 0.6457 0.6472 0.5710 0.8227 0.6453

f1 66.0612 39.6294 30.4633 21.5784 18.6131 16.8627 15.5792 15.5599 14.0726 13.7921

f2 0.5839 0.6344 0.6018 0.6273 0.4149 0.4072 0.3858 0.3539 0.3547 0.3093

f1 66.0460 39.8979 36.1268 24.2572 22.3118 20.6487 20.1535 21.2064 16.0364 19.171

DDM

DCM f2 0.8173 0.6876 0.5062 0.5940 0.5760 0.4507 0.6923 0.6818 0.8068 0.9085

f1 65.4387 30.7244 18.3740 15.3118 13.7229 11.6789 10.6922 12.0700 11.4055 11.6097

f2 0.6272 0.5801 0.3953 0.3565 0.1528 0.4249 0.2859 0.2658 0.5672 0.4522

f1 69.9349 49.1908 21.2191 18.5432 17.5339 17.5122 19.0308 16.2029 15.2826 15.0461

DBM

f2 0.7880 0.6496 0.8159 0.9141 0.9934 0.9302 0.9041 0.8728 0.9195 0.7439

DCM

0.8

f2 0.7081 0.5215 0.4112 0.4905 0.3615 0.5961 0.6020 0.4612 0.4493 0.4329

1

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.8

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.8

0.2

(a)

45

0 15

50

20

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.6 0.4

25 f1 30

35

0 10

40

15

20 f1

0.8

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.6

f2

0 5

10

15

f1 (f)

1

20

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.8

0 0

f2

f2

f2

0.2

0.2

(h)

15

20

0 0

5

f1 10 (g)

20

25

15

20

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.8

0.4

f1 10

15 f1

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

1

0.4

5

10

0.6

0.6

0.2

0 5

0.4

0.4

0 0

0.6

30

0.6

f2

20

25

0.8

f2

15

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.6

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

(d)

0.2

0.8

f2 0.1657±0.0390 0.0715±0.0196 0.0552±0.0319 0.0357±0.0210 0.0519±0.0196 0.0694±0.0286 0.0801±0.0294 0.0958±0.0247 0.0913±0.0276 0.0894±0.0353

0.8

(c)

0.2

f1 (e)

NSJaDE

0.2

0.2

10

f1 45.499±4.9795 28.5864±4.5944 21.4522±2.1600 18.388±2.2422 17.552±2.7333 14.4369±2.8500 13.9041±1.8878 13.3921±1.3020 12.1063±1.3223 13.8334±2.4948

0.2

0.4

0 5

f2 0.6230 0.5594 0.5001 0.5530 0.3343 0.3466 0.3267 0.2870 0.4987 0.5091

ACM f1 92.4371 48.5626 39.6693 39.2133 36.3879 33.1575 29.0617 19.8043 20.0053 19.7398

f2 0.5161±0.1274 0.5337±0.1631 0.5717±0.0959 0.5732±0.0974 0.6228±0.1063 0.6176±0.1521 0.6474±0.1330 0.6591±0.0766 0.7342±0.1490 0.6758±0.1846

0.4

(b)

0.8

f1 35.7263±0.5095 19.5515±0.5312 14.6406±0.4917 12.0156±0.4366 10.2899±0.1988 9.0319±0.2712 8.1508±0.2659 7.5523±0.3178 6.9948±0.1352 6.5065±0.1634

0.4

0.2

40 f1

NSJaDE f2 0.7548 0.6711 0.6425 0.6401 0.5770 0.5340 0.6969 0.7510 0.5716 0.7770

NSJaDE: δ = 0 NSGA-II: δ = 1 NSCDE: δ = 1 NSJaDE: δ = 1

0.6

0.4

35

ACM f1 74.8161 43.5417 36.7076 33.3747 28.2171 27.0773 23.7672 16.2829 15.2549 15.1200

0.6

f2

0.4

ABM f2 f1 f2 0.6590 43.7986 0.8400 0.8302 25.5939 0.8243 0.9905 18.8590 0.8135 1.0104 17.4694 0.9759 0.9599 17.2157 0.7947 0.8393 13.8326 0.8301 0.9162 13.8880 0.7809 0.8068 13.3623 0.7678 1.0164 11.0589 1.0407 1.0483 12.4255 1.0047 Minimizing f2 ADM ABM f1 f2 f1 f2 67.6944 0.4718 57.4223 0.6211 40.8429 0.4900 33.6755 0.5260 27.6654 0.4952 21.7675 0.4672 20.7836 0.5048 21.2178 0.4794 18.2905 0.4383 22.4871 0.4244 14.8751 0.454 20.1307 0.3848 14.1671 0.4555 16.3540 0.4301 17.2003 0.5048 18.1430 0.4702 17.2438 0.4614 17.0392 0.3728 16.4081 0.4769 15.1279 0.4632 f1 56.5904 27.7872 17.7841 14.7988 14.1449 12.1591 11.8183 11.1105 11.9458 11.0197

f2

0.6

0 30

ADM

f2

l = 5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

DBM

f2

DDM

5

f110

(i)

15

20

0 0

5

f1 10

15

20

(j)

Fig. 3. Comparisons of mean effective fronts (MEFs) obtained by NSGA-II, NSCDE, and NSJaDE (with or without uncertainty, i. e., δ = 0, 1) under different l. (a) MEFs under l = 5; (b) MEFs under l = 10; (c) MEFs under l = 15; (d) MEFs under l = 20; (e) MEFs under l = 25; (f) MEFs under l = 30; (g) MEFs under l = 35; (h) MEFs under l = 40; (i) MEFs under l = 45; (j) MEFs under l = 50. 1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

4.1

8

The graph a cat’s brain

The dataset of a cat’s brain network comes from [40]. The coupling matrix of a cat’s brain network is obtained from several subtle steps including cortical parcellation, thalamic parcellation, collation of connection data, and translation from database to connection matrix. The dataset contains the corticocortical and cortico-thalamical projections between areas of one brain hemisphere in a cat. The neuronal network can be divided into 53 cortical regions (N = 53) with 830 fibres of different densities. In addition, the network is composed of four topological clusters coinciding with four thalamocortical systems (functional cortical communities): visual cortex (16 areas), auditory (7 areas), somato-motor (16 areas), and fronto-limbic (14 areas). For more details regarding the properties of the neuronal network of a cat’s brain, one can refer to [25], [40], [41], [43] and the references therein. 4.2

Computer experimental information

The encoding scheme follows [19], [25], which is made up of two components with identical dimension size l: each dimension of the first part is a positive integer search space to represent the number of a controlling region; the other part is a continuous positive real number space to indicate the control gain of driver nodes. The dimension size is D = 2 ∗ l. It should be mentioned that the search spaces in the former and the latter part are (0.5, N + 0.5) and (0, N ), respectively. For the former part, we take the round operator to determine the locations of the driver nodes and the values of the latter part are directly used for the design of the control gains. Apparently, the allowed search range for each dimension is exactly N . The population size of NSGA-II, NSJaDE and NSCDE is SP = 100. The generation number is tmax = D ∗ γ ∗ α, where D is the dimension size, α is to control the number of fitness evaluation for each MOEA being nearly the same at each generation and γ is a predefined parameter. Each MOEA is carried out for 20 runs for eliminating discrepancy. 4.3

The effect of γ on search performance

The function of γ is to make a balance between complexity and search accuracy. As expected, a large γ is conductive to the enhancement of search accuracy and causes heavy computation complexity. A small γ is able to spare computation burdens and may bring about unsatisfactory search accuracy. For a detailed illustration of the parameter γ, we refer to Table 2, which clarifies that the results of γ = 150, 200, 250 do not show significant differences. In order to achieve a balance between search accuracy and computational complexity, we choose γ = 200 for all MOEAs, i. e., NSGA-II, NSJaDE and NSCDE. 4.4 Comparison of NSJaDE with statistical strategies in complex network theory In complex network theory, degree, between-centrality and closeness are important statistical information, which reflect different aspects of complex networks [21]. Here, in order to make a comparison with NSJaDE, driver nodes of the neuronal network are selected according to the statistical strategies by sorting degree, between-centrality (BC) and closeness in either ascending or descending way. It is worth

mentioning that in [10], [43], these kinds of statistical information are also important to characterize structural or functional aspects in brain networks for human or nonhuman. In [19], [25], these statistical information has been applied to the detection of controlling regions of the neuronal network of a cat’s brain. The summary of statistical strategies in the complex network theory is stated in Table 1. Due to the limitation of statistical strategies, control gains are simply supposed to be equal in each driver node and then tuned gradually from an interval [0, N ] with a step-size 0.01. Table 3 summarizes the comparisons of optimizing the mean effective objectives f1 and f2 by means of statistical methods of complex network theory and NSJaDE, when H = 5 and δ = 1. It should be noted that NSJaDE is a MOEA and has a set of nondominated solutions when achieving the termination condition. In order to make a comparison between the statistical methods and NSJaDE, the solution with the minimum f1 is selected among a set of nondominated solutions and record its corresponding f2 in each run. Then, after 20 runs, the mean values of the minimum values f1 and f2 at each run can be computed. For minimization of f2 , the same procedure can be carried out. It is also worth mentioning that NSJaDE can produce a set of solutions at one run by presenting robust Pareto fronts, which is an advantage over statistical methods in complex network theory. When minimizing f1 in (11), it is clear that NSJaDE has the best performance among the seven methods. All six statistical methods perform much worse than NSJaDE, which shows the satisfactory performance of NSJaDE. In addition, ADM, ABM and ACM usually work better than DDM, DCM and DBM, which means that a node in the neuronal network with a small degree plays an important role in controlling the entire network. This phenomenon shows that by selection of driver nodes with a small degree, the control information can then spread among the entire network without neglecting some less connected areas. As l increases, the mean effective function f1 gradually becomes smaller, which is similar to the case without uncertainty [19]. In the latter case of minimizing f2 in (11), NSJaDE is also the most powerful method among the seven methods for various l. Different from the case without uncertainties in selection of driver nodes and control gains [19], [25], NSJaDE cannot achieve zero for f2 as the calculation of the mean effective function f2 relies on the vicinity fitness and density of H. Nonetheless, NSJaDE performs better than the six statistical methods from complex network theory. To summarize, for minimization of f1 and f2 in (11), NSJaDE delivers better results than the statistical methods in robust multiobjective controllability of neuronal networks, as NSJaDE can efficiently exploit and explore robust solutions. 4.4.1 Comparison of NSGA-II, NSCDE and NSJaDE The comparison results of the three MOEAs: NSGA-II, NSCDE and NSJaDE for handling robust multiobjective controllability of neuronal networks are shown in Fig. 3 and Table 4, in which δ = 1 and H = 5 are considered. In addition, the results of the case without uncertainty (δ = 0) of NSJaDE is also provided in Fig. 3 and Table 4. The mean values and minimum values of robust nondominated solutions with minimum f1 (f2 ) and their corresponding f2 (f1 ) are given for NSGA-II, NSCDE and NSJaDE obtained from 20 runs. As expected, NSJaDE for the case without

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

9

0.8

δ=0 δ = 1, H = 5 δ = 1, H = 15

f2

0.4

0.4

0.4

0.2

0.2

0.2

0.8

δ=0 δ = 1, H = 5 δ = 1, H = 15

0.6

f2

0.6

f2

0.6

0.8

δ=0 δ = 1, H = 5 δ = 1, H = 15

δ=0 δ = 1, H = 5 δ = 1, H = 15

0.7 0.6 0.5

f2

0.8

0.4 0.3 0.2 0.1

0 30

f1

35

40

0 10

45

f1

15

(a)

20

0 5

25

f1

10

(b)

15

0 0

20

5

(c)

10

f1

15

(d)

Fig. 4. The comparison of MEFs by NSJaDE under different H. (a) MEFs obtained by NSJaDE with δ = 0, 1 and l = 5; (b) MEFs obtained by NSJaDE with δ = 0, 1 and l = 15; (c) MEFs obtained by NSJaDE with δ = 0, 1 and l = 25; (d) MEFs obtained by NSJaDE with δ = 0, 1 and l = 35.

1

0.8

δ=0 δ = 0.5 δ=1 δ=2

δ=0 δ = 0.5 δ=1 δ=2

0.8

0.8

0.6 0.5 0.4

0.4

0.3

0.4

0.3

0.2

f2

0.6

0.4

f2

0.6

0.5

0.2

0.2

0.2

0.1

0 15

20

25 f1 30

35

40

0 8

δ=0 δ = 0.5 δ=1 δ=2

0.7

f2

0.6

1

f2

0.8

δ=0 δ = 0.5 δ=1 δ=2

0.7

0.1

10

(a)

12 f1 14

16

0 5

18

f1

10

15

0 0

20

5

(c)

(b)

10

f1

15

(d)

Fig. 5. MEFs obtained by NSJaDE under δ = 0, 0.5, 1, 2. (a) MEFs under l = 10; (b) MEFs under l = 20; (c) MEFs under l = 30; (d) MEFs under l = 40.

0.8

δ = 0.5 δ¯1 = 0, δ¯2 = 1 δ¯1 = 1, δ¯2 = 0 δ=1

0.7

δ = 0.5 δ¯1 = 0, δ¯2 = 1 δ¯1 = 1, δ¯2 = 0 δ=1

0.7 0.6

0.8

δ = 0.5 δ¯1 = 0, δ¯2 = 1 δ¯1 = 1, δ¯2 = 0 δ=1

0.7 0.6

0.8

0.6

0.5

0.5

0.4

0.4

0.4

0.4

0.3

0.3

0.3

0.3

0.2

0.2

0.2

0.2

0.1

0.1

0.1

0.1

0 30

35

f1

40

45

f2

0.5

f2

f2

0.5

0 10

15

(a)

f1

20

25

0 5

(b)

δ = 0.5 δ¯1 = 0, δ¯2 = 1 δ¯1 = 1, δ¯2 = 0 δ=1

0.7

f2

0.6

0.8

10

f1

15

20

0 0

(c)

5

f1

10

15

(d)

Fig. 6. Comparisons of robust pareto fronts obtained by NSJaDE for illustrating the importance of driver nodes and control gains. (a) MEFs obtained with l = 5; (b) MEFs obtained with l = 15; (c) MEFs obtained with l = 25; (d) MEFs obtained with l = 35.

0.8

l = 5, δ = 0 l = 5, δ = 1 l = 10, δ = 0 l = 10, δ = 1

0.7

l = 15, δ = 0 l = 15, δ = 1 l = 20, δ = 0 l = 20, δ = 1

0.7

l = 25, δ = 0 l = 25, δ = 1 l = 30, δ = 0 l = 30, δ = 1

0.7 0.6

0.8

0.6

0.5

0.5

0.4

0.4

0.4

0.4

0.3

0.3

0.3

0.3

0.2

0.2

0.2

0.2

0.1

0.1

0.1

0.1

30f1

20

(a)

40

50

0 5

10

15 f1

(b)

20

25

0 5

δ=0 δ=1 δ=0 δ=1

f2

0.5

f2

0.5

0 10

l = 35, l = 35, l = 40, l = 40,

0.7

f2

0.6

0.8

f2

0.6

0.8

10

(c)

f1

15

20

0 0

5

f1

10

15

(d)

Fig. 7. Comparison of MEFs by NSJaDE under different l with δ = 0 and 1. (a) MEFs obtained by NSJaDE with l = 5, 10; (b) MEFs obtained by NSJaDE with l = 15, 20; (c) MEFs obtained by NSJaDE with l = 25, 30; (d) MEFs obtained by NSJaDE with l = 35, 40. 1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

10

TABLE 4 Comparisons of NSGA-II [29], [30], NSCDE [19] and NSJaDE (with out uncertainty and with uncertainty, i. e., δ = 0, 1) for robust multiobjective identification of controlling regions of the neuronal network of a cat brain when l = 5, 10, 15, 20, 25, 30, 35, 40, 45, and 50. NSGA-II minimization of f1

l=5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

l=5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

l=5 l = 10 l = 15 l = 20 l = 25 l = 30 l = 35 l = 40 l = 45 l = 50

min 36.3174 20.4134 14.8658 13.0702 10.5673 9.9763 8.9597 8.0467 7.9250 7.4251

f1 mean(std) 37.3361±0.6422 21.4044±0.4914 16.0194±0.6058 13.6908±0.5184 11.3782±0.4808 10.6048±0.5559 9.4198±0.4768 8.8390±0.5266 8.2876±0.2730 7.9144±0.3428

min 35.3123 18.6761 14.1437 11.8928 9.7754 9.0962 8.2649 7.5499 6.8917 6.3839

minimization of f1 f1 mean(std) min 35.7113±0.3557 0.4219 19.7261±0.7391 0.3736 14.6618±0.4385 0.2463 12.2824±0.2583 0.3029 10.4371±0.2945 0.4594 9.5423±0.3168 0.4254 8.6436±0.2813 0.5411 7.8690±0.2571 0.5239 7.2120±0.2494 0.4869 6.7463±0.2036 0.5054

min 35.2865 18.6793 13.7654 11.4364 10.0462 8.4736 7.6653 7.0858 6.7787 6.3091

minimization of f1 f1 mean(std) min 35.7263±0.5095 0.304 19.5515±0.5312 0.3361 14.6406±0.4917 0.3886 12.0156±0.4366 0.4352 10.2899±0.1988 0.4586 9.0319±0.2712 0.4498 8.1508±0.2659 0.4524 7.5523±0.3178 0.5713 6.9948±0.1352 0.4908 6.5065±0.1634 0.4197

min 0.3148 0.3936 0.3196 0.3368 0.4063 0.4581 0.4445 0.4621 0.4610 0.4439

f2 mean(std) min 0.4409±0.0921 38.1373 0.5389±0.0929 23.5122 0.5053±0.1525 18.2579 0.5303±0.0899 15.6893 0.5516±0.1253 14.1195 0.6032±0.1066 13.1264 0.5896±0.0997 12.5647 0.6301±0.1025 11.9614 0.6524±0.1378 13.7105 0.5741±0.0843 13.1944 NSCDE f2 mean(std) min 0.5019±0.0544 43.2933 0.6175±0.1606 24.0317 0.5372±0.1537 17.8599 0.5334±0.1297 17.5032 0.6867±0.0995 15.4534 0.6232±0.1272 13.2501 0.6486±0.0774 13.839 0.6730±0.1170 12.0002 0.7160±0.1399 11.874 0.6913±0.1221 10.9941 NSJaDE f2 mean(std) 0.5161±0.1274 0.5337±0.1631 0.5717±0.0959 0.5732±0.0974 0.6228±0.1063 0.6176±0.1521 0.6474±0.1330 0.6591±0.0766 0.7342±0.1490 0.6758±0.1846

uncertainty, i. e., δ = 0, works best in all cases. This phenomenon is intuitive, since one can accurately input right control gains in proper areas. By taking into account uncertainty under δ = 1 and H = 5, one can observe from Fig. 3 and Table 4 that NSJaDE performs best in most of the cases among the three MOEAs in dealing with robust multiobjective controllability when minimizing the mean effective function f1 and f2 under δ = 1. The mean effective fronts found by NSJaDE dominate those found by NSGA-II and NSCDE in most cases. Owing to the adaptive mechanism in JaDE, NSJaDE offers the best results, especially when l is large. NSCDE is better than the traditional NSGA-II, since composite differential evolution [47] is adopted. 4.4.2

Effect of H on robust multiobjective controllability

Due to the performance of NSJaDE in robust multiobjective controllability in the previous subsection, we utilize NSJaDE to carry out the following analysis. According to the method described in Fig. 2, it is natural to guess if more neighboring points are chosen for computing the mean effective objectives fi (i = 1, 2) in (11), the objective values will be closer to the true average values. However, the computation time will inevitably increase. Fig. 4 depicts the effect of H on robust multiobjective controllability of the neuronal

min 37.5511 22.8786 18.6849 16.1062 13.7504 10.8796 11.5995 11.5999 9.9770 10.6834

minimization of f2 f1 mean(std) min 44.5211±4.1592 0.1080 28.6677±3.6168 0.0879 21.4530±2.4363 0.03360 18.8167±2.3575 0.0639 16.6919±1.9868 0.0443 15.2309±1.6079 0.0272 15.5854±2.6728 0.0328 15.0002±1.6671 0.0324 15.4628±1.4467 0.0391 14.9346±1.3754 0.0306

f2 mean(std) 0.1829±0.0425 0.1303±0.0488 0.0896±0.0299 0.1129±0.0338 0.0903±0.0258 0.0786±0.0310 0.0769±0.0416 0.0840±0.0368 0.0780±0.0287 0.0737±0.0204

minimization of f2 f1 mean(std) min 51.9164±6.6717 0.1034 32.4252±6.6548 0.0254 23.7491±4.3319 0 19.6136±1.3881 0.0340 16.6246±1.1617 0.0363 17.0716±2.7555 0.0203 16.8841±2.2275 0.0428 15.6621±2.0074 0.0348 14.1512±2.1172 0.0461 14.6357±2.8069 0.0783

f2 mean(std) 0.1322±0.0202 0.0688±0.0310 0.0610±0.0330 0.0570±0.0185 0.0736±0.0264 0.0806±0.0348 0.1094±0.0321 0.0857±0.0274 0.0942±0.0317 0.1087±0.0188

minimization of f2 f1 mean(std) min 45.4990±4.9795 0.0986 28.5864±4.5944 0.0235 21.4522±2.1600 0.0049 18.3880±2.2422 0.013 17.5520±2.7333 0.0255 14.4369±2.8500 0.0345 13.9041±1.8878 0.0163 13.3921±1.3020 0.0615 12.1063±1.3223 0.0543 13.8334±2.4948 0.0231

f2 mean(std) 0.1657±0.0390 0.0715±0.0196 0.0552±0.0319 0.0357±0.0210 0.0519±0.0196 0.0694±0.0286 0.0801±0.0294 0.0958±0.0247 0.0913±0.0276 0.0894±0.0353

network. As shown from Fig. 4, in spite of having a small computational complexity, the mean effective front with a small H overestimates the true robust front. However, the mean effective front with H = 5 can also achieve a satisfactory approximation, compared with H = 15. In applications, one can choose H according to practical situations. 4.4.3 Effect of δ on robust multiobjective controllability In this subsection, we will investigate the impact of neighborhood size δ (or the level of uncertainty δ) on robust multiobjective controllability. The setting of the parameter H is H = 5. First, we consider δj = δ (j = 1, ..., D) for all the dimensions. The results are shown in Fig. 5 for δ = 0, 0.5, 1 and 2. It can be observed from Fig. 5 that as δ increases, the shift in the mean effective front moves away from the original Pareto front. In particular, when l is small, the mean effective front with δ = 0.5 is quite close to the mean effective front with δ = 0. The phenomenon of being away from the original Pareto front (δ = 0) as the increase of δ is natural, since the robust controllability of the neuronal network depends on driver nodes’ and their neighbors’ controllability. In order to further identify the importance of uncertainties in driver nodes and control gains, the following computer

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

11

TABLE 5 Comparisons of the times of selection as driver nodes under different δ by summarizing nondominated solutions under l = 5, 10, 15, 20, 25, 30, 35, 40, 45, 50. δ=0 Node Name AMLS Hipp VPc Sb ALLS AAF PS 21b 21a 3a PLLS 2 P SIV PSb 19 Tem 4 Enr 1 SII 36 VLS PFCMiI RS 17 DLS 20b AII 18 7 AI 3b SSAo 61 Ig PMLS PFCI SSAi AES 5Bm 6m 35 20a 4g Cga 5BI PFCMd Ia EPp 5Am 5AI CGp

δ = 0.5 βi 971 963 929 922 914 914 901 894 852 837 832 832 817 817 743 717 717 717 717 714 694 647 617 617 519 517 517 517 517 460 430 417 417 417 360 360 357 357 317 300 260 257 257 217 200 160 100 100 100 0 0 0 0

Node Name VPc Sb AAF PS 21b ALLS 2 Tem 21a Hipp AMLS Enr SII PLLS P 3a PSb DLS 18 19 1 36 SIV 17 3b 4 AI AES AII 4g 35 PMLS 20b VLS Ig RS 7 61 PFCMiI EPp 5BI SSAo 5Am 5Bm 20a SSAi Ia Cga 6m CGp 5AI PFCMd PFCI

δ = 1.5 βi 414 356 333 318 287 286 286 277 268 267 263 255 241 209 209 205 199 198 184 174 173 167 165 151 148 145 144 143 142 135 134 132 131 130 121 109 94 94 94 85 81 69 57 57 51 50 43 42 40 27 18 14 10

Node Name Sb AAF Hipp VPc P PLLS ALLS AII AMLS 21a 36 Tem Enr 2 18 1 PS 3a 17 Ig 19 DLS AES PSb SII SIV 4 21b EPp 4g 35 3b PMLS AI 7 VLS PFCMiI 20b Cga 61 20a SSAo Ia 5Bm 5AI SSAi PFCI 6m 5Am RS PFCMd 5BI CGp

Combined four cases δ = 0, 0.5, 1, 2

δ=2 βi 317 310 288 271 263 189 178 174 172 161 156 154 141 139 138 138 133 132 131 129 127 126 126 121 119 119 118 117 117 117 112 111 109 108 107 96 84 80 60 49 40 39 39 37 35 33 25 23 18 17 12 10 10

experiment is carried out. For the part of driver nodes in the encoding scheme, δj = δ¯1 , j = 1, ..., D , while for the part of 2 control gains in the encoding scheme, δj = δ¯2 , j = 1, ..., D . 2 The mean effective fronts are shown in Fig. 6. It can be seen that the average allowed variation in all dimensions are equal for the three cases δ = 0.5, δ¯1 = 0, δ¯2 = 1 and δ¯1 = 1, δ¯2 = 0. However, the distinctions among these three cases are evident. Fig. 6 shows that the mean effective front of δ¯1 = 0, δ¯2 = 1 dominates the other mean effective fronts when l = 15, 25, 35. For l = 5, the solutions obtained by δ¯1 = 0, δ¯2 = 1 are close to those of δ = 0.5. The phenomena explicitly show that the uncertainty in the selection of driver nodes heavily deteriorates the controllability of the neuronal network. Variations of the control gains mi do not have great impacts on controllability. The observation is consistent with the result in [33] and [25]. As shown in

Node Name Hipp P AAF Enr VPc AII PMLS EPp 17 1 PLLS DLS PSb 3b AMLS AI VLS SIV 36 19 SII Tem 18 Sb AES ALLS 7 4g 2 21a PS 3a 20b RS 20a 4 Cga 21b Ig CGp 35 Ia 61 5Am PFCMiI PFCI SSAo PFCMd SSAi 5AI 6m 5BI 5Bm

βi 590 365 265 236 218 203 186 186 172 170 165 165 164 160 153 151 149 148 143 136 136 125 124 124 119 118 118 115 107 98 88 88 86 86 78 72 71 70 65 64 62 60 51 42 31 30 29 26 25 24 21 15 12

Node Name Hipp VPc AAF Sb P AMLS ALLS PS PLLS 21a 21b 2 Enr Tem 3a SIV PSb 1 SII 19 36 4 AII DLS VLS 17 18 3b PFCMiI AI 20b PMLS 7 RS AES Ig 4g 35 61 SSAo SSAi PFCI EPp 20a 5Bm 6m Cga Ia 5BI PFCMd 5Am CGp 5AI

β 2108 1832 1822 1719 1654 1559 1496 1440 1395 1379 1368 1364 1349 1273 1262 1249 1227 1195 1190 1154 1113 1052 1036 1006 992 971 906 836 826 820 814 784 749 731 688 675 567 565 554 554 425 422 388 386 366 341 333 242 206 152 117 101 77

community Frontolimbic Auditory Auditory Frontolimbic Auditory Visual Visual Visual Visual Visual Visual Somato-motor Frontolimbic Auditory Somato-motor Somato-motor Frontolimbic Somato-motor Somato-motor Visual Frontolimbic Somato-motor Auditory Visual Visual Visual Visual Somato-motor Frontolimbic Auditory Visual Visual Visual Frontolimbic Visual Frontolimbic Somato-motor Frontolimbic Somato-motor Somato-motor Somato-motor Frontolimbic Auditory Visual Somato-motor Somato-motor Frontolimbic Frontolimbic Somato-motor Frontolimbic Somato-motor Frontolimbic Somato-motor

△k 2 4 3 8 3 7 4 7 5 5 4 7 -1 2 2 5 3 5 3 3 9 3 1 1 -2 1 2 1 -3 -1 0 2 -1 -2 -1 5 -1 7 0 -5 -5 -10 -6 -6 -6 -4 -13 -3 -10 -6 -8 -10 -10

[25], [33], control gains within a certain interval will deliver a similar controllability. The misplacement of the controller in a nearby area could degenerate the control performance. Hence, most of the current existing results focus on the selection of driver nodes instead of designing control gains. Fig. 7 plots the effect of the number of driver nodes l on controllability, together with δ. Intuitively, the increase of l will possess a better controllability. In addition, although more nodes are allowed to control, the solutions of l = 30, δ = 1 are dominated by l = 25, δ = 0. That is, the control performance under a certain level of uncertainty may not be compensated by a certain increase of the number of driver nodes. Hence, the uncertainties in selection of driver nodes and control gains have great impacts on robust multiobjective controllability.

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

12

4.5 Controlling regions identified by NSJaDE under different δ Here, we utilize NSJaDE to detect controlling regions under different levels of δ. Denote βi (i = 1, 2, ..., 53) as the sum of each area to serve as driver nodes in all the obtained nondominated solutions by NSJaDE when l = 5, 10, 15, 20, 25, 30, 35, 40, 45 and 50. Table 5 summarizes the respective times of each area/region/node to be selected as driver nodes for δ = 0, 0.5, 1 and δ = 2. In addition, Table 5 also presents the total summary of times for different δ and the communities of the nodes belonging to. The larger βi is, the more important of the area/region is to be as a driver node. From Table 5, it can be observed that the variations of the driver nodes are apparent. For the case without uncertainty, i. e., δ = 0, there are four areas/regions like EPp, 5Am, 5AI and CGp are never selected as driver nodes. When δ = 0, EPp, 5Am, 5AI and CGp become more important. The reason for the differences of driver nodes among different δ are that the driver nodes should be selected from the area having good controllability and its neighbor also having this property. For the last four columns of Table 5, it summarizes the total times of nodes selected as driver nodes. One can find that the results differ from the case without uncertainty. Although uncertainties are taken into account, the driver nodes are also spread uniformly in four communities, which is similar to the case without uncertainty as shown in [19]. In addition, the △k(i) = kin (i)−kout (i) of each node in the neuronal network is also an important index to measure the importance of nodes as driver nodes like the case without uncertainty δ = 0. The controlling areas are usually selected from the nodes with a large kin (i) and a small kout (i).

5

D ISCUSSION

AND CONCLUSION

This paper has tackled the robust multiobjective controllability of the neuronal network of a cat’s brain. By introducing uncertainties in the selection of driver nodes and the design of control gains, driver nodes are identified under different levels of uncertainties, in which a nondominated sorting adaptive differential evolution (NSJaDE) is proposed and utilized. The performance of the proposed NSJaDE has been manifested through a set of simulations by comparing with statistical methods from complex network theory and two multiobjecitve optimization evolutionary algorithms. The observations from computer experiments show that uncertainties have profound impacts on controllability of neuronal networks. The corresponding robust Pareto fronts are provided under various settings of parameters including the levels of uncertainties and the number of driver nodes. It is also found that the selection of driver nodes is more significant than designing controlling gains. Finally, some future research topics are provided here. Firstly, it is important to include the constraint like control cost [5], [34] to strengthen the study of robust multiobjective controllability. Secondly, we are also interested in making further extensions to leader selection in multi-agent systems [5], [48].

6

Industry (Northeastern University), IRTG 1740 (DFG), and the Alexander von Humboldt Foundation of Germany.

R EFERENCES [1] [2] [3]

[4]

[5]

[6]

[7]

[8] [9]

[10] [11] [12] [13] [14] [15]

[16] [17] [18]

[19]

[20]

ACKNOWLEDGEMENTS

This research is supported by the National Natural Science Foundation of China (61203235, 61333012) and the Key Laboratory of Integrated Automation for the Process

[21]

N. Markov, M. Ercsey-Ravasz, D. V. Essen, K. Knoblauch, Z. Toroczkai, and H. Kennedy, “Cortical high-density counterstream architectures,” Science, vol. 342, p. 1238406, 2013. E. Izhikevich, “Neural excitability, spiking and bursting,” Int. J. Bifur. & Chaos, vol. 10, no. 6, pp. 1171–1266, 2000. M. Ercsey-Ravasz, N. Markov, C. Lamy, D. V. Essen, K. Knoblauch, Z. Toroczkai, and H. Kennedy, “A predictive network model of cerebral cortical connectivity based on a distance rule,” Neuron, vol. 80, pp. 184–197, 2013. A. Goulas, M. Bastiani, G. Bezgin, H. Uylings, A. Roebroeck, and P. Stiers, “Comparative analysis of the macroscale structural connectivity in the macaque and human brain,” PLoS Comput Biol, vol. 10, p. e1003529, 2014. Y. Tang, F. Qian, H. Gao, and J. Kurths, “Synchronization in complex networks and its application-a survey of recent advances and challenges,” Annual Reviews in Control, vol. 38, no. 2, pp. 184–198, 2014. T. Wang, H. Gao, and J. Qiu, “A combined adaptive neural network and nonlinear model predictive control for multirate networked industrial process control,” IEEE Trans. Neural Networks and Learning Systems, accepted, 2015. Y. Tang, H. Gao, J. Lu, and J. Kurths, “Pinning distributed synchronization of stochastic dynamical networks: a mixed optimization approach,” IEEE Trans. Neural Networks and Learning Systems, vol. 25, no. 10, pp. 1804–1815, 2014. J. Lu, J. Zhong, Y. Tang, T. Huang, J. Cao, and J. Kurths, “Synchronization in output-coupled temporal boolean networks,” Scientific Reports, vol. 4, p. 6292, 2014. W. Zou, D. V. Senthilkumar, R. Nagao, I. Z. Kiss, Y. Tang, A. Koseska, J. Duan, and J. Kurths, “Restoration of rhythmicity in diffusively coupled dynamical networks,” Nature Comm., vol. 6, p. 7709, 2015. E. Bullmore and O. Sporns, “The economy of brain network organization,” Nature Review Neuroscince, vol. 13, pp. 336–349, 2012. V. Klyachko and C. Stevens, “Connectivity optimization and the positioning of cortical areas,” Proc Natl Acad Sci, vol. 100, pp. 7937–7941, 2003. E. Bullmore and O. Sporns, “Complex brain networks: graph theoretical analysis of structural and functional systems.” Nature Review Neuroscince, vol. 10, p. 1, 2009. D. Modha and R. Singh, “Network architecture of the longdistance pathways in the macaque brain,” Proc Natl Acad Sci, vol. 107, pp. 13 485–13 490, 2010. M. Young, “Objective analysis of the topological organization of the primate cortical visual system,” Nature, vol. 358, pp. 152–155, 1992. D. Zhang, L. Yu, Q. Wang, and C. J. Ong, “Estimator design for discrete-time switched neural networks with asynchronous switching and time-varying delay,” IEEE Trans. Neural Netw. Learning Syst, vol. 23, pp. 827–834, 2012. A. Engel, P. Fries, and W. Singer, “Rapid feature selective neuronal synchronization through correlated latency shifting,” Nat Rev Neurosc, vol. 2, pp. 704–716, 2001. P. Uhlhaas and W. Singer, “Neural synchrony in brain disorders: Relevance for cognitive dysfunctions and pathophysiology,” Neuron, vol. 52, pp. 155–168, 2006. N. Kopell and B. Ermentrout, “Chemical and electrical synapses perform complementary roles in the synchronization of interneuronal networks,” Proc. Natl Acad. Sci. USA, vol. 101, pp. 15 482–15 487, 2004. Y. Tang, H. Gao, and J. Kurths, “Multiobjective identification of controlling areas in neuronal networks,” IEEE/ACM Transactions On Computational Biology and Bioinformatics, vol. 10, no. 3, pp. 708–720, 2013. C. Zhou, L. Zemanov´a, G. Zamora, C. Hilgetag, and J. Kurths, “Hierarchical organization unveiled by functional connectivity in complex brain networks,” Phys. Rev. Lett., vol. 97, p. 238103, 2006. A. Arenas, A. Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Physics Reports, vol. 469, pp. 93–153, 2008.

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

13

[22] W. Zhang, Y. Tang, X. Wu, and J. Fang, “Synchronization of nonlinear dynamical networks with heterogeneous impulses,” IEEE Transactions on Circuits and Systems-I: Regular Papers, vol. 61, no. 4, pp. 1220–1228, 2014. [23] J. M. Palva, S. Monto, S. Kulashekhar, and S. Palva, “Neuronal synchrony reveals working memory networks and predicts individual memory capacity,” Proc. Natl Acad. Sci. USA, vol. 107, pp. 7580–7585, 2010. [24] Y. Liu, J. Slotine, and A. Barabasi, “Controllability of complex networks,” Nature, vol. 473, pp. 167–173, 2011. [25] Y. Tang, Z. Wang, H. Gao, S. Swift, and J. Kurths, “A constrained evolutionary computation method for detecting controlling regions of cortical networks,” IEEE/ACM Trans. On Computational Biology and Bioinformatics, vol. 9, pp. 1569–1581, 2012. [26] J. Gomez, G. Zamora, Y. Moreno, and A. Arenas, “From modular to centralized organization of synchronization in functional areas of the cat cerebral cortex,” PLoS ONE, vol. 5, p. e12313, 2010. [27] S. Cornelius, W. Kath, and A. Motter, “Realistic control of network dynamics,” Nature Communications, vol. 4, p. 1942, 2013. [28] K. Deb, Multi-Objective Optimization using Evolutionary Algorithms. Wiley., 2002. [29] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” IEEE Trans. On Evolutionary Computation, vol. 6, pp. 182–197, 2002. [30] K. Deb and H. Gupta, “Introducing robustness in multiobjective optimization,” Evolutionary Computation, vol. 14, no. 4, 2006. [31] K. Zhou, J. Doyle, and K. Glover, Robust and optimal control. Prentice Hall, USA, 1996. [32] I. Petersen and R. Tempo, “Robust control of uncertain systems: classical results and recent developments,” Automatica, vol. 50, no. 5, pp. 1315–1335, 2014. [33] F. Sorrentino, M. Bernardo, F. Garofalo, and G. Chen, “Controllability of complex networks via pinning,” Phys. Rev. E, vol. 75, p. 046103, 2007. [34] Y. Tang, Z. Wang, H. Gao, H. Qiao, and J. Kurths, “On controllability of neuronal networks with constraints on the average of control gains,” IEEE Trans. Cybernetics, vol. 44, no. 12, pp. 2670–2681, 2014. [35] Y. Tang, H. Gao, and J. Kurths, “Distributed robust synchronization of dynamical networks with stochastic coupling,” IEEE Transactions On Circuits and Systems-I: Regular Papers, vol. 61, no. 5, pp. 1508–1519, 2014. [36] S. Yin, X. Zhu, and O. Kaynak, “Improved PLS focused on key-performance-indicator-related fault diagnosis,” IEEE Trans. Ind. Electron., vol. 62, pp. 1651–1658, 2015. [37] P. Shi, Y. Zhang, and R. Agarwal, “Stochastic finite-time state estimation for discrete time-delay neural networks with markovian jumps,” Neurocomputing, vol. 151, pp. 168–74, 2015. [38] Y. Li, J. Li, and M. Hua, “New results of h∞ filtering for neural network with time-varying delay,” Int. J. Innovative Computing, Information and Control, vol. 10, pp. 2309–2323, 2014. [39] C. Xu, M. Liao, and Q. Zhang, “On the mean square exponential stability for a stochastic fuzzy cellular neural network with distributed delays and time-varying delays,” Int. J. Innovative Computing, Information and Control, vol. 11, pp. 247–256, 2015. [40] J. Scannell, G. Burns, C. Hilgetag, M. O’Neill, and M. Young, “The connectional organization of the cortico-thalamic system of the cat,” Cereb Cortex, vol. 9, pp. 277–299, 1999. [41] C. Hilgetag, G. Burns, M. O’Neill, J. Scannell, and M. Young, “Anatomical connectivity defines the organization of clusters of cortical areas in the macaque monkey and the cat,” Phil Trans R Soc London B, vol. 355, pp. 91–110, 2000. [42] G. Zamora-Lopez, C. S. Zhou, and J. Kurths, “Exploring brain function from anatomical connectivity,” Front Neurosci, vol. 5, p. 83, 2011. [43] O. Sporns, C. Honey, and R. Kotter, “Identification and classification of hubs in brain networks,” PLoS ONE, vol. 10, p. e1049, 2007. [44] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, “Spontaneous synchrony in power-grid networks,” Nature Physics, vol. 9, pp. 191–197, 2013.

[45] L. Ghaoui, F. Oustry, and H. Lebret, “Robust solutions to uncertain semidefinite programs,” SIAM Journal on Optimization, vol. 9, p. 1, 1998. [46] J. Q. Zhang and A. C. Sanderson, “Jade: adaptive differential evolution with optional external archive,” IEEE Trans. on Evolutionary Computation, vol. 13, pp. 945–958, 2009. [47] Y. Wang, Z. Cai, and Q. Zhang, “Differential evolution with composite trial vector generation strategies and control parameters,” IEEE Trans. On Evolutionary Computation, vol. 15, pp. 55–66, 2011. [48] Y. Tang, H. Gao, W. Zhang, and J. Kurths, “Leader-following consensus of a class of stochastic delayed multi-agent systems with partial mixed impulses,” Automatica, vol. 53, pp. 346–354, 2015.

Yang Tang (M’11) received the B.S. and the Ph.D. degrees in electrical engineering from Donghua University, Shanghai, China, in 2006 and 2010, respectively. He was a Research Associate with The Hong Kong Polytechnic University, Kowloon, Hong Kong, from 2008 to 2010. From 2011-2015, he conducted his postdoctoral research with the Humboldt University of Berlin, Berlin, Germany and the Potsdam Institute for Climate Impact Research, Potsdam, Germany. Since 2015, he has been a Professor at East China University of Science and Technology, Shanghai, China. He has published more than 50 refereed papers in international journals. His current research interests include synchronization/consensus, networked control systems, evolutionary computation, and bioinformatics and their applications. Dr. Tang was the recipient of the Alexander von Humboldt Fellowship in 2011. He is an Associate Editor of the Journal of the Franklin Institute, Neurocomputing and a Leading Guest Editor of the Journal of the Franklin Institute.

Huijun Gao (SM’09-F’14) received the Ph.D. degree in Control Science and Engineering from the Harbin Institute of Technology, Harbin, China, in 2005. From 2005 to 2007, he conducted his postdoctoral research with the Department of Electrical and Computer Engineering, University of Alberta, Edmonton, AB, Canada. Since November 2004, he has been with the Harbin Institute of Technology, where he is currently a Professor and the Director of the Research Institute of Intelligent Control and Systems. His current research interests include network-based control, robust control/filter theory, and timedelay systems and their engineering applications. He is Co-EiC of the IEEE Transactions on Industrial Electronics, an Associate Editor of Automatica, the IEEE Transactions on Cybernetics, the IEEE Transactions on Fuzzy Systems, the IEEE/ASME Transactions on Mechatronics, and the IEEE Transactions on Control Systems Technology. He is an Administrative Committee Member of the IEEE Industrial Electronics Society.

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2015.2485226, IEEE/ACM Transactions on Computational Biology and Bioinformatics FINAL VERSION

14

Wei Du received the B.S. and M.S. degrees in electrical engineering from Donghua University, Shanghai, China, in 2009 and 2012, respectively. He is currently pursuing the Ph.D. degree with The Hong Kong Polytechnic University, Hong Kong. His current research interests include evolutionary computation, especially differential evolution, multiobjective evolutionary algorithm, and their applications.

Jurgen ¨ Kurths studied mathematics at the University of Rostock, Rostock, Germany, and received the Ph.D. degree from the GDR Academy of Sciences, Berlin, Germany, in 1983. He was a Full Professor at the University of Potsdam, Potsdam, Germany, from 1994 to 2008. Since 2008, he has been a Professor of nonlinear dynamics at the Humboldt University, Berlin, and the Chair of the research domain Transdisciplinary Concepts of the Potsdam Institute for Climate Impact Research, Potsdam. Since 2009, he has been the 6th Century Chair of Aberdeen University, Aberdeen, U.K. He has published more than 650 papers (H-factor: 68). His current research interests include synchronization, complex networks, and time series analysis and their applications. Dr. Kurths is a fellow of the American Physical Society. He received the Alexander von Humboldt Research Award from CSIR, India, in 2005, and the Honorary Doctorate from the Lobachevsky University Nizhny Novgorod, in 2008, and from the State University Saratov, in 2012. He became a member of the Academia Europaea in 2010 and of the Macedonian Academy of Sciences and Arts in 2012. He is an Editor of PLoS ONE, the Journal of Nonlinear Science, CHAOS, Europ. Physics J. Special Topics etc.

Jianquan Lu (M’13) received the B.S. degree in mathematics from Zhejiang Normal University, Zhejiang, China, in 2003, the M.S. degree in mathematics from Southeast University, Nanjing, China, in 2006, and the Ph.D. degree from City University of Hong Kong, Hong Kong, in 2009. He is currently a professor at the Department of Mathematics, Southeast University, Nanjing, China. His current research interests include collective behavior in complex dynamical networks and multi-agent systems, and Boolean networks. He has published over 30 papers in refereed international journals. Dr. Lu is an associate editor of Neural Processing Letters (Springer) and a guest editor of Mathematics and Computers in Simulation (Elsevier). He is the recipient of an Alexander von Humboldt Fellowship in 2010, Program for New Century Excellent Talents in University by The Ministry of Education, China in 2010, and The First Award of Jiangsu Provincial Progress in Science and Technology in 2010 as the Second Project Member.

Athanasios V. Vasilakos (SM’11) is recently Professor with the Lulea University of Technology, Sweden. He served or is serving as an Editor for many technical journals, such as the IEEE TRANSACTIONS ON NETWORK AND SERVICE MANAGEMENT; IEEE TRANSACTIONS ON CLOUD COMPUTING, IEEE TRANSACTIONS ON INFORMATION FORENSICS AND SECURITY, IEEE TRANSACTIONS ON CYBERNETICS; IEEE TRANSACTIONS ON NANOBIOSCIENCE; IEEE TRANSACTIONS ON INFORMATION TECHNOLOGY IN BIOMEDICINE; ACM Transactions on Autonomous and Adaptive Systems; IEEE JOURNAL ON SELECTED AREAS IN COMMUNICATIONS. He is also General Chair of the European Alliances for Innovation (www.eai.eu).

1545-5963 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

Robust Multiobjective Controllability of Complex Neuronal Networks.

This paper addresses robust multiobjective identification of driver nodes in the neuronal network of a cat's brain, in which uncertainties in determin...
1KB Sizes 1 Downloads 6 Views