PHYSICAL REVIEW E 78, 046210 共2008兲

Optimal weighted networks of phase oscillators for synchronization
Takuma Tanaka*
Department of Morphological Brain Science, Graduate School of Medicine, Kyoto University, Kyoto, Japan

Toshio Aoyagi
Department of Applied Analysis and Complex Dynamical Systems, Graduate School of Informatics, Kyoto University, Kyoto, Japanand CREST, JST, Saitama, Japan
共Received 30 March 2008; revised manuscript received 12 August 2008; published 20 October 2008兲
The phase order parameter of oscillators on a network is optimized using two different sets of constraints.
First, the maximization is achieved by adjusting the coupling strengths among the oscillators without changing
the total coupling strength and the natural frequencies of the oscillators. This optimization reveals that astronger weight tends to be assigned to a connection between two oscillators with greatly different natural
frequencies. Second, we vary both coupling strengths and natural frequencies while maximizing the phaseorder and minimizing the penalty function which prevents the natural frequencies of the oscillators from taking
the same value. This optimization reveals that a large total coupling strength makes oscillators take two natural
frequencies 共two-group state兲, whereas a small total coupling strength facilitates the convergence of natural
frequencies to one single value 共one-group state兲. Small and large penalty parameters make the optimized
network take the one- and two-group states, respectively. This phase transition is observed in all-to-all, lattice,and scale-free networks although the clustering coefficient of the strongest links in the optimized network
reflects the difference of the underlying network topologies.
DOI: 10.1103/PhysRevE.78.046210

PACS number共s兲: 05.45.Xt, 64.60.aq, 87.55.de

I. INTRODUCTION

Collective phenomena of oscillatory elements underlie alarge number of fundamental problems in physical, biological, chemical, and social systems. The method of phase description is effective in analyzing these systems, and the resultant systems of coupled phase oscillators yield powerfultheoretical tools applicable to understanding synchronization
observed in various real-world systems 关1兴. In its simplest
solvable form, the evolution of the phase ␾i共t兲 of oscillator i
共i = 1 , . . . , N兲 is governed by the equationNd␾i
K= ␻i + 兺 sin共␾ j − ␾i兲,
dt
N j=1

共1兲

where ␻i is the natural frequency of oscillator i and K is thecoupling strength. Because each oscillator in Eq. 共1兲 is uniformly coupled to all oscillators in the system, this system iscalled a globally coupled phase oscillator system with meanfield interaction. The degree of synchronization is measured
by the order parameter defined asR=1N冑冉兺 cos␾ 冊 + 冉兺 sin␾ 冊 .
2

2

i

i

i

i

Kuramoto showed that the globally coupled phase oscillator
system exhibits a transition from an incoherent state to asynchronized state at a critical coupling strength Kc. In otherwords, the order parameter R takes nonzero values above Kc,
whereas R vanishes below Kc. This type of transition waslater reported for globally coupled phase oscillator systems
with more general coupling functions. Systems with cou-

*ttakuma@mbs.med.kyoto-u.ac.jp1539-3755/2008/78共4兲/046210共10兲

pling functions containing only odd harmonics exhibit transitions from an incoherent state to a synchronized state at a
critical coupling strength 关2兴.
Since the classical work by Kuramoto, the system ofphase oscillators has been a widely studied topic in nonlinear
physics 关3兴, and recently, several groups have moved attention to the investigation of oscillator systems on more complex network structures, such as small world networks andscale-free networks 关4–9兴. Moreno et al. 关4兴 numericallysimulated the behavior of oscillators on scale-free networks
whose dynamics are governed by the evolution of the equationd␾i= ␻i + K 兺 sin共␾ j − ␾i兲,dt
j苸I共i兲
where I共i兲 is the set of neighbors of oscillator i. They reported that the transition to the synchronized state occurred
at a smaller value of the coupling strength K than in a systemof globally coupled oscillators. Theoretical work by Ichinomiya 关6兴 showed that the critical coupling constant at
which transition occurs becomes 0 in random scale-free networks with degree distribution P共k兲 ⬀ k−␥, if 2 ⬍ ␥ 艋 3. A
more recent report compared the synchronizability of scalefree networks and Erdös-Renyi networks 关9兴, and concluded
that oscillators on scale-free networks are more readily synchronized than on Erdös-Renyi networks for small values ofK, whereas for large values of K, the reverse is true.
In these previous studies, the networks had a uniform coupling strength K. However, in real-world networks, the coupling strengths 共called connection weights Wij in the terminology of networks兲, are generally not uniform but
potentially take any positive value depending on the link

046210-1©2008 The American Physical Society

PHYSICAL REVIEW E 78, 046210 共2008兲

TAKUMA TANAKA AND TOSHIO AOYAGI

关10–12兴. Investigations of synchronization on these generalweighted networks have been limited to oscillator systems in
the formNẋi = F共xi兲 − ␴ 兺 GijH关x j兴,
j=1

where ẋ = F共x兲 governs the dynamics of oscillator i, H关x兴 is alinear function, ␴ is the coupling strength, and G is the coupling matrix 关13,14兴. Phase oscillator systems, which cannot
be written in this form because the oscillators are coupled bya nonlinear function, have not been studied on these generalweighted networks.
In this paper, we will investigate the optimal weighted
network of phase oscillators for synchronization by focusing
on the order parameter R. This endeavor requires many complicated questions to be answered, among them the questionsof how large the coupling strength assigned to each linkshould be for better synchronization. We attempt to answer
this question by estimating how each coupling strength andnatural frequency contribute to the facilitation of synchronization of an oscillator for a given network structure. Changing the weight of links in a way that the order parameterincreases, we can obtain the weighted network optimized forsynchronization. Examining the weight obtained through optimization will allow us to understand the essential conditions that must be met for the synchronization of oscillators
on weighted networks. In this paper, we first derive the optimization algorithm to maximize the order parameter R by
using the steepest gradient method and then present numerical results showing the successful optimization of network
weights in the case of all-to-all, lattice-type, and scale-freeconnection.

