Feature

XIAOFAN LIU

Synchronization and Control of Complex Networks via Contraction, Adaptation and Evolution Pietro DeLellis, Mario di Bernardo, Thomas E. Gorochowski, and Giovanni Russo

Digital Object Identifier 10.1109/MCAS.2010.937884

64

IEEE CIRCUITS AND SYSTEMS MAGAZINE

1. Introduction omplex networked systems abound in Nature and Technology. They consist of a multitude of interacting agents communicating with each other over a web of complex interconnections. Flocks of birds, platoon of cooperating robots, swirling fishes in the Ocean are all examples whose intricate dynamics can be modeled in terms of three essential ingredients: (i) a mathematical description of the dynamical behavior of each of the agents in the network; (ii) an interaction (or coupling) protocol used by agents to communicate with each other and (iii) a graph describing the network of interconnections between neighboring agents. These three elements are actually mapped onto the mathematical model usually considered in the literature to describe a complex network which uses appropriate equations to describe the node dynamics, the coupling protocol and the

C

1531-636X/10/$26.00©2010 IEEE

THIRD QUARTER 2010

network topology (see Sec. 2.2; e.g., [16] and references therein for further details). The most striking feature of complex networked systems observed in Nature and Technology is their ability to show emergent behavior that cannot be explained in terms of the individual dynamics of each single agent in the network. Synchronization and coordination of motion are possibly the simplest and most frequent of these types of phenomena. Think for example of the incredible patterns formed by flocks of birds or schools of fishes [24], the synchronous flashing of fireflies in Amazonia [19] and the many other situations where the emergence of coordinated behavior has been observed. Over the past few years, much research effort has been focused on understanding and explaining the onset of synchronization in complex networks, with the aim of uncovering the very mechanism causing its occurrence. Also, the dependence has been studied of synchronization on specific features of the network such as the form of the coupling strategy or the topology of the network. It has been found that synchronization emerges even in the presence of simple, diffusive coupling where neighboring nodes adjust their dynamics proportionally to the mismatch between some output function of their states [6, 14, 15, 22, 23, 53, 55, 71, 74, 80, 108]. Also, it was observed that synchronization depends on the network topology so that it is simpler to synchronize networks with certain topological features. Synchronizability of small world, scale-free, disassortative or assortative networks was studied showing that the topology has indeed an influence on the dynamics of the network (see for example [9, 79, 105, 106, 110]). In all cases, when synchronization occurs all nodes tend towards the same asymptotic evolution which is unknown a priori and emerges from the negotiation process taking place on the network between neighboring nodes. A different problem is when the aim is to control the network onto some desired asymptotic solution. In this case, a possible solution would be to control every node in the network so as to steer the dynamics of each of them towards the desired evolution. In practice, this is not viable as the number of nodes in the network is typically very large and the amount of control effort available bounded. Pinning control was introduced as a viable alternative, where only a (small) fraction of the network nodes are directly controlled by means of some external input, with the control action being propagated to the rest of the nodes through their interconnections [21, 45, 46, 56, 54, 77, 78, 92, 96, 107, 109]. Controllabil-

ity of networks was also defined in [96] and further investigated in [77]. One of the key features of most of the techniques for the control and synchronization of networks presented so far in the literature is the time-invariant nature of the coupling strength between nodes and the topology of their interconnections. Specifically, it is often supposed that the coupling gain determining the strength of the interaction between nodes is constant as is the adjacency matrix describing the topology of the interconnections between them. This is not in agreement with much evidence from the natural world, showing that instead in such examples as flocks of birds, swarms of fireflies or schools of fishes, individuals in the network are often able to form or suppress interconnections between themselves and adjust (or adapt) the strength and structure of their coupling. Also, often, rather than being simply diffusive, the coupling functions between agents in the network are typically nonlinear with a functional form that has been engineered by Nature to best fulfill a certain need. Over the past few years, strategies to cope with the case of time varying topologies have been proposed but often with the aim of showing that synchronization can be achieved even in the presence of variations in the network topology, e.g. [11, 12]. Adaptation and evolution in Nature are instead desired phenomena, engineered over the years in order for complex networks of dynamical agents to adapt the strength of their interactions and evolve an appropriate network of interconnections so as to perform some function of interest. The aim of our work is to attempt at mimicking these features of natural networks and propose strategies for the adaptation and evolution of complex networks of dynamical systems that guarantee the emergence of some asymptotic behavior of interest, namely a synchronous solution. Such strategies will be, in general, local and decentralized so that it is up to pairs or groups of agents to self-determine the strength and topology of their interconnections. We start with the generalization of the concept of dynamic graph first stated by Šiljak in [91] by incorporating some of the features of a complex adaptive system as defined by Holland in [47]. Then, we propose two alternative strategies to design the coupling functions between nodes. One is based on finding appropriate adaptation laws to evolve the coupling gain between neighboring nodes, and the other on the ad hoc design of a nonlinear coupling protocol based on contraction theory [62]. We will discuss how in both cases it is possible to prove convergence of all network nodes onto a

Pietro DeLellis, Mario di Bernardo, and Giovanni Russo are with the Department of Systems and Computer Engineering, University of Naples Federico II, Italy. Mario di Bernardo and Thomas E. Gorochowski are with the Bristol Centre for Complexity Sciences, University of Bristol, U.K.

THIRD QUARTER 2010

IEEE CIRCUITS AND SYSTEMS MAGAZINE

65

common solution and briefly present a set of representative examples. We then move to the case where the topology of the interconnections between nodes in the network evolves. Firstly, we will consider a novel strategy, termed as edgesnapping in [29]. Specifically, we will make neighboring nodes able to decide whether to activate or suppress a link between themselves according to the strength of their mismatch. We will show that in this manner networks with a certain final topology emerge with interesting features that will be discussed in Sec. 4. Next, we take a different approach and consider the case where evolution of the network must be such as to optimize some performance index defined on the network and its dynamics. A novel computational tool, NetEvo, will be presented and used to evolve networks in order to maximize their synchronizability, i.e. the range of values of the coupling gain guaranteeing synchronization. We will compare the topologies of the networks evolved using the Laplacian eigenratio as a static measure of the network synchronizability (e.g. [9, 75]) with those obtained by considering dynamical synchronizability measures. We will observe that in both cases networks can present some striking features in terms of their topologies and focus on the presence of certain characteristics such as the frequencies of motifs [51, 73] and other macroscopic observables. The possibility of having topological bifurcations will also be discussed. Finally, we discuss some of those that we believe are the most pressing open problems and challenges in the area of complex dynamical networks suggesting that much work is still needed to achieve the ultimate goal of engineering fully adaptive, self-evolving networks able to show the incredible features, we are able to marvel about in all the natural complex systems surrounding us. 2. Background 2.1. Describing a Complex Network Complex networks provide a flexible tool for modeling many real-world systems and consist of three main attributes: topology, dynamics and an evolutionary process. In order to fully define a complex network and study the influence each of these attributes has, a common mathematical framework is required. This is to not only to provide a means to analyzing various properties a particular system possesses, but also to act as a shared language in which to describe the systems themselves. Complex networks are a purposefully general concept and so being able to relate problems from different fields provides an important first step to understanding underlying principles they may share. Several attempts have 66

IEEE CIRCUITS AND SYSTEMS MAGAZINE

