This work is licensed under a Creative Commons Attribution-NonCommercial-
NoDerivatives 4.0 International License.
This is a postprint version of the following published document:
Sevilla-Escoboza, R., Buldú, J. M., Boccaletti, S., Papo, D.,
Hwang, D.-U., Huerta-Cuellar, G., Gutiérrez, R. (2016).
Experimental implementation of maximally synchronizable
networks. Physica A: Statistical Mechanics and Its
Applications, 448, 113–121
DOI: https://doi.org/10.1016/j.physa.2015.12.086
© Elsevier, 2016
Experimental implementation of maximally synchronizable networks
R. Sevilla-Escoboza,
1
J. M. Buld
´
u,
2, 3
S. Boccaletti,
4, 5
D. Papo,
2
D.-U. Hwang,
6
G. Huerta-Cuellar,
1
and R. Guti
´
errez
7, 8
1
Centro Universitario de los Lagos, Universidad de Guadalajara,
Enrique D
´
ıaz de Leon, Paseos de la Monta
˜
na, Lagos de Moreno, Jalisco 47460, Mexico
2
Laboratory of Biological Networks, Center for Biomedical Technology,
Technical University of Madrid, Pozuelo de Alarc
´
on, 28223 Madrid, Spain
3
Complex Systems Group & GISC, Universidad Rey Juan Carlos, 28933 M
´
ostoles, Madrid, Spain
4
CNR-Istituto dei Sistemi Complessi, Via Madonna del Piano, 10, 50019 Sesto Fiorentino, Italy
5
The Italian Embassy in Israel, 25 Hamered Street, 68125 Tel Aviv, Israel
6
National Institute for Mathematical Sciences, Daejeon 305-811, South Korea
7
Department of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel
8
School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, United Kingdom
Maximally synchronizable networks (MSNs) are acyclic directed networks that maximize synchronizability.
In this paper, we investigate the feasibility of transforming networks of coupled oscillators into their correspond-
ing MSNs. By tuning the weights of any given network so as to reach the lowest possible eigenratio λ
N
/λ
2
,
the synchronized state is guaranteed to be maintained across the longest possible range of coupling strengths.
We check the robustness of the resulting MSNs with an experimental implementation of a network of nonlinear
electronic oscillators and study the propagation of the synchronization errors through the network. Importantly,
a method to study the effects of topological uncertainties on the synchronizability is proposed and explored both
theoretically and experimentally.
PACS: 89.75.Hc, 89.75.Fb
INTRODUCTION
Synchronization is a paradigmatic example of collective
behavior in physical and biological sciences. The scientific
study of synchronization started with the pioneering obser-
vations on the dynamics of pendulum clocks hanging from a
common beam by Christiaan Huygens in the 17th century [1],
and continue these days with studies showing different forms
of synchronization in networks of both regular and chaotic
self-oscillating dynamics. Throughout the years, evidence of
synchronized behavior has been found in mechanical, elec-
tronic and neuronal systems, as well as in chemical reactions
and biological rhythms, to name but a few examples. A flurry
of interest in synchronization started in the mid 1990s, with
some earlier breakthroughs such as the Kuramoto model [2]
paving the ground for it. The results of those fruitful years of
research are reviewed in [1, 3].
Since the beginning of the previous decade, the focus has
been moving from the study of the synchronized dynamics of
just two or a few dynamical system, towards the study or large
ensembles of oscillators with complex coupling arrangements.
In this regard, the development of network theory [4, 5] al-
lowed for a fruitful interaction between a dynamical point of
view and a more topological perspective that has led to the
present-day field of the synchronization of complex networks
[4, 6].
Among these new developments, the master stability func-
tion (MSF) approach [7, 8] is a cornerstone of network syn-
chronization research, and it provides the language in which
the ideas we develop are written. The MSF describes the sta-
bility of the synchronized state of a set of N coupled dynam-
ical units as a function of the coupling parameter σ and the
topology of the network under study. Given a network of N
dynamical systems, the evolution (in isolation) of the state of
each node is described by a set of n-dimensional differential
equations ˙x
i
= F(x
i
), where x
i
∈ R
n
is the dynamic state
vector of node i. If nodes have a certain interaction between
them, the evolution of the coupled systems is given by:
˙x
i
= F(x
i
) − σ
N
X
j=1
L
ij
H(x
j
), i = 1, . . . , N (1)
where σ is the coupling strength, H(x) : R
n
→ R
n
is a
vectorial coupling function and L
ij
are the elements of the
network Laplacian matrix. We here consider weighted net-
works: the connectivity is given by a weighted adjacency
matrix whose elements W
ij
are real numbers giving the link
weights between nodes i and j if they are connected and
zero otherwise. Equivalently, the coupling can be represented
(as in Eq. 1) by a weighted Laplacian matrix defined as
L
ij
= δ
ij
P
k
W
ik
− W
ij
. Note that Eq. (1) is equivalent
to a diffusive coupling between the nodes of the network, the
weight of the coupling contained in W
ij
. Due to the zero-
row sum property of the Laplacian matrix (
P
j
L
ij
= 0) the
synchronization manifold x
1
= x
2
= ... = x
N
≡ x
s
, with
˙x
s
= F (x
s
), is an invariant set of the dynamics. Under these
conditions, the MSF approach provides a framework for the
study of the stability of synchronization in which the topology
and the dynamics are in some sense uncoupled. The stability
of the synchronized dynamics or synchronizability is estab-
lished by computing the maximum Lyapunov exponent of a
suitably modified kernel that depends on a parameter ν, which
is proportional to the coupling strength in the network ν = σλ
[7]. This Lyapunov exponent can be seen as the parameter
giving the exponential divergence/convergence of perturba-
tions orthogonal to the synchronization manifold, and when
parameterized in terms of ν gives the MSF curve, which we
2
denote as Λ(ν). The proportionality constant λ represents one
of the nonzero eigenvalues of the graph Laplacian matrix. The
graph Laplacian matrices of the networks considered in this
work have a real and non-negative spectrum, with eigenvalues
0 = λ
1
< λ
2
≤ λ
3
≤ ··· ≤ λ
N
, where N is the number
of nodes (dynamical units) of the network. Moreover, as the
networks are connected, only one eigenvalue is zero, and thus
the ν corresponding to the different oscillation eigenmodes
will be positive. For those values of ν for which the MSF
is negative, perturbations transversal to the synchronization
manifold damp out exponentially fast. For an uncoupled net-
work, ν = 0, the MSF corresponds to that of the autonomous
dynamics at each node, and is therefore either zero (for regu-
lar dynamics) or positive (for chaotic dynamics). As ν is var-
ied, the MSF may become negative in some regions, and it is
of specially interest to study the boundaries of those regions.
When Λ(ν) is negative for all σλ
i
with i = 2, 3, . . . , N, the
synchronization manifold is stable and the network is said to
be synchronizable [7].
For a given dynamical system and coupling function, three
categories can be defined according to the classification pro-
posed in Ref. [4]: Class I systems have a non-negative
MSF (and are therefore not synchronizable), class II systems
are those that have an unbounded synchronizability region
(Λ(ν) < 0 for ν > ν
c
, where ν
c
represents the only zero of
the MSF) and class III systems are those that have a bounded
synchronizability region (Λ(ν) < 0 for ν ∈ (ν
1
, ν
2
), where ν
1
and ν
2
are the two zeros of the MSF). In practice, only a finite
range of ν that is expected to cover all the cases of interest is
considered, and the classification is applied in this restricted
sense. If one considers larger ν ranges, MSFs with more than
two zero crossings are indeed possible, and some complicated
MSF curves have been reported for example in [9].
Of these three main classes, class I and class II are, in some
sense, trivial. Class I is not synchronizable, while class II sys-
tems are always synchronizable for large enough σ, specifi-
cally for σ > ν
c
/λ
2
. However, the conditions for class III
systems to be synchronizable, λ
2
σ > ν
1
and λ
N
σ < ν
2
, im-
ply that a topology such that λ
N
/λ
2
> ν
2
/ν
1
is not synchro-
nizable for any value of σ. In fact, whatever the system, if it
is class III, a topology with a smaller eigenratio R ≡ λ
N
/λ
2
is easier to synchronize as the different ν
2
≤ ··· ≤ ν
N
that
have to be accommodated within the synchronization region
are bunched together more closely. Following this argument,
the optimal case is that for which λ
2
= λ
3
= ··· = λ
N
,
in which all the relevant values of ν for the given topology,
namely σλ
i
for i = 2, 3, . . . , N , become equal, and the eigen-
ratio R reaches its minimum R = 1, leading to an optimal
synchronizability (i.e., the one being stable for a larger range
of σ). Networks with this particular structure are known as
maximally synchronizable networks (MSN) [10].
In this paper, we investigate the feasibility of transforming a
given network into its MSN and test its robustness to noise and
parameter mismatch in real systems. While some publications
have shown the possibility of enhancing the synchronizability
of networks by rearranging the links [11, 12], in many cases
the creation or deletion of links is not available due to exper-
imental restriction. Here, we consider the architecture under-
lying the network topology to be fixed (nodes are connected
or disconnected once and for all), and try to achieve the graph
that optimizes the stability of the synchronized state by tuning
the link weights. We design an experimental implementation,
by means of electronic circuits, of the MSN and investigate
its stability as a function of the coupling strength. Interest-
ingly, we observe how the synchronization error is propagated
through the network when the system is close to the synchro-
nization boundaries. Next, we analyze the effects of the devia-
tions from the optimal topology on the synchronization of the
whole system, and show the interplay between the topological
noise with the coupling strength of the whole network.
MAXIMALLY SYNCHRONIZABLE NETWORKS:
EXPERIMENTAL IMPLEMENTATION
In this section we discuss the experimental implementation
of MSNs and the corresponding results. It is divided into three
subsections: in the first one, we describe the algorithm that
produces MSNs out of arbitrary topologies; then we describe
the dynamical system of use, and the experimental setup based
on electronic circuits; finally, we discuss the experimental re-
sults.
Maximally synchronizable network algorithm
Any undirected (connected) network can be converted into
a MSN by the following procedure [10]:
1. Select any node of the network as the initial node (from
now on, node 1).
2. Number the k
1
neighbors of node 1 sequentially (i.e.,
give the numers 2, 3, . . . , k
1
+ 1 to them).
3. Repeat the process with the second neighbors (i.e. the
neighbors of 2, 3, . . . , k
1
+ 1), third neighbors (i.e. the
neighbors of the neighbors of 2, 3, . . . , k
1
+ 1). and
so on, until we have numbered to all the nodes of the
network.
4. Transform the links into unidirectional links pointing
from the node that has been assigned the lower number
to the node with the greater number.
5. Give a weight of 1/k
in
i
to all incoming links of nodes
i = 2, . . . , N, k
in
i
being the in-degree of node i (i.e. the
number of neighbors whose links point to it), so that the
total incoming weight of all nodes will be one, with the
exception of node 1 (which is zero).
By following these steps, we obtain a directed weighted
network with a weighted adjacency matrix such that W
ij
=
3
FIG. 1. Construction of Maximally Synchronizable Network
(MSN). (a) Starting from a weighted undirected network, (b) the first
step is the numbering process, which consists in selecting an initial
node and sequentially numbering its first neighbours. The process is
repeated for successive layers of neighbors until the whole network
has been numbered. (c) Next, directions are given to the links with
arrows pointing from the node with the lower number to the node
with the higher number. (d) Finally, the weights of the links pointing
to node i are set to 1/k
in
i
.
1/k
in
i
if there is a link (j → i) or zero otherwise. The pro-
cedure is equally applicable whether the original network is
weighted or not, as only the basic structure as given by the un-
weighted (binary) adjacency matrix is used. The correspond-
ing elements of the resulting Laplacian matrix are
L
ij
=
1 if i = j > 1
−1/k
in
i
if j < i and i and j are connected
0 otherwise
(2)
The full Laplacian matrix therefore is a lower triangular ma-
trix
L =
0 0 0 ··· 0
−1 1 0 ··· 0
−1/k
in
3
(0) −1/k
in
3
(0) 1 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
−1/k
in
N
(0) −1/k
in
N
(0) −1/k
in
N
(0) ··· 1
(3)
where the element L
ij
for j < i being denoted as −1/k
in
i
(0)
means that it can be either −1/k
in
i
if there is a link from j to
i or zero if there is not. As the network is connected, node
1 must have at least one neighbor, and therefore there must
be a link from 1 to 2, and that is why we do not apply this
notation to L
21
, which is known to be −1. The procedure and
the resulting Laplacian matrix are illustrated in Fig. 1. As the
FIG. 2. Experimental setup. Electronic implementation of a net-
work of R
¨
ossler oscillators. The three variables of a R
¨
ossler chaotic
circuit are recoded by means of an analog-to-digital card (ADC) and
recorded in a computer. The same circuit is used to simulate all nodes
of the network, whose coupling matrix (and weights) is sent from the
computer to the circuit through a digital-to-analog card (DAC). A
digital line (DO), controls the coupling strength and the gain product
of perturbations by means of digital potentiometers. See the Ap-
pendix for a detailed description of the circuit parameters.
spectrum of a triangular matrix is given by the elements along
the main diagonal, it is obvious that the matrix corresponds to
a MSN.
Experimental realization with non-linear electronic circuits
The procedure described above can be applied experimen-
tally to arbitrary large networks of electronic circuits by adapt-
ing the methodology first developed in Ref. [13]. In simple
terms, a large unidirectional network of nonlinear circuits is
obtained by the sequential recording of the time series of suc-
cessive layers of neighbors and the weighted reinjection of
the data from previous layers using just one electronic circuit.
The technical details are discussed in the Appendix, and an
illustration of the experimental setup is provided in Fig. 2.
Our experimental system is a network of piecewise linear
R
¨
ossler-like electronic circuits, as illustrated in Fig. 2. The
dynamics of node i is given by the following equations [14]:
˙x = −α
1
i
x
i
+ βy
i
+ Γz
i
+ σψ
P
N
j=1
L
ij
x
j
,
˙y = −α
2
i
(−γx
i
+ [1 − δ]y
i
) ,
˙z = −α
3
i
(−g(x
i
) + z
i
) ,
(4)
where x, y and z are the oscillator state variables. The piece-
wise linear function g(x) defined as
g(x
i
) =
0 x
i
≤ 3,
µ (x
i
− 3) x
i
> 3,
(5)
introduces the nonlinearity in the system that leads to a chaotic
behavior. The parameter values are α
1
= 500, α
2
= 200,
α
3
= 10000, β = 10, Γ = 20, γ = 50, δ = 10.0402,
ψ = 20 and µ = 15. The parameter σ is the coupling strength,
which can be adjusted. For a detailed study of this dynamical
system, whose attractor is sketched in Fig. 2, see Refs. [14–
16]. Concerning the initial topologies of the networks, we
next report results obtained with scale-free networks, as this
4
is a very relevant case from an experimental point of view.
The network size considered is N = 200.
Experimental results and a comparison with MSF predictions
We capture the dynamics of the x(t) variable of each circuit
and compute the overall synchronization error in the network,
hei =
P
i,j
D
x
i
x
j
/N
2
, where D
x
i
x
j
= h|x
i
(t) − x
j
(t)|i and
the angular brackets stand for a time averaging. Times series
have a length of 50000 points after a transient of 5000 is disre-
garded. The synchronization error hei is shown in Fig. 3 (a),
where points corresponding to six different dynamical realiza-
tions of the same MSN are shown, and the continuous line is
the average. If we focus on the lower values of e to study when
the system is the closest to complete synchronization, we can
see that the system is synchronized from, roughly, σ = 0.40
to σ = 2.40. The fact that the system becomes unsynchro-
nizable for low and high values of σ indicates that the R
¨
ossler
system coupled through the x variable is a class III system.
This fact is confirmed when the MSF corresponding to Eq. 4
is calculated numerically, as shown in Fig. 3 (b). This cal-
culation is performed by linearizing the equations of motion
according to the reasoning described in [7]. The calculation
of the maximum Lyapunov exponent, which is based on the
time evolution of a vector in the tangent space that is periodi-
cally renormalized, follows the method proposed in [17]. The
justification for the use of such a method to calculate the Lya-
punov exponent of a piecewise-linear system representing the
macroscopic behavior of electronic systems, such as the one
given by Eq. 4, is provided in [18].
Since the MSN has λ
2
= ··· = λ
N
= 1, one can identify ν
with σ. We can see that while the second zero ν
2
= 2.337 is
close to the upper boundary of the synchronization region in
the experimental results of Fig. 3 (a), the first zero ν
1
= 0.137
is significantly smaller than the lower boundary in relative
terms. To further investigate this discrepancy we have ob-
tained the subgraph synchronization error hei
sub
around the
values of ν
1
and ν
2
, computed as the synchronization error
within the subgraph given by a node and the next M nodes in
increasing order as given by the node labelling. Fig. 4 (a) and
(b) show the values of hei
sub
as a function of the subgraph size
M and the coupling strengths surrounding the corresponding
values of ν
1
and ν
2
. We can observe how at the boundaries of
the synchronization region, hei
sub
increases with M, which
indicates that the experimental error is propagating through
the network.
We conclude that the system size dependence of the syn-
chronization region that is found experimentally can be at-
tributed to the accumulation of synchronization error across
the network layers. While one could be tempted to attribute
this effects to the 12-bit resolution of the analog-to-digital ac-
quisition card (ADC) combined with a sampling of 100 kS/s
(kilosamples per second), which leads to truncation errors
which are amplified due to the chaotic dynamics of the sys-
tem, or to drifts in the electronic component properties that
FIG. 3. Synchronization error and MSF of the system. (a) Syn-
chronization error hei of experimental MSN for 6 different realiza-
tions (black dots) and the average across realizations (red continuous
line). (b) Numerically obtained MSF for the system described in Eq.
4. Dashed lines indicate the values of ν
1
= 0.137 and ν
2
= 2.337.
FIG. 4. Subgraph synchronization error hei
sub
around the values
of ν
1
and ν
2
. (a) Subgraph synchronization error hei
sub
around the
first zero ν
1
= 0.137 of the MSF for different subgraph sizes M. (b)
The same as in (a) but around the second zero ν
2
= 2.337.
![FIG. 6. Synchronization error around the first MSF zero for different topological noise strengths. Experimental study of the synchronization error for Σ = 0.000, 0.015, 0.030 with N = 50. The results here reported are averages of 20 independent realizations. Inset: same results for the full coupling range that has been experimentally explored (σ ∈ [0, 0.2]).](/figures/fig-6-synchronization-error-around-the-first-msf-zero-for-1qxzrexj.webp)