because without this constraint the coupling strength Kij cangrow without limit. We assume that the connection is symmetric, that is, aij = a ji and Kij = K ji. In the following, we assume that N oscillators in the system are phase locked, andhence the frequency of an oscillator is equal to the average ofthe frequencies of the oscillators,
1
1
1
d␾i 1
= 兺 ␻i + 兺 aijKij sin共␾ j − ␾i兲 = 兺 ␻i = ⍀,兺N i dt N i
N ij
N i
where ⍀ is the average of the natural frequencies. Therefore,
the equation

␻i + 兺 aijKij sin共␾ j − ␾i兲 = ⍀jholds for the system considered in this paper.
First, we derive the steepest gradient of Q1 with respect toKij. When the coupling strength Kij is changed to Kij + ⌬Kij,
the difference of a phases ␾ j − ␾i between oscillators i and jchanges to ␾ j + ⌬␾ j − ␾i − ⌬␾i, which is a solution of

␻i + 兺 aij共Kij + ⌬Kij兲sin共␾ j + ⌬␾ j − ␾i − ⌬␾i兲 = ⍀.
jAssuming ⌬Kij to be infinitesimally small, the change ⌬␾ j
− ⌬␾i is described by

␻i − ⍀ + 兺 aij共Kij + ⌬Kij兲 ⫻ sin共␾ j + ⌬␾ j − ␾i − ⌬␾i兲j

= ␻i − ⍀ + 兺 aij共Kij + ⌬Kij兲关sin共␾ j − ␾i兲cos共⌬␾ j − ⌬␾i兲j

+ cos共␾ j − ␾i兲sin共⌬␾ j − ⌬␾i兲兴⬇ ␻i − ⍀ + 兺 aij共Kij + ⌬Kij兲 ⫻ 关sin共␾ j − ␾i兲j+ cos共␾ j − ␾i兲共⌬␾ j − ⌬␾i兲兴

II. ALGORITHM⬇ 兺 aij⌬Kij sin共␾ j − ␾i兲 + aijKij cos共␾ j − ␾i兲

A. Optimization with respect to coupling strengthj

We derive here the algorithm to obtain the couplingstrengths so that the optimization function
Q1 = R共2兲

2

is maximized for coupled Kuramoto oscillators. We use theoptimization function Q1 instead of R to simplify the calculation. The dynamics of the oscillators are described byd␾i= ␻i + 兺 aijKij sin共␾ j − ␾i兲dtj

共i = 1, . . . ,N兲,

Note that this system of linear equations with variables ⌬␾i
共i = 1 , . . . , N兲 is degenerate because these variables appear
only in the form of ⌬␾ j − ⌬␾i. Although, in general, degenerate system of linear equations cannot be solved, this system can be solved as follows. The equation reads

共3兲

where ␻i is the natural frequency of oscillator i, aij is theadjacency matrix 共1 if oscillators i and j are connected and 0
otherwise兲, and Kij is the coupling strength between oscillators i and j. We maximize Q1 under the constraint␣,兺ij K␣ij = Ktotal⫻共⌬␾ j − ⌬␾i兲 = 0.

共4兲where ␣ and Ktotal are constant parameters. Each summation
for Kij is assumed to be taken over all 共i , j兲, where 1 艋 i
艋 N and 1 艋 j 艋 N. Thus, the L␣ norm of the couplingstrength is fixed to Ktotal. The constraint Eq. 共4兲 is essential

共5兲

Mx = y,where
x = 关⌬␾i兴,
y=

冋兺
j

and

冉

册

aij⌬Kij sin共␾ j − ␾i兲 ,冊

M = ␦ij 兺 ailKil cos共␾l − ␾i兲 − aijKij cos共␾ j − ␾i兲 .
lBecause Mx = 0 if ⌬␾1 = ⌬␾2 = ¯ = ⌬␾N, the rank of matrix
M is N − 1. Thus, its eigenvector corresponding to zero ei-

046210-2

PHYSICAL REVIEW E 78, 046210 共2008兲

OPTIMAL WEIGHTED NETWORKS OF PHASE…

f(ωi-ωj)genvalue is e = 共1 , 1 , . . . , 1兲, and because M is a symmetricmatrix, the other eigenvectors are orthogonal to e. Hence, thesolution x exists if y is orthogonal to e. This condition issatisfied since the coupling matrix is symmetric. Thus, using
the specific value decomposition
M = U⌺UT ,we obtain-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
ωi-ωj

x = U⌺+UTy,
where ⌺+ is the matrix ⌺ with every nonzero entry replacedby its reciprocal. Therefore, using Kij = K ji, we obtain⳵␾k ⳵␾k ⳵␾k ⳵ y i ⳵␾k ⳵ y j==+
⳵Kij ⳵K ji ⳵y i ⳵Kij ⳵y j ⳵K ji

FIG. 1. Penalty function f共␻i − ␻ j兲 for A = 1 introduces repulsion
and attraction among the natural frequencies of the oscillators. Thisfunction takes its minima at ␻i − ␻ j = ⫾ 1.

= Wkiaij sin共␾ j − ␾i兲 + Wkja ji sin共␾i − ␾ j兲