been made with this in mind, including coupled cell networks, dynamic graphs and complex adaptive systems. Coupled cell networks were introduced by Golubitsky and Stewart [41] as a means to investigating how network topology and system-level dynamics are linked. These networks consist of a number of cells connected by a fixed set of edges. Equivalence relations are defined over each of these sets, assigning a type to every element. These relations can represent the form of dynamical system at each cell and the type of coupling across an edge, although the exact forms are not necessarily required. By studying the structure of these systems with particular dynamics, it has been shown that certain symmetries within the network topology and equivalence relations can predict robust synchronous states and other related phenomena. Although this formalism can accommodate any form of fixed topology and dynamics, it unfortunately does not facilitate evolution of the structure itself. In contrast, the dynamic graphs of Šiljak [100] focus on how the structure of a network evolves over time. A weighted graph of a fixed number of nodes is placed within a linear vector space, and a dynamic graph then defined as a one-parameter group of transformations of this space into itself. Put more simply, a dynamic graph can be viewed as one where edge weights vary in time due to some differential equation. A major benefit of this formalism is that much of the existing dynamical systems theory can be directly applied using adjusted forms of stability. This approach allows for the evolution of edge structure to be incorporated into a model, but does have limitations. Edge states can only take a single state value (their weight) and other aspects of evolution, such as growth, are not possible. To allow for systems that structurally evolve and grow we turn to the work of Holland and his concept of a complex adaptive system [47]. This framework was developed during the 1970’s in an attempt to understand the process of adaptation. It comprises of four main elements: 1. A – A set of possible structures; 2. V – A set of operators for modifying structures; 3. I – A set of inputs from the environment; 4. t : I 3 A S V – The adaptive plan. Evolution takes place in two stages, first the current structure and input from the environment are mapped to an operator v [ V using the adaptive plan t. Then, the chosen operator is applied to return the evolved structure. It should be noted that both the adaptive plan and operators for modifying structures can be probabilistic, leading to a stochastic evolution of the structure. The main advantage of this framework is its flexibility, allowing for the set of admissible structures, operators THIRD QUARTER 2010

and adaptive plan to be defined any way we choose. This provides endless possibilities, but unfortunately gives no indication as to how network topology and dynamics can be best incorporated. As we have seen, several attempts have been made to develop a framework for the study of complex networks. In each case, trade-offs were made that resulted in none fully addressing the three attributes discussed earlier. With this in mind and taking inspiration from these previous ideas, we now introduce the concept of an evolving dynamical network (EDN) which aims to meet this goal. We begin by first introducing the concept of a generalized dynamic graph. This has the role of providing a structure in which network topology and dynamics can be fully defined. Taking a directed graph G 5 1 V, E 2 of a fixed number of nodes and edges, we associate a state to every node and edge 1 V, E 2 . This state can be of any form and dimension required. Furthermore, these states embody the dynamics of the system, with their motion described by a mapping F which is also defined over some set of times T. Given the current set of states and any external inputs U, this mapping returns the new state of the system. It is worth noting that we consider network dynamics to cover any aspects related to changes in system states. To incorporate evolution of the structure and dynamics, i.e. G and F, we consider the set of all possible generalized dynamic graphs D. We see evolution of a system as exploration of this space, taking place either online, as the network dynamics occur, or off-line, having the network state reset after each evolutionary step. To embody this feature we view an evolving dynamical network as a complex adaptive system where the set of admissible structures is D. We consider a set of structural operators V that can alter a given structure, a set of inputs I that can be from the environment or for control purposes, and the evolutionary plan t, which given an input and structure decides the structural operator to apply. With these components we are now in the position to fully describe the three main attribute of a complex network. Formal definitions and an illustration of the connection between the two aspects of an EDN, network dynamics and evolution, can be found in Box 1. In following sections EDNs will hold a central role, helping us relate similarities and differences between various approaches to synchronization.

collective behaviors, like synchronization and consensus, in networks of interconnected dynamical systems. This class of network, often named complex networks, is an instance of a type of EDN described in Box 1. In fact, the classic complex network model can be viewed as a generalized dynamic graph in which: ■ V 5 5 x1, x2, c, xn 6 , where xi [ Rn, 4i : vi [ V. The network nodes are dynamical systems, and therefore the nodes assignments coincide with the states xi of the nodes belonging to the real coordinate space Rn. ■ E 5 5 eij 6 , 4 ij : 1 vi, vj 2 [ E, with eij 5 s [ R. The constant s represents the overall coupling strength between the network nodes. ■ U 5 [. There are no external inputs. ■ F 5 V 3 E 3 T S V 3 E can be split into two operators Fv and Fe, where Fv 5 V 3 E 3 T S V is a continuous differential operator that describes the nodes’ evolution, while Fe 5 E S E is the identity operator describing the constant edge states. In terms of ODEs, the evolution of the nodes’ states is described by the following set of differential equations

2.2. Problem Statements: Synchronization and Control As explained in Section 1, scientists from different communities, including physicists, biologists, control engineers and mathematicians, are making an effort to understand the dynamics leading to the emergence of

where degi is the cardinality of the set Ni, that is, the N set of incident nodes in i. Notice that ,ii 5 2 a j 51 ,ij, 4i 5 1, c, N. Throughout the paper, we typically consider undirected networks, that is, networks in which if 1 vi, vj 2 [ E, then also 1 vj, vi 2 [ E. Therefore, L is symmetric and zero-column sum.

THIRD QUARTER 2010

# xi 5 f 1 xi, t 2 1 ci 1 x, t 2 ,

(1)

with i 5 1, c, N. In (1) the function f : Rn 3 R 1 S Rn is the intrinsic dynamics of the i-th node, while the function ci : RnN 3 R 1 S Rn, termed as interaction or coupling function, describes the interaction of the i-th node with the other nodes composing the interconnected system. Typically, ci depends only on the dynamics of those nodes directly linked to node i selected by the Laplacian matrix describing the network topology, see e.g. [49]. Specifically, often the coupling protocol is chosen in (1) as ci 1 x, t 2 5 2s a ,ij h 1 xj 2 , N

i [ V.

(2)

j51

Here, h : Rn S Rn is the output function through which the network nodes are coupled, ,ij is the ij th element of the zero-row sum Laplacian matrix L associated to the graph G 5 5 V, E 6 describing the network topology, defined as: degi ,ij 5 • 21 0

i5j 1 vi, vj 2 [ E, otherwise

(3)

IEEE CIRCUITS AND SYSTEMS MAGAZINE

67

BOX 1. EVOLVING DYNAMICAL NETWORKS

I

llustration of how topology, dynamics and evolution are Network Dynamics

integrated within an evolving dynamical network (EDN). It

can be seen that there are two main dimensions, with the op-

a vital step towards more general theories spanning multiple application domains. Definition 4 (Generalized Dynamic Graph) We define a generalized dynamic graph for a fixed number of nodes N, network structure E, and opera-

tor F, as the collection D 5 1 V, E, V, E, U, T, F 2 [ D

Network Evolution

Having a framework like this for studying complex networks is

Φ

D1

erators F and v [ V determining which direction to take.

ω1

Update Structure, Rewire an Edge

Φ

D2

Φ

Φ

Grow Structure, Add Node, and Connecting Edges

ω2 D3

Φ

Φ

where,

• V 5 5 v1, v2, c, vN 6 — finite set of nodes.

• E 5 5 e1, e2, c, eM 6 where ei [ V 3 V — finite

Definition 5 (Evolving Dynamical Network) Using generalized dynamic graphs as an underlying

set of edges.

• V 5 5 v1, v2, c, vN 6 where vi [ Vi — set of

structure, we define an evolving dynamical network (EDN) as the collection E^ 5 1 D, V, I, t 2 [ E^ where, • D 5 5 D1, D2, c 6 — set of admissible gener-

node state assignments.

• E 5 5 e1, e2, c, eM 6 where ei [ Ei — set of

alized dynamic graphs.

• V 5 5 v 1, v 2, c 6 — set of structural operators such that v [ V is a mapping v : D S P

edge state assignments.

• U 5 5 u1, u2, c, uP 6 where ui [ Ui — set of external input assignments.

where P is a probability distribution over D.

• T – set of times relating to dynamics.

• I — inputs from the environment or for

• F : V 3 E 3 U 3 T S V 3 E — mapping from state to state for nodes, edges and control input where the mapping F can take any form,

determines the operator v [ V to be used

e.g. deterministic, probabilistic, etc.

when transitioning to a new structure.

Now we are ready to give a formal definition of the concept of synchronization in a complex network. Definition 1 Network (1) is said to achieve asymptotic synchronization if and only if lim 1 xi 1 t 2 2xj 1 t 2 2 5 0,

t S 1`

i, j : 1 vi, vj 2 [ E.

(4)

Asymptotic synchronization is the most intuitive case of synchronization. However, other forms of synchronization are also possible, e.g. lag synchronization as defined in [17, 81]. The emergence of the synchronous behavior depends on both the individual nodes’ dynamics and the structure of the interconnections. Different analytical tools have been used to give conditions for synchronization, such as Lyapunov-based techniques, the Master Stability Func68

control purposes. • t : I 3 D S V — evolutionary plan which

IEEE CIRCUITS AND SYSTEMS MAGAZINE

tion approach and contraction theory (see Box 2 for further details). In many real networks, synchronization is not the only goal that one may want to achieve. In formation control [7, 10, 25, 36, 66], for instance, the aim is not only to synchronize the motion of the set of nodes (the agents) of the network, but we want to tame the dynamics of the whole network to a desired reference trajectory. To solve this kind of problem, it is necessary to give external inputs ui 1 x 2 [ Rm, m # n to the network nodes. In the EDN framework, the only difference is that the set U is not anymore empty and the domain of the operator F is V 3 E 3 U 3 T. A possible solution to this control problem is to add a controller to each of the network nodes. This approach is not feasible when the aim is to control large networks. This is why in the recent literature the so-called pinning control scheme for controlling complex networks has been proposed [21, THIRD QUARTER 2010

45, 54, 77, 107, 109]. In this scheme, an extra-node vs is added in the network, the pinner. This extra-node is identical to the other nodes of the network and defines the desired trajectory by being only unidirectionally coupled to the rest of the network. Only a small subset of nodes, say Npin V N, are directly controlled by the pinner. These nodes are called pinned nodes. The equation describing the closed-loop network are the following: # xi 5f 1 xi, t 2 2 s a ,ijh 1 xj 2 2 diqi 1 xi 2 xs 2 , N

# xs 5 f 1 xs, t 2 di 5 e

i 5 1, c, N (5)

i51

(6)

1, for i 5 1, c, Npin, 0, for i 5 1 Npin 1 1 2 , c, N,

(7)

Definition 2 Network (5)–(7) is said to achieve controlled network synchronization if and only if lim 1 xi 1 t 2 2 xs 1 t 2 2 5 0,

i [ V.

(8)

To study the local controllability of the network nodes to the desired trajectory, the concept of pinning controllability was introduced in [96], applying the MSF approach (see box 2) to an extended network, which includes the external pinner. With this approach, theoretical bounds on the control gains qi and on the coupling gain s were given. In the following sections, after highlighting the limitations of the previous models (2) and (5)–(7), we will propose novel strategies for designing the nodes’ coupling (Section 3) and for evolving the network topology (Section 4). 3. Designing Novel Coupling Strategies 3.1. Synchronization via Decentralized Gain Adaptation Many real phenomena are characterized by the presence of adaptive mechanisms; for example, wireless networks of sensors that gather and communicate data to a central base station [68], or networks of robots when environmental conditions change unexpectedly (i.e. a robot loses a sensor) [98]. Moreover, examples of adaptive networks can be found in biology and the natural world [24, 38, 50]. In all these cases, it is realistic to assume that the strength of the interactions among nodes, characterized mathematically by s, is not identical for every node and timeinvariant. Real-world networks are often characterized instead by evolving, adapting couplings which vary in time THIRD QUARTER 2010

N # xi 5 f 1 xi, t 2 2 a lij 1 t 2 h 1 xj 2 ,

i 5 1, c, N

(9)

1 vi, vj 2 [ E,

(10)

j51

where xs denotes the desired trajectory defined by the pinner, ui 5 2 diqi 1 xi 2 xs 2 is the control input, with qi [ R being the i th control gain.

t S 1`

according to different environmental conditions. Adaptive gain strategies have been recently presented in a number of papers, see for example, [20, 27, 28, 31, 57, 97, 103, 112, 113]. In these strategies, the network gains evolve according to integro-differential adaptation laws driven by global or local synchronization errors. For instance, a common adaptive gain for all the network nodes is studied in [20, 57] while decentralized gain adaptation is proposed, using different gains for each node in [113] or edge in [28]. Specifically, the coupling protocol is chosen as a timevarying function so that the network model becomes:

# sij 5 g 1 t, sij, xi, xj 2

where sij [ R is the time-varying coupling gain associated to edge 1 vi, vj 2 , and lij 1 t 2 is ij th element of the timevarying laplacian matrix L 1 t 2 , defined as 2sij 1 t 2 lij 1 t 2 5 µ a sik k:vk [Ni

0

1 vi, vj 2 [ E i5j otherwise,

(11)

where Ni ( V is the set of neighbors of node vi, that is, the set of nodes incident to vi. In designing the adaptive approach, a key point is the appropriate selection of the function g. The design of the adaptive law is clearly dependent on the nature of the individual system, described by f, and on the output function h. In what follows, we assume that the output function is h 1 z 2 5 Gz, where G is a non-null diagonal positive semidefinite matrix, commonly named as the inner coupling matrix. This matrix defines the subset of state variables through which the network nodes are coupled. As for the node dynamics, we give the following definition: Definition 3 We say that the vector field f 1 z, t 2 or, equivalently, the dynam# ical system z 5 f 1 z, t 2 is QUAD 1 D, v 2 if and only if 1 x2y 2 T 1 f 1 x, t 2 2 f 1 y, t 22 2 1 x2y 2 T D 1 x2y 2

# 2v 1 x2y 2 T 1 x2y 2 , 4x, y [ Rn,

(12)

where D is an arbitrary diagonal matrix and v is a positive scalar. This condition gives an upper bound on the rate of variation of the vector field f. This is an assumption which is often assumed to prove asymptotic synchronization in complex networks, see for instance [13, 27, 32, 60, 61, 64]. Notice that the assumption is quite mild and is satisfied by all linear systems and by many chaotic systems, as for instance the Chua’s circuit and the Lorenz oscillator IEEE CIRCUITS AND SYSTEMS MAGAZINE

69

[28, 111]. The relationship between QUAD assumption and other assumptions usually made to prove synchronization, such as the well-known Lipschitz condition and the contraction theory, are expounded in [30]. Under the assumption that the QUAD condition is satisfied, analytical conditions guaranteeing asymptotic synchronization and convergence of the coupling gains were derived using a Lyapunov-based approach. (Here, we give some of the more significant results, referring the reader to [28, 111] for the proofs.) Let us assume that g depends only on the difference between the states of the endpoints of the corresponding edges, that is g 1 t, sij, xi, xj 2 5 k 1 0 xj 2 xi 0 2 ,

(13)

where k is a function belonging to one of the following classes: 1. Class 1. k 1 z 2 5 azp, with 0 , p # 2, a . 0. 2. Class 2. k is a bounded k function1. Theorem 1 Assume that f is a continuously differentiable QUAD vector field, with D 2 vI # 0, and G 5 I. If the scalar function k belongs to class 1 or 2, then network (9), (10), (13) reaches asymptotic synchronization and each coupling gain sij converges to some finite value. It is worth emphasizing that, under the assumptions of the theorem, the network synchronizes for any connected topology. This is a quite strong result, but is limited to the case in which the network nodes are connected through all their states and the vector field f is required to belong to a subset of QUAD systems, satisfying D 2 vI # 0. In order to relax these assumption, in the following theorem a specific adaptive law is selected. Theorem 2 Assume that f is QUAD and G . 0. If we chose the adaptive law # (14) sij 5 1 xi 2 xj 2 TG 1 xi 2 xj 2 , then network (9) reaches asymptotic synchronization and each coupling gain sij converges to some finite value. Therefore, with this choice of g, we can show that any QUAD system, coupled on all its state variables through G, can be asymptotically synchronized through a decentralized adaptive modulation of the coupling gains. It is also possible to generalize the above Theorem and give sufficient conditions for asymptotic synchronization,

even when the coupling is only on a subset of the state variables. A fully decentralized pinning control strategy is presented in [32]. 3.1. Numerical Example In what follows, we consider, as a testbed example, a network of Lorenz systems [63, 72]. This nonlinear model is described by the set of the following three differential equations: # p 5 a1 1 q 2 p 2 # q 5 a2p 2 q 2 pr # r 5 pq 2 a3r,

(15) (16) (17)

where x 5 3 p, q, r 4 T is the state vector and a1, a2 and a3 are three positive parameters. In our simulation, we consider a network of N 5 500 Lorenz systems, coupled only on the first two state variables, that is 1 G 5 £0 0

0 1 0

0 0§. 0

The network topology is generated using the BarabàsiAlbert scale free model (we address the reader to [8] for further details). As for the node dynamics, we choose the parameters of the model guaranteeing a chaotic behavior. Namely, we set a1 5 10, a2 5 28 and a3 5 8/3. Moreover, the initial conditions on the states are taken from a normal distribution with 0-mean and standard deviation equal to 40. From trivial manipulations of the results achieved in [60], we have that the Lorenz system satisfies the assumptions. Therefore, we can select the adaptive law (14) to synchronize the network, with null initial conditions on the edge states, that is, sij 1 0 2 5 0. As we can see from Figure 1, synchronization is asymptotically achieved and all the coupling gains converge to finite steady state values. 3.2. Designing Protocols via Contraction Theory Typically, network coordination is reduced to a stability problem of some invariant set in the network phase space. The analysis/design of decentralized coordination strategies is then addressed by means of Lyapunov based techniques (Box 2), and by making some hypotheses on the vector field of each of the nodes (see e.g. the QUAD condition above). An alternative approach is that of using contraction theory, instead of Lyapunov tools, for proving network coordination (Box 2). The

A function F : I S R is positive definite if F 1x 2 . 0, 4x [ I, x 2 0 and F 10 2 5 0. A function f : R $0 S R $0 is of class k if it is continuous, positive defi nite and strictly increasing.

1

70

IEEE CIRCUITS AND SYSTEMS MAGAZINE

THIRD QUARTER 2010

xi (t )

50 0 −50 0

0.5

1

1.5 t (a)

2

2.5

0

0.5

1

1.5 t (b)

2

2.5

σij (t )

3 2 1 0

Figure 1. Network of 500 Lorenz systems coupled via edgebased strategy (14). (a) Evolution of the states and (b) the coupling gains.

main advantage of this approach lies in the fact that stability is defined as an incremental property of the trajectories themselves. Notice that the invariant set towards which they all converge is not known a priory. This, in turn, allows to separate the analysis/control problem in two well-defined steps: i) determine some invariant subspace for the network evolution; ii) prove/ impose contraction. The general principles used for the proofs of the results presented in what follows are based on two key ideas. We prove contraction of only some directions in phase space (the directions transverse to the synchronization manifold). Notice that contraction is a property of network trajectories, which is metric-dependent. To prove this directional contraction we make use of non-Euclidean norms. We illustrate the design of decentralized network coordination strategies and turn our attention to the problem of imposing some desired coordinated behavior for a set

where rix,y and vix,y represent the position and velocity vectors of agent i in a fixed reference frame. As in (1), the agents are coupled via a coupling function ci 1 x, t 2 5 3 0, 0, aix, aiy 4 T, with aix and aiy being acceleration inputs on the agent dynamics, chosen as: aix 5 s1 a 1 xj 2 xi 2 2 1 s1 1 s2 2 vix j[Ni

aiy 5 s3 a 1 yj 2 yi 2 2 1 s3 1 s4 2 viy .

1

0.4

0.8

0.2

0.6 viy

0

−0.2

2

2.5

3 (a)

3.5

4

0.4 0.2

−0.4

0

−0.6

−0.2

−0.8

(19)

j[Ni

0.6

vx

2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

of N . 1 nodes. The general principle that we use in our proofs can be summarized as follows. We differentiate the network dynamics so as to obtain network Jacobian, J. Now, let M be an invariant subspace for the network and define as V the orthonormal matrix spanning the null of M. Then, all network trajectories globally exponentially converge towards M if the matrix V J V T is contracting as defined in Box 2. An equivalent approach, which has no straightforward geometric interpretation, makes use of some virtual system, see [85] and [83] for further details. In Figure 2 an application to the rendezvous-problem is shown. Rendezvous is a typical example of coordination, that has also been used as a testbed problem in the literature. Namely, the problem is that of designing a distributed control strategy driving all the agents towards a common point in space, where they all have zero velocities (see for example [34, 58] and references therein). The rendezvous problem can be then formalized as the problem of finding an appropriate (distributed) control strategy which guides all the agents towards a common position, say r, i.e. such that limt S ` ri 1 t 2 5 r, 4i 5 1, c, N. In accordance with the existing literature [52, 101], we choose the nodes’ dynamics as: # # rix 5 vix, riy 5 viy, # # vix 5 aix, viy 5 aiy, (18)

0

5

10 Time (b)

15

20

−0.4

0

5

10 Time (c)

15

20

Figure 2. Designing communication protocols solving the rendezvous problem. (a) Evolution of agents’ positions in the x 2 y plane. (b) and (c) Shows time evolution of velocities vix, vij.

THIRD QUARTER 2010

IEEE CIRCUITS AND SYSTEMS MAGAZINE

71

2

1

6

5

3

4 (a) 5

xi (t )

0

−5

−10

0

5

10

15 Time (b)

20

25

30

Figure 3. (a) Network used to test the multi-scale approach and (b) temporal behavior of (20) with R 5 1, G 5 K 5 0.9.

We assume that the sensing region of each agent is bounded in space and is represented by a disk of radius R. In practice, this implies that each agent communicates only with those agents inside its sensing disk, i.e. the i-th agent communicates only with its neighbors, Ni. From the topological viewpoint, the assumption of limited sensing region makes the network of interest switched as the neighbors of the i-th agent can enter/exit its sensing disk. Clearly, the switches of network topology (due to the interactions among agents) occurs with a finite frequency. As shown in Figure 2, using contraction theory it is possible to show convergence of all trajectories towards each other (for more details see [86]). Sometimes, designing the coupling protocol alone is either not feasible (e.g. quorum sensing in biological networks) or insufficient to guarantee synchronization. In this case, it is possible to engineer the network topology and/ or the node dynamics to guarantee synchronization. Contraction theory can again be used as an invaluable tool for this purpose. The idea is to break down the analysis and design into two independent steps: (a) At a global level, properties of the network or interconnection graph are imposed so as to guarantee a desired behavior for the full interconnected system. 72

IEEE CIRCUITS AND SYSTEMS MAGAZINE

In this analysis, subsystems may be characterized as “black boxes” with assumed input-output characteristics, but detailed knowledge of their internal structure is not required. (b) At a local level of analysis, one imposes constraints on the structure and behavior of individual subsystems (components), so as to fit the requirements of the global approach. These requirements are verified independently of the overall network structure. This multi-scale or hierarchical methodology is robust in so far as a large degree of uncertainty can be tolerated in the components, only constrained by meeting appropriate behavioral requirements. There are many examples of such approaches in systems/control theory, including among others (1) the use of small-gain theorems to guarantee stability of a negative feedback loop provided that the components are individually stable (qualitative property of components) and the overall loop H` gain is less than one, as well as nonlinear generalizations based on input to state stability [26, 93], (2) input/output monotone systems theory [2, 3, 94], in which input-output characteristics are the only required “quantitative” data, and (3) the use of passivity-based tools [5, 69, 102]. We now present some examples where the problem is solved using contraction theory, see [88] for further details. Specifically, it is possible to prove that a contraction property on a reduced-order matrix that quantifies the interconnection structure, coupled with contractivity/expansion estimates on the individual component subsystems, suffices to ensure that trajectories of the overall system converge towards each other. See [88] for further details. Note that the construction of such reduced-order matrix uses only norm and matrix-measure estimates, but no precise knowledge, of components (i.e. of the nodes composing the network). As a first example, we show the result of the design for a network of Hopfield neural models in Figure 3. Here the nodes are coupled via the nonlinear coupling functions in Table 1 engineered using contraction to guarantee synchronization (see [88] for further details).

Table 1. Coupling function for network (20). Edge 142, 6 S 2, 2 S 3, 5 S 6 3 S 1, 5 4 2 344

hij 1 x 2

G arctan 1 x 2 Kx 1 1 2 e 2x 2 / 1 1 1 e 2x 2

THIRD QUARTER 2010

Figure 3 shows one of the application of such a multiscale approach. Each of the nodes in such a figure, is a Hopfield neural model [48]: cl

LuxR

Luxl

Al

u (t ) Injection of u (t)

4. Evolving the Network Topology In all the network models presented above, the topology of the connections is considered fixed. Conversely, as mentioned in the introduction, real networks are characterized by evolving topologies. Similar mechanisms appear also in technological applications, such as flocking and rendez-vous. Here, we propose different possible ways of evolving the network dynamics. In Sections 4.1, we propose to evolve the network on-line, on the basis of the dynamical evolution of the network nodes. The edge evolution is modeled as a second order dynamical system, stimulated by an external input, that is the norm of the mismatch between the endpoints of the edge. In Section 4.2, a novel computational tool is presented for the off-line evolution of network topologies that optimise some user-specified performance measure. 4.1. Edge-Snapping The network model (9), (10) allows the decentralized adaptation of the coupling gains associated to each edge in the network. Nonetheless, the network structure is fixed a priori, while we aim at a network that is able

Cell 1

u (t ) Cell 2

u (t ) Cell 3

u (t )

Cell 4

u (t ) Cell 5

u (t ) Cell n

u (t )

(a) cI (arbitray units)

Another notable example is that of synchronization in biological systems where the coupling between nodes is determined by physical characteristics of the systems and of the external environment, which are hard or even impossible to modify with the technologies currently available. To overcome the above technological problem, our approach is that of engineering the intrinsic dynamics of the nodes by tuning their internal parameters. Those parameters are related to biochemical characteristics of the systems of interest, like e.g. protein degradations and promoter strengths which can be easily modified. From a methodological viewpoint, as shown in [82, 89, 90], the above biochemical parameters can be tuned so as to guarantee contraction of the node. Then, network synchronization is attained, since networks of contracting nodes synchronize in the presence of diffusive coupling. To this aim, we use a graphical algorithm, which allows to easily prove/impose contraction [82, 84, 86, 87]. Figure 4 shows an application to the synchronization of a network of Repressilators, coupled by means of a quorum sensing mechanism: the shared quantity, representing the environment where nodes live is in turn affected by an exogenous input, u 1 t 2 .

THIRD QUARTER 2010

Lacl

Lacl

(20)

u(t) (arbitray units)

xi # xi 5 2 1 a 1 hij 1 xj 2 2 hij 1 xi 2 2 1 u 1 t 2 . R j[Ni

TetR

200 150 100 50 0

0

50

100

0

50

100

150 (b)

200

250

300

150 200 Time (min) (c)

250

300

0.8 0.6 0.4 0.2 0

Figure 4. The Repressilator is a synthetic biological circuit that consists of three genes that inhibit each other in a cyclic way [37]. Specifically, gene lacI (associated to the state variable c in the model) expresses protein LacI, which inhibits the transcription of gene tetR. This translates into protein TetR, which inhibits transcription of gene cI. Finally, the protein CI, translated from cI, inhibits expression of lacI, completing the cycle. A modular addition to the three-genes circuit is shown in (a). The module was first presented in [40] and makes the Repressilator circuit sensitive to the concentration of the auto-inducer which is a small molecule that can pass through the cell membrane. Specifically, the module makes use of two proteins: (i) LuxI, which synthesizes the auto-inducer; (ii) LuxR, with which the auto-inducer synthesized by LuxI forms a complex that activates the transcription of various genes. In this figure, the behavior of one protein, representing the output of the Represilator circuit, cI, is shown. Notice that at steady state synchronization is attained onto a periodic orbit having the same period as u 1 t 2 5 0.4 1 0.4 sin 1 0.5t 2 .

to self-determine its topology. With this in mind and according to the model proposed in [29], we associate a dynamical system to each edge in the network, with each coupling gain being the one-dimensional output of the dynamical system: # sij 5 g 1 sij, xj 2 xi, t 2 ,

sij 5 c 1 sij 2 .

(21) (22)

IEEE CIRCUITS AND SYSTEMS MAGAZINE

73

where k is a class k` function2 of the norm of the mismatch between the states of the endpoints. In what follows we choose k 1 z 2 5 az2 while the potential V is such that system (23), (24) has two stable equilibria in 0 and 1, representing respectively the absence and the activation of the corresponding edge. Namely, we select the potential V 1 z 2 5 bz2 1 z 2 1 2 2 which is depicted in Figure 5(a). In this model, the set E of the edges can be viewed as the set of network edges that can be potentially activated. In fact, for physical or economical reasons, there are surely some edges that cannot be activated. Nonetheless, we leave the network nodes free to negotiate what edges have to be activated at steady state among all the feasible links. We think of the coupling gain associated to each edge as a mass which is at rest at the coordinate origin, and then is pushed to the right by an external force given by the squared mismatch between the states of the endpoints of the edge. If this force is strong enough, then the gain snaps to the second equilibrium point (link on), otherwise it asymptotically comes back to the origin, so that the corresponding edge is respectively activated or switched off.

V(σij)

b /2

b /4

b /16 0

0

0.5 σij (a)

1

1.5

xi (t)

20 0 −20 0

2

4

6

8

t 1 σij (t )

σij (t )

20 10 0

0

0.5

0

5

10

15

t

t (b)

Figure 5. Network of 50 Lorenz oscillators coupled through edge snapping: (a) bistable potential driving the evolution of each sij; (b) evolutions of the node states (top) and of the coupling gains (bottom).

Here, sij [ Rm is the state associated to the edge 1 vi, vj 2 , g : Rm 3 Rn 3 R 1 S Rm is the vector field defining the evolution of the edges, and c 1 t 2 : Rm S R is the output function of the edge dynamics, defining the coupling gains. Notice that the network (9), (21), (22) is an example of a generalized dynamic graph without external inputs. It is worth remarking here that an appropriate selection of the function g is fundamental to achieve our goal. According to [29], we consider the so-called snapping dynamics for the edge evolution, given by: # sij112 5 sij122,

(23)

d # sij122 5 2zsij122 2 112 V 1 sij112 2 1 k 1 00 xj 2 xi 00 2 , dsij

(24)

sij 5 sij112 ,

(25)

2

where z is a damping parameter. Therefore, the coupling gain is modeled as a unitary mass in a double-well potential V subjected to an external force k 1 7 xj 2 xi 7 2 ,

Numerical Example As a numerical example, we consider a network of 50 Lorenz systems, coupled on all the state variables. We set the damping z 5 20 while the parameter b of the potential, determining the height of the barrier between the two equilibria, is chosen equal to 16. The set of potential edges, represented by the set E, is generated using a Barabàsi-Albert scale free model. The edge states are assumed to be initially in the coordinates origin 1 sij112 1 0 2 5 s2ij 1 0 2 5 0 2 , implying that at the onset of the evolution the network topology is disconnected. As for the nodes’ initial conditions, they are taken from a random normal distribution with zero mean and standard deviation of 15. As we can see from Figure 5(b), synchronization is asymptotically achieved, while at steady state a topology emerges, which is depicted in Figure 8. The features of this topology will be described later in Section 4.3. 4.2. NetEvo: A Computational Approach An alternative approach is to compute some optimal network topology offline to minimize some performance measure of interest. NetEvo [42] is a computational framework developed to study the evolution of dynamical complex networks3. It provides tools to simulate, evolve and analyze a wide range of systems in the hope of understanding possible links between topology, dynamics and evolution. A central concept within NetEvo is that of a supervisor.

A function F : I S R is positive definite if F 1x 2 . 0, 4x [ I, x 2 0 and F 10 2 5 0. A function f : R $0 S R $0 is of class k if it is continuous, positive defi nite and strictly increasing. An unbounded function of class k belongs to class k`. 3 For further information see the project website at http://www.netevo.org. 2

74

IEEE CIRCUITS AND SYSTEMS MAGAZINE

THIRD QUARTER 2010

This entity has the role of making changes to a systems underlying topology and dynamical parameters, with the aim of minimizing a user-specified performance measure Q. This acts as a guide to the evolutionary process and can be calculated using both topological information and simulated dynamics from the system. Figure 6 illustrates the main components of a supervised network. To use NetEvo a user must provide an initial topology, the form of node and edge dynamics, a method to alter the network topology and dynamical parameters, and a performance measure. Together with a set of internal modules (see Fig. 7) these features are combined and evolution of the system performed. To aid with the analysis of this process, NetEvo has been built using several well established open-source libraries. These include igraph to hold and analyze network topologies, and the GNU Scientific Library for numerical simulation and spectral methods. At present, NetEvo allows for dynamics in the form of ordinary differential equations (ODEs) and contains a built-in supervisor based on a simulated annealing metaheuristic. This supervisor assumes no knowledge of structure within the system of interest and will perform an unbiased probabilistic search of system configurations. An interesting future direction is to improve this search process to make use of either expert knowledge, i.e. it has been proven that certain topological features yield higher performance, or by machine-based learning techniques to find structural and dynamical relationships on-the-fly. In regards to synchronization, NetEvo has been used to study the effect that allowing a systems simulated dynamics to influence the evolutionary process has on enhanced topologies [44, 43]. Previously, the study of evolving topologies for optimal synchronization has focused on topological features of a system [35]. This is in part due to the synchronizability of a network— range of coupling strengths for with full synchronization is achieved—being linked by Pecora and Carrol to minimization of the network Laplacian eigenratio (lN/l2) [9, 75]. Although this type of approach provides insight

Node Dynamics

Edge Dynamics

Dynamics Simulator ODE Solver System Structure Graph and Dynamics

Data Output Numerical and Visual

Updates Parameters

Node/Coupling Dynamics dxi /dt = F (xi ) + coupling

Alters Coupling

Supervisor Computerized Agent

Q

Performance Measure Network Topology Laplacian/Adjacency Matrix

Updates Topology

Figure 6. Flow diagram of a supervised network. The supervisor attempts to minimize a performance measure Q by making changes to the underlying topology and parameters of the dynamics.

into proven optimal topologies, it does not necessarily tell us anything about synchronous systems found in Nature, where evolution is instead directed by the dynamical attributes (phenotype) of an organism. To address this, we attempted to take an alternative approach and used an order parameter performance measure [110] (calculated from system-level dynamics) to guide the evolutionary process. A number of separate evolutions were performed for systems of identical Rössler oscillators with diffusive coupling, and a variety of fixed coupling strengths s^ during the process. Varying the fixed coupling strength altered the stability of the synchronization manifold and allowed us to investigate how the resultant topologies adapted to alternative scenarios, i.e. when full synchronization was not possible. Furthermore, we also considered systems of alternative node dynamics (Lorenz and Chua) and levels of heterogeneity in the dynamical parameters of each node. 4.3. Features of the Evolved Topologies In this section, we discuss some of the emerging features found when using NetEvo and the edge snapping techniques to carry out network evolution. Furthermore, network topologies are analyzed with

Performance Measures

Network Mutations

Supervisor Simulated Annealing Performance Multithreading

Analysis Dynamical and Topological

Figure 7. The main modules of NetEvo split into two main types, core and interface. Core modules, shown in orange, are designed to carry out analysis, simulation and evolution. Interface modules, shown in blue, are to describe the particular system of interest and should be provided by an end-user.

THIRD QUARTER 2010

IEEE CIRCUITS AND SYSTEMS MAGAZINE

75

BOX 2. PROVING STABILITY

O

The Master Stability Function Approach

well as an EDN are well described by a set of differential

The Master Stability Function (or MSF) is an approach for

equations. In terms of the notation introduced in Box 1, this

studying synchronization of diffusively coupled networks (2)

means that both the operators F and t satisfy some non-

mainly used within the Physics community. Specifically, net-

triviality, restriction, semigroup, identity axioms (see e.g.

work dynamics are first linearized around the synchronization

[95]). The problem of analyzing the steady state behavior

manifold so as to obtain a variational equation from (2), de-

of (1) has been addressed in the literature using several ap-

scribing small variations, jk, of network trajectories from the synchronous evolution, say xs 1 t 2 . This equation is then block

ften, in applications a generalized dynamic graph as

proaches, briefly itemized in what follows. Lyapunov Based Proofs

The key idea for the proofs of synchronization based on the use of Lyapunov functions is that of recasting a synchronization problem onto a stability problem of some fixed point. In this view, the typical approach is that of defining an error variable, ei J xi 2 xs, where xs 1 t 2 is the synchronous soT T lution (which is assumed to exist). Let e J 3 e1 , c, eTN 4

and notice that network global synchronization is attained if e 5 0 is a globally asymptotically stable fixed point for the # dynamics of e. Now, stability of such a fixed point can be analyzed by means of some Lyapunov Function, V 1 e 2 . A typical choice for V 1 e 2 is

V 1 e 2 J eTMe. The matrix M is assumed to be positive definite. Depending

'h 1 xs 2 'f 1 xs 2 # 2 slk jk 5 c d jk, k 5 1, c, N. (26) 'x 'x Each of the above equations, or modes, represents the variations (around the synchronization manifold) along a direction of the network phase space. The mode associated to k 5 1 corresponds to the variations along the synchronization manifold, while all other k’s correspond to the transverse eigenmodes to such a manifold. Thus, local synchronization occurs if all such transverse modes are stable. To this aim, the maximum Lyapunov exponent is computed as a function of s, for k . 0. It immediately follows that stability is then guaranteed for all the s’s such that all the maximum Lyapunov exponents are negative. See [75, 74, 9] for further details.

on the choice of such a matrix and on the hypotheses made # on the network nodes which make V , 0 along network

Contraction Based Proofs

trajectories, some sufficient condition for synchronization

lies in the fact that synchronization is seen as a property of

can be then obtained. See e.g. [60, 61, 64, 55, 28, 59] for

some invariant set. The viewpoint of contraction theory is,

further details.

instead, completely different. In the contraction framework,

the aim of better understanding if properties of the emerging topologies can be tuned by parameters of the evolutionary process. Edge Snapping It is straightforward that the main constraint influencing the properties of the emerging topology is the fundamental edge snapping topology E, representing the set of potential network edges. Therefore, to avoid any bias in the analysis of the emerging topology, we refer to the case in which E 5 V 3 V. This is the case in which, at the onset of the evolution, each couple of nodes can communicate, and then decide whether activating or not the corresponding edge. Moreover, in the analysis, we exclude the case of very low potential barriers b, in which almost all the potential nodes are activated. 76

diagonalized to give a set of decoupled variational equations:

IEEE CIRCUITS AND SYSTEMS MAGAZINE

The common feature of the above presented approaches

In [29] it was observed that, for many diverse node dynamics, the initial conditions x 1 0 2 of the network nodes plays a key role in determining the main topological features of the emerging topology, such as the degree distribution or the betweenness centrality. We recall here that the degree of a node i is the number of incident edges, while the betweenness centrality is a measure of importance of a node in the network [4, 39]. Therefore, a node starting from a more scattered initial condition is more likely to activate more connections and have a higher betweenness centrality. Moreover, it was also observed that the edge snapping dynamics induce a topology in which the maximum eigenvector (that is, the eigenvector associated to the maximum eigenvalue) of the Laplacian matrix L is strongly related to the initial conditions of the nodes. Therefore, in the case of the classical consensus problem (see for THIRD QUARTER 2010

stability is intended as a property of trajectories and not as

• a geometric approach, which allows to study also

a property of some invariant set. From a qualitative view-

cluster synchronization, is presented in [76]. Here,

point, a contraction region is a convex open set in phase

contraction is proved only for some directions in

space, where all the trajectories converge towards each

phase space. Specifically, consider the linear invariant

other. A system is said to be contracting if the contraction

subspace of (1), M. Let V be an orthonormal matrix

region coincides with the whole state space. Formally, a

spanning the null of M. Then, the trajectories trans-

system is said to be contracting if there exist some matrix

versal to the invariant subspace globally exponentially

measure for system Jacobian which is uniformly negative

contract towards M if the matrix VJV T is contracting.

definite, see [62, 89]. From the methodological viewpoint,

In network synchronization problems, the subspace M

two main approaches are used to study network conver-

can be chosen as x1 5 c 5 xN.

gence via contraction. Specifically: • a first idea is to use the concept of partial contraction discussed in [104]. Basically, given a smooth nonlinear # n-dimensional system of the form x 5 f 1 x, x, t 2 , assume # that the auxiliary system y 5 f 1 y, x, t 2 is contracting with respect to y. If a particular solution of the auxiliary y-system verifies a smooth specific property, then all

Linking Stability Tools

In [30] it is shown that when the nodes’ dynamics are QUAD within a certain range of parameters, then it is equivalent to contraction. Similarly, in [83] it is shown that if network dynamics are contracting towards the synchronization manifold, then the MSF is negative.

trajectories of the original x-system verify this property C

exponentially. The original system is said to be partially

contracting. Indeed, the virtual y-system has two par-

ticular solutions, namely y 1 t 2 5 x 1 t 2 for all t $ 0 and

x1(t0)

the particular solution with the specific property. Since all trajectories of the y-system converge exponentially to a single trajectory, this implies that x 1 t 2 verifies the specific property exponentially. The specific property above may be e.g. a relationship between state variables, or simply a particular trajectory. In the case of network synchronization, such a relationship is y1 5 c 5 yN.

instance [70]), in which f 1 x, t 2 5 0, edge snapping can be used as a viable method to build topologies wellsuited to quickly synchronize integrators with given initial conditions. These features are not significantly affected by the height of the potential barrier b or by the damping parameter b. Nonetheless, these parameters can be used to tune the other topological properties. For instance, given the network initial conditions, the barrier b can be used to increase the sparseness of the topology, while preserving network synchronizability. As an example, we show the effect of the increase of the barrier in a network of 30 Lorenz oscillators. In this case we consider a random network with average degree 5. As illustrated in Figure 8, as the barrier is increased from 1 to 10, the percentage of active nodes reduces from the 87.5% (Fig. 8(a)) to the 78% (Fig. 8(c)). THIRD QUARTER 2010

x2(t0)

x1(t ) x2(t )

Figure 10. A schematic representation of trajectories inside the contraction region, C. In such a region, any two trajectories, i.e. x1 1 t 2 , x2 1 t 2 , having initial conditions, i.e. x1 1 t0 2 , x2 1 t0 2 , belonging to C converge towards each other.

NetEvo Using NetEvo and the dynamical order parameter performance measure, we see yet more features emerging in the resultant topologies. Some representative examples can be found in Fig. 9. For comparison, Fig. 9(a) shows a network evolved using the topological eigenratio performance measure. This has what is termed an “entangled” structure with a low diameter and narrow betweenness and degree distributions. Of these evolved topologies it is clear to see two main groups. The first consists of our eigenratio evolved network and the order parameter evolved network for Rössler dynamics and a fixed coupling of s^ 5 0.6, shown in Fig. 9(a) and 9(b) respectively. Both these topologies show a convergence in statistical properties towards very homogeneous features and display high levels of synchronization. The fixed coupling strength s^ 5 0.6 falls IEEE CIRCUITS AND SYSTEMS MAGAZINE

77

(a)

(b)

(c)

Figure 8. Network of 30 Lorenz oscillators coupled through the edge snapping topology. Increasing the barrier b between the wells of the potential, the sparseness of the network is increased. In orange are represented the edges belonging to the emerging topology, while in gray are represented the edges that are asymptotically switched off. In (a), the barrier b is set to 1 and the percentage pa of activated edges represents is the 87.5%. In (b), b is increased to 5, with a subsequent decrease of pa to 82.3%. In (c), b 5 10 implies pa 5 78%.

within the stable region of the synchronization manifold. Therefore, the features we see here are likely to be representative of those necessary for optimal synchronization. We will refer to these types of networks as Type 1. The second group, containing Fig. 9(c)–9(e), show very different topological properties. Rather than the homogeneous appearance of Type 1 networks, we see the formation of hubs containing low degree node that surround a highly interconnected central region. We broadly refer to networks with these properties as Type 2. This change is caused by selecting a fixed coupling strength s^ outside the stable region for synchronization. By imposing this constraint, the dynamical performance measure has adapted the resultant topology for the new scenario. This is further shown in the amount of synchronization exhibited by these networks. Although they are unable to see full synchronization, they do see partial synchronization for a much wider range of fixed coupling strengths and at strengths where “entangled” or Type 1 networks see none.

(a)

(b)

(c)

A fascinating aspect of this change is the speed at which it happens. As the fixed coupling strength is slowly varied towards the unstable region, there is not a smooth transition in topological features from Type 1 to Type 2. Instead, the system jumps from one form to the other in what we call a topological bifurcation [44]. It is worth noting that topological measures do not consider dynamics at all and so cannot display adaptive behavior due to changes in a systems dynamical parameters. Another interesting feature shown to play an important role in many complex networks is that of motifs; small sub-graphs found at greater frequencies than would be expected by random chance [67, 1]. To understand if these played a part in our Type 1 and 2 structures, average motif distributions were calculated. Results show (see Table 2) that networks evolved using a topological measure (eigenratio) display relatively few types of motif statistically over or under expressed. In contrast, a more striking pattern is seen when a dynamical performance

(d)

(e)

Figure 9. Topologies evolved using NetEvo for various performance measures, node dynamics and fixed coupling strengths s^ . (a) Eigenratio, (b) Rössler order parameter s^ = 0.6, (c) Rössler order parameter s^ = 0.9, (d) Lorenz order parameter s^ = 4.5, (e) Chua order parameter s^ = 0.85. Two main forms of topology are evident with (a) and (b) displaying “entangled” Type 1, and (c) to (e) hub-like Type 2 features.

78

IEEE CIRCUITS AND SYSTEMS MAGAZINE

THIRD QUARTER 2010

Table 2. Average motif frequencies for topologies evolved using NetEvo. Arrows represent statistically significant over c and under T expression in comparison to a randomized network (p-value * 0.01). Networks evolved using a dynamical performance measure show an over expression of closed loop feedback motifs 2, 4 and 6. Dynamics

Performance measure

Topology

Eigenratio

Type 1

99.90

0.10

Order Parameter s^ = 0.6 Order Parameter s^ = 0.9

Type 1 Type 2

96.74 ≈ 97.18 ≈

3.26 Ω 2.82 Ω

Lorenz

Order Parameter s^ = 4.5

Type 2

98.39 ≈

Chua

Order Parameter s^ = 0.85

Type 2

93.98 ≈

– Rössler Rössler

1

measure (order parameter) is used, regardless of coupling strength or node dynamics. In all cases, these distributions exhibited an over expression of three node closed loop feedback motifs 2, 4 and 6. In [44] it was conjectured that these may help improve localized stability of a shared dynamic and would therefore be an artifact of dynamics influencing topology through the evolutionary process. Furthermore, and in support of this idea, evolution of networks with differing levels of variability in node dynamics showed a correlation up to a threshold between motif expression frequency and the percentage of variation in dynamical parameters [43]. Further investigation into these features is still required, but points towards the interesting possibility of using localized structures to help shape system-level behaviors. 5. Conclusions and Future Work We have shown that synchronization is indeed possible in evolving dynamical networks in a number of different ways. Starting from networks with constant diffusive linear coupling, we showed that it is possible to consider decentralized adaptive gain evolutions. In this case, the strength of the coupling between neighboring nodes is determined by a set of integro-differential equations describing how each of the time-varying gains in the network should be adjusted as a function of the mismatch between the states of those nodes connected to each other. We then showed that an alternative approach is that of using contraction theory so as to construct nonlinear coupling protocols that make trajectories of the nodes in the network converge towards each other. This, we believe, is a particularly timely approach in all those problems where the interest is to prove the network synchronizes rather than showing stability of some asymptotic solution of interest (which is often unknown a priori in synchronization problems). To further mimic the adaptation and time-varying nature of networks of agents in the natural world, we then THIRD QUARTER 2010

2

3

4

5

6

7

19.71 ≈

0.20

15.95 ≈ 40.25 Ω

5.83 Ω 5.53 Ω

1.61 Ω

25.90 ≈

6.02 Ω

36.74 Ω

8

80.05 Ω

–

0.04 ≈

–

77.77 ≈ 52.50 ≈

0.17 Ω 0.33 Ω

0.28 ≈ 1.09 ≈

0.002 Ω 0.01 Ω

3.58 Ω

69.56 ≈

0.11 Ω

0.86 ≈

–

8.60 Ω

51.93 ≈

1.36 Ω

1.28 ≈

0.10 Ω

proposed that the structure of the network can itself be evolved either to guarantee synchronization (edge-snapping) or in order to optimize some cost function of interest. NetEvo was presented as an effective computational framework to deal with this latter type of problems. Interestingly, we noticed that when the network structure is evolved, its topological features depend on the dynamics at the nodes, the cost function and even the structure and parameters of the coupling function. Feedback motifs were for example found to be more likely to occur when dynamic performance measures are used rather than static ones. Despite progress in the area of complex network theory, there are many pressing open problems that still need to be addressed. First of all, synchronization itself is a rather basic type of emergent phenomenon. One might want to induce more sophisticated types of complex behavior in networks. For instance, concurrent or cluster synchronization has been proposed as an important phenomenon in neural networks [76] and other biological networks [65, 109] where different groups of interacting agents synchronize onto different evolutions. An open problem is to understand if it is possible to devise fully decentralized strategies to achieve this aim and, maybe, even evolve network topologies guaranteeing the emergence of cluster synchronization. Another interesting challenge is to further characterize the topological features of evolving dynamical networks and understand what the mechanisms are that induce the presence of modularity, motifs and other topological features in the evolved networks. For example, in [51], it is suggested that modularity and motifs can spontaneously emerge in the presence of switching network goals (cost functions). Here we noticed that, when dynamics is considered, certain motifs can become more frequent than others but there are certainly other possible causes, still unaccounted for, that might explain some of the features observed in natural networks. Also, it is currently assumed that the network nodes share the same identical dynamics, often, modeled in term IEEE CIRCUITS AND SYSTEMS MAGAZINE

79

of some smooth vector field. In practice, nodes will not be identical and the case of networks of near-identical or nonidentical nodes is seldom studied in the literature (some results can be found in [99]) and is waiting to be properly addressed. A striking open problem is when the nodes or coupling protocols are described by piecewise-smooth models [33] as is the case for example of networks of vibroimpacting mechanical devices, gene regulatory networks and switching converters. Finally, the grand challenge of optimally evolving dynamical network structures to perform a desired function remain unsolved. Can we use contraction theory for instance to evolve optimal network topologies? Can we always reduce the problem to some convex optimization problem e.g see [18]? These are only some of the many challenges and open problems that remain to be solved in this new exciting area of research in circuits and systems. The task ahead of us and the next generation of scientists is gigantic. It is likely new analytical and computational tools will be needed and, possibly, even a completely new theoretical framework. But, if successful, the potential impact of these new tools on our understanding of complex systems would be immense and the possible applications so promising that attempting to solve these problems is of utmost importance. Pietro De Lellis is currently a postdoctoral researcher at the University of Naples Federico II. He got his Master Degree in Management Engineering in September 2006 from the same university. His Masters Thesis was dealing with the analysis of the dynamics in complex networks of agents. In November 2006 he won the competition for entering the Ph.D. course in Informatics and Automation Engineering. From January to July 2009 he was a visiting Ph.D. student at the Department of Aerospace and Mechanical Engineering of the Polytechnic Institute of NYU. He received his Ph.D. on December 2009 from the University of Naples Federico II. His current research activity is focused on the analysis and control of complex networks with special emphasis on the novel adaptive strategies for the network synchronization. Mario di Bernardo is currently Associate Professor of Automatic Control at the University of Naples Federico II in Naples, Italy. In 1998 he obtained a Ph.D. in Engineering Mathematics from the University of Bristol, U.K. He was appointed to a Lectureship at the Department of Engineering Mathematics of the same University in 1997 and then promoted to a Readership and a Full Professorship. From 2001 till 2004, he was Assistant Professor of Automatic Control at the University of Sannio, Italy. His research interests 80

IEEE CIRCUITS AND SYSTEMS MAGAZINE

are within the broad area of nonlinear systems, on both dynamics and control. He authored and co-authored more than 150 international scientific publications. From July 1999 to December 2001, he served as Associate Editor of the IEEE Transactions on Circuits and Systems Part I, and from January 2002 to 2006 he was Associate Editor of the IEEE Transactions on Circuits and Systems Express Letters. He is currently AE of the IEEE Transactions on Circuits and Systems Part I. He is a member of the organizing committees of the IEEE Symposium on Circuits and Systems and has been chair or co-chair of many scientific events. In 2004 he was elected to the governing board of the Italian Society for Chaos and Complexity and in 2006 and 2009 to the Board of Governors of the IEEE Circuits and Systems Society. He was elected President of the Italian Society for Chaos and Complexity for the term 2010–2014. He received funding from major funding bodies and industries including the EPSRC, the European Union, the Italian Ministry of Research and University, Jaguar Engineering Centre, QinetiQ. Together with Bernard Brogliato (INRIA, France), he was the organizer and scientific coordinator of the EUR 2.8M EU Project SICONOS on the simulation and control of nonsmooth dynamical systems. On the 28th February 2007, he was honoured with the title of ‘Cavaliere della Repubblica Italiana’ (equivalent to a British OBE) for scientific merits by the President of the Italian Republic. Thomas E. Gorochowski received his Computer Science M.Eng. degree from the University of Warwick, UK in 2004. From 2004 until 2007 he was with Accenture Ltd. working as a technology consultant in the domain of business intelligence and information management. Following this, he joined the Centre for Complexity Sciences at the University of Bristol, UK, receiving a Complexity Sciences M.Res. degree in 2009 and has continued studying towards a Ph.D. in Complexity Sciences at this University. His research interests include the study of structure, dynamics and evolution of complex networks with possible applications to engineering (specifically synthetic biology), novel computing architectures and in the analysis of biological systems. Giovanni Russo was born in Naples on March 28, 1984. He received from the University of Naples Federico II, the Master Degree in Automation Engineering on October 30, 2007 with the thesis “Modelling, Analysis, Identification and Control of Biological Networks”. In November 2007 he won the competition for entering the Ph.D. course at the Department of Systems and THIRD QUARTER 2010

Computer Engineering of the University of Naples Federico II. His research interests are focused on the analysis and control of biochemical and complex networked systems. He spent the period from October 1, 2009 to March 31, 2010 at the Nonlinear Systems Laboratory of the Massachusetts Institute of Technology. References [1] U. Alon, “Network motifs: Theory and experimental approaches,” Nat. Rev. Genet., vol. 8, pp. 450–461, 2007. [2] D. Angeli, J. E. Ferrell, and E. Sontag, “Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems,” Proc. Nat. Acad. Sci. USA, vol. 101, pp. 1822–1827, 2004. [3] D. Angeli and E. Sontag, “Monotone control systems,” IEEE Trans. Automat. Contr., vol. 48, pp. 1684–1698, 2003. [4] J. M. Anthonisse, “The rush in a directed graph,” Stichting Mathematisch Centrum, Amsterdam, The Netherlands, Tech. Rep. BN 9/71, 1971. [5] M. Arcak and E. Sontag, “A passivity-based stability criterion for a class of interconnected systems and applications to biochemical reaction networks,” Math. Biosci. Eng., vol. 5, pp. 1–19, 2008. [6] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Phys. Rep., vol. 469, pp. 93–153, 2008. [7] T. Balch and R. C. Arkin, “Behavior-based formation control for multirobot teams,” IEEE Trans. Robot. Automat., vol. 14, pp. 926–939, 1998. [8] A. L. Barabàsi and L. Albert, “Emergence of scaling in random networks,” Science, vol. 286, pp. 509–512, 1999. [9] M. Barahona and L. M. Pecora, “Synchronization in small-world systems,” Phys. Rev. Lett., vol. 89, p. 054101, 2002. [10] R. W. Beard, A. W. Beard, J. Lawton, and F. Y. Hadaegh, “A coordination architecture for spacecraft formation control,” IEEE Trans. Contr. Syst. Technol., vol. 9, pp. 777–790, 1999. [11] I. Belykh, V. Belykh, and M. Hasler, “Synchronization in complex networks with blinking interactions,” in Proc. 2005 Int. Conf. Physics and Control, 2005, pp. 86–91. [12] I. V. Belykh, V. N. Belykh, and M. Hasler, “Blinking model and synchronization in small-world networks with a time-varying coupling,” Physica D, vol. 195, pp. 188–206, 2004. [13] V. N. Belykh, I. V. Belykh, and M. Hasler, “Connection graph stability method for synchronize coupled chaotic systems,” Physica D, vol. 195, pp. 159–187, 2004. [14] V. N. Belykh, N. N. Verichev, L. J. Kocarev, and L. O. Chua, “On chaotic synchronization in a linear array of Chua’s circuits,” J. Circuits Syst. Comput., vol. 3, pp. 579–589, 1993. [15] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, and C. S. Zhou, “The synchronization of chaotic systems,” Phys. Rep., vol. 366, pp. 1–101, 2002. [16] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. U. Hwang, “Complex networks: Structure and dynamics,” Phys. Rep., vol. 424, pp. 175–308, 2006. [17] S. Boccaletti and D. L. Valladares, “Characterization of intermittent lag synchronization,” Phys. Rev. E, vol. 62, pp. 7497–7500, 2000. [18] S. Boyd, “Convex optimization of graph laplacian eigenvalues,” in Proc. Int. Congr. Mathematicians, 2006, pp. 1311–1319. [19] J. Buck, “Synchronous rhythmic flashing of fireflies II,” Quart. Rev. Biol., vol. 63, pp. 265–289, 1988. [20] T. Chen and X. Liu. Network synchronization with an adaptive coupling strength. 2006. [21] T. Chen, X. Liu, and W. Lu, “Pinning complex networks by a single controller,” IEEE Trans. Circuits Syst. I, vol. 54, pp. 1317–1326, 2007. [22] L. O. Chua, M. Itoh, L. Kocarev, and K. Eckert, “Chaos synchronization in Chua’s circuit,” Univ. California, Berkeley, Tech. Rep., 1992. [23] L. O. Chua, L. Kocarev, K. Eckert, and M. Itoh, “Experimental chaos synchronization in Chua’s circuit,” Int. J. Bifurcat. Chaos, vol. 2, pp. 705– 708, 1992. [24] I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin, “Effective leadership and decision-making in animal groups on the move,” Nature, vol. 433, pp. 513–516, 2005. [25] A. Das, R. Fierro, V. Kumar, J. Ostrowski, J. Spletzer, and C. Taylor, “A vision-based formation control framework,” IEEE Trans. Robot. Automat., vol. 18, pp. 813–825, 2002. [26] S. Dashkovskiy, B. Rüffer, and F. Wirth, “An ISS small-gain theorem for general networks,” Math. Control Signal Syst., vol. 19, pp. 93–122, 2007. THIRD QUARTER 2010

[27] P. DeLellis, M. diBernardo, and F. Garofalo, “Synchronization of complex networks through local adaptive coupling,” Chaos, vol. 18, p. 037110, 2008. [28] P. DeLellis, M. diBernardo, and F. Garofalo, “Novel decentralized adaptive strategies for the synchronization of complex networks,” Automatica, vol. 45, pp. 1312–1318, 2009. [29] P. DeLellis, M. diBernardo, F. Garofalo, and M. Porfiri, “Evolution of complex networks via edge snapping,” IEEE Trans. Circuits Syst. I, to be published. [30] P. DeLellis, M. diBernardo, and G. Russo, “On quad, lipschitz and contracting vector fields for consensus and synchronization of networks,” IEEE Trans. Circuits Syst., to be published. [31] P. DeLellis, M. diBernardo, F. Sorrentino, and A. Tierno, “Adaptive synchronization of complex networks,” Int. J. Comput. Math., vol. 85, pp. 1189–1218, 2008. [32] P. DeLellis, M. diBernardo, and L. F. Turci, “Pinning control of complex networked systems via a fully adaptive decentralized strategy,” IEEE Trans. Automat. Contr., submitted for publication. [33] M. diBernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk, Piecewise-Smooth Dynamical Systems. London: Springer-Verlag, 2008. [34] D. V. Dimarogonas and K. J. Kyriakopoulos, “On the rendezvous problem for multiple nonholonomic agents,” IEEE Trans. Automat. Contr., vol. 52, pp. 916–922, 2007. [35] L. Donetti, P. Hurtado, and M. Muñoz, “Entangled networks, synchronization, and optimal network topology,” Phys. Rev. Lett., vol. 95, p. 188701, 2005. [36] M. Egerstedt and H. Xiaoming, “Formation constrained multi-agent control,” IEEE Trans. Robot. Automat., vol. 17, pp. 947–951, 2001. [37] M. B. Elowitz and S. Leibler, “A synthetic oscillatory network of transcriptional regulators,” Nature, vol. 403, pp. 335–338, 2000. [38] J. H. Fewell, “Social insect networks,” Science, vol. 301, pp. 1867– 1870, 2003. [39] L. Freeman, “A set of measures of centrality based upon betweenness,” Sociometry, vol. 40, pp. 35–41, 1977. [40] J. Garcia-Ojalvo, M. B. Elowitz, and S. H. Strogatz, “Modeling a synthetic multicellular clock: Repressilators coupled by quorum sensing,” Proc. Nat. Acad. Sci., vol. 101, pp. 10955–10960, 2004. [41] M. Golubitsky and I. Stewart, “Nonlinear dynamics of networks: The groupoid formalism,” Bull. Amer. Math. Soc., vol. 43, p. 305, 2006. [42] T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, “Netevo: A computational framework for the evolution of complex dynamical networks,” to be published. [43] T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, “A dynamic approach to the evolution of complex networks,” in Proc. 19th Int. Symp. Mathematical Theory of Networks and Systems, 2010, pp. 1–5. [44] T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, “Evolving enhanced topologies for the synchronization of dynamical complex networks,” Phys. Rev. E, vol. 81, no. 5, p. 056212, 2010. [45] R. O. Grigoriev, M. C. Cross, and H. G. Schuster, “Pinning control of spatiotemporal chaos,” Phys. Rev. Lett., vol. 79, pp. 2795–2798, 1997. [46] W. Guo, F. Austin, S. Chen, and W. Chen, “Pinning synchronization of the complex networks with non-delayed and delayed coupling,” Phys. Lett. A, vol. 373, pp. 1565–1572, 2009. [47] J. H. Holland, Adaptation in Natural and Artificial Systems. Cambridge, MA: MIT Press, 1992. [48] J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proc. Nat. Acad. Sci. USA, vol. 79, pp. 2554–2558, 1982. [49] R. A. Horn and C. R. Johnson, Matrix Analysis. Cambridge, U.K.: Cambridge Univ. Press, 1985. [50] A. Kashiwagi, I. Urabe, K. Kaneko, and T. Yomo, “Adaptive response of a gene network to environmental changes by fitness-induced attractor selection,” PLoS One, vol. 1, p. e49, 2006. [51] N. Kashtan and U. Alon, “Spontaneous evolution of modularity and network motifs,” Proc. Nat. Acad. Sci. USA, vol. 102, p. 13733, 2005. [52] N. E. Leonard and E. Fiorelli, “Virtual leaders, artificial potentials and coordinated control of groups,” in Proc. 40th Conf. Decision and Control, 2001. [53] C. Li, L. Chen, and K. Aihara, “Synchronization of coupled nonidentical genetic oscillators,” Phys. Biol., vol. 3, pp. 37–44, 2006. [54] X. Li, X. Wang, and G. Chen, “Pinning a complex dynamical network to its equilibrium,” IEEE Trans. Circuits Syst. I, vol. 51, pp. 2074–2087, 2004. [55] Z. Li and G. Chen, “Global synchronization and asymptotic stability of complex dynamical networks,” IEEE Trans. Circuits Syst. II, vol. 53, pp. 28–33, 2006. IEEE CIRCUITS AND SYSTEMS MAGAZINE

81

[56] Z. Li, Z. Duan, and G. Chen, “Cost and effects of pinning control for network synchronization,” Chin. Phys., vol. 18, pp. 106–118, 2009. [57] Z. Li, L. Jiao, and J.-J. Lee, “Robust adaptive global synchronization of complex dynamical networks by adjusting time-varying coupling strength,” Phys. A: Stat. Mech. Applicat., vol. 387, pp. 1369–1380, 2008. [58] J. Lin, A. Morse, and B. Anderson, “The multi-agent rendezvous problem,” in Proc. 42nd IEEE Conf. Decision and Control, 2003. [59] X. Liu and T. Chen, “Exponential synchronization of nonlinear coupled dynamical networks with a delayed coupling,” Physica A, vol. 381, pp. 82–92, 2007. [60] X. Liu and T. Chen, “Boundedness and synchronization of y-coupled lorenz systems with or without controller,” Physica D, vol. 237, pp. 630–639, 2008. [61] X. Liu and T. Chen, “Synchronization analysis for nonlinearly-coupled complex networks with an asymmetrical coupling matrix,” Physica A, vol. 387, pp. 4429–4439, 2008. [62] W. Lohmiller and J. Slotine, “On contraction analysis for non-linear systems,” Automatica, vol. 34, pp. 683–696, 1998. [63] E. N. Lorenz, “Deterministic nonperiodic flow,” J. Atmos. Sci., vol. 20, pp. 130–141, 1963. [64] J. Lu and D. W. Ho, “Local and global synchronization in general complex dynamical networks with delay coupling,” Chaos Solitons Fractals, vol. 37, pp. 1497–1510, 2008. [65] X. B. Lu and B. Z. Qin, “Adaptive cluster synchronization in complex dynamical networks,” Phys. Lett. A, vol. 373, pp. 3650–3658, 2009. [66] M. Mesbahi and F. Y. Hadaegh, “Formation flying control of multiple spacecraft via graphs, matrix inequalities, and switching,” J. Guid. Control Dyn., vol. 24, pp. 369–377, 2001. [67] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon, “Network motifs: Simple building blocks of complex networks,” Science, vol. 298, pp. 824–827, 2002. [68] C. C. Moallemi and B. van Roy, “Distributed optimization in adaptive networks,” Adv. Neural Inform. Process. Syst., vol. 16, 2004. [69] P. Moylan and D. Hill, “Stability criteria for large-scale systems,” IEEE Trans. Automat. Contr., vol. 23, pp. 143–149, 1978. [70] R. Olfati-Saber and R. Murray, “Consensus problems in networks of agents with switching topology and time-delays,” IEEE Trans. Automat. Contr., vol. 49, pp. 1520–1533, 2004. [71] G. V. Osipov, J. Kurths, and C. Zhou, Synchronization in Oscillatory Networks. Berlin, Germany: Springer-Verlag, 2007. [72] E. Ott, Chaos in Dynamical Systems. Cambridge, U.K.: Cambridge Univ. Press, 1993. [73] M. Parter, N. Kashtan, and U. Alon, “Facilitated variation: How evolution learns from past environments to generalize to new environments,” PLoS Comput. Biol., vol. 4, p. e1000206, 2008. [74] L. M. Pecora and T. L. Carroll, “Synchronization in chaotic systems,” Phys. Rev. Lett., vol. 64, pp. 821–824, 1990. [75] L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Phys. Rev. Lett., vol. 80, pp. 2109–2112, 1998. [76] Q. C. Pham and J. J. E. Slotine, “Stable concurrent synchronization in dynamic system networks,” Neural Netw., vol. 20, pp. 62–77, 2007. [77] M. Porfiri and M. di Bernardo, “Criteria for global pinning-controllability of complex networks,” Automatica, vol. 44, pp. 3100–3106, 2008. [78] M. Porfiri and F. Fiorilli, “Node-to-node pinning control of complex networks,” Chaos, vol. 19, p. 013122, 2009. [79] M. Porfiri, D. Stilwell, and E. Bollt, “Synchronization in random weighted directed networks,” IEEE Trans. Circuit Syst. I, vol. 55, pp. 3170–3177, 2008. [80] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “Phase synchronization of chaotic oscillators,” Phys. Rev. Lett., vol. 76, pp. 1804–1807, 1996. [81] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “From phase to lag synchronization in coupled chaotic oscillators,” Phys. Rev. Lett., vol. 78, pp. 4193–4196, 1997. [82] G. Russo and M. di Bernardo, “An algorithm for the construction of synthetic self synchronizing biological circuits,” in Proc. Int. Symp. Circuits and Systems, 2009, pp. 305–308. [83] G. Russo and M. di Bernardo, “Contraction theory and the master stability function: Linking two approaches to study synchronization in complex networks,” IEEE Trans. Circuit Syst. II, vol. 56, pp. 177–181, 2009. [84] G. Russo and M. di Bernardo, “How to synchronize biological clocks,” J. Comput. Biol., vol. 16, pp. 379–393, 2009. [85] G. Russo and M. di Bernardo, “Solving the rendezvous problem for multi-agent systems using contraction theory,” in Proc. Int. Conf. Decision and Control, 2009.

82

IEEE CIRCUITS AND SYSTEMS MAGAZINE

[86] G. Russo, M. di Bernardo, and J. J. E. Slotine, “A graphical algorithm to prove contraction of nonlinear circuits and systems,” IEEE Trans. Circuits Syst. I, submitted for publication. [87] G. Russo, M. di Bernardo, and J. J. E. Slotine, “An algorithm to prove contraction, consensus, and network synchronization,” in Proc. 1st IFAC Workshop on Estimation and Control of Networked Systems, 2009. [88] G. Russo, M. di Bernardo, and E. Sontag, “Stability of networked systems: A multi-scale approach using contraction,” in Proc. IEEE Conf. Decision and Control, submitted for publication. [89] G. Russo, M. di Bernardo, and E. Sontag, “Global entrainment of transcriptional systems to periodic inputs,” PLoS Comput. Biol., vol. 6, p. e1000739, 2010. [90] G. Russo and J. Slotine, “Global convergence of quorum-sensing networks,” J. R. Soc. Interface, submitted for publication. [91] D. D. Siljak, Large-Scale Dynamic Systems: Stability and Structure. New York: North Holland, 1978. [92] Q. Song and J. Cao, “On pinning synchronization of directed and undirected complex dynamical networks,” IEEE Trans. Circuits Syst. I, to be published. [93] E. Sontag, “Input to state stability: Basic concepts and results,” in Nonlinear and Optimal Control Theory, P. Nistri and G. Stefani, Eds. Berlin: Springer-Verlag, 2007, pp. 163–220. [94] E. Sontag, “Monotone and near-monotone biochemical networks,” Syst. Synth. Biol., vol. 1, pp. 59–87, 2007. [95] E. D. Sontag, Mathematical Control Theory: Deterministic Finite-Dimensional Systems. New York: Springer-Verlag, 1998. [96] F. Sorrentino, M. diBernardo, F. Garofalo, and G. Chen, “Controllability of complex networks via pinning,” Phys. Rev. E, Stat. Nonlin. Soft Matter Phys., vol. 75, p. 046103, 2007. [97] F. Sorrentino and E. Ott, “Adaptive synchronization of dynamics on evolving complex networks,” Phys. Rev. Lett., vol. 100, p. 114101, 2008. [98] K. O. Stanley, B. D. Bryant, and R. Mikkulainen, “Evolving adaptive neural networks with and without adaptive synapses,” in Proc. 2003 Congr. Evolutionary Computation (CEC’03). [99] J. Sun, E. Bollt, and T. Nishikawa, “Master stability function for coupled nearly identical dynamical systems,” Europhys. Lett., vol. 85, p. 60011, 2009. [100] D. Syliak, “Dynamic graphs,” Nonlin. Anal.: Hybrid Syst., vol. 2, pp. 544–567, 2008. [101] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, “Stale flocking of mobile agents, part II: Dynamic topology,” in Proc. 42nd Conf. Decision and Control, 2003. [102] M. Vidyasagar, Input-Output Analysis of Large Scale Interconnected Systems. New York, Berlin: Springer-Verlag, 1981. [103] L. Wang, H. P. Dai, H. Dong, Y. Y. Cao, and Y. X. Sun, “Adaptive synchronization of weighted complex dynamical networks through pinning,” Eur. Phys. J. B, vol. 61, pp. 335–342, 2008. [104] W. Wang and J. J. E. Slotine, “On partial contraction analysis for coupled nonlinear oscillators,” Biol. Cybern., vol. 92, pp. 38–53, 2005. [105] X. F. Wang and G. Chen, “Synchronization in scale-free dynamical networks: Robustness and fragility,” IEEE Trans. Circuits Syst. I, vol. 49, pp. 54–62, 2002. [106] X. F. Wang and G. Chen, “Synchronization in small-world dynamical networks,” Int. J. Bifurcat. Chaos, vol. 12, pp. 187–192, 2002. [107] C. W. Wu, “Localization of effective pinning control in complex networks of dynamical systems,” in Proc. IEEE Int. Symp. Circuits and Systems, 2008, pp. 2530–2533. [108] C. W. Wu and G. Chen, “Synchronization in an array of linearly coupled dynamical systems,” IEEE Trans. Circuits Syst., vol. 42, pp. 430–447, 1995. [109] W. Wu, W. Zhou, and T. Chen, “Cluster synchronization of linearly coupled complex networks under pinning control,” IEEE Trans. Circuits Syst. I, vol. 56, pp. 829–839, 2009. [110] S.-H. Yook and H. Meyer-Ortmanns, “Synchronization of rössler oscillators on scale-free topologies,” Phys. A: Statist. Theor. Phys., vol. 371, pp. 781–789, 2006. [111] W. Yu, P. DeLellis, G. Chen, M. diBernardo, and J. Kurths, “Distributed adaptive control of synchronization in complex networks,” IEEE Trans. Automat. Cont., submitted for publication. [112] Q. Zhang, J. Lu, J. Lu, and C. Tse, “Adaptive feedback synchronization of a general complex dynamical network with delayed nodes,” IEEE Trans. Circuits Syst. II, vol. 55, pp. 183–187, 2008. [113] C. Zhou and J. Kurths, “Dynamical weights and enhanced synchronization in adaptive complex networks,” Phys. Rev. Lett., vol. 96, p. 164102, 2006. THIRD QUARTER 2010

XIAOFAN LIU

Synchronization and Control of Complex Networks via Contraction, Adaptation and Evolution Pietro DeLellis, Mario di Bernardo, Thomas E. Gorochowski, and Giovanni Russo

Digital Object Identifier 10.1109/MCAS.2010.937884

64

IEEE CIRCUITS AND SYSTEMS MAGAZINE

1. Introduction omplex networked systems abound in Nature and Technology. They consist of a multitude of interacting agents communicating with each other over a web of complex interconnections. Flocks of birds, platoon of cooperating robots, swirling fishes in the Ocean are all examples whose intricate dynamics can be modeled in terms of three essential ingredients: (i) a mathematical description of the dynamical behavior of each of the agents in the network; (ii) an interaction (or coupling) protocol used by agents to communicate with each other and (iii) a graph describing the network of interconnections between neighboring agents. These three elements are actually mapped onto the mathematical model usually considered in the literature to describe a complex network which uses appropriate equations to describe the node dynamics, the coupling protocol and the

C

1531-636X/10/$26.00©2010 IEEE

THIRD QUARTER 2010

network topology (see Sec. 2.2; e.g., [16] and references therein for further details). The most striking feature of complex networked systems observed in Nature and Technology is their ability to show emergent behavior that cannot be explained in terms of the individual dynamics of each single agent in the network. Synchronization and coordination of motion are possibly the simplest and most frequent of these types of phenomena. Think for example of the incredible patterns formed by flocks of birds or schools of fishes [24], the synchronous flashing of fireflies in Amazonia [19] and the many other situations where the emergence of coordinated behavior has been observed. Over the past few years, much research effort has been focused on understanding and explaining the onset of synchronization in complex networks, with the aim of uncovering the very mechanism causing its occurrence. Also, the dependence has been studied of synchronization on specific features of the network such as the form of the coupling strategy or the topology of the network. It has been found that synchronization emerges even in the presence of simple, diffusive coupling where neighboring nodes adjust their dynamics proportionally to the mismatch between some output function of their states [6, 14, 15, 22, 23, 53, 55, 71, 74, 80, 108]. Also, it was observed that synchronization depends on the network topology so that it is simpler to synchronize networks with certain topological features. Synchronizability of small world, scale-free, disassortative or assortative networks was studied showing that the topology has indeed an influence on the dynamics of the network (see for example [9, 79, 105, 106, 110]). In all cases, when synchronization occurs all nodes tend towards the same asymptotic evolution which is unknown a priori and emerges from the negotiation process taking place on the network between neighboring nodes. A different problem is when the aim is to control the network onto some desired asymptotic solution. In this case, a possible solution would be to control every node in the network so as to steer the dynamics of each of them towards the desired evolution. In practice, this is not viable as the number of nodes in the network is typically very large and the amount of control effort available bounded. Pinning control was introduced as a viable alternative, where only a (small) fraction of the network nodes are directly controlled by means of some external input, with the control action being propagated to the rest of the nodes through their interconnections [21, 45, 46, 56, 54, 77, 78, 92, 96, 107, 109]. Controllabil-

ity of networks was also defined in [96] and further investigated in [77]. One of the key features of most of the techniques for the control and synchronization of networks presented so far in the literature is the time-invariant nature of the coupling strength between nodes and the topology of their interconnections. Specifically, it is often supposed that the coupling gain determining the strength of the interaction between nodes is constant as is the adjacency matrix describing the topology of the interconnections between them. This is not in agreement with much evidence from the natural world, showing that instead in such examples as flocks of birds, swarms of fireflies or schools of fishes, individuals in the network are often able to form or suppress interconnections between themselves and adjust (or adapt) the strength and structure of their coupling. Also, often, rather than being simply diffusive, the coupling functions between agents in the network are typically nonlinear with a functional form that has been engineered by Nature to best fulfill a certain need. Over the past few years, strategies to cope with the case of time varying topologies have been proposed but often with the aim of showing that synchronization can be achieved even in the presence of variations in the network topology, e.g. [11, 12]. Adaptation and evolution in Nature are instead desired phenomena, engineered over the years in order for complex networks of dynamical agents to adapt the strength of their interactions and evolve an appropriate network of interconnections so as to perform some function of interest. The aim of our work is to attempt at mimicking these features of natural networks and propose strategies for the adaptation and evolution of complex networks of dynamical systems that guarantee the emergence of some asymptotic behavior of interest, namely a synchronous solution. Such strategies will be, in general, local and decentralized so that it is up to pairs or groups of agents to self-determine the strength and topology of their interconnections. We start with the generalization of the concept of dynamic graph first stated by Šiljak in [91] by incorporating some of the features of a complex adaptive system as defined by Holland in [47]. Then, we propose two alternative strategies to design the coupling functions between nodes. One is based on finding appropriate adaptation laws to evolve the coupling gain between neighboring nodes, and the other on the ad hoc design of a nonlinear coupling protocol based on contraction theory [62]. We will discuss how in both cases it is possible to prove convergence of all network nodes onto a

Pietro DeLellis, Mario di Bernardo, and Giovanni Russo are with the Department of Systems and Computer Engineering, University of Naples Federico II, Italy. Mario di Bernardo and Thomas E. Gorochowski are with the Bristol Centre for Complexity Sciences, University of Bristol, U.K.

THIRD QUARTER 2010

IEEE CIRCUITS AND SYSTEMS MAGAZINE

65

common solution and briefly present a set of representative examples. We then move to the case where the topology of the interconnections between nodes in the network evolves. Firstly, we will consider a novel strategy, termed as edgesnapping in [29]. Specifically, we will make neighboring nodes able to decide whether to activate or suppress a link between themselves according to the strength of their mismatch. We will show that in this manner networks with a certain final topology emerge with interesting features that will be discussed in Sec. 4. Next, we take a different approach and consider the case where evolution of the network must be such as to optimize some performance index defined on the network and its dynamics. A novel computational tool, NetEvo, will be presented and used to evolve networks in order to maximize their synchronizability, i.e. the range of values of the coupling gain guaranteeing synchronization. We will compare the topologies of the networks evolved using the Laplacian eigenratio as a static measure of the network synchronizability (e.g. [9, 75]) with those obtained by considering dynamical synchronizability measures. We will observe that in both cases networks can present some striking features in terms of their topologies and focus on the presence of certain characteristics such as the frequencies of motifs [51, 73] and other macroscopic observables. The possibility of having topological bifurcations will also be discussed. Finally, we discuss some of those that we believe are the most pressing open problems and challenges in the area of complex dynamical networks suggesting that much work is still needed to achieve the ultimate goal of engineering fully adaptive, self-evolving networks able to show the incredible features, we are able to marvel about in all the natural complex systems surrounding us. 2. Background 2.1. Describing a Complex Network Complex networks provide a flexible tool for modeling many real-world systems and consist of three main attributes: topology, dynamics and an evolutionary process. In order to fully define a complex network and study the influence each of these attributes has, a common mathematical framework is required. This is to not only to provide a means to analyzing various properties a particular system possesses, but also to act as a shared language in which to describe the systems themselves. Complex networks are a purposefully general concept and so being able to relate problems from different fields provides an important first step to understanding underlying principles they may share. Several attempts have 66

IEEE CIRCUITS AND SYSTEMS MAGAZINE

been made with this in mind, including coupled cell networks, dynamic graphs and complex adaptive systems. Coupled cell networks were introduced by Golubitsky and Stewart [41] as a means to investigating how network topology and system-level dynamics are linked. These networks consist of a number of cells connected by a fixed set of edges. Equivalence relations are defined over each of these sets, assigning a type to every element. These relations can represent the form of dynamical system at each cell and the type of coupling across an edge, although the exact forms are not necessarily required. By studying the structure of these systems with particular dynamics, it has been shown that certain symmetries within the network topology and equivalence relations can predict robust synchronous states and other related phenomena. Although this formalism can accommodate any form of fixed topology and dynamics, it unfortunately does not facilitate evolution of the structure itself. In contrast, the dynamic graphs of Šiljak [100] focus on how the structure of a network evolves over time. A weighted graph of a fixed number of nodes is placed within a linear vector space, and a dynamic graph then defined as a one-parameter group of transformations of this space into itself. Put more simply, a dynamic graph can be viewed as one where edge weights vary in time due to some differential equation. A major benefit of this formalism is that much of the existing dynamical systems theory can be directly applied using adjusted forms of stability. This approach allows for the evolution of edge structure to be incorporated into a model, but does have limitations. Edge states can only take a single state value (their weight) and other aspects of evolution, such as growth, are not possible. To allow for systems that structurally evolve and grow we turn to the work of Holland and his concept of a complex adaptive system [47]. This framework was developed during the 1970’s in an attempt to understand the process of adaptation. It comprises of four main elements: 1. A – A set of possible structures; 2. V – A set of operators for modifying structures; 3. I – A set of inputs from the environment; 4. t : I 3 A S V – The adaptive plan. Evolution takes place in two stages, first the current structure and input from the environment are mapped to an operator v [ V using the adaptive plan t. Then, the chosen operator is applied to return the evolved structure. It should be noted that both the adaptive plan and operators for modifying structures can be probabilistic, leading to a stochastic evolution of the structure. The main advantage of this framework is its flexibility, allowing for the set of admissible structures, operators THIRD QUARTER 2010

and adaptive plan to be defined any way we choose. This provides endless possibilities, but unfortunately gives no indication as to how network topology and dynamics can be best incorporated. As we have seen, several attempts have been made to develop a framework for the study of complex networks. In each case, trade-offs were made that resulted in none fully addressing the three attributes discussed earlier. With this in mind and taking inspiration from these previous ideas, we now introduce the concept of an evolving dynamical network (EDN) which aims to meet this goal. We begin by first introducing the concept of a generalized dynamic graph. This has the role of providing a structure in which network topology and dynamics can be fully defined. Taking a directed graph G 5 1 V, E 2 of a fixed number of nodes and edges, we associate a state to every node and edge 1 V, E 2 . This state can be of any form and dimension required. Furthermore, these states embody the dynamics of the system, with their motion described by a mapping F which is also defined over some set of times T. Given the current set of states and any external inputs U, this mapping returns the new state of the system. It is worth noting that we consider network dynamics to cover any aspects related to changes in system states. To incorporate evolution of the structure and dynamics, i.e. G and F, we consider the set of all possible generalized dynamic graphs D. We see evolution of a system as exploration of this space, taking place either online, as the network dynamics occur, or off-line, having the network state reset after each evolutionary step. To embody this feature we view an evolving dynamical network as a complex adaptive system where the set of admissible structures is D. We consider a set of structural operators V that can alter a given structure, a set of inputs I that can be from the environment or for control purposes, and the evolutionary plan t, which given an input and structure decides the structural operator to apply. With these components we are now in the position to fully describe the three main attribute of a complex network. Formal definitions and an illustration of the connection between the two aspects of an EDN, network dynamics and evolution, can be found in Box 1. In following sections EDNs will hold a central role, helping us relate similarities and differences between various approaches to synchronization.

collective behaviors, like synchronization and consensus, in networks of interconnected dynamical systems. This class of network, often named complex networks, is an instance of a type of EDN described in Box 1. In fact, the classic complex network model can be viewed as a generalized dynamic graph in which: ■ V 5 5 x1, x2, c, xn 6 , where xi [ Rn, 4i : vi [ V. The network nodes are dynamical systems, and therefore the nodes assignments coincide with the states xi of the nodes belonging to the real coordinate space Rn. ■ E 5 5 eij 6 , 4 ij : 1 vi, vj 2 [ E, with eij 5 s [ R. The constant s represents the overall coupling strength between the network nodes. ■ U 5 [. There are no external inputs. ■ F 5 V 3 E 3 T S V 3 E can be split into two operators Fv and Fe, where Fv 5 V 3 E 3 T S V is a continuous differential operator that describes the nodes’ evolution, while Fe 5 E S E is the identity operator describing the constant edge states. In terms of ODEs, the evolution of the nodes’ states is described by the following set of differential equations

2.2. Problem Statements: Synchronization and Control As explained in Section 1, scientists from different communities, including physicists, biologists, control engineers and mathematicians, are making an effort to understand the dynamics leading to the emergence of

where degi is the cardinality of the set Ni, that is, the N set of incident nodes in i. Notice that ,ii 5 2 a j 51 ,ij, 4i 5 1, c, N. Throughout the paper, we typically consider undirected networks, that is, networks in which if 1 vi, vj 2 [ E, then also 1 vj, vi 2 [ E. Therefore, L is symmetric and zero-column sum.

THIRD QUARTER 2010

# xi 5 f 1 xi, t 2 1 ci 1 x, t 2 ,

(1)

with i 5 1, c, N. In (1) the function f : Rn 3 R 1 S Rn is the intrinsic dynamics of the i-th node, while the function ci : RnN 3 R 1 S Rn, termed as interaction or coupling function, describes the interaction of the i-th node with the other nodes composing the interconnected system. Typically, ci depends only on the dynamics of those nodes directly linked to node i selected by the Laplacian matrix describing the network topology, see e.g. [49]. Specifically, often the coupling protocol is chosen in (1) as ci 1 x, t 2 5 2s a ,ij h 1 xj 2 , N

i [ V.

(2)

j51

Here, h : Rn S Rn is the output function through which the network nodes are coupled, ,ij is the ij th element of the zero-row sum Laplacian matrix L associated to the graph G 5 5 V, E 6 describing the network topology, defined as: degi ,ij 5 • 21 0

i5j 1 vi, vj 2 [ E, otherwise

(3)

IEEE CIRCUITS AND SYSTEMS MAGAZINE

67

BOX 1. EVOLVING DYNAMICAL NETWORKS

I

llustration of how topology, dynamics and evolution are Network Dynamics

integrated within an evolving dynamical network (EDN). It

can be seen that there are two main dimensions, with the op-

a vital step towards more general theories spanning multiple application domains. Definition 4 (Generalized Dynamic Graph) We define a generalized dynamic graph for a fixed number of nodes N, network structure E, and opera-

tor F, as the collection D 5 1 V, E, V, E, U, T, F 2 [ D

Network Evolution

Having a framework like this for studying complex networks is

Φ

D1

erators F and v [ V determining which direction to take.

ω1

Update Structure, Rewire an Edge

Φ

D2

Φ

Φ

Grow Structure, Add Node, and Connecting Edges

ω2 D3

Φ

Φ

where,

• V 5 5 v1, v2, c, vN 6 — finite set of nodes.

• E 5 5 e1, e2, c, eM 6 where ei [ V 3 V — finite

Definition 5 (Evolving Dynamical Network) Using generalized dynamic graphs as an underlying

set of edges.

• V 5 5 v1, v2, c, vN 6 where vi [ Vi — set of

structure, we define an evolving dynamical network (EDN) as the collection E^ 5 1 D, V, I, t 2 [ E^ where, • D 5 5 D1, D2, c 6 — set of admissible gener-

node state assignments.

• E 5 5 e1, e2, c, eM 6 where ei [ Ei — set of

alized dynamic graphs.

• V 5 5 v 1, v 2, c 6 — set of structural operators such that v [ V is a mapping v : D S P

edge state assignments.

• U 5 5 u1, u2, c, uP 6 where ui [ Ui — set of external input assignments.

where P is a probability distribution over D.

• T – set of times relating to dynamics.

• I — inputs from the environment or for

• F : V 3 E 3 U 3 T S V 3 E — mapping from state to state for nodes, edges and control input where the mapping F can take any form,

determines the operator v [ V to be used

e.g. deterministic, probabilistic, etc.

when transitioning to a new structure.

Now we are ready to give a formal definition of the concept of synchronization in a complex network. Definition 1 Network (1) is said to achieve asymptotic synchronization if and only if lim 1 xi 1 t 2 2xj 1 t 2 2 5 0,

t S 1`

i, j : 1 vi, vj 2 [ E.

(4)

Asymptotic synchronization is the most intuitive case of synchronization. However, other forms of synchronization are also possible, e.g. lag synchronization as defined in [17, 81]. The emergence of the synchronous behavior depends on both the individual nodes’ dynamics and the structure of the interconnections. Different analytical tools have been used to give conditions for synchronization, such as Lyapunov-based techniques, the Master Stability Func68

control purposes. • t : I 3 D S V — evolutionary plan which

IEEE CIRCUITS AND SYSTEMS MAGAZINE

tion approach and contraction theory (see Box 2 for further details). In many real networks, synchronization is not the only goal that one may want to achieve. In formation control [7, 10, 25, 36, 66], for instance, the aim is not only to synchronize the motion of the set of nodes (the agents) of the network, but we want to tame the dynamics of the whole network to a desired reference trajectory. To solve this kind of problem, it is necessary to give external inputs ui 1 x 2 [ Rm, m # n to the network nodes. In the EDN framework, the only difference is that the set U is not anymore empty and the domain of the operator F is V 3 E 3 U 3 T. A possible solution to this control problem is to add a controller to each of the network nodes. This approach is not feasible when the aim is to control large networks. This is why in the recent literature the so-called pinning control scheme for controlling complex networks has been proposed [21, THIRD QUARTER 2010

45, 54, 77, 107, 109]. In this scheme, an extra-node vs is added in the network, the pinner. This extra-node is identical to the other nodes of the network and defines the desired trajectory by being only unidirectionally coupled to the rest of the network. Only a small subset of nodes, say Npin V N, are directly controlled by the pinner. These nodes are called pinned nodes. The equation describing the closed-loop network are the following: # xi 5f 1 xi, t 2 2 s a ,ijh 1 xj 2 2 diqi 1 xi 2 xs 2 , N

# xs 5 f 1 xs, t 2 di 5 e

i 5 1, c, N (5)

i51

(6)

1, for i 5 1, c, Npin, 0, for i 5 1 Npin 1 1 2 , c, N,

(7)

Definition 2 Network (5)–(7) is said to achieve controlled network synchronization if and only if lim 1 xi 1 t 2 2 xs 1 t 2 2 5 0,

i [ V.

(8)

To study the local controllability of the network nodes to the desired trajectory, the concept of pinning controllability was introduced in [96], applying the MSF approach (see box 2) to an extended network, which includes the external pinner. With this approach, theoretical bounds on the control gains qi and on the coupling gain s were given. In the following sections, after highlighting the limitations of the previous models (2) and (5)–(7), we will propose novel strategies for designing the nodes’ coupling (Section 3) and for evolving the network topology (Section 4). 3. Designing Novel Coupling Strategies 3.1. Synchronization via Decentralized Gain Adaptation Many real phenomena are characterized by the presence of adaptive mechanisms; for example, wireless networks of sensors that gather and communicate data to a central base station [68], or networks of robots when environmental conditions change unexpectedly (i.e. a robot loses a sensor) [98]. Moreover, examples of adaptive networks can be found in biology and the natural world [24, 38, 50]. In all these cases, it is realistic to assume that the strength of the interactions among nodes, characterized mathematically by s, is not identical for every node and timeinvariant. Real-world networks are often characterized instead by evolving, adapting couplings which vary in time THIRD QUARTER 2010

N # xi 5 f 1 xi, t 2 2 a lij 1 t 2 h 1 xj 2 ,

i 5 1, c, N

(9)

1 vi, vj 2 [ E,

(10)

j51

where xs denotes the desired trajectory defined by the pinner, ui 5 2 diqi 1 xi 2 xs 2 is the control input, with qi [ R being the i th control gain.

t S 1`

according to different environmental conditions. Adaptive gain strategies have been recently presented in a number of papers, see for example, [20, 27, 28, 31, 57, 97, 103, 112, 113]. In these strategies, the network gains evolve according to integro-differential adaptation laws driven by global or local synchronization errors. For instance, a common adaptive gain for all the network nodes is studied in [20, 57] while decentralized gain adaptation is proposed, using different gains for each node in [113] or edge in [28]. Specifically, the coupling protocol is chosen as a timevarying function so that the network model becomes:

# sij 5 g 1 t, sij, xi, xj 2

where sij [ R is the time-varying coupling gain associated to edge 1 vi, vj 2 , and lij 1 t 2 is ij th element of the timevarying laplacian matrix L 1 t 2 , defined as 2sij 1 t 2 lij 1 t 2 5 µ a sik k:vk [Ni

0

1 vi, vj 2 [ E i5j otherwise,

(11)

where Ni ( V is the set of neighbors of node vi, that is, the set of nodes incident to vi. In designing the adaptive approach, a key point is the appropriate selection of the function g. The design of the adaptive law is clearly dependent on the nature of the individual system, described by f, and on the output function h. In what follows, we assume that the output function is h 1 z 2 5 Gz, where G is a non-null diagonal positive semidefinite matrix, commonly named as the inner coupling matrix. This matrix defines the subset of state variables through which the network nodes are coupled. As for the node dynamics, we give the following definition: Definition 3 We say that the vector field f 1 z, t 2 or, equivalently, the dynam# ical system z 5 f 1 z, t 2 is QUAD 1 D, v 2 if and only if 1 x2y 2 T 1 f 1 x, t 2 2 f 1 y, t 22 2 1 x2y 2 T D 1 x2y 2

# 2v 1 x2y 2 T 1 x2y 2 , 4x, y [ Rn,

(12)

where D is an arbitrary diagonal matrix and v is a positive scalar. This condition gives an upper bound on the rate of variation of the vector field f. This is an assumption which is often assumed to prove asymptotic synchronization in complex networks, see for instance [13, 27, 32, 60, 61, 64]. Notice that the assumption is quite mild and is satisfied by all linear systems and by many chaotic systems, as for instance the Chua’s circuit and the Lorenz oscillator IEEE CIRCUITS AND SYSTEMS MAGAZINE

69

[28, 111]. The relationship between QUAD assumption and other assumptions usually made to prove synchronization, such as the well-known Lipschitz condition and the contraction theory, are expounded in [30]. Under the assumption that the QUAD condition is satisfied, analytical conditions guaranteeing asymptotic synchronization and convergence of the coupling gains were derived using a Lyapunov-based approach. (Here, we give some of the more significant results, referring the reader to [28, 111] for the proofs.) Let us assume that g depends only on the difference between the states of the endpoints of the corresponding edges, that is g 1 t, sij, xi, xj 2 5 k 1 0 xj 2 xi 0 2 ,

(13)

where k is a function belonging to one of the following classes: 1. Class 1. k 1 z 2 5 azp, with 0 , p # 2, a . 0. 2. Class 2. k is a bounded k function1. Theorem 1 Assume that f is a continuously differentiable QUAD vector field, with D 2 vI # 0, and G 5 I. If the scalar function k belongs to class 1 or 2, then network (9), (10), (13) reaches asymptotic synchronization and each coupling gain sij converges to some finite value. It is worth emphasizing that, under the assumptions of the theorem, the network synchronizes for any connected topology. This is a quite strong result, but is limited to the case in which the network nodes are connected through all their states and the vector field f is required to belong to a subset of QUAD systems, satisfying D 2 vI # 0. In order to relax these assumption, in the following theorem a specific adaptive law is selected. Theorem 2 Assume that f is QUAD and G . 0. If we chose the adaptive law # (14) sij 5 1 xi 2 xj 2 TG 1 xi 2 xj 2 , then network (9) reaches asymptotic synchronization and each coupling gain sij converges to some finite value. Therefore, with this choice of g, we can show that any QUAD system, coupled on all its state variables through G, can be asymptotically synchronized through a decentralized adaptive modulation of the coupling gains. It is also possible to generalize the above Theorem and give sufficient conditions for asymptotic synchronization,

even when the coupling is only on a subset of the state variables. A fully decentralized pinning control strategy is presented in [32]. 3.1. Numerical Example In what follows, we consider, as a testbed example, a network of Lorenz systems [63, 72]. This nonlinear model is described by the set of the following three differential equations: # p 5 a1 1 q 2 p 2 # q 5 a2p 2 q 2 pr # r 5 pq 2 a3r,

(15) (16) (17)

where x 5 3 p, q, r 4 T is the state vector and a1, a2 and a3 are three positive parameters. In our simulation, we consider a network of N 5 500 Lorenz systems, coupled only on the first two state variables, that is 1 G 5 £0 0

0 1 0

0 0§. 0

The network topology is generated using the BarabàsiAlbert scale free model (we address the reader to [8] for further details). As for the node dynamics, we choose the parameters of the model guaranteeing a chaotic behavior. Namely, we set a1 5 10, a2 5 28 and a3 5 8/3. Moreover, the initial conditions on the states are taken from a normal distribution with 0-mean and standard deviation equal to 40. From trivial manipulations of the results achieved in [60], we have that the Lorenz system satisfies the assumptions. Therefore, we can select the adaptive law (14) to synchronize the network, with null initial conditions on the edge states, that is, sij 1 0 2 5 0. As we can see from Figure 1, synchronization is asymptotically achieved and all the coupling gains converge to finite steady state values. 3.2. Designing Protocols via Contraction Theory Typically, network coordination is reduced to a stability problem of some invariant set in the network phase space. The analysis/design of decentralized coordination strategies is then addressed by means of Lyapunov based techniques (Box 2), and by making some hypotheses on the vector field of each of the nodes (see e.g. the QUAD condition above). An alternative approach is that of using contraction theory, instead of Lyapunov tools, for proving network coordination (Box 2). The

A function F : I S R is positive definite if F 1x 2 . 0, 4x [ I, x 2 0 and F 10 2 5 0. A function f : R $0 S R $0 is of class k if it is continuous, positive defi nite and strictly increasing.

1

70

IEEE CIRCUITS AND SYSTEMS MAGAZINE

THIRD QUARTER 2010

xi (t )

50 0 −50 0

0.5

1

1.5 t (a)

2

2.5

0

0.5

1

1.5 t (b)

2

2.5

σij (t )

3 2 1 0

Figure 1. Network of 500 Lorenz systems coupled via edgebased strategy (14). (a) Evolution of the states and (b) the coupling gains.

main advantage of this approach lies in the fact that stability is defined as an incremental property of the trajectories themselves. Notice that the invariant set towards which they all converge is not known a priory. This, in turn, allows to separate the analysis/control problem in two well-defined steps: i) determine some invariant subspace for the network evolution; ii) prove/ impose contraction. The general principles used for the proofs of the results presented in what follows are based on two key ideas. We prove contraction of only some directions in phase space (the directions transverse to the synchronization manifold). Notice that contraction is a property of network trajectories, which is metric-dependent. To prove this directional contraction we make use of non-Euclidean norms. We illustrate the design of decentralized network coordination strategies and turn our attention to the problem of imposing some desired coordinated behavior for a set

where rix,y and vix,y represent the position and velocity vectors of agent i in a fixed reference frame. As in (1), the agents are coupled via a coupling function ci 1 x, t 2 5 3 0, 0, aix, aiy 4 T, with aix and aiy being acceleration inputs on the agent dynamics, chosen as: aix 5 s1 a 1 xj 2 xi 2 2 1 s1 1 s2 2 vix j[Ni

aiy 5 s3 a 1 yj 2 yi 2 2 1 s3 1 s4 2 viy .

1

0.4

0.8

0.2

0.6 viy

0

−0.2

2

2.5

3 (a)

3.5

4

0.4 0.2

−0.4

0

−0.6

−0.2

−0.8

(19)

j[Ni

0.6

vx

2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

of N . 1 nodes. The general principle that we use in our proofs can be summarized as follows. We differentiate the network dynamics so as to obtain network Jacobian, J. Now, let M be an invariant subspace for the network and define as V the orthonormal matrix spanning the null of M. Then, all network trajectories globally exponentially converge towards M if the matrix V J V T is contracting as defined in Box 2. An equivalent approach, which has no straightforward geometric interpretation, makes use of some virtual system, see [85] and [83] for further details. In Figure 2 an application to the rendezvous-problem is shown. Rendezvous is a typical example of coordination, that has also been used as a testbed problem in the literature. Namely, the problem is that of designing a distributed control strategy driving all the agents towards a common point in space, where they all have zero velocities (see for example [34, 58] and references therein). The rendezvous problem can be then formalized as the problem of finding an appropriate (distributed) control strategy which guides all the agents towards a common position, say r, i.e. such that limt S ` ri 1 t 2 5 r, 4i 5 1, c, N. In accordance with the existing literature [52, 101], we choose the nodes’ dynamics as: # # rix 5 vix, riy 5 viy, # # vix 5 aix, viy 5 aiy, (18)

0

5

10 Time (b)

15

20

−0.4

0

5

10 Time (c)

15

20

Figure 2. Designing communication protocols solving the rendezvous problem. (a) Evolution of agents’ positions in the x 2 y plane. (b) and (c) Shows time evolution of velocities vix, vij.

THIRD QUARTER 2010

IEEE CIRCUITS AND SYSTEMS MAGAZINE

71

2

1

6

5

3

4 (a) 5

xi (t )

0

−5

−10

0

5

10

15 Time (b)

20

25

30

Figure 3. (a) Network used to test the multi-scale approach and (b) temporal behavior of (20) with R 5 1, G 5 K 5 0.9.

We assume that the sensing region of each agent is bounded in space and is represented by a disk of radius R. In practice, this implies that each agent communicates only with those agents inside its sensing disk, i.e. the i-th agent communicates only with its neighbors, Ni. From the topological viewpoint, the assumption of limited sensing region makes the network of interest switched as the neighbors of the i-th agent can enter/exit its sensing disk. Clearly, the switches of network topology (due to the interactions among agents) occurs with a finite frequency. As shown in Figure 2, using contraction theory it is possible to show convergence of all trajectories towards each other (for more details see [86]). Sometimes, designing the coupling protocol alone is either not feasible (e.g. quorum sensing in biological networks) or insufficient to guarantee synchronization. In this case, it is possible to engineer the network topology and/ or the node dynamics to guarantee synchronization. Contraction theory can again be used as an invaluable tool for this purpose. The idea is to break down the analysis and design into two independent steps: (a) At a global level, properties of the network or interconnection graph are imposed so as to guarantee a desired behavior for the full interconnected system. 72

IEEE CIRCUITS AND SYSTEMS MAGAZINE

In this analysis, subsystems may be characterized as “black boxes” with assumed input-output characteristics, but detailed knowledge of their internal structure is not required. (b) At a local level of analysis, one imposes constraints on the structure and behavior of individual subsystems (components), so as to fit the requirements of the global approach. These requirements are verified independently of the overall network structure. This multi-scale or hierarchical methodology is robust in so far as a large degree of uncertainty can be tolerated in the components, only constrained by meeting appropriate behavioral requirements. There are many examples of such approaches in systems/control theory, including among others (1) the use of small-gain theorems to guarantee stability of a negative feedback loop provided that the components are individually stable (qualitative property of components) and the overall loop H` gain is less than one, as well as nonlinear generalizations based on input to state stability [26, 93], (2) input/output monotone systems theory [2, 3, 94], in which input-output characteristics are the only required “quantitative” data, and (3) the use of passivity-based tools [5, 69, 102]. We now present some examples where the problem is solved using contraction theory, see [88] for further details. Specifically, it is possible to prove that a contraction property on a reduced-order matrix that quantifies the interconnection structure, coupled with contractivity/expansion estimates on the individual component subsystems, suffices to ensure that trajectories of the overall system converge towards each other. See [88] for further details. Note that the construction of such reduced-order matrix uses only norm and matrix-measure estimates, but no precise knowledge, of components (i.e. of the nodes composing the network). As a first example, we show the result of the design for a network of Hopfield neural models in Figure 3. Here the nodes are coupled via the nonlinear coupling functions in Table 1 engineered using contraction to guarantee synchronization (see [88] for further details).

Table 1. Coupling function for network (20). Edge 142, 6 S 2, 2 S 3, 5 S 6 3 S 1, 5 4 2 344

hij 1 x 2

G arctan 1 x 2 Kx 1 1 2 e 2x 2 / 1 1 1 e 2x 2

THIRD QUARTER 2010

Figure 3 shows one of the application of such a multiscale approach. Each of the nodes in such a figure, is a Hopfield neural model [48]: cl

LuxR

Luxl

Al

u (t ) Injection of u (t)

4. Evolving the Network Topology In all the network models presented above, the topology of the connections is considered fixed. Conversely, as mentioned in the introduction, real networks are characterized by evolving topologies. Similar mechanisms appear also in technological applications, such as flocking and rendez-vous. Here, we propose different possible ways of evolving the network dynamics. In Sections 4.1, we propose to evolve the network on-line, on the basis of the dynamical evolution of the network nodes. The edge evolution is modeled as a second order dynamical system, stimulated by an external input, that is the norm of the mismatch between the endpoints of the edge. In Section 4.2, a novel computational tool is presented for the off-line evolution of network topologies that optimise some user-specified performance measure. 4.1. Edge-Snapping The network model (9), (10) allows the decentralized adaptation of the coupling gains associated to each edge in the network. Nonetheless, the network structure is fixed a priori, while we aim at a network that is able

Cell 1

u (t ) Cell 2

u (t ) Cell 3

u (t )

Cell 4

u (t ) Cell 5

u (t ) Cell n

u (t )

(a) cI (arbitray units)

Another notable example is that of synchronization in biological systems where the coupling between nodes is determined by physical characteristics of the systems and of the external environment, which are hard or even impossible to modify with the technologies currently available. To overcome the above technological problem, our approach is that of engineering the intrinsic dynamics of the nodes by tuning their internal parameters. Those parameters are related to biochemical characteristics of the systems of interest, like e.g. protein degradations and promoter strengths which can be easily modified. From a methodological viewpoint, as shown in [82, 89, 90], the above biochemical parameters can be tuned so as to guarantee contraction of the node. Then, network synchronization is attained, since networks of contracting nodes synchronize in the presence of diffusive coupling. To this aim, we use a graphical algorithm, which allows to easily prove/impose contraction [82, 84, 86, 87]. Figure 4 shows an application to the synchronization of a network of Repressilators, coupled by means of a quorum sensing mechanism: the shared quantity, representing the environment where nodes live is in turn affected by an exogenous input, u 1 t 2 .

THIRD QUARTER 2010

Lacl

Lacl

(20)

u(t) (arbitray units)

xi # xi 5 2 1 a 1 hij 1 xj 2 2 hij 1 xi 2 2 1 u 1 t 2 . R j[Ni

TetR

200 150 100 50 0

0

50

100

0

50

100

150 (b)

200

250

300

150 200 Time (min) (c)

250

300

0.8 0.6 0.4 0.2 0

Figure 4. The Repressilator is a synthetic biological circuit that consists of three genes that inhibit each other in a cyclic way [37]. Specifically, gene lacI (associated to the state variable c in the model) expresses protein LacI, which inhibits the transcription of gene tetR. This translates into protein TetR, which inhibits transcription of gene cI. Finally, the protein CI, translated from cI, inhibits expression of lacI, completing the cycle. A modular addition to the three-genes circuit is shown in (a). The module was first presented in [40] and makes the Repressilator circuit sensitive to the concentration of the auto-inducer which is a small molecule that can pass through the cell membrane. Specifically, the module makes use of two proteins: (i) LuxI, which synthesizes the auto-inducer; (ii) LuxR, with which the auto-inducer synthesized by LuxI forms a complex that activates the transcription of various genes. In this figure, the behavior of one protein, representing the output of the Represilator circuit, cI, is shown. Notice that at steady state synchronization is attained onto a periodic orbit having the same period as u 1 t 2 5 0.4 1 0.4 sin 1 0.5t 2 .

to self-determine its topology. With this in mind and according to the model proposed in [29], we associate a dynamical system to each edge in the network, with each coupling gain being the one-dimensional output of the dynamical system: # sij 5 g 1 sij, xj 2 xi, t 2 ,

sij 5 c 1 sij 2 .

(21) (22)

IEEE CIRCUITS AND SYSTEMS MAGAZINE

73

where k is a class k` function2 of the norm of the mismatch between the states of the endpoints. In what follows we choose k 1 z 2 5 az2 while the potential V is such that system (23), (24) has two stable equilibria in 0 and 1, representing respectively the absence and the activation of the corresponding edge. Namely, we select the potential V 1 z 2 5 bz2 1 z 2 1 2 2 which is depicted in Figure 5(a). In this model, the set E of the edges can be viewed as the set of network edges that can be potentially activated. In fact, for physical or economical reasons, there are surely some edges that cannot be activated. Nonetheless, we leave the network nodes free to negotiate what edges have to be activated at steady state among all the feasible links. We think of the coupling gain associated to each edge as a mass which is at rest at the coordinate origin, and then is pushed to the right by an external force given by the squared mismatch between the states of the endpoints of the edge. If this force is strong enough, then the gain snaps to the second equilibrium point (link on), otherwise it asymptotically comes back to the origin, so that the corresponding edge is respectively activated or switched off.

V(σij)

b /2

b /4

b /16 0

0

0.5 σij (a)

1

1.5

xi (t)

20 0 −20 0

2

4

6

8

t 1 σij (t )

σij (t )

20 10 0

0

0.5

0

5

10

15

t

t (b)

Figure 5. Network of 50 Lorenz oscillators coupled through edge snapping: (a) bistable potential driving the evolution of each sij; (b) evolutions of the node states (top) and of the coupling gains (bottom).

Here, sij [ Rm is the state associated to the edge 1 vi, vj 2 , g : Rm 3 Rn 3 R 1 S Rm is the vector field defining the evolution of the edges, and c 1 t 2 : Rm S R is the output function of the edge dynamics, defining the coupling gains. Notice that the network (9), (21), (22) is an example of a generalized dynamic graph without external inputs. It is worth remarking here that an appropriate selection of the function g is fundamental to achieve our goal. According to [29], we consider the so-called snapping dynamics for the edge evolution, given by: # sij112 5 sij122,

(23)

d # sij122 5 2zsij122 2 112 V 1 sij112 2 1 k 1 00 xj 2 xi 00 2 , dsij

(24)

sij 5 sij112 ,

(25)

2

where z is a damping parameter. Therefore, the coupling gain is modeled as a unitary mass in a double-well potential V subjected to an external force k 1 7 xj 2 xi 7 2 ,

Numerical Example As a numerical example, we consider a network of 50 Lorenz systems, coupled on all the state variables. We set the damping z 5 20 while the parameter b of the potential, determining the height of the barrier between the two equilibria, is chosen equal to 16. The set of potential edges, represented by the set E, is generated using a Barabàsi-Albert scale free model. The edge states are assumed to be initially in the coordinates origin 1 sij112 1 0 2 5 s2ij 1 0 2 5 0 2 , implying that at the onset of the evolution the network topology is disconnected. As for the nodes’ initial conditions, they are taken from a random normal distribution with zero mean and standard deviation of 15. As we can see from Figure 5(b), synchronization is asymptotically achieved, while at steady state a topology emerges, which is depicted in Figure 8. The features of this topology will be described later in Section 4.3. 4.2. NetEvo: A Computational Approach An alternative approach is to compute some optimal network topology offline to minimize some performance measure of interest. NetEvo [42] is a computational framework developed to study the evolution of dynamical complex networks3. It provides tools to simulate, evolve and analyze a wide range of systems in the hope of understanding possible links between topology, dynamics and evolution. A central concept within NetEvo is that of a supervisor.

A function F : I S R is positive definite if F 1x 2 . 0, 4x [ I, x 2 0 and F 10 2 5 0. A function f : R $0 S R $0 is of class k if it is continuous, positive defi nite and strictly increasing. An unbounded function of class k belongs to class k`. 3 For further information see the project website at http://www.netevo.org. 2

74

IEEE CIRCUITS AND SYSTEMS MAGAZINE

THIRD QUARTER 2010

This entity has the role of making changes to a systems underlying topology and dynamical parameters, with the aim of minimizing a user-specified performance measure Q. This acts as a guide to the evolutionary process and can be calculated using both topological information and simulated dynamics from the system. Figure 6 illustrates the main components of a supervised network. To use NetEvo a user must provide an initial topology, the form of node and edge dynamics, a method to alter the network topology and dynamical parameters, and a performance measure. Together with a set of internal modules (see Fig. 7) these features are combined and evolution of the system performed. To aid with the analysis of this process, NetEvo has been built using several well established open-source libraries. These include igraph to hold and analyze network topologies, and the GNU Scientific Library for numerical simulation and spectral methods. At present, NetEvo allows for dynamics in the form of ordinary differential equations (ODEs) and contains a built-in supervisor based on a simulated annealing metaheuristic. This supervisor assumes no knowledge of structure within the system of interest and will perform an unbiased probabilistic search of system configurations. An interesting future direction is to improve this search process to make use of either expert knowledge, i.e. it has been proven that certain topological features yield higher performance, or by machine-based learning techniques to find structural and dynamical relationships on-the-fly. In regards to synchronization, NetEvo has been used to study the effect that allowing a systems simulated dynamics to influence the evolutionary process has on enhanced topologies [44, 43]. Previously, the study of evolving topologies for optimal synchronization has focused on topological features of a system [35]. This is in part due to the synchronizability of a network— range of coupling strengths for with full synchronization is achieved—being linked by Pecora and Carrol to minimization of the network Laplacian eigenratio (lN/l2) [9, 75]. Although this type of approach provides insight

Node Dynamics

Edge Dynamics

Dynamics Simulator ODE Solver System Structure Graph and Dynamics

Data Output Numerical and Visual

Updates Parameters

Node/Coupling Dynamics dxi /dt = F (xi ) + coupling

Alters Coupling

Supervisor Computerized Agent

Q

Performance Measure Network Topology Laplacian/Adjacency Matrix

Updates Topology

Figure 6. Flow diagram of a supervised network. The supervisor attempts to minimize a performance measure Q by making changes to the underlying topology and parameters of the dynamics.

into proven optimal topologies, it does not necessarily tell us anything about synchronous systems found in Nature, where evolution is instead directed by the dynamical attributes (phenotype) of an organism. To address this, we attempted to take an alternative approach and used an order parameter performance measure [110] (calculated from system-level dynamics) to guide the evolutionary process. A number of separate evolutions were performed for systems of identical Rössler oscillators with diffusive coupling, and a variety of fixed coupling strengths s^ during the process. Varying the fixed coupling strength altered the stability of the synchronization manifold and allowed us to investigate how the resultant topologies adapted to alternative scenarios, i.e. when full synchronization was not possible. Furthermore, we also considered systems of alternative node dynamics (Lorenz and Chua) and levels of heterogeneity in the dynamical parameters of each node. 4.3. Features of the Evolved Topologies In this section, we discuss some of the emerging features found when using NetEvo and the edge snapping techniques to carry out network evolution. Furthermore, network topologies are analyzed with

Performance Measures

Network Mutations

Supervisor Simulated Annealing Performance Multithreading

Analysis Dynamical and Topological

Figure 7. The main modules of NetEvo split into two main types, core and interface. Core modules, shown in orange, are designed to carry out analysis, simulation and evolution. Interface modules, shown in blue, are to describe the particular system of interest and should be provided by an end-user.

THIRD QUARTER 2010

IEEE CIRCUITS AND SYSTEMS MAGAZINE

75

BOX 2. PROVING STABILITY

O

The Master Stability Function Approach

well as an EDN are well described by a set of differential

The Master Stability Function (or MSF) is an approach for

equations. In terms of the notation introduced in Box 1, this

studying synchronization of diffusively coupled networks (2)

means that both the operators F and t satisfy some non-

mainly used within the Physics community. Specifically, net-

triviality, restriction, semigroup, identity axioms (see e.g.

work dynamics are first linearized around the synchronization

[95]). The problem of analyzing the steady state behavior

manifold so as to obtain a variational equation from (2), de-

of (1) has been addressed in the literature using several ap-

scribing small variations, jk, of network trajectories from the synchronous evolution, say xs 1 t 2 . This equation is then block

ften, in applications a generalized dynamic graph as

proaches, briefly itemized in what follows. Lyapunov Based Proofs

The key idea for the proofs of synchronization based on the use of Lyapunov functions is that of recasting a synchronization problem onto a stability problem of some fixed point. In this view, the typical approach is that of defining an error variable, ei J xi 2 xs, where xs 1 t 2 is the synchronous soT T lution (which is assumed to exist). Let e J 3 e1 , c, eTN 4

and notice that network global synchronization is attained if e 5 0 is a globally asymptotically stable fixed point for the # dynamics of e. Now, stability of such a fixed point can be analyzed by means of some Lyapunov Function, V 1 e 2 . A typical choice for V 1 e 2 is

V 1 e 2 J eTMe. The matrix M is assumed to be positive definite. Depending

'h 1 xs 2 'f 1 xs 2 # 2 slk jk 5 c d jk, k 5 1, c, N. (26) 'x 'x Each of the above equations, or modes, represents the variations (around the synchronization manifold) along a direction of the network phase space. The mode associated to k 5 1 corresponds to the variations along the synchronization manifold, while all other k’s correspond to the transverse eigenmodes to such a manifold. Thus, local synchronization occurs if all such transverse modes are stable. To this aim, the maximum Lyapunov exponent is computed as a function of s, for k . 0. It immediately follows that stability is then guaranteed for all the s’s such that all the maximum Lyapunov exponents are negative. See [75, 74, 9] for further details.

on the choice of such a matrix and on the hypotheses made # on the network nodes which make V , 0 along network

Contraction Based Proofs

trajectories, some sufficient condition for synchronization

lies in the fact that synchronization is seen as a property of

can be then obtained. See e.g. [60, 61, 64, 55, 28, 59] for

some invariant set. The viewpoint of contraction theory is,

further details.

instead, completely different. In the contraction framework,

the aim of better understanding if properties of the emerging topologies can be tuned by parameters of the evolutionary process. Edge Snapping It is straightforward that the main constraint influencing the properties of the emerging topology is the fundamental edge snapping topology E, representing the set of potential network edges. Therefore, to avoid any bias in the analysis of the emerging topology, we refer to the case in which E 5 V 3 V. This is the case in which, at the onset of the evolution, each couple of nodes can communicate, and then decide whether activating or not the corresponding edge. Moreover, in the analysis, we exclude the case of very low potential barriers b, in which almost all the potential nodes are activated. 76

diagonalized to give a set of decoupled variational equations:

IEEE CIRCUITS AND SYSTEMS MAGAZINE

The common feature of the above presented approaches

In [29] it was observed that, for many diverse node dynamics, the initial conditions x 1 0 2 of the network nodes plays a key role in determining the main topological features of the emerging topology, such as the degree distribution or the betweenness centrality. We recall here that the degree of a node i is the number of incident edges, while the betweenness centrality is a measure of importance of a node in the network [4, 39]. Therefore, a node starting from a more scattered initial condition is more likely to activate more connections and have a higher betweenness centrality. Moreover, it was also observed that the edge snapping dynamics induce a topology in which the maximum eigenvector (that is, the eigenvector associated to the maximum eigenvalue) of the Laplacian matrix L is strongly related to the initial conditions of the nodes. Therefore, in the case of the classical consensus problem (see for THIRD QUARTER 2010

stability is intended as a property of trajectories and not as

• a geometric approach, which allows to study also

a property of some invariant set. From a qualitative view-

cluster synchronization, is presented in [76]. Here,

point, a contraction region is a convex open set in phase

contraction is proved only for some directions in

space, where all the trajectories converge towards each

phase space. Specifically, consider the linear invariant

other. A system is said to be contracting if the contraction

subspace of (1), M. Let V be an orthonormal matrix

region coincides with the whole state space. Formally, a

spanning the null of M. Then, the trajectories trans-

system is said to be contracting if there exist some matrix

versal to the invariant subspace globally exponentially

measure for system Jacobian which is uniformly negative

contract towards M if the matrix VJV T is contracting.

definite, see [62, 89]. From the methodological viewpoint,

In network synchronization problems, the subspace M

two main approaches are used to study network conver-

can be chosen as x1 5 c 5 xN.

gence via contraction. Specifically: • a first idea is to use the concept of partial contraction discussed in [104]. Basically, given a smooth nonlinear # n-dimensional system of the form x 5 f 1 x, x, t 2 , assume # that the auxiliary system y 5 f 1 y, x, t 2 is contracting with respect to y. If a particular solution of the auxiliary y-system verifies a smooth specific property, then all

Linking Stability Tools

In [30] it is shown that when the nodes’ dynamics are QUAD within a certain range of parameters, then it is equivalent to contraction. Similarly, in [83] it is shown that if network dynamics are contracting towards the synchronization manifold, then the MSF is negative.

trajectories of the original x-system verify this property C

exponentially. The original system is said to be partially

contracting. Indeed, the virtual y-system has two par-

ticular solutions, namely y 1 t 2 5 x 1 t 2 for all t $ 0 and

x1(t0)

the particular solution with the specific property. Since all trajectories of the y-system converge exponentially to a single trajectory, this implies that x 1 t 2 verifies the specific property exponentially. The specific property above may be e.g. a relationship between state variables, or simply a particular trajectory. In the case of network synchronization, such a relationship is y1 5 c 5 yN.

instance [70]), in which f 1 x, t 2 5 0, edge snapping can be used as a viable method to build topologies wellsuited to quickly synchronize integrators with given initial conditions. These features are not significantly affected by the height of the potential barrier b or by the damping parameter b. Nonetheless, these parameters can be used to tune the other topological properties. For instance, given the network initial conditions, the barrier b can be used to increase the sparseness of the topology, while preserving network synchronizability. As an example, we show the effect of the increase of the barrier in a network of 30 Lorenz oscillators. In this case we consider a random network with average degree 5. As illustrated in Figure 8, as the barrier is increased from 1 to 10, the percentage of active nodes reduces from the 87.5% (Fig. 8(a)) to the 78% (Fig. 8(c)). THIRD QUARTER 2010

x2(t0)

x1(t ) x2(t )

Figure 10. A schematic representation of trajectories inside the contraction region, C. In such a region, any two trajectories, i.e. x1 1 t 2 , x2 1 t 2 , having initial conditions, i.e. x1 1 t0 2 , x2 1 t0 2 , belonging to C converge towards each other.

NetEvo Using NetEvo and the dynamical order parameter performance measure, we see yet more features emerging in the resultant topologies. Some representative examples can be found in Fig. 9. For comparison, Fig. 9(a) shows a network evolved using the topological eigenratio performance measure. This has what is termed an “entangled” structure with a low diameter and narrow betweenness and degree distributions. Of these evolved topologies it is clear to see two main groups. The first consists of our eigenratio evolved network and the order parameter evolved network for Rössler dynamics and a fixed coupling of s^ 5 0.6, shown in Fig. 9(a) and 9(b) respectively. Both these topologies show a convergence in statistical properties towards very homogeneous features and display high levels of synchronization. The fixed coupling strength s^ 5 0.6 falls IEEE CIRCUITS AND SYSTEMS MAGAZINE

77

(a)

(b)

(c)

Figure 8. Network of 30 Lorenz oscillators coupled through the edge snapping topology. Increasing the barrier b between the wells of the potential, the sparseness of the network is increased. In orange are represented the edges belonging to the emerging topology, while in gray are represented the edges that are asymptotically switched off. In (a), the barrier b is set to 1 and the percentage pa of activated edges represents is the 87.5%. In (b), b is increased to 5, with a subsequent decrease of pa to 82.3%. In (c), b 5 10 implies pa 5 78%.

within the stable region of the synchronization manifold. Therefore, the features we see here are likely to be representative of those necessary for optimal synchronization. We will refer to these types of networks as Type 1. The second group, containing Fig. 9(c)–9(e), show very different topological properties. Rather than the homogeneous appearance of Type 1 networks, we see the formation of hubs containing low degree node that surround a highly interconnected central region. We broadly refer to networks with these properties as Type 2. This change is caused by selecting a fixed coupling strength s^ outside the stable region for synchronization. By imposing this constraint, the dynamical performance measure has adapted the resultant topology for the new scenario. This is further shown in the amount of synchronization exhibited by these networks. Although they are unable to see full synchronization, they do see partial synchronization for a much wider range of fixed coupling strengths and at strengths where “entangled” or Type 1 networks see none.

(a)

(b)

(c)

A fascinating aspect of this change is the speed at which it happens. As the fixed coupling strength is slowly varied towards the unstable region, there is not a smooth transition in topological features from Type 1 to Type 2. Instead, the system jumps from one form to the other in what we call a topological bifurcation [44]. It is worth noting that topological measures do not consider dynamics at all and so cannot display adaptive behavior due to changes in a systems dynamical parameters. Another interesting feature shown to play an important role in many complex networks is that of motifs; small sub-graphs found at greater frequencies than would be expected by random chance [67, 1]. To understand if these played a part in our Type 1 and 2 structures, average motif distributions were calculated. Results show (see Table 2) that networks evolved using a topological measure (eigenratio) display relatively few types of motif statistically over or under expressed. In contrast, a more striking pattern is seen when a dynamical performance

(d)

(e)

Figure 9. Topologies evolved using NetEvo for various performance measures, node dynamics and fixed coupling strengths s^ . (a) Eigenratio, (b) Rössler order parameter s^ = 0.6, (c) Rössler order parameter s^ = 0.9, (d) Lorenz order parameter s^ = 4.5, (e) Chua order parameter s^ = 0.85. Two main forms of topology are evident with (a) and (b) displaying “entangled” Type 1, and (c) to (e) hub-like Type 2 features.

78

IEEE CIRCUITS AND SYSTEMS MAGAZINE

THIRD QUARTER 2010

Table 2. Average motif frequencies for topologies evolved using NetEvo. Arrows represent statistically significant over c and under T expression in comparison to a randomized network (p-value * 0.01). Networks evolved using a dynamical performance measure show an over expression of closed loop feedback motifs 2, 4 and 6. Dynamics

Performance measure

Topology

Eigenratio

Type 1

99.90

0.10

Order Parameter s^ = 0.6 Order Parameter s^ = 0.9

Type 1 Type 2

96.74 ≈ 97.18 ≈

3.26 Ω 2.82 Ω

Lorenz

Order Parameter s^ = 4.5

Type 2

98.39 ≈

Chua

Order Parameter s^ = 0.85

Type 2

93.98 ≈

– Rössler Rössler

1

measure (order parameter) is used, regardless of coupling strength or node dynamics. In all cases, these distributions exhibited an over expression of three node closed loop feedback motifs 2, 4 and 6. In [44] it was conjectured that these may help improve localized stability of a shared dynamic and would therefore be an artifact of dynamics influencing topology through the evolutionary process. Furthermore, and in support of this idea, evolution of networks with differing levels of variability in node dynamics showed a correlation up to a threshold between motif expression frequency and the percentage of variation in dynamical parameters [43]. Further investigation into these features is still required, but points towards the interesting possibility of using localized structures to help shape system-level behaviors. 5. Conclusions and Future Work We have shown that synchronization is indeed possible in evolving dynamical networks in a number of different ways. Starting from networks with constant diffusive linear coupling, we showed that it is possible to consider decentralized adaptive gain evolutions. In this case, the strength of the coupling between neighboring nodes is determined by a set of integro-differential equations describing how each of the time-varying gains in the network should be adjusted as a function of the mismatch between the states of those nodes connected to each other. We then showed that an alternative approach is that of using contraction theory so as to construct nonlinear coupling protocols that make trajectories of the nodes in the network converge towards each other. This, we believe, is a particularly timely approach in all those problems where the interest is to prove the network synchronizes rather than showing stability of some asymptotic solution of interest (which is often unknown a priori in synchronization problems). To further mimic the adaptation and time-varying nature of networks of agents in the natural world, we then THIRD QUARTER 2010

2

3

4

5

6

7

19.71 ≈

0.20

15.95 ≈ 40.25 Ω

5.83 Ω 5.53 Ω

1.61 Ω

25.90 ≈

6.02 Ω

36.74 Ω

8

80.05 Ω

–

0.04 ≈

–

77.77 ≈ 52.50 ≈

0.17 Ω 0.33 Ω

0.28 ≈ 1.09 ≈

0.002 Ω 0.01 Ω

3.58 Ω

69.56 ≈

0.11 Ω

0.86 ≈

–

8.60 Ω

51.93 ≈

1.36 Ω

1.28 ≈

0.10 Ω

proposed that the structure of the network can itself be evolved either to guarantee synchronization (edge-snapping) or in order to optimize some cost function of interest. NetEvo was presented as an effective computational framework to deal with this latter type of problems. Interestingly, we noticed that when the network structure is evolved, its topological features depend on the dynamics at the nodes, the cost function and even the structure and parameters of the coupling function. Feedback motifs were for example found to be more likely to occur when dynamic performance measures are used rather than static ones. Despite progress in the area of complex network theory, there are many pressing open problems that still need to be addressed. First of all, synchronization itself is a rather basic type of emergent phenomenon. One might want to induce more sophisticated types of complex behavior in networks. For instance, concurrent or cluster synchronization has been proposed as an important phenomenon in neural networks [76] and other biological networks [65, 109] where different groups of interacting agents synchronize onto different evolutions. An open problem is to understand if it is possible to devise fully decentralized strategies to achieve this aim and, maybe, even evolve network topologies guaranteeing the emergence of cluster synchronization. Another interesting challenge is to further characterize the topological features of evolving dynamical networks and understand what the mechanisms are that induce the presence of modularity, motifs and other topological features in the evolved networks. For example, in [51], it is suggested that modularity and motifs can spontaneously emerge in the presence of switching network goals (cost functions). Here we noticed that, when dynamics is considered, certain motifs can become more frequent than others but there are certainly other possible causes, still unaccounted for, that might explain some of the features observed in natural networks. Also, it is currently assumed that the network nodes share the same identical dynamics, often, modeled in term IEEE CIRCUITS AND SYSTEMS MAGAZINE

79

of some smooth vector field. In practice, nodes will not be identical and the case of networks of near-identical or nonidentical nodes is seldom studied in the literature (some results can be found in [99]) and is waiting to be properly addressed. A striking open problem is when the nodes or coupling protocols are described by piecewise-smooth models [33] as is the case for example of networks of vibroimpacting mechanical devices, gene regulatory networks and switching converters. Finally, the grand challenge of optimally evolving dynamical network structures to perform a desired function remain unsolved. Can we use contraction theory for instance to evolve optimal network topologies? Can we always reduce the problem to some convex optimization problem e.g see [18]? These are only some of the many challenges and open problems that remain to be solved in this new exciting area of research in circuits and systems. The task ahead of us and the next generation of scientists is gigantic. It is likely new analytical and computational tools will be needed and, possibly, even a completely new theoretical framework. But, if successful, the potential impact of these new tools on our understanding of complex systems would be immense and the possible applications so promising that attempting to solve these problems is of utmost importance. Pietro De Lellis is currently a postdoctoral researcher at the University of Naples Federico II. He got his Master Degree in Management Engineering in September 2006 from the same university. His Masters Thesis was dealing with the analysis of the dynamics in complex networks of agents. In November 2006 he won the competition for entering the Ph.D. course in Informatics and Automation Engineering. From January to July 2009 he was a visiting Ph.D. student at the Department of Aerospace and Mechanical Engineering of the Polytechnic Institute of NYU. He received his Ph.D. on December 2009 from the University of Naples Federico II. His current research activity is focused on the analysis and control of complex networks with special emphasis on the novel adaptive strategies for the network synchronization. Mario di Bernardo is currently Associate Professor of Automatic Control at the University of Naples Federico II in Naples, Italy. In 1998 he obtained a Ph.D. in Engineering Mathematics from the University of Bristol, U.K. He was appointed to a Lectureship at the Department of Engineering Mathematics of the same University in 1997 and then promoted to a Readership and a Full Professorship. From 2001 till 2004, he was Assistant Professor of Automatic Control at the University of Sannio, Italy. His research interests 80

IEEE CIRCUITS AND SYSTEMS MAGAZINE

are within the broad area of nonlinear systems, on both dynamics and control. He authored and co-authored more than 150 international scientific publications. From July 1999 to December 2001, he served as Associate Editor of the IEEE Transactions on Circuits and Systems Part I, and from January 2002 to 2006 he was Associate Editor of the IEEE Transactions on Circuits and Systems Express Letters. He is currently AE of the IEEE Transactions on Circuits and Systems Part I. He is a member of the organizing committees of the IEEE Symposium on Circuits and Systems and has been chair or co-chair of many scientific events. In 2004 he was elected to the governing board of the Italian Society for Chaos and Complexity and in 2006 and 2009 to the Board of Governors of the IEEE Circuits and Systems Society. He was elected President of the Italian Society for Chaos and Complexity for the term 2010–2014. He received funding from major funding bodies and industries including the EPSRC, the European Union, the Italian Ministry of Research and University, Jaguar Engineering Centre, QinetiQ. Together with Bernard Brogliato (INRIA, France), he was the organizer and scientific coordinator of the EUR 2.8M EU Project SICONOS on the simulation and control of nonsmooth dynamical systems. On the 28th February 2007, he was honoured with the title of ‘Cavaliere della Repubblica Italiana’ (equivalent to a British OBE) for scientific merits by the President of the Italian Republic. Thomas E. Gorochowski received his Computer Science M.Eng. degree from the University of Warwick, UK in 2004. From 2004 until 2007 he was with Accenture Ltd. working as a technology consultant in the domain of business intelligence and information management. Following this, he joined the Centre for Complexity Sciences at the University of Bristol, UK, receiving a Complexity Sciences M.Res. degree in 2009 and has continued studying towards a Ph.D. in Complexity Sciences at this University. His research interests include the study of structure, dynamics and evolution of complex networks with possible applications to engineering (specifically synthetic biology), novel computing architectures and in the analysis of biological systems. Giovanni Russo was born in Naples on March 28, 1984. He received from the University of Naples Federico II, the Master Degree in Automation Engineering on October 30, 2007 with the thesis “Modelling, Analysis, Identification and Control of Biological Networks”. In November 2007 he won the competition for entering the Ph.D. course at the Department of Systems and THIRD QUARTER 2010

Computer Engineering of the University of Naples Federico II. His research interests are focused on the analysis and control of biochemical and complex networked systems. He spent the period from October 1, 2009 to March 31, 2010 at the Nonlinear Systems Laboratory of the Massachusetts Institute of Technology. References [1] U. Alon, “Network motifs: Theory and experimental approaches,” Nat. Rev. Genet., vol. 8, pp. 450–461, 2007. [2] D. Angeli, J. E. Ferrell, and E. Sontag, “Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems,” Proc. Nat. Acad. Sci. USA, vol. 101, pp. 1822–1827, 2004. [3] D. Angeli and E. Sontag, “Monotone control systems,” IEEE Trans. Automat. Contr., vol. 48, pp. 1684–1698, 2003. [4] J. M. Anthonisse, “The rush in a directed graph,” Stichting Mathematisch Centrum, Amsterdam, The Netherlands, Tech. Rep. BN 9/71, 1971. [5] M. Arcak and E. Sontag, “A passivity-based stability criterion for a class of interconnected systems and applications to biochemical reaction networks,” Math. Biosci. Eng., vol. 5, pp. 1–19, 2008. [6] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Phys. Rep., vol. 469, pp. 93–153, 2008. [7] T. Balch and R. C. Arkin, “Behavior-based formation control for multirobot teams,” IEEE Trans. Robot. Automat., vol. 14, pp. 926–939, 1998. [8] A. L. Barabàsi and L. Albert, “Emergence of scaling in random networks,” Science, vol. 286, pp. 509–512, 1999. [9] M. Barahona and L. M. Pecora, “Synchronization in small-world systems,” Phys. Rev. Lett., vol. 89, p. 054101, 2002. [10] R. W. Beard, A. W. Beard, J. Lawton, and F. Y. Hadaegh, “A coordination architecture for spacecraft formation control,” IEEE Trans. Contr. Syst. Technol., vol. 9, pp. 777–790, 1999. [11] I. Belykh, V. Belykh, and M. Hasler, “Synchronization in complex networks with blinking interactions,” in Proc. 2005 Int. Conf. Physics and Control, 2005, pp. 86–91. [12] I. V. Belykh, V. N. Belykh, and M. Hasler, “Blinking model and synchronization in small-world networks with a time-varying coupling,” Physica D, vol. 195, pp. 188–206, 2004. [13] V. N. Belykh, I. V. Belykh, and M. Hasler, “Connection graph stability method for synchronize coupled chaotic systems,” Physica D, vol. 195, pp. 159–187, 2004. [14] V. N. Belykh, N. N. Verichev, L. J. Kocarev, and L. O. Chua, “On chaotic synchronization in a linear array of Chua’s circuits,” J. Circuits Syst. Comput., vol. 3, pp. 579–589, 1993. [15] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, and C. S. Zhou, “The synchronization of chaotic systems,” Phys. Rep., vol. 366, pp. 1–101, 2002. [16] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. U. Hwang, “Complex networks: Structure and dynamics,” Phys. Rep., vol. 424, pp. 175–308, 2006. [17] S. Boccaletti and D. L. Valladares, “Characterization of intermittent lag synchronization,” Phys. Rev. E, vol. 62, pp. 7497–7500, 2000. [18] S. Boyd, “Convex optimization of graph laplacian eigenvalues,” in Proc. Int. Congr. Mathematicians, 2006, pp. 1311–1319. [19] J. Buck, “Synchronous rhythmic flashing of fireflies II,” Quart. Rev. Biol., vol. 63, pp. 265–289, 1988. [20] T. Chen and X. Liu. Network synchronization with an adaptive coupling strength. 2006. [21] T. Chen, X. Liu, and W. Lu, “Pinning complex networks by a single controller,” IEEE Trans. Circuits Syst. I, vol. 54, pp. 1317–1326, 2007. [22] L. O. Chua, M. Itoh, L. Kocarev, and K. Eckert, “Chaos synchronization in Chua’s circuit,” Univ. California, Berkeley, Tech. Rep., 1992. [23] L. O. Chua, L. Kocarev, K. Eckert, and M. Itoh, “Experimental chaos synchronization in Chua’s circuit,” Int. J. Bifurcat. Chaos, vol. 2, pp. 705– 708, 1992. [24] I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin, “Effective leadership and decision-making in animal groups on the move,” Nature, vol. 433, pp. 513–516, 2005. [25] A. Das, R. Fierro, V. Kumar, J. Ostrowski, J. Spletzer, and C. Taylor, “A vision-based formation control framework,” IEEE Trans. Robot. Automat., vol. 18, pp. 813–825, 2002. [26] S. Dashkovskiy, B. Rüffer, and F. Wirth, “An ISS small-gain theorem for general networks,” Math. Control Signal Syst., vol. 19, pp. 93–122, 2007. THIRD QUARTER 2010

[27] P. DeLellis, M. diBernardo, and F. Garofalo, “Synchronization of complex networks through local adaptive coupling,” Chaos, vol. 18, p. 037110, 2008. [28] P. DeLellis, M. diBernardo, and F. Garofalo, “Novel decentralized adaptive strategies for the synchronization of complex networks,” Automatica, vol. 45, pp. 1312–1318, 2009. [29] P. DeLellis, M. diBernardo, F. Garofalo, and M. Porfiri, “Evolution of complex networks via edge snapping,” IEEE Trans. Circuits Syst. I, to be published. [30] P. DeLellis, M. diBernardo, and G. Russo, “On quad, lipschitz and contracting vector fields for consensus and synchronization of networks,” IEEE Trans. Circuits Syst., to be published. [31] P. DeLellis, M. diBernardo, F. Sorrentino, and A. Tierno, “Adaptive synchronization of complex networks,” Int. J. Comput. Math., vol. 85, pp. 1189–1218, 2008. [32] P. DeLellis, M. diBernardo, and L. F. Turci, “Pinning control of complex networked systems via a fully adaptive decentralized strategy,” IEEE Trans. Automat. Contr., submitted for publication. [33] M. diBernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk, Piecewise-Smooth Dynamical Systems. London: Springer-Verlag, 2008. [34] D. V. Dimarogonas and K. J. Kyriakopoulos, “On the rendezvous problem for multiple nonholonomic agents,” IEEE Trans. Automat. Contr., vol. 52, pp. 916–922, 2007. [35] L. Donetti, P. Hurtado, and M. Muñoz, “Entangled networks, synchronization, and optimal network topology,” Phys. Rev. Lett., vol. 95, p. 188701, 2005. [36] M. Egerstedt and H. Xiaoming, “Formation constrained multi-agent control,” IEEE Trans. Robot. Automat., vol. 17, pp. 947–951, 2001. [37] M. B. Elowitz and S. Leibler, “A synthetic oscillatory network of transcriptional regulators,” Nature, vol. 403, pp. 335–338, 2000. [38] J. H. Fewell, “Social insect networks,” Science, vol. 301, pp. 1867– 1870, 2003. [39] L. Freeman, “A set of measures of centrality based upon betweenness,” Sociometry, vol. 40, pp. 35–41, 1977. [40] J. Garcia-Ojalvo, M. B. Elowitz, and S. H. Strogatz, “Modeling a synthetic multicellular clock: Repressilators coupled by quorum sensing,” Proc. Nat. Acad. Sci., vol. 101, pp. 10955–10960, 2004. [41] M. Golubitsky and I. Stewart, “Nonlinear dynamics of networks: The groupoid formalism,” Bull. Amer. Math. Soc., vol. 43, p. 305, 2006. [42] T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, “Netevo: A computational framework for the evolution of complex dynamical networks,” to be published. [43] T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, “A dynamic approach to the evolution of complex networks,” in Proc. 19th Int. Symp. Mathematical Theory of Networks and Systems, 2010, pp. 1–5. [44] T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, “Evolving enhanced topologies for the synchronization of dynamical complex networks,” Phys. Rev. E, vol. 81, no. 5, p. 056212, 2010. [45] R. O. Grigoriev, M. C. Cross, and H. G. Schuster, “Pinning control of spatiotemporal chaos,” Phys. Rev. Lett., vol. 79, pp. 2795–2798, 1997. [46] W. Guo, F. Austin, S. Chen, and W. Chen, “Pinning synchronization of the complex networks with non-delayed and delayed coupling,” Phys. Lett. A, vol. 373, pp. 1565–1572, 2009. [47] J. H. Holland, Adaptation in Natural and Artificial Systems. Cambridge, MA: MIT Press, 1992. [48] J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proc. Nat. Acad. Sci. USA, vol. 79, pp. 2554–2558, 1982. [49] R. A. Horn and C. R. Johnson, Matrix Analysis. Cambridge, U.K.: Cambridge Univ. Press, 1985. [50] A. Kashiwagi, I. Urabe, K. Kaneko, and T. Yomo, “Adaptive response of a gene network to environmental changes by fitness-induced attractor selection,” PLoS One, vol. 1, p. e49, 2006. [51] N. Kashtan and U. Alon, “Spontaneous evolution of modularity and network motifs,” Proc. Nat. Acad. Sci. USA, vol. 102, p. 13733, 2005. [52] N. E. Leonard and E. Fiorelli, “Virtual leaders, artificial potentials and coordinated control of groups,” in Proc. 40th Conf. Decision and Control, 2001. [53] C. Li, L. Chen, and K. Aihara, “Synchronization of coupled nonidentical genetic oscillators,” Phys. Biol., vol. 3, pp. 37–44, 2006. [54] X. Li, X. Wang, and G. Chen, “Pinning a complex dynamical network to its equilibrium,” IEEE Trans. Circuits Syst. I, vol. 51, pp. 2074–2087, 2004. [55] Z. Li and G. Chen, “Global synchronization and asymptotic stability of complex dynamical networks,” IEEE Trans. Circuits Syst. II, vol. 53, pp. 28–33, 2006. IEEE CIRCUITS AND SYSTEMS MAGAZINE

81

[56] Z. Li, Z. Duan, and G. Chen, “Cost and effects of pinning control for network synchronization,” Chin. Phys., vol. 18, pp. 106–118, 2009. [57] Z. Li, L. Jiao, and J.-J. Lee, “Robust adaptive global synchronization of complex dynamical networks by adjusting time-varying coupling strength,” Phys. A: Stat. Mech. Applicat., vol. 387, pp. 1369–1380, 2008. [58] J. Lin, A. Morse, and B. Anderson, “The multi-agent rendezvous problem,” in Proc. 42nd IEEE Conf. Decision and Control, 2003. [59] X. Liu and T. Chen, “Exponential synchronization of nonlinear coupled dynamical networks with a delayed coupling,” Physica A, vol. 381, pp. 82–92, 2007. [60] X. Liu and T. Chen, “Boundedness and synchronization of y-coupled lorenz systems with or without controller,” Physica D, vol. 237, pp. 630–639, 2008. [61] X. Liu and T. Chen, “Synchronization analysis for nonlinearly-coupled complex networks with an asymmetrical coupling matrix,” Physica A, vol. 387, pp. 4429–4439, 2008. [62] W. Lohmiller and J. Slotine, “On contraction analysis for non-linear systems,” Automatica, vol. 34, pp. 683–696, 1998. [63] E. N. Lorenz, “Deterministic nonperiodic flow,” J. Atmos. Sci., vol. 20, pp. 130–141, 1963. [64] J. Lu and D. W. Ho, “Local and global synchronization in general complex dynamical networks with delay coupling,” Chaos Solitons Fractals, vol. 37, pp. 1497–1510, 2008. [65] X. B. Lu and B. Z. Qin, “Adaptive cluster synchronization in complex dynamical networks,” Phys. Lett. A, vol. 373, pp. 3650–3658, 2009. [66] M. Mesbahi and F. Y. Hadaegh, “Formation flying control of multiple spacecraft via graphs, matrix inequalities, and switching,” J. Guid. Control Dyn., vol. 24, pp. 369–377, 2001. [67] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon, “Network motifs: Simple building blocks of complex networks,” Science, vol. 298, pp. 824–827, 2002. [68] C. C. Moallemi and B. van Roy, “Distributed optimization in adaptive networks,” Adv. Neural Inform. Process. Syst., vol. 16, 2004. [69] P. Moylan and D. Hill, “Stability criteria for large-scale systems,” IEEE Trans. Automat. Contr., vol. 23, pp. 143–149, 1978. [70] R. Olfati-Saber and R. Murray, “Consensus problems in networks of agents with switching topology and time-delays,” IEEE Trans. Automat. Contr., vol. 49, pp. 1520–1533, 2004. [71] G. V. Osipov, J. Kurths, and C. Zhou, Synchronization in Oscillatory Networks. Berlin, Germany: Springer-Verlag, 2007. [72] E. Ott, Chaos in Dynamical Systems. Cambridge, U.K.: Cambridge Univ. Press, 1993. [73] M. Parter, N. Kashtan, and U. Alon, “Facilitated variation: How evolution learns from past environments to generalize to new environments,” PLoS Comput. Biol., vol. 4, p. e1000206, 2008. [74] L. M. Pecora and T. L. Carroll, “Synchronization in chaotic systems,” Phys. Rev. Lett., vol. 64, pp. 821–824, 1990. [75] L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Phys. Rev. Lett., vol. 80, pp. 2109–2112, 1998. [76] Q. C. Pham and J. J. E. Slotine, “Stable concurrent synchronization in dynamic system networks,” Neural Netw., vol. 20, pp. 62–77, 2007. [77] M. Porfiri and M. di Bernardo, “Criteria for global pinning-controllability of complex networks,” Automatica, vol. 44, pp. 3100–3106, 2008. [78] M. Porfiri and F. Fiorilli, “Node-to-node pinning control of complex networks,” Chaos, vol. 19, p. 013122, 2009. [79] M. Porfiri, D. Stilwell, and E. Bollt, “Synchronization in random weighted directed networks,” IEEE Trans. Circuit Syst. I, vol. 55, pp. 3170–3177, 2008. [80] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “Phase synchronization of chaotic oscillators,” Phys. Rev. Lett., vol. 76, pp. 1804–1807, 1996. [81] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “From phase to lag synchronization in coupled chaotic oscillators,” Phys. Rev. Lett., vol. 78, pp. 4193–4196, 1997. [82] G. Russo and M. di Bernardo, “An algorithm for the construction of synthetic self synchronizing biological circuits,” in Proc. Int. Symp. Circuits and Systems, 2009, pp. 305–308. [83] G. Russo and M. di Bernardo, “Contraction theory and the master stability function: Linking two approaches to study synchronization in complex networks,” IEEE Trans. Circuit Syst. II, vol. 56, pp. 177–181, 2009. [84] G. Russo and M. di Bernardo, “How to synchronize biological clocks,” J. Comput. Biol., vol. 16, pp. 379–393, 2009. [85] G. Russo and M. di Bernardo, “Solving the rendezvous problem for multi-agent systems using contraction theory,” in Proc. Int. Conf. Decision and Control, 2009.

82

IEEE CIRCUITS AND SYSTEMS MAGAZINE

[86] G. Russo, M. di Bernardo, and J. J. E. Slotine, “A graphical algorithm to prove contraction of nonlinear circuits and systems,” IEEE Trans. Circuits Syst. I, submitted for publication. [87] G. Russo, M. di Bernardo, and J. J. E. Slotine, “An algorithm to prove contraction, consensus, and network synchronization,” in Proc. 1st IFAC Workshop on Estimation and Control of Networked Systems, 2009. [88] G. Russo, M. di Bernardo, and E. Sontag, “Stability of networked systems: A multi-scale approach using contraction,” in Proc. IEEE Conf. Decision and Control, submitted for publication. [89] G. Russo, M. di Bernardo, and E. Sontag, “Global entrainment of transcriptional systems to periodic inputs,” PLoS Comput. Biol., vol. 6, p. e1000739, 2010. [90] G. Russo and J. Slotine, “Global convergence of quorum-sensing networks,” J. R. Soc. Interface, submitted for publication. [91] D. D. Siljak, Large-Scale Dynamic Systems: Stability and Structure. New York: North Holland, 1978. [92] Q. Song and J. Cao, “On pinning synchronization of directed and undirected complex dynamical networks,” IEEE Trans. Circuits Syst. I, to be published. [93] E. Sontag, “Input to state stability: Basic concepts and results,” in Nonlinear and Optimal Control Theory, P. Nistri and G. Stefani, Eds. Berlin: Springer-Verlag, 2007, pp. 163–220. [94] E. Sontag, “Monotone and near-monotone biochemical networks,” Syst. Synth. Biol., vol. 1, pp. 59–87, 2007. [95] E. D. Sontag, Mathematical Control Theory: Deterministic Finite-Dimensional Systems. New York: Springer-Verlag, 1998. [96] F. Sorrentino, M. diBernardo, F. Garofalo, and G. Chen, “Controllability of complex networks via pinning,” Phys. Rev. E, Stat. Nonlin. Soft Matter Phys., vol. 75, p. 046103, 2007. [97] F. Sorrentino and E. Ott, “Adaptive synchronization of dynamics on evolving complex networks,” Phys. Rev. Lett., vol. 100, p. 114101, 2008. [98] K. O. Stanley, B. D. Bryant, and R. Mikkulainen, “Evolving adaptive neural networks with and without adaptive synapses,” in Proc. 2003 Congr. Evolutionary Computation (CEC’03). [99] J. Sun, E. Bollt, and T. Nishikawa, “Master stability function for coupled nearly identical dynamical systems,” Europhys. Lett., vol. 85, p. 60011, 2009. [100] D. Syliak, “Dynamic graphs,” Nonlin. Anal.: Hybrid Syst., vol. 2, pp. 544–567, 2008. [101] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, “Stale flocking of mobile agents, part II: Dynamic topology,” in Proc. 42nd Conf. Decision and Control, 2003. [102] M. Vidyasagar, Input-Output Analysis of Large Scale Interconnected Systems. New York, Berlin: Springer-Verlag, 1981. [103] L. Wang, H. P. Dai, H. Dong, Y. Y. Cao, and Y. X. Sun, “Adaptive synchronization of weighted complex dynamical networks through pinning,” Eur. Phys. J. B, vol. 61, pp. 335–342, 2008. [104] W. Wang and J. J. E. Slotine, “On partial contraction analysis for coupled nonlinear oscillators,” Biol. Cybern., vol. 92, pp. 38–53, 2005. [105] X. F. Wang and G. Chen, “Synchronization in scale-free dynamical networks: Robustness and fragility,” IEEE Trans. Circuits Syst. I, vol. 49, pp. 54–62, 2002. [106] X. F. Wang and G. Chen, “Synchronization in small-world dynamical networks,” Int. J. Bifurcat. Chaos, vol. 12, pp. 187–192, 2002. [107] C. W. Wu, “Localization of effective pinning control in complex networks of dynamical systems,” in Proc. IEEE Int. Symp. Circuits and Systems, 2008, pp. 2530–2533. [108] C. W. Wu and G. Chen, “Synchronization in an array of linearly coupled dynamical systems,” IEEE Trans. Circuits Syst., vol. 42, pp. 430–447, 1995. [109] W. Wu, W. Zhou, and T. Chen, “Cluster synchronization of linearly coupled complex networks under pinning control,” IEEE Trans. Circuits Syst. I, vol. 56, pp. 829–839, 2009. [110] S.-H. Yook and H. Meyer-Ortmanns, “Synchronization of rössler oscillators on scale-free topologies,” Phys. A: Statist. Theor. Phys., vol. 371, pp. 781–789, 2006. [111] W. Yu, P. DeLellis, G. Chen, M. diBernardo, and J. Kurths, “Distributed adaptive control of synchronization in complex networks,” IEEE Trans. Automat. Cont., submitted for publication. [112] Q. Zhang, J. Lu, J. Lu, and C. Tse, “Adaptive feedback synchronization of a general complex dynamical network with delayed nodes,” IEEE Trans. Circuits Syst. II, vol. 55, pp. 183–187, 2008. [113] C. Zhou and J. Kurths, “Dynamical weights and enhanced synchronization in adaptive complex networks,” Phys. Rev. Lett., vol. 96, p. 164102, 2006. THIRD QUARTER 2010