Kij ← Kij + ⑀= 共Wki − Wkj兲aij sin共␾ j − ␾i兲,where Wkj is the 共k , j兲 component of the matrix U⌺+UT. Differentiating Q1 with respect to Kij, we obtain冉⳵␾k⳵Q1 ⳵Q1 1==− 2 兺 cos ␾k 兺 sin ␾k⳵Kij⳵Kij ⳵K ji N2kk+ 2 兺 sin ␾k 兺 cos ␾k
k

k

冊

⳵␾k
= Fij .
⳵Kij

␣共␣ − 1兲 ␣−2 2
Kij uij ,2

␣ − 1 ␣−2 2
Kij uij = 0.
2

Differentiation of the function
L ⬇ Q1 + ␮␣ 兺 K␣ij−1uij +ij

␣ − 1 ␣−2 2
Kij uij
2with respect to uij gives the equation
Fij + ␮␣关K␣ij−1 + 共␣ − 1兲K␣ij−2uij兴 = 0.
Rearranging, we obtain the optimal direction⌬Kij = uij = −

冉

冊

1
Fij
+ Kij .
␣ − 1 ␮␣K␣ij−2

Because the second term does not affect the direction, we
can replace it with
⌬Kij = −␭兺 f共␻i − ␻ j兲,
N2 ij

共8兲

where f共x兲 = 共x − A兲2共x + A兲2 and ␭ is a constant parameter. In
the following, we set A = 1. The first term is the square of theorder parameter and the second term is a penalty term, which
prevents oscillators from having identical natural frequencies
共Fig. 1兲. We add the second term to the optimization function
because maximizing the order parameter without the penalty
term while allowing the natural frequencies to vary justmakes the oscillators converge to the same natural frequency.
The second term introduces repulsion and attraction among
the natural frequencies of the oscillators.
Because maximization of the optimization function Eq.
共8兲 can be done by changing Kij and ␻i independently, we
can derive the steepest gradient with respect to Kij and ␻i
separately. Since the steepest gradient of the optimization
function with respect to Kij is given by Eq. 共7兲, it is sufficientfor us to derive the steepest gradient with respect to ␻i. Proceeding in the same way as in the derivation of the gradientwith respect to Kij, we obtain the differentiation of phase ␾ jwith respect to ␻i,

⳵␾ j= V ji ,⳵␻iFij␮␣共␣ − 1兲K␣ij−2

共7兲

In order to investigate the effect of natural frequencies inthe synchronization of phase oscillators, we examine the optimization with respect to not only coupling strength but also
the natural frequencies of the oscillators in the network to
maximize the optimization function
Q2 = Q1 −and assuming uij Ⰶ Kij, we can replace the constraint Eq. 共4兲
with

ij

,

and renormalize it with the constraint Eq. 共4兲. Note that,since the approximation in Eq. 共6兲 does not hold for ␣ = 1,
the optimal direction Eq. 共7兲 cannot be used in this case.

共6兲␣ 兺 K␣ij−1uij +FijK␣ij−2

B. Optimization with respect to coupling strengthand natural frequency

The optimal direction ⌬Kij = uij under the constraint Eq.
共4兲 is given by the Lagrange method. Using兺ij 共Kij + uij兲␣ ⬇ 兺ij K␣ij + ␣K␣ij−1uij +

9
8
7
6
5
4
3
2
1
0

where
Pz = w,and thereby we update Kij with
046210-3

PHYSICAL REVIEW E 78, 046210 共2008兲

TAKUMA TANAKA AND TOSHIO AOYAGIz = 关⌬␾i兴,
w = 关⌬␻i兴,

冉

冊

k

R

P = ␦ij 兺 aikKik cos共␾k − ␾i兲 − aijKij cos共␾ j − ␾i兲 ,z = Y⌶+Y Tw = Vw,

⌬␻i =
=

冉
冉 兺

⳵
␭
Q1 − 2 兺 f共␻ j − ␻k兲
⳵␻i
N jk

冊

1
⳵␾ j
cos ␾ j 兺 sin ␾ j
2 −2
⳵␻i
N
j
j
+ 2 兺 sin ␾ j 兺 cos ␾ jjj

The natural frequency ␻i is updated by

冉

冊

⳵␾ j
− 2␭ 兺 f ⬘共␻ j − ␻i兲 .
⳵␻ij␻i ← ␻i + ␩ ⌬␻i −冊

1兺 ⌬␻ j .
N j

共9兲The last summation is subtracted from Eq. 共9兲 in order to set
the average of the natural frequencies to 0.
III. RESULTS
A. Optimization with respect to coupling strength

Since the present algorithm requires all the oscillators tobe phase locked, Ktotal must be large enough at the beginningof the optimization process to allow the oscillators connectedby the initial random coupling matrix to be phase locked. To
obtain the optimized coupling strength in the network with asmall Ktotal, we adopted the following simulation schedule:
start
and run the simulation for period
共i兲 set Ktotal to a large Ktotal
T1 according to Eq. 共3兲 without running the optimization
algorithm; 共ii兲 run the simulation and optimization algorithm,
start
end
to Ktotal
linearly for period T2,
decreasing Ktotal from Ktotal
start
end
start
兲 / T2; and 共iii兲 run thethat is, Ktotal = Ktotal + 共t − T1兲共Ktotal − Ktotal
simulation and optimization algorithm without changing
Ktotal for period T3. We updated the coupling strength using
Eq. 共7兲 at t = 10n, where n is a positive integer, to allow thephase of oscillators to converge to a new equilibrium state in
the intervals of updates. In the simulations presented here,
the interval of 10 was long enough for the phase of oscillators to converge.
First, we optimized the coupling strength without changing the natural frequency. The all-to-all coupling strengths
for the N = 20 oscillators were drawn from an uniform distribution and normalized using Ktotal = 1 and ␣ = 2. Natural frequency of oscillator i 共i = 1 , . . . , 20兲 was set to 0.05共i − 10.5兲.

Ktotal,end=1
Ktotal,end=0.801000 2000 3000 4000 5000 6000
Time

FIG. 2. Order parameter of the oscillators on an all-to-all network improved through optimization. Here the following parametervalues were used: N = 20, ␣ = 2, ⑀ = 0.1, ␩ = 0, Kstart
total = 1, T1 = 1000,
T2 = 4000, and T3 = 1000.

After the phases of oscillators and the order parameter converged in first T1 = 1000, we started to update the couplingstrength using Eq. 共7兲 and ⑀ = 0.1. We stopped updating thetotal coupling strength at t = T1 + T2 = 5000. Figure 2 showsendthat the order parameter reached 0.98 in the case of Ktotal
end= 1, and 0.97 in the case of Ktotal = 0.8 after optimization.
After t = 1500, the order parameter decreased in the case ofend= 0.8 since Ktotal started to decrease at t = 1000.
Ktotalend
Figure 3 shows the final state of the network for Ktotal= 1. The coupling strength of a link became large only when
the difference of natural frequencies of the oscillators connected by this link was large, and thus the oscillators withnatural frequencies far from the average gained strong connections. This is a natural result due to the following reasons.
Assigning a large weight to the coupling between two oscillators with close natural frequencies results in little increasein the order parameter, whereas assigning a large weight tothe link between oscillators with very different natural frequencies improves the order parameter more effectively. In
other words, the cost-efficient way is to concentrate the coupling between two oscillators whose natural frequencies arevery different from the average in opposite directions. We(a)(b)
Index of the oscillator

and Y⌶Y T is a specific value decomposition of P. In a similar way as for M in the previous derivation, P is degenerateand symmetric, and thus w must be orthogonal to the eigenvector e. We hence set the average of the natural frequencies
to 0. Differentiating the optimization function Eq. 共8兲 withrespect to ␻i yields the steepest direction,

1
0.99
0.98
0.97
0.96
0.95
0.94
0.93
0.92
0.91
0.9

0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0

5
10
15
20
5
10 15 20
Index of the oscillator

FIG. 3. Final coupling strength of the network with Kend
total = 1 inFig. 2. 共a兲 The oscillators are indicated by circles and their natural
frequencies are indicated by gray levels. Black and white correspond to lowest and highest natural frequencies, respectively. Oscillators with high and low natural frequencies have strong connections indicated by thick lines. 共b兲 Large coupling strength is foundonly on links connecting two oscillators with large deviation of thenatural frequency from the mean. The coupling strength between
oscillators with close natural frequencies is weak 共diagonal whitebelt兲.

046210-4

PHYSICAL REVIEW E 78, 046210 共2008兲

OPTIMAL WEIGHTED NETWORKS OF PHASE…

Large ␭⬇A

Small ␭⬇0

兩 ␻ 1 − ␻ 2兩Large Ktotal⬇A

Small Ktotal
⬇0

3

(b)

1

0.05

B. Optimization with respect to coupling strength and natural
frequency

Next, we examined the outcome of optimizing the function Eq. 共8兲 by changing both coupling strength Kij and natural frequency ␻i. We updated the coupling strength and natural frequency using Eqs. 共7兲 and 共9兲 at t = 10n 共T1 ⬍ t兲, where
n is a positive integer. To understand the essential factors
affecting the resultant network structure, let us examine how
parameters ␭ and Ktotal change the natural frequencies ofoscillators in the optimized network of two oscillators before
presenting the simulation result. If ␭ is large and thereby the
first term of Eq. 共8兲 can be ignored, the optimization function
Eq. 共8兲 is maximized by setting the natural frequencies of theoscillators to −A / 2 and A / 2, because the penalty function
f共x兲 of the natural frequency has minima at x = ⫾ A. In contrast, if ␭ is small and we can ignore the second term, thenatural frequencies of all oscillators tend to converge to thesame value in order to facilitate synchronization. Large Ktotal
allows the system of two oscillators with greatly differing
natural frequencies to attain R2 ⬇ 1, and thus the optimization
function Eq. 共8兲 can be approximated by␭兺 f共␻i − ␻ j兲.
N2 ij

共10兲

Under this condition, maximizing Eq. 共10兲 sets the natural
frequencies of the oscillators to −A / 2 and A / 2. On the other
hand, if Ktotal is small, the natural frequencies of oscillatorsconverge to a single value in order to improve the first termof Eq. 共8兲 because otherwise a small Ktotal may not allow theoscillators to be phase locked. Thus, a large Ktotal and a large
␭ make the oscillators take two natural frequencies, and asmall Ktotal and a small ␭ facilitate the convergence of natural
frequencies to one single value 共Table I兲.
Next, let us examine whether or not the qualitative results
described above are applicable to systems with more than
two oscillators. Figure 4共a兲 shows the phase diagram for the
optimized network of 50 oscillators coupled in an all-to-allmanner. In this diagram, the total coupling strength Ktotal atendranges from 1 to 5 and ␭
the end of the simulation Ktotal

0

10

20

0

10

20

i

30

40

50

30

40

50

(d)
0.6
0.4
0.2
0
-0.2
-0.4
-0.6

ωi

ωi

found that this result held for systems with the distribution ofnatural frequencies other than a uniform distribution, such asa Gaussian or a Lorentzian.
When we started the optimization process from different
initial coupling matrices, the resultant coupling converged toalmost the same matrices. The present optimization algorithm hence is expected to be able to reach the global maximum.

0.1

λ

(c)

Q2 ⬇ 1 −

*

(d)

2

0.6
0.4
0.2
0
-0.2
-0.4
-0.6

0.4
0.2
0
-0.2
-0.4
-0.6

(c)

4

ωi

兩 ␻ 1 − ␻ 2兩(b)0.6(a) 5
end
KtotalTABLE I. Natural frequencies of the optimized two-oscillator
system.

0

10

20

i

30

40

50

i

FIG. 4. Optimization of the oscillators for an all-to-all networkwith respect to coupling strength and natural frequency. 共a兲 Phase
diagram of the optimized network. Simulation results are classified
into the 共almost兲 symmetric two-group state 共dark gray兲, the onegroup state 共white兲, and the other states 共light gray兲, which include
the asymmetric two-group state and the multigroup state. The blackline shows the theoretical border of one- and two-group states.
More detailed simulation results for Kend
total = 2 共dashed line兲 areshown in Fig. 5. 共b兲–共d兲 Examples of the natural frequencies ofoscillators after optimization. We present the one-group state 共␭= 0.01 and Ktotal = 1兲 in 共b兲, the two-group state 共␭ = 0.1 and Ktotal= 4兲 in 共c兲, and a state that is classified in neither group 共␭ = 0.05 andKtotal = 2兲 in 共d兲. Here the following parameter values were used:
N = 50, ␣ = 2, ⑀ = 10, ␩ = 5, A = 1, Kstart
total = 20, T1 = 100, T2 = 1900, andT3 = 8000.

ranges from 0.01 to 0.1. The resultant networks are classified
into three categories, one-group state, two-group state, andmiscellaneous. The natural frequencies of the oscillators in
the one-group state 关Fig. 4共b兲兴 converged to the same value,
and in the two-group state 关Fig. 4共c兲兴 the natural frequenciesof the oscillators converged to values close to −A / 2 or A / 2.
In the two-group state, couplings linking two oscillators inthe same group converged to zero, and couplings linking two
oscillators from different groups converged to the same value
共not shown兲. Note that ␭ in the region where two-groupstates were observed is larger than that in the one-group region, and Ktotal in the region where two-group states wereobserved is larger than that in the one-group region, as in the
two-oscillator networks considered above. In addition, there
were many networks which were classified as neither onegroup nor two-group states 关Fig. 4共d兲兴.
In contrast to the optimization with respect to only Kij, theoptimization with respect to both Kij and ␻i did not necessarily converge to unique Kijs and ␻is. In the two-groupstate, whether the natural frequency of oscillator i converged
to A / 2 or −A / 2 depended on the initial connections andnatural frequencies.
Whether the two-group state is more favorable than the
one-group state can be discriminated by analytical calculation. We derive the optimization function of Eq. 共8兲 in twogroup and one-group states and compare them to estimate

046210-5

PHYSICAL REVIEW E 78, 046210 共2008兲

TAKUMA TANAKA AND TOSHIO AOYAGI␣,= Ktotal

1
N2
−

i

cos ␾i冊 冉兺
2

sin ␾i+iN␭N− 2 N f共2␻兲 + f共0兲
N
2
2

冊册
2

冊

␭
= cos2 ␾ − 共f共2␻兲 + f共0兲兲
2

再冋 冉 冊

K␣ 1/␣
␭
f 2n totalsin 2␾
= cos2 ␾ −22m

冎

0.02 0.04 0.06 0.08
λ

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

0

10 20 30 40 50
i

0

10

0.1

(d) 0.5

0

10

20

30

40

50

0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4

20

30

40

50

i

FIG. 5. Number of connections in 0.8-core, clustering coefficient 共CC兲, and strength Si of the nodes exhibiting critical behavior
at ␭ ⬇ 0.02 in a network with all-to-all connection optimized withKend
total = 2 共a兲. The optimized network structures with ␭ = 0.01, 0.02,
and 0.1 are shown in 共b兲, 共c兲, and 共d兲, respectively. Here the same
parameter values as in Fig. 4 were used.

1 2
␭2
2 兺 f共␻i − ␻ j兲 = 2 共N cos ␾兲
N ij
N

冉

2

Σi σi

i

␭
兺 f共␻i − ␻ j兲
N2 ij

冋冉 兺

(c)
ωi

where we assume a total of m connections have the same
coupling strength K, the oscillators have n connections eachto the oscillators in the other group, and the two groups havethe same number of oscillators, N / 2. In the case of all-to-all
coupling and two-dimensional 共2D兲 lattice coupling, theseparameters take the values n = N / 2 and m = Nn / 2 = N2 / 4, andn = 4 and m = 2N, respectively. If the maximum of the optimization function of the two-group state,=

CC in core1
0
0

3×10-13
2×10-13
1×10-130-1×10-13-2×10-13

Connections in core

0.2
0.1
0

␻ + nK sin共− 2␾兲 = 0,

2Qtwo
2 共␾兲 = R −

(b)

400
300
200
100

ωi

2mK

␣

(a)

ωi

which is more favorable. To obtain the optimization functionof a two-group state in which two groups have natural frequencies ␻ and −␻, and phases ␾ and −␾, we can derive theequationsthat most of coupling strength is gained by a relatively small
number of connections in networks optimized with ␭ = 0.02.
To quantify this observation, we define the p-core C p of aweighted network as the smallest set of strongest connections which satisfies

册␣
2K␣ij 艌 pKtotal,兺共i,j兲苸C ,i⬍jp

+ f共0兲 ,is larger than that of the one-group state, Qone
2 = 1 − ␭f共0兲, thetwo-group state becomes stable. Because Qtwo
2 共␾兲 is an even
one
共0兲
=
Q
,
the
sign
of
the
secondderivative
function and Qtwo
22of Qtwo
2 共␾兲 at ␾ = 0 determines the stability of the two-group2state. Specifically, if 兩d2Qtwo
2 / d␾ 兩␾=0 ⬎ 0, the two-groupstate is more favorable than the one-group state, and if
2is less favorable. Be兩d2Qtwo
2 / d␾ 兩␾=0 ⬍ 0, two-group state
␣
2 two
2
2 2 Ktotal 2/␣cause d Q2 / d␾ = −2 + 32␭A n 共 2m 兲 at ␾ = 0, this is theborder of the two- and one-group states. When the oscillators
are coupled in an all-to-all manner, in the two-group state, an
oscillator in a group has n = N / 2 connections to the oscillators in the other group, and the connection strength to theoscillators in the same group becomes 0. Thus, the total number of connections is m = N2 / 4, and the border of the one2
= 1 for ␣ = 2 andand two-group states is given by 8␭Ktotal
A = 1 关solid line in Fig. 4共a兲兴. This theoretical border is consistent with the simulation results. Large ␭ and Ktotal favor
the two-group state and small ␭ and Ktotal favor the onegroup state, as discussed above.
Next, we examined what takes place at the border of theend=2one- and two-group states in the network with Ktotal
关dashed line with asterisk on the side in Fig. 4共a兲兴. We foundwhere ␣ is set to 2. In Fig. 5共a兲, the natural frequencies ofoscillators are indicated by the y coordinate of the dots, andconnections in the 0.8-core are indicated by lines, which arefewer in the networks optimized with ␭ = 0.02 than in thosewith ␭ = 0.01 and ␭ = 0.1. Optimizing M = 20 networks withall-to-all connectivity starting from random connection
strengths and natural frequencies, we found that the numberof connections in the 0.8-core was a minimum at ␭ = 0.02
关Fig. 5共a兲兴. The middle panel of Fig. 5共a兲 shows that theclustering coefficient of the 0.8-core peaks at this point, andstarts to decrease until ␭ = 0.04 where it vanishes. The clustering coefficient in the core vanishes at large ␭ because the
network takes the two-group state, and in this state connections among the oscillators in the same group are weakened
关Fig. 5共d兲兴. The network optimized with ␭ = 0.01 关Fig. 5共b兲兴
also exhibits a small clustering coefficient because strong
connections are formed only between the oscillators whose
natural frequencies are of opposite signs, although the natural frequencies of all oscillators are close to zero and thereby
this network is classified into the one-group state. In contrast,
the clustering coefficient of the network optimized with ␭= 0.02 can be larger. With this parameter, in some networks
the natural frequencies of oscillators i and j converge to finite values ␻i and ␻ j 共␻i ⬍ 0 ⬍ ␻ j兲, and the natural frequencies of other oscillators vanish 关Fig. 5共c兲兴. In this case, theconnection between oscillators i and j is strengthened, and046210-6

PHYSICAL REVIEW E 78, 046210 共2008兲

OPTIMAL WEIGHTED NETWORKS OF PHASE…

15

5

10
x

15

20

(d)

y

0

10
x

15

20

(g)

0.4
0
-0.2

15

5

10
x

15

20

0
-0.2
-0.4

20

-0.6

5

10
x

15

20

-0.6

0.4
0.2

10

0
-0.2
-0.4
10
x

15

20

(i)

-0.6

0.6
0.4

5

0.2

-0.6

0.6

5

0.4

10

20

20

-0.6

0.6

15

-0.4

20

20

5

0.2

10

15

(h)

0.6

5

10
x

15

15

-0.4
5

10
x

5

-0.2

20

-0.6

-0.4

(f)

0.2
0

0
-0.2

5

0.4

10

0.2

10

20

-0.6

0.6

15

-0.4
5

20

5

-0.2

20

15

(e)

0.2

15

10
x

0.4

15

-0.4
5

0.4

10

-0.2

20

-0.6

0.6

5

0

y

10

0.6

5

0.2

15

-0.4

20

y

y

-0.2

y

y

0

(c)

0.4

5

0.2

10

0.6

y

0.4

5

y

(b)

0.6

y

(a)

0.2

10

0
-0.2

15

-0.4

20
5

10
x

15

20

-0.6

FIG. 6. Natural frequencies of the oscillators on a lattice structure. Here the following parameter values were used: N = 400, ␣ = 2, ⑀
end
= 1, ␩ = 4, Kstart
total = 2000, Ktotal = 32, T1 = 100, T2 = 400, and T3 = 99 500. Penalty parameter ␭ was set to 共a兲 0.001, 共b兲 0.002, 共c兲 0.003, 共d兲 0.004,
共e兲 0.005, 共f兲 0.006, 共g兲 0.007, 共h兲 0.01, and 共i兲 0.02. Theoretical border of the one-group and two-group states is ␭ = 0.0061. All the systems
presented in this figure started from the same initial condition.

oscillators k with natural frequencies close to zero are connected by strong links to both oscillators i and j. This structure contains a large number of triangles composed of thethree oscillators i, j, k. This means that the clustering coefficient grows larger in this network structure. This abruptchange in the clustering coefficient at ␭ = 0.02 suggests that aphase transition occurs at ␭ ⬇ 0.02. We now examine thechange of the network structure further, and define K共m兲ij as
the optimized connection strength between nodes i and j in
the network starting from the mth initial condition and S共m兲i
as the strength of node i, in the mth network as
共m兲 ␣
S共m兲
i = 兺 共Kij 兲 .
jMeasuring the variance of the strength of node i,

␴2i =

冉

冊

2
1
1
2
共S共m兲
S共m兲
,
兺
兺
i 兲 −
i
M m
M m

where the summation is taken over M initial conditions, wefound that the summation of the variance 兺i␴2i has a peak at␭ ⬇ 0.02 关Fig. 5共a兲兴. The peaks in the variance and clustering
coefficient, and the drop in the number of connections in thecore, indicate a phase transition from a one-group to a twogroup state in the optimized network structure.
In the previous paragraphs, we described the results for an
all-to-all network. Next, we examined whether or not we cangeneralize this relationship between the parameter ␭ and thegroup formation of natural frequencies to other network
structures. We take as an example a 2D lattice of 20⫻ 20
oscillators, each of which is connected to oscillators on thetop, bottom, left, and right. We use periodic boundary conditions under which the topmost oscillators are connected tothe bottommost oscillators and the leftmost and rightmost
oscillators are also connected. Figure 6 shows the natural
frequencies of the oscillators in this system at the end of the
optimization process. In this simulation, Ktotal was fixed to 32
and ␭ ranged from 0.001 to 0.02. The natural frequencies
converged to zero for small ␭, and they showed frustrated
checker patterns, which correspond to a two-group state, forlarge ␭. Thus, as in an all-to-all network, the network takes
the two-group state for large ␭ whereas a small ␭ favors the
one-group state. The transition from the one-to two-groupstate takes place around the parameter ␭ = 0.0061 predictedby the above analysis using n = 4 and m = 2N. Hence, therelationship between the parameter ␭ and the group forma-

046210-7

PHYSICAL REVIEW E 78, 046210 共2008兲

TAKUMA TANAKA AND TOSHIO AOYAGI

= 

> 

? 

@ 

FIG. 7. 共Color online兲 Optimized connections and natural frequencies of the oscillators on a Barabási-Albert network. The oscillators areindicated by circles in which their degree numbers are shown. Their natural frequencies are indicated by colors. Positive, negative, and zeronatural frequencies are indicated by red, blue, and white respectively. Connections in 0.8-core 共see definition in the main text兲 are indicatedend
by green lines. Here the following parameter values were used: N = 100, ␣ = 2, ⑀ = 10, ␩ = 5. Kstart
total = 1000, and Ktotal = 10. Penalty parameter ␭
was set to 共a兲 0.001, 共b兲 0.01, 共c兲 0.1, and 共d兲 0.1.

tion of natural frequencies holds for the 2D lattice.
Finally, we examined the structure emerging through the
optimization of synchronization on a scale-free network. We
optimized the synchronizability of a Barabási-Albert network with 100 nodes 共Fig. 7兲. In this simulation, we adopted
start
the following simulation schedule: 共i兲 set Ktotal to Ktotal= 1000 and run the simulation for period T1 = 100 according
to Eq. 共3兲 without running the optimization algorithm; 共ii兲run the simulation and optimization algorithm, decreasing
middle= 100 linearly for period T2 = 100; 共iii兲 run the
Ktotal to Ktotal
simulation and optimization algorithm, decreasing Ktotal toend= 10 linearly for period T2⬘ = 1900; and 共iv兲 run the simuKtotallation and optimization algorithm without changing Ktotal for
period T3 = 98 000. We found that all natural frequencies converged to zero if ␭ = 0.001, and that almost all natural frequencies converged to around ⫾A / 2 if ␭ = 1, similarly to the
behavior observed in networks with all-to-all coupling andlattice structures. However, the networks optimized using ␭= 0.01 and ␭ = 0.1 cannot be classified into the one-groupstate or the two-group state. In these networks, the natural
frequencies of oscillators with large degree numbers tended
to converge to a value far removed from the mean, whereas
those of oscillators with small degree numbers vanished 共Fig.
8兲. This is not the case with the one-group and two-group
networks resulting from the optimization with ␭ = 0.001 and␭ = 1, which exhibited a rather uniform degree-natural frequency relationship. In addition, we found that a relatively
small number of connections gained most of the couplingstrength in the networks optimized with ␭ = 0.01 关Fig. 7共b兲兴
and 0.1 关Fig. 7共c兲兴. In Fig. 7, connections in the 0.8-core,
indicated by green lines, are fewer in the networks optimized
with ␭ = 0.01 关Fig. 7共b兲兴 than in those with ␭ = 0.001 关Fig.
7共a兲兴 and 1 关Fig. 7共d兲兴. Optimizing M = 20 networks with thesame topological structure starting from random connection
strengths and natural frequencies, we found that the numberof connections in the 0.8-core was minimal at ␭ ⬇ 0.008 共Fig.

046210-8

PHYSICAL REVIEW E 78, 046210 共2008兲

Average of |ωi|

OPTIMAL WEIGHTED NETWORKS OF PHASE…
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

80
70
60
50
0.016
0.012
0.008
0.004
0
200
150
100
50
0
0.001

λ=0.001
λ=0.01
λ=0.1
λ=1
0

5

10

15

20

25

ki

Connections in core

CC in core

2

Σi σi

0.01

0.1

1

λ

FIG. 8. Average of 兩␻i兩 as a function of ki tends to be monotonically increasing. In this simulation, the same parameter values, network topology, and optimization schedule as in Fig. 7 were used.

9, top兲. The middle panel of Fig. 9 shows a steep decrease in,and subsequent vanishing around ␭ = 0.01, of the clustering
coefficient of the 0.8-core. This is in sharp contrast to thecase with an all-to-all network, because in an all-to-all network a peak of the clustering coefficient was observed at thevalue of ␭ where the number of connections in the core isminimal. However, in the same way as in the all-to-all network, this suggests that an abrupt structural change occurs at
␭ ⬇ 0.008. The summation of the variance 兺i␴2i has a peak at␭ ⬇ 0.008 共Fig. 9, bottom兲. We confirmed that this peak coincided with the sharp declines of the number of connectionsin the core and the clustering coefficient in the networks withlarger and smaller number of nodes. These peaks and declines suggest that there takes place a phase transition of the
network structure obtained through the optimization of synchronization on a scale-free network.

IV. DISCUSSION

The optimization of the degree of synchronization of oscillators on a network with a given structure was discussed inthis paper. One of the advantages of our approach over previous ones is that the present algorithm allows us to understand how the weights should be distributed in the networkin order to maximize synchronization. In our study, the relationship between connection strength and natural frequency
in networks favorable for synchronization was revealed byour optimization algorithm. More specifically, strongly connecting oscillators with largely differing natural frequencies
optimizes the synchronization. This property, which would
be useful in application for real-world systems, has not been
paid much attention in the previous works considering large
random scale-free and Erdös-Renyi networks.
In addition, we observed the transition of the optimal network structure between the one- and two-group states. Thistransition was observed in networks with all-to-all, lattice,and scale-free network topologies. We confirmed that in both
all-to-all and scale-free networks, the number of connectionsin the core decreased at the border of these two states and thevariability of the optimized network structure has a peak at
this border. In contrast, at this border, the clustering coefficient has a peak in the all-to-all network, whereas it showed

FIG. 9. Number of connections in 0.8-core, clustering coefficient 共CC兲, and strength Si of the nodes exhibiting a critical behavior at ␭ ⬇ 0.008 in a scale-free network. Clustering coefficient of thewhole network is 0.049. In this simulation, the same parametervalues, network topology, and optimization schedule as shown in
Fig. 7 were used.

a sharp drop in the scale-free network. This difference needs
further investigation.
Our algorithm has some limitations. First, our optimization algorithm is not applicable to systems with ␣ = 1. This isbecause the approximation Eq. 共6兲 cannot be used in thiscase. For the system with ␣ ⬎ 1, we optimized the couplingstrength of the network using Eq. 共7兲 and found that theorder parameter was maximized when links between two oscillators with greatly different natural frequencies had strong
connections, as we saw in the system with ␣ = 2 共data notshown兲.
Second, our algorithm assumes all oscillators in the system are phase locked. Thus, optimization through this algorithm does not necessarily generate the network structure inwhich oscillators with random phases rapidly converge to aphase-locked state. Further work is required to characterize
the network structure which facilitates fast synchronization.
Third, in the optimization with respect to natural frequency, we used the optimization function Eq. 共8兲. Since thegenerated patterns depend not only on the parameters ␭ andKtotal as we saw in the previous section, but also on the
penalty function f, further work will be required to characterize the synchronization of the oscillators on the weighted
network in a more general situation.
In this paper, we derived an optimization algorithm for theoscillators with the coupling function sin共␾ j − ␾i兲 because itis the simplest and most intensively studied coupling function for oscillator systems, and not because the optimization
functions for the other coupling functions cannot be derived.
Optimization algorithms for systems with more complex
coupling functions can be derived. These algorithms will answer the question of how to assign strong connections inindividual physical, chemical, and biological systems with agiven structure and coupling function. In addition, our algorithm will be useful not only in constructing highly synchronized networks, but also in examining to what extent theconnection strength of a given network is optimized for synchronization.

046210-9

PHYSICAL REVIEW E 78, 046210 共2008兲

TAKUMA TANAKA AND TOSHIO AOYAGI关1兴 Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence
共Springer-Verlag, Berlin, 1984兲.
关2兴 H. Daido, Prog. Theor. Phys. 88, 1213 共1992兲.
关3兴 J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R.
Spigler, Rev. Mod. Phys. 77, 137 共2005兲.
关4兴 Y. Moreno and A. F. Pacheco, Europhys. Lett. 68, 603 共2004兲.
关5兴 Y. Moreno Vega, M. Vázquez-Prada, and A. F. Pacheco,
Physica A 343, 279 共2004兲.
关6兴 T. Ichinomiya, Phys. Rev. E 70, 026116 共2004兲.
关7兴 D.-S. Lee, Phys. Rev. E 72, 026208 共2005兲.
关8兴 J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. E 71,
036151 共2005兲.

关9兴 J. Gomez-Gardenes, Y. Moreno, and A. Arenas, Phys. Rev.
Lett. 98, 034101 共2007兲.
关10兴 M. E. J. Newman, Phys. Rev. E 64, 016131 共2001兲.
关11兴 A. Barrat, M. Barthélemy, R. Pastor-Satorras, and A. Vespignani, Proc. Natl. Acad. Sci. U.S.A. 101, 3747 共2004兲.
关12兴 R. Ferrer i Cancho and R. V. Solé, Santa Fe Institute Technical
Report No. 01-03-016, 2001 共unpublished兲, http://ideas.repec.org/p/wop/safiwp/01-03-016.html
关13兴 M. Chavez, D.-U. Hwang, A. Amann, H. G. E. Hentschel, and
S. Boccaletti, Phys. Rev. Lett. 94, 218701 共2005兲.
关14兴 F. Sorrentino, M. di Bernardoa, G. Huerta Cuéllar, and S. Boccaletti, Physica D 224, 123 共2006兲.

046210-10