# Node differentiation dynamics along the route to synchronization in complex
networks

Christophe Letellier http://www.atomosyd.net/spip.php?article1  Rouen
Normandie University — CORIA, Campus Universitaire du Madrillet, F-76800
Saint-Etienne du Rouvray, France<EMAIL_ADDRESS>Irene Sendiña-
Nadal  Complex Systems Group & GISC, Universidad Rey Juan Carlos, 28933
Móstoles, Madrid, Spain<EMAIL_ADDRESS>Center for Biomedical
Technology, Universidad Politécnica de Madrid, 28223 Pozuelo de Alarcón,
Madrid, Spain  Ludovico Minati  Center for Mind/Brain Sciences (CIMeC),
University of Trento, 38123 Trento, Italy  Institute of Innovative Research,
Tokyo Institute of Technology, Yokohama, 226-8503, Japan<EMAIL_ADDRESS>I.
Leyva  Complex Systems Group & GISC, Universidad Rey Juan Carlos, 28933
Móstoles, Madrid, Spain  Center for Biomedical Technology, Universidad
Politécnica de Madrid, 28223 Pozuelo de Alarcón, Madrid, Spain
<EMAIL_ADDRESS>
###### Abstract

Synchronization has been the subject of intense research during decades mainly
focused on determining the structural and dynamical conditions driving a set
of interacting units to a coherent state globally stable. However, little
attention has been paid to the description of the dynamical development of
each individual networked unit in the process towards the synchronization of
the whole ensemble. In this paper, we show how in a network of identical
dynamical systems, nodes belonging to the same degree class differentiate in
the same manner visiting a sequence of states of diverse complexity along the
route to synchronization independently on the global network structure. In
particular, we observe, just after interaction starts pulling orbits from the
initially uncoupled attractor, a general reduction of the complexity of the
dynamics of all units being more pronounced in those with higher connectivity.
In the weak coupling regime, when synchronization starts to build up, there is
an increase in the dynamical complexity whose maximum is achieved, in general,
first in the hubs due to their earlier synchronization with the mean field.
For very strong coupling, just before complete synchronization, we found a
hierarchical dynamical differentiation with lower degree nodes being the ones
exhibiting the largest complexity departure. We unveil how this
differentiation route holds for several models of nonlinear dynamics including
toroidal chaos and how it depends on the coupling function. This study
provides new insights to understand better strategies for network
identification and control or to devise effective methods for network
inference.

Complex networks — Dynamical complexity — Route to synchronization — Chaos

##  I Introduction

The description of networks can be structural, based on the characterization
and modelling of their topological properties,Newman (2003); Boccaletti  _et
al._ (2006, 2014) or dynamical, usually referring to the collective behavior
emerging from the interaction between node dynamics and the network
architecture. Synchronization is the most thoroughly investigated ensemble
dynamics,Arenas, Díaz-Guilera, and Pérez-Vicente (2006); Boccaletti  _et al._
(2006); Rodrigues  _et al._ (2016) including a rich variety of related
behaviours such as chimera states,Kuramoto and Battogtokh (2002); Abrams and
Strogatz (2004); Hagerstrom  _et al._ (2012) cluster synchronization or Pecora
_et al._ (2014) explosive synchronization.Gómez-Gardenes  _et al._ (2011);
Leyva  _et al._ (2013)

The route to the synchronous state has been approached from different points
of view. Microscopically, it has been analyzed the different ways local
synchronization grows as the coupling increases depending on the connectivity
structure.Arenas, Díaz-Guilera, and Pérez-Vicente (2006); Restrepo, Ott, and
Hunt (2006); Gómez-Gardeñes, Moreno, and Arenas (2007); Pereira, van Strien,
and Tanzi (2020) Globally, the synchronizability of a population of identical
networked oscillators, that is, the stability of the coherent state, has been
tackled through the master stability function and foresees the emergence of
patterns when the stability of the network uniform state is lost.Pecora and
Carroll (1998); Barahona and Pecora (2002); Pecora (2008)

While these approaches help to understand how synchronization clusters grow
and eventually merge into a macroscopic coherent state or predict its
stability for a particular network structure and coupling function, few works
pay attention to the description of the nodal dynamics along the route to
synchronization. Recently, it has been shown that the degree of complexity of
the node dynamics is strongly dependent on the connectivity; in particular,
the local connectedness appears to confer considerable node differentiation
hallmarks, albeit following mechanisms that are intricate and delineate a non-
monotonic effect.Tlaie  _et al._ (2019a, b); Minati  _et al._ (2019) These
ensemble fingerprints in the node dynamics can be used to infer the network
statistical description, or even the detailed connectivity in certain
conditions.Eroglu  _et al._ (2020) By node differentiation it is here meant as
the dynamical departure of an isolated node due to the interaction with its
neighborhood, being more pronounced in networks with nodes having degree
heterogeneity. In particular, a relatively large range of coupling strengths
over which nodes with a higher degree have a less complex dynamics has been
reported,Tlaie  _et al._ (2019b) explaining the low complexity observed in
hubs of functional brain networks.Martínez  _et al._ (2018) However, the
reverse situation, the existence of an even lower range of coupling strengths
where high degree nodes are more complex than lower degree ones, remains under
study.Minati  _et al._ (2019) The purpose of the present work is to elucidate
the route of the node differentiation dynamics in a complex network before
reaching a synchronous state, and, in particular, how it depends on the
dynamical system, the coupling function and on the topological properties of
the network structure.

In order to characterize and quantify this route to node differentiation
dynamics as the coupling strength increases, we used the complexity measure
recently introduced by Letellier et al.Letellier, Leyva, and Sendiña Nadal
(2020) able to distinguish organized from disorganized chaotic behavior
together with maps of the dynamical patterns associated to relevant degree
classes. Section II is devoted to introducing terminology and providing a
brief description of the dynamical complexity measure. We then investigated in
Section III the dynamical differences among Rössler oscillators in a star
network according to their degree, and how those differences dependent upon
the coupling function, the nominal dynamics, and the size of the star. In
Section IV, we extended the analysis to star networks of dynamical systems
featuring higher nominal complexity, namely the symmetric Lorenz system, the
high-dimensional Mackey-Glass equation, and the toroidal Saito system, to
determine whether a general scenario can be outlined. In Section V, we
considered Rössler systems in larger networks to explore the node
differentiation induced by embedding the oscillators in a much more complex
topological environment. Finally, Section VI offers general conclusions.

##  II Dynamical complexity measure

Let us start by considering the general description of a network comprising
$N$ diffusively-coupled $m$-dimensional identical dynamical systems, whose
state vector $\mbox{x}_{i}$ evolves as

$\dot{\mbox{x}}_{i}=\mbox{f}(\mbox{x}_{i})-d\sum_{j=1}^{N}\mathcal{L}_{ij}\,\mbox{h}(\mbox{x}_{j})$
(1)

where $\mbox{f}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ and
$\mbox{h}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ are the nominal dynamics —
here, nominal refers to the vector field of the $i$-th nodal dynamics when
isolated — and the coupling function, respectively, and $d$ is the coupling
strength. $\mathcal{L}_{ij}$ are the elements of the Laplacian matrix encoding
the network’s connectivity, with $\mathcal{L}_{ii}=k_{i}$ the node degree,
$\mathcal{L}_{ij}=-1$ if nodes $i$ and $j$ are connected, and
$\mathcal{L}_{ij}=0$ otherwise.

To characterize the nodal dynamics in a more reliable way, we computed a
Poincaré section to the trajectory $\mbox{x}_{i}(t)$ in the state space to
rule out the local linear component and focus on the nonlinear signatures
governing the dynamics.Letellier (2006) Thus, for each oscillator, in the
following, we computed a two-dimensional Poincaré section of the attractor
when the dynamics has a toroidal structure, or a first-return map built with
one of the coordinates of the Poincaré section in case it is non-toroidal.
Then, each oscillator’s map is used to obtain a complexity coefficient
$C_{\rm{D}}$ introduced in Ref. Letellier, Leyva, and Sendiña Nadal (2020) and
defined as

$C_{{\rm D}}=S_{{\rm p}}+\Delta,$ (2)

being $S_{\rm p}$ a permutation entropy and $\Delta$ a structurality marker.
We chose this complexity definition because of its ability to distinguish
organized unpredictable behavior (dissipative chaos) from disorganized
unpredictable phenomena (noise or conservative chaos). While the entropy
quantifies the unpredictability of the dynamics, the structurality accounts
for its undescribability, two very different aspects contributing to the
complexity of a dynamics.

The permutation entropy $S_{\rm p}$ is defined as in Ref. Bandt and Pompe
(2002):

$S_{\rm p}=-\frac{1}{\log(N_{\rm s}!)}\sum_{\pi}^{N_{\rm
s}!}p_{\pi}\log{p_{\pi}}\in[0,1]\,.$ (3)

It is based on the probability distribution $p_{\pi}$ of the $N_{\rm s}!$
possible ordinal patterns constructed from the order relations of $N_{\rm s}$
successive data-points in the Poincaré section. The structurality is computed
by dividing the first-return map (or the two-dimensional projection of the
Poincaré section) into a $N_{q}\times N_{q}$ boxes and

$\Delta=\sum_{i,j=1}^{N_{q}}\frac{q_{ij}}{N_{q}^{2}}\in[0,1]\,,$ (4)

where $q_{ij}=1$ if the box $(i,j)$ is visited at least once, and $q_{ij}=0$
otherwise. The structurality $\Delta$ is typically low ($\Delta<0.2$) for
organized dynamics, and large ($\Delta>0.80$) for disorganized dynamics. We
thus have $C_{\rm D}\approx 0$ for a limit cycle ($S_{\rm p}\approx 0$ and
$\Delta\approx 0$), $0.5<C_{\rm D}<1.0$ for dissipative chaotic behavior
(organized but unpredictable), and $C_{\rm D}>1.5$ for weakly dissipative
chaos or stochastic processes (both are unpredictable and disorganized).

In all our computations we used $N_{s}=6$ as the length of the ordinal
patterns. The number $N_{\rm d}$ of data points in the Poincaré section to
compute $S_{\rm p}$ and $\Delta$ was $N_{\rm d}=10,000$, largely above the
requirement $N_{\rm d}>5N_{\rm s}!$ recommended by Riedl et al.Riedl, Müller,
and Wessel (2013) The number of boxes was $N_{q}=5\log(N_{\rm d})$ chosen such
that $N_{\rm d}\gg N_{q}$, as prescribed in Ref. Letellier, Leyva, and Sendiña
Nadal (2020). The box width $\delta_{\rm p}$ is determined by the largest
range visited in our simulations for a given system as $\delta_{\rm p}=(x_{\rm
max}-x_{\rm min})/N_{q}$, where $x_{\rm max}$ and $x_{\rm min}$ are the
maximal and minimal values recorded along one axis of the first-return map. To
avoid excessively promoting noise contamination or small fluctuations around
period-1 limit cycle, we introduced a “noise filter” to interpret the ordered
patterns of the $N_{\rm s}$ data points to allow permutation only if
$|x_{i}-x_{j}|>\delta_{\rm p}$.

Figure 1: Dynamical characterization of a single Rössler system upon variation
of the $a$ parameter and fixing $b=0.2$ and $c=5.7$. Bifurcation diagram using
the Poincaré section ${\cal P}_{x}$ (a) and dynamical complexity $C_{\rm D}$
(b). The dashed vertical lines in (a) and (b) are located at $a_{1}=0.26$,
$a_{2}=0.31$, and $a_{3}=0.38$. The corresponding first-return maps from
$y_{n}$ are shown in (c)-(e) and their complexity values at the bottom of each
panel.

Let us consider the paradigmatic Rössler systemRössler (1976) to exemplify the
ability of Eq.(2) to discriminate between different dynamics. The vector flow
in Eq. (1) is

$\mbox{f}(\mbox{x})=\left[-y-z,x+ay,b+z(x-c)\right]$ (5)

taking $b=0.2$ and $c=5.7$ as fixed parameters and $a\in[0.08,0.38]$ as the
bifurcation parameter. Equation (1) was integrated using a fourth-order Runge-
Kutta scheme Press, Teukolsky, and Flannery (1992) with a time step $\delta
t=0.01$. Initial conditions are randomly selected within a small neighborhood
of radius $0.5$ centered at the origin of the state space. Figure 1(a) shows
the bifurcation diagram as a function of the parameter $a$ for a single
Rössler system ($d=0$) computed from the $y_{n}$ coordinate of the Poincaré
section

${\cal
P}_{x}\equiv\left\\{\displaystyle(y_{n},z_{n})\in\mathbb{R}^{2}~{}|~{}x_{n}=x_{-},\dot{x}_{n}>0\right\\}\,,$
(6)

where $x_{-}=(c-\sqrt{c^{2}-4ab})/2$ is the coordinate of the inner singular
point.Letellier, Dutertre, and Maheu (1995) The classical period-doubling
cascade as a route to chaos is quantitatively characterized by the dynamical
complexity $C_{\rm D}$ in Fig. 1(b). It perfectly captures the periodic
windows ($C_{\rm D}<0.2$) intermingled with chaotic behavior. Panels (c)-(e)
show the first-return maps for the three $a$ values marked in (a)-(b) with
vertical lines. The first value, $a=0.26$ is located between a period-3 and a
period-2 window and its first-return map is bimodal [Fig. 1(c)] but with a
slightly developed third branch appearing just after the period-3 window. For
$a=0.31$, Fig. 1(d), the bimodal map has the three branches. A more developed
chaos, characterized with three monotone branches, is observed for $a=0.38$ in
Fig. 1(e). In particular, the right end of the third branch approaches the
vicinity of the bisecting line, indicating that a new period-1 orbit is about
to be created in the population of unstable periodic orbits. This is the most
developed chaos that can be observed along this line of the parameter space,
as indicated by the dynamical complexity ($C_{\rm D}=0.89$), in agreement with
the most developed first-return map.

Figure 2: Colormaps of the complexity values for a star of $N=16$ coupled
Rössler oscillators in the $d$-$a$ parameter plane. Actual complexity change
of the hub $C_{\rm hub}-C_{0}$ (a)-(b) and the difference $C_{\rm hub}-\langle
C_{\rm leaves}\rangle$ (c)-(d) when oscillators are coupled through the $y$
(a,c), and $x$ (b,d) variables, respectively. Bottom panels (e)-(f) show the
synchronization error $E$. The color scales are provided on the right sides.
In (a)-(f) panels, white curves are the corresponding null isolines, and
horizontal dashed lines in (a)-(b) are plotted at $a=0.26,0.31$ and $a=0.38$.
White domains correspond to the ejection of the trajectory to infinity. Other
parameter values are set as in Fig. 1.

##  III A star network of Rössler systems

Let us now consider a star network of $N$ Rössler systems coupled according to
Eq. (1), with $N-1$ peripheral nodes and one central node acting as the hub,
aiming at exploring the influence of the coupling function in the parameter
space $d-a$. Top panels (a) and (b) of Fig. 2 show the actual change in
complexity $C_{\rm hub}-C_{0}$ of the hub in a $N=16$ star with respect to the
reference complexity value $C_{0}$ (which is a function of the parameter $a$)
of an uncoupled node, for coupling schemes through variables $y$ and $x$,
respectively. Middle panels (c) and (d) provide the complexity difference
between the hub and the leaves $C_{\rm hub}-\langle C_{\rm leaves}\rangle$.

From the colormaps, the first clear observation is that the coupling function
has a role in the node differentiation dynamics, as reflected by the different
distribution of complexity change on the $a$-$d$ plane, which strongly depends
on the network synchronizability. Panels (e,f) in Fig.2 show the time-averaged
synchronization error computed as

$E=\frac{2}{N(N-1)}\displaystyle\sum_{i\neq
j}\|\mbox{x}_{i}-\mbox{x}_{j}\|\,.$ (7)

which is in agreement with the prediction given by the master stability
function (MSF) for Rössler systems coupled through the $y$ variable, type I
synchronizability class -the MSF becomes negative above a critical value, and
$x$ variable, type II class -the MSF is negative in a bounded interval.Pecora
and Carroll (1998); Huang  _et al._ (2009); Sendiña-Nadal, Boccaletti, and
Letellier (2016) For the particular network structure of a star of size
$N=16$, complete synchronization is indeed reached for a larger number of
pairs ($a$,$d$) when Rössler systems are $y$-coupled than when $x$-coupled. In
the latter case, full synchronization cannot be reached for $a\gtrsim 0.25$
before a boundary crisis ejects the trajectory to infinity (white domains in
Fig. 2(f)). For this configuration, the location of $d\lambda_{i}$, being
$\lambda_{i}$ the $i$-th eigenvalue of the Laplacian matrix, is outside the
range where the MSF is negative for any value of $d$, and the synchronous
solution is always unstable.

Roughly, in the regions where the synchronization error is not null, for both
$x$ and $y$ coupling, there is a large region of the parameter space with
yelowish areas delimited by the white null isolines in panels in Figs. 2(a)
and 2(b) where the effect of coupling is to render the dynamics of the hub
more complex with respect to its uncoupled regime ($C_{\rm hub}>C_{0}$).
Outside this region, there are smaller islands (in blue) where the dynamics
turns out to be less complex. However, when comparing this departure from the
uncoupled dynamics between the hub and the leaves, Figs. 2(c) and 2(d), there
are regions where the hub dynamics appears to be more complex than that of the
leaves — for intermediate values of the coupling and of the parameter $a$ —
and regions where the opposite is realized. For example, for the $y$ coupling,
$C_{\rm hub}<\langle C_{\rm leaves}\rangle$ at the left and right side of the
yellow region.

Figure 3: Evolution of the relative complexity $C_{j0}$ as a function of the
coupling strength $d$ for the hub (in black) and a leaf (in red) of a star
network of $y$ coupled (a)-(d) and $x$ coupled (e)-(h) Rössler oscillators.
(a,e) $a=0.26$, (b,f) $a=0.31$, and (c,g) $a=0.38$ and $N=16$. (d,h) $a=0.31$
and $N=32$. Insets in all panels show the phase synchronization $S_{j}$ of the
hub (black) and the leaf (red) with the mean field. Other parameter values:
$b=0.2$, and $c=5.7$. One run per $d$-value.

To analyse with further detail the dynamical differentiation of hub and leaves
along the route to synchronization, we show in Fig. 3 cuts of the complexity
difference $C_{j0}=C_{j}-C_{0}$, with $j=\\{{\rm hub,leaf}\\}$ (black for the
hub and red for one leaf) along the three horizontal dashed lines shown in
Fig. 2(a). These curves clearly show that hub and leaf experience unalike
differentiation routes before either the star synchronizes at $d_{\rm c}$ to
the same dynamical state as in $d=0$ for $y$-coupling [Fig. 3(a)-(d)] or it
reaches a boundary crisis at $d_{\rm b}$ for the $x$ coupling [Fig. 3(e)-(h)].
In all cases, the maximal differentiation occurs always before for the hub
than for the leaf, that is, the maximum of the $C_{j0}$ curves is located at a
coupling that is lower for the hub than for the leaf. As it is shown at the
inset of each panel, this shift coincides with the earlier synchronization of
the hub to the mean field, here measured as

$S_{j}=\displaystyle\langle\operatorname{Re}(e^{{\rm
i}(\theta_{j}-\Phi)})\rangle_{t}\,,$ (8)

being $\langle\dots\rangle_{t}$ the time average, $\operatorname{Re}$ stands
for the real part, $\theta_{j}=\arctan\left(y_{j}/x_{j}\right)$ the phase of
the $j={\rm hub/leaf}$, and $\Phi$ the mean-field global phase. By definition,
$S_{j}\in[0,1]$, and approaches $1$ as the phase of the $j$-th oscillator is
locked to the phase of the mean field.

Figure 4: First-return maps to the Poincaré section ${\cal P}$ for the hub
(top row) and one of the leaves (bottom row) in a star network of $N=16$
$y$-coupled Rössler oscillators for $a=0.31$ and for different $d$ values. The
corresponding complexities $C_{\rm D}$ are reported. The nominal complexity is
$C_{0}=0.82$.

Another observed systematic behavior is that, for increasing values of $a$, as
the uncoupled dynamics is more developed (see Fig. 1), the relative complexity
of the hub diminishes while increases for the leaf. Therefore, the relative
position between the hub and leaf curves changes such that for $a=0.26$ the
hub’s curve is almost always above the leaf’s one, $C_{\rm hub,0}>C_{\rm
leaf,0}$, while for $a=0.38$, when the isolated dynamics is the most
developed, the hub’s curve is always below, $C_{\rm hub,0}<C_{\rm leaf,0}$.
Another salient feature is that in the very weak coupling regime ($d<0.05$),
just after oscillators start interacting, there is a marked fall of the
complexity with respect to $C_{0}$ for both hub and leaf. Finally, the impact
of the size of the star is analysed in Fig. 3(d,h) for $N=32$ and $a=0.31$
which has to be compared with the panels (b) and (f) of the same figure for
$N=16$. Precisely, the hub complexity is slightly affected and only in the
region of low coupling. As predicted by the MSF, for the $y$ coupling, the
threshold for synchronization $d_{c}$ does only depend on the smallest non-
zero eigenvalue $\lambda_{2}$, which equals $\lambda_{2}=1$ in both
configurations. On the other hand, for the $x$ coupling, the boundary crisis
is anticipated to smaller coupling values as the largest eigenvalue
$\lambda_{N}=N$ increases with the size.

To picture the different dynamical states the hub and the leaf are visiting,
we monitored their first-return maps at a specific points along the route to
synchronization as shown in Fig. 3(b) with vertical dashed lines. Having in
mind the map characterizing the uncoupled dynamics for $a=0.31$ and shown in
Fig. 1(d), for very low coupling [Fig. 4(A)], the leaf dynamics is nearly
unaffected while the hub is exhibiting the presence of incoherent small
perturbations and a shorter third branch. A slight increase in the coupling
strength [Fig. 4(B)] amplifies this effect turning the hub’s attractor into a
very thick unimodal map and leaving the leaf’s map almost unaltered. The next
scenario [Fig. 4(C)] corresponds to the largest deviation from the uncoupled
dynamics with both types of nodes displaying a much less complex dynamics. The
hub is constrained to a small neighborhood of the inner singular point as
revealed by the feeble cloud of points in the bisecting line. The leaf is
locked on a period-2 limit cycle. Beyond this drop of complexity, the hub
dynamics develops into a single huge cloud of points, displaced along the
bisecting line towards the period-1 limit cycle [Fig. 4(D)]. The large
complexity value attained by the hub $C_{\rm hub}=1.18>1$ indicates a very
disorganized dynamics far from the typical noisy limit cycle. At the same
time, the leaf moves into a first-return map with three thick branches due to
the noisy feedback from the hub. Pursuing along the route to synchronization
both hub and leaf each recover three branched maps slightly distorted until a
synchronous state is recovered for $d>d_{c}$ in Fig. 4(F) and all nodes share
the same dynamics as they had when uncoupled.

##  IV Star networks made of other dynamical systems

###  IV.1 Lorenz systems

Figure 5: (a)-(b) Relative complexity for the hub (black) and one of the
leaves (red) as a function of the coupling strength $d$ for a star network
with $N=32$ $y$-coupled Lorenz systems for two $R$-values. First-return maps
of the hub (top) and of one leaf (bottom) for the cases A, B and C indicated
by the blue dashed lines in the top panel. Other parameter values:
$\sigma=10$, $b=\frac{8}{3}$. A magnified view of the maps is plotted for the
case B: the range used is not the same between the hub and leaf.

Let us now consider the Lorenz 63 system Lorenz (1963) whose vector flow in
Eq. (1) is

$\mbox{f}(\mbox{x})=\left[\sigma(y-x),Rx-y-xz,-bz+xy\right]$ (9)

which is equivariant under a rotation symmetry around the $z$-axisLetellier
and Gilmore (2001). This global property is the main difference with respect
to the Rössler dynamics since, when the symmetry is modded out, the Lorenz
attractor is topologically equivalent to the Rössler one.Letellier, Dutertre,
and Gouesbet (1994); Letellier (1994); Letellier and Gilmore (2001) Lorenz
systems coupled through variable $y$ present a type-I synchronizability class
and, therefore, complete synchronization is stable above a critical coupling
threshold.

To compute the first-return map we proceeded as follows. The common Lorenz
attractor for $R=28$, resembling the two wings of a butterfly, is bounded by a
genus-3 torus and the Poincaré section is the union of the two components

${\cal P}_{\pm}\equiv\left\\{\left(\displaystyle
y_{n},z_{n}\right)\in\mathbb{R}^{2}~{}|~{}x_{n}=\pm
x_{\pm},\dot{x}_{n}\lessgtr 0\right\\}$ (10)

where $x_{\pm}=\pm\sqrt{b(R-1)}$.Letellier, Dutertre, and Gouesbet (1994);
Tsankov and Gilmore (2004); Letellier, Aguirre, and Maquet (2005) The interval
visited by each one of these components is normalized to the unit interval,
and the component ${\cal P}_{-}$ is shifted by $-1$, leading to variable
$\rho$. A quite close example of the first-return map for $R=28$ is shown at
the bottom of Fig. 5A featuring four branches paired due to the rotation
symmetry.Byrne, Gilmore, and Letellier (2004) The two increasing branches
correspond to the reinjection of the trajectory into the wing from which it is
issued, while the two decreasing branches are associated with the transition
from one wing to the other. The two left (right) branches correspond to the
$n$th intersection in the left (right) wing and to the ($n+1$)th intersection
in the left (right) or right (left) wing depending on the sign of the slope.

Figure 5 shows the relative complexity $C_{j0}$ of the hub and a leaf of a
$N=32$ star network $y$-coupled Lorenz systems for two values of the parameter
$R$. For $R=28$, full synchronization is reached for $d_{\rm c}\approx 2.22$.
As with the Rössler systems, the complexity of the hub and leaves drop below
$C_{0}$ in the weak coupling regime ($d\lesssim 1$) and beyond that point the
hub becomes more complex than the nominal dynamics before full synchronization
is eventually reached. For this parameter setting, the Lorenz hub experiences
a sudden drop in complexity as illustrated in the first-return map [Fig.
5(A)]: while the hub dynamics is characterized by small fluctuations around
the singular point in the centre of one of the wings, the leaves are nearly
unperturbed with their four branches. When $d$ is further increased [Fig.
5(B)], the leaf dynamics also collapses, with all nodes surprisingly
exhibiting a quasi-periodic dynamics, a regime not observed in an isolated
Lorenz system. Notice that the size of the two tori are different. When the
nodes’ dynamics have topologically equivalent attractors (here, tori) but with
a scaling factor, the phenomenon is known as amplitude enveloppe
synchronization.Gonzalez-Miranda (2002) Finally, beyond this point of reduced
dynamics and just before full synchronization [Fig. 5(C)], the interaction
between hub and leaves moves the dynamics of the leaves again to the four
paired branches although slightly more noisy than the uncoupled one, while the
hub is strongly perturbed with a very disorganized dynamics. This scenario
resembles the one observed for the Rössler system depicted in panel D of Fig.
4. Finally, we explored a second regime with an even more developed uncoupled
dynamics for $R=175$ whose first-return map has six branches (not shown). As
in Fig. 3(c) for the Rössler system, the switch between the relative
complexities of hub and leaf is no longer observed in Fig. 5 and the leaf
curve is always above $C_{\rm leaf,0}>C_{\rm hub,0}$.

###  IV.2 Mackey-Glass delay differential equation

Figure 6: (a) Bifurcation diagram for the Mackey-Glass delay differential
equation (11) versus the time delay $\tau$. (b-c) Bifurcation diagrams
computed versus the coupling strength $d$ for the hub and one of the leaves of
a star network of $N=32$ coupled Mackey-Glass equations (11) for $\tau=4.8$
[see red-dashed vertical line in panel (a)]. For avoiding the multi-stability
which is observed in such a network, the bifurcation diagram is computed
without a reset of the initial condition at each new value of the coupling
strength. Other parameter values are $\mu=1.2$ and $p=18.5$.

It is possible to produce a chaotic attractor characterized by a smooth
unimodal map via one-dimensional delay differential equations, for instance,
the Mackey-Glass (MG) equationMackey and Glass (1977); Glass and Mackey (1979)

$\dot{x}=\mu\frac{x_{\tau}}{1+x_{\tau}^{p}}-x_{t}$ (11)

where $x$ is the population of blood cells, $x_{t}=x(t)$ and
$x_{\tau}=x(t-\tau)$. This equation was initially proposed for the control of
hematopoiesis (the production of blood cells). Typically, the delay $\tau$ is
the time-scale for proliferation and maturation of these blood cells. The
dimension of the effective state space associated with a delay differential
equation is dependent on the delay.Gumowski (1974); Farmer (1982) Parameter
values for parameters $\mu$ and $p$ are such that a period-doubling cascade as
a route to chaos is observed [Fig. 6(a)]. We chose a delay $\tau=4.8$ (red
dashed line in Fig. 6(a) whose attractor is characterized by a smooth unimodal
map with equivalent topological properties to the Rössler system. Increasing
$\tau$ further leads to a much more complex behavior and computing a reliable
Poincaré section is rather tricky.

Our goal here is to investigate the route to synchronization for unimodal
dynamics produced by a potentially high-dimensional system. A network of MG
systems can be fully synchronized using bidirectional coupling.Pyragas (1998)

The attractor is bounded by a genus-1 torus and the single-component Poincaré
section can be defined as

${\cal P}_{\rm
MG}\equiv\left\\{x_{n},\in\mathbb{R}~{}|~{}\dot{x}_{n}=0.025,\ddot{x}_{n}>0,x_{n}<0.9\right\\}\,.$
(12)

The first-return map is one-dimensional and slightly foliated and quite
similar to the one shown for the leaf in Fig.7(A) for $N=16$ coupled MG
systems. Figure 7(a) shows the relative complexity of the hub and leaf as a
function of the coupling strength. As observed in the Rössler and Lorenz
systems, for low $d$-values, the hub experiences a sudden complexity drop
while the leaves keep almost their nominal dynamics. As the coupling is
increased, the hub starts synchronizing with the mean field and its relative
complexity rises above $0$. Beyond point B, all nodes reduce their complexity
down to a minimal value (point C) before finally increasing up to full
synchronization is reached. The noisy curves in that region is due to the
extreme sensitivity to initial conditions found in this system. The maps for
the hub and leaf [Figs. 7(A)-7(C)] corresponding to the setting points A,B,
and C marked in Fig. 7(a) sketch the dynamical differentiation route.

Figure 7: (a) Relative complexity as a function of the coupling strength $d$
for a star network of $N=16$ coupled Mackey-Glass delay differential equations
(11) for $\tau=4.8$. (A)-(C) First-return maps of the hub (top) and of one
leaf (bottom) for three different values of $d$. Other parameter values:
$\mu=1.2$ and $p=18.5$.

The bifurcation diagrams for the hub and leaf respectively are computed as a
function of the coupling $d$ [Figs. 6(b)-6(c)]. In both cases, an inverse
period-doubling cascade is observed up to a period-2 limit cycle is settled at
$d=0.18$. Curiously, the size of both limit cycles is different which, again,
is an example of an amplitude envelope synchronization. The crisis leading to
a larger chaotic attractor is strongly dependent on the initial conditions. In
summary, the prevalent lines of the route to synchronization previously
sketched are visible observed, albeit with some differences. Further
investigations are necessary to fully understand their origin.

###  IV.3 Four-dimensional Saito model

Figure 8: Relative complexity $C_{k0}$ as a function of the coupling strength
$d$ of a star network with $N=16$ $y$-coupled Saito models (13) for (a)
$\delta=0.58$ and (b) $\delta=0.82$. Poincaré sections for the hub and a leaf
for (A) $\delta=0.58$ and for (B) $\delta=0.82$, both with $d=0.07$. Other
parameter values: $\rho=14$, $\eta=1$, and $\epsilon=0.01$.

In this section, to confirm the generality of the nodal dynamical
differentiation route, we will investigate a completely different model of
nonlinear dynamics, the Saito model,Saito (1990) which is able to generate a
very rich dynamical behavior including toroidal chaos. The vector flow in Eq.
(1) reads

$\mbox{f}(\mbox{x})=\left[\rho(-y+z),x+2\delta
y,-x-w,(z-\Phi(w))/\epsilon\right]$ (13)

where $\mbox{x}=(x,y,z,w)$ and

$\Phi(w)=\left|\begin{array}[]{ccl}w-(1+\eta)&&w\geq\eta\\\\[2.84544pt]
\displaystyle-\frac{w}{\eta}&\mbox{ if }&|w|<\eta\\\\[8.5359pt]
w+(1+\eta)&&w\leq-\eta\,.\end{array}\right.$ (14)

This is a four-dimensional system involving a linear piecewise function as a
switch mechanism. The dynamics produced by this model can be chaotic, quasi-
periodic, toroidal chaotic, or even hyperchaotic, also structured around a
torus. Parameters $\rho$, $\eta$, and $\epsilon$ are fixed and $\delta$ is
used as a bifurcation parameter. Here we chose $\delta=0.58$ for toroidal
chaos, and $\delta=82$ to produce hyperchaotic toroidal chaos. Close examples
of the Poincaré sections of these attractors can be grasped, respectively, in
Fig. 8(A) and 8(B) for a leaf of a star of $N=16$ $y$-coupled Saito models
whose dynamics are almost identical to the uncoupled scenarios. The main
difference between these two dynamics is the “thickness” of the Poincaré
section being much thiner for $\delta=0.58$ than for $\delta=0.82$. The
hyperchaotic nature is revealed by the overlapping structures of the Poincaré
section as observed in the folded-towel map introduced by Rössler.Rössler
(1979) This is partly confirmed with the Lyapunov exponents which for
$\delta=0.58$ are

$\lambda_{1}=0.047>\lambda_{2}=0.012\approx\lambda_{3}=-0.022>\lambda_{4}=-94.79\,,$

with two null exponents as expected for toroidal chaos structured around a
torus T2,Letellier and Rössler (2020) and for $\delta=0.82$ are

$\lambda_{1}=0.164>\lambda_{2}=0.069\approx\lambda_{3}=-0.033>\lambda_{4}=-94.70\,,$

with two positive and one null exponents as needed for toroidal hyperchaos.
Note that in the latter case it is still unclear whether the second positive
exponent is merged or not with the third null exponent and this is an issue
currently under study. Investigations, which are out of the scope of the
present work, would allow determining whether the fourth dimension is required
for embedding the toroidal chaos ($\delta=0.58$). Being hyperchaotic for
$\delta=0.82$, the dynamics are necessarily four-dimensional.

When these models are coupled through the variable $y$, they synchronize as
shown in Figs. 8(a) and 8(b) with the convergence of their relative complexity
to zero, being the critical coupling $d_{\rm c}$ larger for the hyperchaotic
dynamics. Nevertheless, in both cases, what is preserved along the route to
synchronization with these toroidal chaotic and hyperchaotic dynamics is the
node differentiation mainly of the hub whose dynamics turns less developed
than the uncoupled one while the leaves sustain it over the whole range of
coupling strengths [compare the first-return maps in Fig. 8(A), and in Fig.
8(B)]. The route to synchronization appears very simple, most likely due to
the constrained toroidal structure of the nominal dynamics.

##  V Networks of Rössler oscillators

Having characterized the relationship between the degree centrality and
dynamical complexity in star networks, we move forward to generalize our
results to networks with a broader degree distribution. This issue was
partially tackled in Refs.Tlaie  _et al._ (2019b); Minati  _et al._ (2019),
where a strong correlation between degree $k$ and complexity $C_{D}$ allowed
establishing a node hierarchy. However, given the present novel results
revealing the important role played by the nominal nodal dynamics, we extend
our present study to larger networks.

Figure 9:  (a,b) Average $k$-class $S_{k}$ phase synchronization parameter for
several values of degree $k$ for $N=300$, $\langle k\rangle=4$ SF networks of
$y$-coupled Rössler systems when $a=0.26$ (left panels) and $a=0.31$ (right
panels), (c,d) averaged relative complexity $C_{k0}$. The vertical dotted
lines mark the couplings analyzed in Fig. 10 the lower panels: $d=0.04$
(black), $d=0.08$ (blue) and $d=0.22$ (magenta). (e,f) Relative complexity for
different $k$-classes as a function of the phase synchronization $S_{k}$.

We maximize the degree heterogeneity by using Barabasi-Albert scale-free (SF)
networks of $N$ identical $y$-coupled Rössler oscillators retaining the
parameter settings provided in Section III. Since we expect that nodes having
the same degree $k$ play equivalent roles in the network, we calculate the
evolution of $C_{k}$ within a degree class $k$ by averaging over the $N_{k}$
nodes having degree $k$, that is,

$C_{k}=\frac{1}{N_{k}}\sum_{[j|k_{j}=k]}C_{j}\,,$ (15)

where $C_{j}$ is the dynamical complexity of the $j$th node. Here, we use the
relative complexity $C_{k0}=C_{k}-C_{0}$ which helps to better assess the
effects of both the coupling and the topology in the complexity of the
$k$-class nodes. In addition, to evaluate the impact of the nodal environment,
we perform our calculations for networks made of $N=300$ nodes, wherein
$\langle k\rangle=4$. All the results are averaged over 10 different networks
realizations.

First, in Fig. 9(a-b) we plot the time averaged phase synchronization $S_{k}$
of $k$-class nodes with respect to the phase of the mean field for $N=300$ SF
networks as a function of the coupling $d$. We define $S_{k}$ as

$S_{k}=\frac{1}{N_{k}}\displaystyle\sum_{{[j|k_{j}=k]}}^{N}S_{j}\,.$ (16)

where $S_{j}$ is the time averaged phase synchronization of node $j$ with the
mean-field global phase, defined in Eq. (8). The result is then ensemble
averaged over 10 network realizations. Along the route to synchronization, the
$k$-classes synchronize hierarchically to the mean field before all of them
lock at the critical coupling ($d_{\rm c}=0.37$ for $a=0.26$ and $d_{\rm
c}=0.52$ for $a=0.31$),Zhou and Kurths (2006); Pereira (2010); Tlaie  _et al._
(2019b) therefore existing node differentiation also in heterogeneous
networks. While the phase synchronization $S_{k}$ is monotonously increasing
for most of the coupling range and classes, the weakly coupled regime presents
anomalous synchronizationBlasius, Montbrió, and Kurths (2003); Boaretto  _et
al._ (2018) over which most of the $S_{k}$ are below the basal, uncoupled
level. This anomalous range, more prominent for $a=0.31$, is associated with
the maximal node differentiation [compare Fig. 9(b) and 9(d)].

As already observed in Section III, the less developed the dynamics, the
smaller the critical value $d_{\rm c}$ [Fig. 9(c)-(d)]. In the weakly coupled
regime, the relative complexity $C_{k0}$ shows a strong node differentiation:
while the hubs present a markedly reduced complexity with respect to the
nominal value (with a clear negative minimal), for the smaller degrees the
complexity increases. This increment in the less connected nodes ($k=2$ and
$k=5$ in the example) is non-monotonous with $d$ and $k$, as observed in the
stars in Sections III and IV.

For stronger coupling all the nodes increase their complexity well above the
nominal value, ordered following the reverse degree ranking [Fig. 9(c)-(d)],
recovering the scenario observed for the hub in Section III. All these
features are qualitatively shared for the two different $a$-values, but the
deviations from the uncoupled value are larger for the more developed chaos,
$a=0.31$.

Plotting the relative complexity $C_{k0}$ as a function of the phase
synchronization $S_{k}$ reveals a dependency which is stronger for nodes with
a larger degree [Fig. 9(e)-(f)]. Typically, low-degree nodes exhibit a
dynamics which is nearly independent of the synchronization while it is the
opposite for large-degree nodes. It also clearly shows that the relative
complexity converges to zero for larger $d$-values when $k$ increases. This
delineates another signature of node differentiation.

Therefore, we conclude that in most of the regimes it is possible to correlate
the node degree with the relative complexity. Furthermore, the degree
centrality is the single structural parameter that affects the node behaviour,
while the rest of environmental topological features has no impact. This is
shown in Fig. 10, where we plot the value of the relative complexity $C_{k0}$
as a function of $k$ for three representative values of $d$, both for SF
networks and ER networks with $\langle k\rangle=4$. The ER and SF curves
overlap and, therefore, the dependence of $C_{k}$ on $a$ and $d$ are the same
regardless of the topology. This is quite remarkable since the ER and SF
networks have a different critical coupling $d_{\rm c}$: for a given
$d$-value, their global dynamics are different, but the nodes of degree $k$
have equivalent dynamics independently of the environment. The same result is
obtained for different sizes of both ER and SF networks (not shown).

Figure 10: Relative complexity $C_{k0}$ as a function of the degree $k$ for
(a) $a=0.26$, and (b) $a=0.31$. In each panel, curves correspond to SF (void
symbols) and ER (full symbols) with the coupling strengths marked in
Fig.10(c,d) with vertical lines: $d=0.04$ (circles), $d=0.08$ (squares), and
$d=0.22$ (triangles). Results are averaged over 10 different network
realizations of $N=300$. Other parameter values as in Fig. 1.

To better illustrate the node differentiation in larger networks, we plotted
the first-return maps for three different $k$-classes of nodes along the route
to synchronization (Fig. 11). Low-degree nodes ($k=2$) in SF networks produce
first-return maps whose thickness increases with the coupling strength, up to
$d<d_{\rm c}$ (top row in Fig. 11): for these nodes, the relative complexity
is always positive, and the first maps resemble the map of the uncoupled
dynamics but thicker (compare with the map for $a=0.31$ in Fig. 1). For large
degree nodes (third row in Fig. 11), once the minimum complexity is reached in
the weakly coupled regime, the complexity increases with the coupling strength
$d$. Around the minimum, the maps comprise a small cloud of points in the
neighborhood of the inner singular point (bottom left of the first-return
map); this is progressively transformed into a “noisy” period-1 limit cycle,
which is characterized by a cloud of points elongated perpendicularly to the
bisecting line and located around the centre of the map ($d=0.08$ with
$C_{27}=0.48$). Before the onset of synchronization, large degree nodes
produce a map which resembles the nominal one but is slightly thicker
($d=0.22$ with $C_{27}=0.88$). The nodes with an intermediary degree ($k=8$ in
the example) produce a map which has features intermediate between the two
extreme cases previously discussed: the node differentiation is, thus,
evidently correlated with the node degree. As previously discussed, in large
networks, node differentiation is mostly a monotonous function of the coupling
strength and of the degree. When SF networks are replaced by ER ones, the
degree range is narrower, but the maps for $k=2$ and $k=8$ are very similar to
their counterparts in SF networks. We conclude that node degree is clearly the
most important factor for the node differentiation along the route to
synchronization.

Figure 11: Different first-return maps to a Poincaré section for three
$k$-classes of nodes along a route to synchronization in a $N=300$ ($\langle
k\rangle=4$) SF network of $y$-coupled Rössler nodes ($a=0.31$ and other
parameters as in Fig. 1). The complexity $C_{\rm D}$ is reported in each case.

##  VI Conclusion

Routes to synchronization need to be elucidated to attain a better
understanding and knowledge of the possible scenarios that may be encountered
depending on the coupling, nominal node dynamics, topology, and network size.
Here, we confirmed and extended previous work depicting a non-trivial effect
of connectedness on node dynamics, particularly the existence of a non-
monotonic relationship between the complexity of node dynamics and coupling
strength. There is, indeed, a node differentiation that evolves with the
coupling strength. Typically, when the coupling function provides type-i
synchronizability, low coupling strengths induce a significant reduction in
the complexity of the large-degree nodes, while those having a small degree
are left nearly unaffected. Consequently, increasing the $d$-value, all the
nodes in small networks or those with a large degree in large networks present
a minimal dynamical complexity. In every network, all nodes have a complexity
that increases, often reaching a greater level than the nominal one, before
the onset of full synchronization. This sketch for the route to
synchronization is clearly observed for the two three-dimensional systems
(Rössler and Lorenz). With more complex node dynamics (Mackey-Glass and
Saito), the decrease towards a minimal complexity is only observed in the hub
of star networks; further studies with other network types are still needed
for attaining a more general view of the latter system. The node
differentiation is not a monotonic function of the coupling nor of the
synchrony. When the coupling provides a type-ii synchronizability, the route
to synchronization is an abridged version of the route observed with a type-i
synchronizability.

## Acknowledgments

ISN and IL acknowledge financial support from the Ministerio de Economía,
Industria y Competitividad of Spain under project FIS2017-84151-P.

## References

  * Newman (2003) M. E. J. Newman, “The structure and function of complex networks,” SIAM Review 45, 167–256 (2003).
  * Boccaletti  _et al._ (2006) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, “Complex networks: Structure and dynamics,” Physics Reports 424, 175–308 (2006).
  * Boccaletti  _et al._ (2014) S. Boccaletti, G. Bianconi, R. Criado, C. del Genio, J. Gómez-Gardñes, M. Romance, I. Sendiña-Nadal, Z. Wang, and M. Zanin, “The structure and dynamics of multilayer networks,” Physics Reports 544, 1–122 (2014).
  * Arenas, Díaz-Guilera, and Pérez-Vicente (2006) A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, “Synchronization processes in complex networks,” Physica D 224, 27–34 (2006).
  * Rodrigues  _et al._ (2016) F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, “The Kuramoto model in complex networks,” Physics Reports 610, 1–98 (2016).
  * Kuramoto and Battogtokh (2002) Y. Kuramoto and D. Battogtokh, “Coexistence of coherence and incoherence in nonlocally coupled phase oscillators,” Nonlinear Phenomena in Complex Systems 5, 380–385 (2002).
  * Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, “Chimera states for coupled oscillators,” Physical Review Letters 93, 174102 (2004).
  * Hagerstrom  _et al._ (2012) A. M. Hagerstrom, T. E. Murphy, R. Roy, P. Hövel, I. Omelchenko, and E. Schöll, “Experimental observation of chimeras in coupled-map lattices,” Nature Physics 8, 658–661 (2012).
  * Pecora  _et al._ (2014) L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy, “Cluster synchronization and isolated desynchronization in complex networks with symmetries,” Nature Communications 5, 4079 (2014).
  * Gómez-Gardenes  _et al._ (2011) J. Gómez-Gardenes, S. Gómez, A. Arenas, and Y. Moreno, “Explosive synchronization transitions in scale-free networks,” Physical Review Letters 106, 128701 (2011).
  * Leyva  _et al._ (2013) I. Leyva, I. Sendiña-Nadal, J. Almendral, A. Navas, S. Olmi, and S. Boccaletti, “Explosive synchronization in weighted complex networks,” Physical Review E 88, 042808 (2013).
  * Restrepo, Ott, and Hunt (2006) J. G. Restrepo, E. Ott, and B. R. Hunt, “Synchronization in large directed networks of coupled phase oscillators,” Chaos 16, 015107 (2006).
  * Gómez-Gardeñes, Moreno, and Arenas (2007) J. Gómez-Gardeñes, Y. Moreno, and A. Arenas, “Paths to synchronization on complex networks,” Physical Review Letters 98, 034101 (2007).
  * Pereira, van Strien, and Tanzi (2020) T. Pereira, S. van Strien, and M. Tanzi, “Heterogeneously coupled maps: hub dynamics and emergence across connectivity layers,” Journal of the European Mathematical Society 22, 2183–2252 (2020).
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Physical Review Letters 80, 2109–2112 (1998).
  * Barahona and Pecora (2002) M. Barahona and L. M. Pecora, “Synchronization in small-world systems,” Physical Review Letters 89, 054101 (2002).
  * Pecora (2008) L. M. Pecora, “Synchronization of oscillators in complex networks,” Pramana 70, 1175–1198 (2008).
  * Tlaie  _et al._ (2019a) A. Tlaie, I. Leyva, R. Sevilla-Escoboza, V. P. Vera-Avila, and I. Sendiña Nadal, “Dynamical complexity as a proxy for the network degree distribution,” Physical Review E 99, 012310 (2019a).
  * Tlaie  _et al._ (2019b) A. Tlaie, L. M. Ballesteros-Esteban, I. Leyva, and I. Sendiña-Nadal, “Statistical complexity and connectivity relationship in cultured neural networks,” Chaos, Solitons & Fractals 119, 284–290 (2019b).
  * Minati  _et al._ (2019) L. Minati, H. Ito, A. Perinelli, L. Ricci, L. Faes, N. Yoshimura, Y. Koike, and M. Frasca, “Connectivity influences on nonlinear dynamics in weakly-synchronized networks: Insights from Rössler systems, electronic chaotic oscillators, model and biological neurons,” IEEE Access 7, 174793–174821 (2019).
  * Eroglu  _et al._ (2020) D. Eroglu, M. Tanzi, S. van Strien, and T. Pereira, “Revealing dynamics, communities, and criticality from data,” Physical Review X 10, 021047 (2020).
  * Martínez  _et al._ (2018) J. H. Martínez, M. E. López, P. Ariza, M. Chavez, J. A. Pineda-Pardo, D. López-Sanz, P. Gil, F. Maestú, and J. M. Buldú, “Functional brain networks reveal the existence of cognitive reserve and the interplay between network topology and dynamics,” Scientific reports 8, 10525 (2018).
  * Letellier, Leyva, and Sendiña Nadal (2020) C. Letellier, I. Leyva, and I. Sendiña Nadal, “Dynamical complexity measure to distinguish organized from disorganized dynamics,” Physical Review E 101, 022204 (2020).
  * Letellier (2006) C. Letellier, “Estimating the Shannon entropy: Recurrence plots versus symbolic dynamics,” Physical Review Letters 96, 254102 (2006).
  * Bandt and Pompe (2002) C. Bandt and B. Pompe, “Permutation entropy: A natural complexity measure for time series,” Physical Review Letters 88, 174102 (2002).
  * Riedl, Müller, and Wessel (2013) M. Riedl, A. Müller, and N. Wessel, “Practical considerations of permutation entropy,” European Physical Journal Specical Topics 222, 249–262 (2013).
  * Rössler (1976) O. E. Rössler, “An equation for continuous chaos,” Physics Letters A 57, 397–398 (1976).
  * Press, Teukolsky, and Flannery (1992) W. H. Press, S. A. Teukolsky, and W. T. V. B. P. Flannery,  _Numerical Recipes in C. The Art of Scientific Computing_ , 2nd ed. (Cambdridge University Press, Cambridge, New York, Port Chester, Melbourne, Sydney, 1992).
  * Letellier, Dutertre, and Maheu (1995) C. Letellier, P. Dutertre, and B. Maheu, “Unstable periodic orbits and templates of the Rössler system: Toward a systematic topological characterization,” Chaos 5, 271–282 (1995).
  * Huang  _et al._ (2009) L. Huang, Q. Chen, Y.-C. Lai, and L. M. Pecora, “Generic behavior of master-stability functions in coupled nonlinear dynamical systems,” Physical Review E 80, 036204 (2009).
  * Sendiña-Nadal, Boccaletti, and Letellier (2016) I. Sendiña-Nadal, S. Boccaletti, and C. Letellier, “Observability coefficients for predicting the class of synchronizability from the algebraic structure of the local oscillators,” Physical Review E 94, 042205 (2016).
  * Lorenz (1963) E. N. Lorenz, “Deterministic nonperiodic flow,” Journal of the Atmospheric Sciences 20, 130–141 (1963).
  * Letellier and Gilmore (2001) C. Letellier and R. Gilmore, “Covering dynamical systems: Two-fold covers,” Physical Review E 63, 016206 (2001).
  * Letellier, Dutertre, and Gouesbet (1994) C. Letellier, P. Dutertre, and G. Gouesbet, “Characterization of the Lorenz system, taking into account the equivariance of the vector field,” Physical Review E 49, 3492–3495 (1994).
  * Letellier (1994) C. Letellier,  _Caractérisation topologique et reconstruction des attracteurs étranges_ , Ph.D. thesis, University of Paris VII, Paris, France (1994).
  * Tsankov and Gilmore (2004) T. D. Tsankov and R. Gilmore, “Topological aspects of the structure of chaotic attractors in $\mathbb{R}^{3}$,” Physical Review E 69, 056206 (2004).
  * Letellier, Aguirre, and Maquet (2005) C. Letellier, L. A. Aguirre, and J. Maquet, “Relation between observability and differential embeddings for nonlinear dynamics,” Physical Review E 71, 066213 (2005).
  * Byrne, Gilmore, and Letellier (2004) G. Byrne, R. Gilmore, and C. Letellier, “Distinguishing between folding and tearing mechanisms in strange attractors,” Physical Review E 70, 056214 (2004).
  * Gonzalez-Miranda (2002) J. M. Gonzalez-Miranda, “Amplitude envelope synchronization in coupled chaotic oscillators,” Physical Review E 65, 036232 (2002).
  * Mackey and Glass (1977) M. C. Mackey and L. Glass, “Oscillation and chaos in physiological control systems,” Science 197, 287–289 (1977).
  * Glass and Mackey (1979) L. Glass and M. C. Mackey, “Pathological conditions resulting from instabilities in physiological control system,” Annals of the New York Academy of Sciences 316, 214–235 (1979).
  * Gumowski (1974) I. Gumowski, “Sensitivity of certain dynamic systems with respect to a small delay,” Automatica 10, 659–674 (1974).
  * Farmer (1982) J. D. Farmer, “Chaotic attractors of an infinite-dimensional dynamical system,” Physica D 4, 366–393 (1982).
  * Pyragas (1998) K. Pyragas, “Synchronization of coupled time-delay systems: Analytical estimations,” Physical Review E 58, 3067–3071 (1998).
  * Saito (1990) T. Saito, “An approach toward higher dimensional hysteresis chaos generators,” IEEE Transactions on Circuits and Systems 37, 399–409 (1990).
  * Rössler (1979) O. E. Rössler, “Chaos,” in  _Structural Stability in Physics_ , edited by G. Güttinger and H. Eikemeier (Springer, Berlin Heidelberg, 1979) pp. 290–309, proceedings of Two International Symposia on Applications of Catastrophe Theory and Topological Concepts in Physics Tübingen, May 2-6 and December 11-14, 1978.
  * Letellier and Rössler (2020) C. Letellier and O. E. Rössler, “Chaos: The world of nonperiodic oscillations,” (Springer, Cham, Switzerland, 2020) Chap. An updated hierarchy of chaos, pp. 181–203.
  * Zhou and Kurths (2006) C. Zhou and J. Kurths, “Hierarchical synchronization in complex networks with heterogeneous degrees,” Chaos 16, 015104 (2006).
  * Pereira (2010) T. Pereira, “Hub synchronization in scale-free networks,” Physical Review E 82, 036201 (2010).
  * Blasius, Montbrió, and Kurths (2003) B. Blasius, E. Montbrió, and J. Kurths, “Anomalous phase synchronization in populations of nonidentical oscillators,” Physical Review E 67, 035204 (2003).
  * Boaretto  _et al._ (2018) B. Boaretto, R. Budzinski, T. Prado, J. Kurths, and S. Lopes, “Neuron dynamics variability and anomalous phase synchronization of neural networks,” Chaos 28, 106304 (2018).


* Pandey and Srivastava (2020) Shalini Pandey and Jaideep Srivastava. 2020.  Rkt: Relation-aware self-attention for knowledge tracing. In _Proceedings of the 29th ACM International Conference on Information & Knowledge Management_. 1205–1214.
  * Pardos et al. (2013) Zachary A. Pardos, Yoav Bergner, Daniel T. Seaton, and David E. Pritchard. 2013.  Adapting Bayesian knowledge tracing to a massive open online course in edX. In _Proceedings of the 6th International Conference on Educational Data Mining, EDM_. International Educational Data Mining Society.
  * Pardos et al. (2014) Zachary A Pardos, Ryan SJD Baker, Maria OCZ San Pedro, Sujith M Gowda, and Supreeth M Gowda. 2014\.  Affective States and State Tests: Investigating How Affect and Engagement during the School Year Predict End-of-Year Learning Outcomes.  _Journal of Learning Analytics_ 1, 1 (2014), 107–128.
  * Pardos and Heffernan (2010a) Zachary A. Pardos and Neil T. Heffernan. 2010a.  Modeling Individualization in a Bayesian Networks Implementation of Knowledge Tracing. In _Proceedings of the 18th International Conference on User Modeling, Adaptation, and Personalization_ _(UMAP)_. Berlin, Heidelberg, 255–266.
  * Pardos and Heffernan (2010b) Zachary A. Pardos and Neil T. Heffernan. 2010b.  Modeling Individualization in a Bayesian Networks Implementation of Knowledge Tracing. In _User Modeling, Adaptation, and Personalization, 18th International Conference, UMAP, Big Island, HI, USA_. 255–266.
  * Pardos and Heffernan (2011) Zachary A. Pardos and Neil T. Heffernan. 2011.  KT-IDEM: Introducing Item Difficulty to the Knowledge Tracing Model. In _User Modeling, Adaption and Personalization - 19th International Conference,UMAP, Girona, Spain_. 243–254.
  * Pavlik et al. (2009) Philip I. Pavlik, Hao Cen, and Kenneth R. Koedinger. 2009\.  Performance Factors Analysis - A New Alternative to Knowledge Tracing. In _Artificial Intelligence in Education: Building Learning Systems that Care: From Knowledge Representation to Affective Modelling, Proceedings of the 14th International Conference on Artificial Intelligence in Education, AIED , Brighton, UK_ _(Frontiers in Artificial Intelligence and Applications)_ , Vol. 200. IOS Press, 531–538.
  * Pelánek (2015a) Radek Pelánek. 2015a.  Modeling Students’ Memory for Application in Adaptive Educational Systems. In _Proceedings of the 8th International Conference on Educational Data Mining, EDM, Madrid, Spain_. International Educational Data Mining Society (IEDMS), 480–483.
  * Pelánek (2015b) Radek Pelánek. 2015b.  Modeling Students’ Memory for Application in Adaptive Educational Systems.  (2015), 480–483.
  * Piech et al. (2015) Chris Piech, Jonathan Bassen, Jonathan Huang, Surya Ganguli, Mehran Sahami, Leonidas J. Guibas, and Jascha Sohl-Dickstein. 2015.  Deep Knowledge Tracing. In _Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems, NIPS,Montreal, Quebec, Canada_. 505–513.
  * Psotka et al. (1988) Joseph Psotka, Leonard Daniel Massey, and Sharon A Mutter. 1988.  _Intelligent tutoring systems: Lessons learned_.  Psychology Press.
  * Qiu et al. (2011) Yumeng Qiu, Yingmei Qi, Hanyuan Lu, Zachary A. Pardos, and Neil T. Heffernan. 2011.  Does Time Matter? Modeling the Effect of Time with Bayesian Knowledge Tracing.. In _Proceedings of the 4th International Conference on Educational Data Mining, Eindhoven, The Netherlands_. www.educationaldatamining.org, 139–148.
  * Rasch (1993) Georg Rasch. 1993\.  _Probabilistic models for some intelligence and attainment tests._ ERIC.
  * Scarselli et al. (2009a) Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. 2009a.  The Graph Neural Network Model.  _IEEE Transactions on Neural Networks_ 20, 1 (2009), 61–80.
  * Scarselli et al. (2009b) Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. 2009b.  The Graph Neural Network Model.  _IEEE Transactions on Neural Networks_ 20, 1 (2009), 61–80.
  * Schmidhuber (2015) Jürgen Schmidhuber. 2015\.  Deep learning in neural networks: An overview.  _Neural Networks_ 61 (2015), 85–117.
  * Schodde et al. (2017) Thorsten Schodde, Kirsten Bergmann, and Stefan Kopp. 2017\.  Adaptive Robot Language Tutoring Based on Bayesian Knowledge Tracing and Predictive Decision-Making. In _Proceedings of the 2017 ACM/IEEE International Conference on Human-Robot Interaction_ _(HRI)_. Association for Computing Machinery, 128–136.
  * Sense et al. (2016) Florian Sense, Friederike Behrens, Rob R. Meijer, and Hedderik van Rijn. 2016.  An Individual’s Rate of Forgetting Is Stable Over Time but Differs Across Materials.  _Topics in cognitive science_ 8, 1 (2016), 305–321.
  * Settles and Meeder (2016) Burr Settles and Brendan Meeder. 2016.  A trainable spaced repetition model for language learning. In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics, ACL ,Berlin, Germany, Volume 1: Long Papers_. The Association for Computer Linguistics, 1848–1858.
  * Shin et al. (2021) Dongmin Shin, Yugeun Shim, Hangyeol Yu, Seewoo Lee, Byungsoo Kim, and Youngduck Choi. 2021\.  SAINT+: Integrating Temporal Features for EdNet Correctness Prediction. In _LAK’21: 11th International Learning Analytics and Knowledge Conference, Irvine, CA, USA_ , Maren Scheffel, Nia Dowell, Srecko Joksimovic, and George Siemens (Eds.). ACM, 490–496.
  * Stamper et al. (2010) J Stamper, A Niculescu-Mizil, S Ritter, GJ Gordon, and KR Koedinger. 2010.  Algebra i 2008–2009. challenge data set from kdd cup 2010 educational data mining challenge.  _Retrieved April_ 25 (2010), 2015.
  * Su et al. (2018) Yu Su, Qingwen Liu, Qi Liu, Zhenya Huang, Yu Yin, Enhong Chen, Chris H. Q. Ding, Si Wei, and Guoping Hu. 2018\.  Exercise-Enhanced Sequential Modeling for Student Performance Prediction. In _Proceedings of the Thirty-Second Conference on Artificial Intelligence, the 30th innovative Applications of Artificial Intelligence, New Orleans, Louisiana, USA_ , Sheila A. McIlraith and Kilian Q. Weinberger (Eds.). AAAI Press, 2435–2443.
  * Sutskever et al. (2014) Ilya Sutskever, Oriol Vinyals, and Quoc V Le. 2014\.  Sequence to Sequence Learning with Neural Networks. In _Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems,NIPS , Montreal, Quebec, Canada_. 3104–3112.
  * Sutton and Barto (2018) Richard S Sutton and Andrew G Barto. 2018.  _Reinforcement learning: An introduction_.  MIT press.
  * Sweeney et al. (2016) Mack Sweeney, Jaime Lester, Huzefa Rangwala, and Aditya Johri. 2016.  Next-Term Student Performance Prediction: A Recommender Systems Approach. In _Proceedings of the 9th International Conference on Educational Data Mining, EDM, Raleigh, North Carolina, USA_. International Educational Data Mining Society (IEDMS), 7.
  * Szabados (2010) Tamás Szabados. 2010\.  An elementary introduction to the Wiener process and stochastic integrals.  _arXiv preprint arXiv:1008.1510_ (2010).
  * Thai-Nghe et al. (2012) Nguyen Thai-Nghe, Lucas Drumond, Tomás Horváth, and Lars Schmidt-Thieme. 2012.  Using factorization machines for student modeling. In _Workshop and Poster Proceedings of the 20th Conference on User Modeling, Adaptation, and Personalization, Montreal, Canada_ _(CEUR Workshop Proceedings)_ , Vol. 872. CEUR-WS.org.
  * Tong et al. (2020b) Hanshuang Tong, Yun Zhou, and Zhen Wang. 2020b.  HGKT : Introducing Problem Schema with Hierarchical Exercise Graph for Knowledge Tracing.  _CoRR_ abs/2006.16915 (2020).
  * Tong et al. (2020a) Shiwei Tong, Qi Liu, Wei Huang, Zhenya Huang, Enhong Chen, Chuanren Liu, Haiping Ma, and Shijin Wang. 2020a.  Structure-based Knowledge Tracing: An Influence Propagation View. In _20th IEEE International Conference on Data Mining, ICDM,Sorrento, Italy_ , Claudia Plant, Haixun Wang, Alfredo Cuzzocrea, Carlo Zaniolo, and Xindong Wu (Eds.). 541–550.
  * Vardi (2012) Moshe Y Vardi. 2012\.  Will MOOCs destroy academia?  _Commun. ACM_ 55, 11 (2012), 5–5.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017\.  Attention is All you Need. In _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems,NIPS, Long Beach, CA, USA_. 6000–6010.
  * Vie and Kashima (2019) Jill-Jênn Vie and Hisashi Kashima. 2019.  Knowledge Tracing Machines: Factorization Machines for Knowledge Tracing. In _The Thirty-Third Conference on Artificial Intelligence, The Thirty-First Innovative Applications of Artificial Intelligence Conference, The Ninth Symposium on Educational Advances in Artificial Intelligence, Honolulu, Hawaii,USA_. AAAI Press, 750–757.
  * Villano (1992a) Michael Villano. 1992a.  Probabilistic Student Models: Bayesian Belief Networks and Knowledge Space Theory. In _Intelligent Tutoring Systems, Second International Conference, ITS, Montréal, Canada, Proceedings_ _(Lecture Notes in Computer Science)_ , Vol. 608. Springer, 491–498.
  * Villano (1992b) Michael Villano. 1992b.  Probabilistic Student Models: Bayesian Belief Networks and Knowledge Space Theory. In _Intelligent Tutoring Systems, Second International Conference, ITS, Montréal, Canada, Proceedings_. 491–498.
  * Wang et al. (2021) Chenyang Wang, Weizhi Ma, Min Zhang, Chuancheng Lv, Fengyuan Wan, Huijie Lin, Taoran Tang, Yiqun Liu, and Shaoping Ma. 2021.  Temporal Cross-Effects in Knowledge Tracing. In _Proceedings of the 14th ACM International Conference on Web Search and Data Mining_ _(WSDM)_. 517–525.
  * Wang and Liu (2021) Feng Wang and Huaping Liu. 2021.  Understanding the Behaviour of Contrastive Loss. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_. 2495–2504.
  * Wilson et al. (2016) Kevin H. Wilson, Yan Karklin, Bojian Han, and Chaitanya Ekanadham. 2016.  Back to the basics: Bayesian extensions of IRT outperform neural networks for proficiency estimation. In _Proceedings of the 9th International Conference on Educational Data Mining, EDM , Raleigh, North Carolina, USA_. International Educational Data Mining Society (IEDMS), 539–544.
  * Wright (1995) Raymond E. Wright. 1995\.  Logistic Regression.  In _Reading and understanding multivariate statistics_. American Psychological Association, Washington, DC, US, 217–244.
  * Wu et al. (2018) Lijun Wu, Fei Tian, Yingce Xia, Yang Fan, Tao Qin, Jian-Huang Lai, and Tie-Yan Liu. 2018.  Learning to Teach with Dynamic Loss Functions. In _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NIPS , Montréal, Canada_. 6467–6478.
  * Xiong et al. (2016) Xiaolu Xiong, Siyuan Zhao, Eric Van Inwegen, and Joseph Beck. 2016.  Going Deeper with Deep Knowledge Tracing. In _Proceedings of the 9th International Conference on Educational Data Mining, EDM, Raleigh, North Carolina, USA_ , Tiffany Barnes, Min Chi, and Mingyu Feng (Eds.). International Educational Data Mining Society (IEDMS),Raleigh, North Carolina, USA, 545–550.
  * Yang et al. (2020) Yang Yang, Jian Shen, Yanru Qu, Yunfei Liu, Kerong Wang, Yaoming Zhu, Weinan Zhang, and Yong Yu. 2020\.  GIKT: A Graph-based Interaction Model for Knowledge Tracing.  12457 (2020), 299–315.
  * Yeung and Yeung (2018) Chun-Kit Yeung and Dit-Yan Yeung. 2018.  Addressing two problems in deep knowledge tracing via prediction-consistent regularization. In _Proceedings of the Fifth Annual Conference on Learning at Scale, London, UK_ , Rose Luckin, Scott Klemmer, and Kenneth R. Koedinger (Eds.). ACM, 5:1–5:10.
  * Yin et al. (2019) Yu Yin, Qi Liu, Zhenya Huang, Enhong Chen, Wei Tong, Shijin Wang, and Yu Su. 2019.  QuesNet: A Unified Representation for Heterogeneous Test Questions. In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD , Anchorage, AK,USA_. 1328–1336.
  * Yudelson et al. (2013) Michael V Yudelson, Kenneth R Koedinger, and Geoffrey J Gordon. 2013.  Individualized bayesian knowledge tracing models. In _Artificial Intelligence in Education - 16th International Conference,(AIED), Memphis, TN, USA, Proceedings_. 171–180.
  * Zhai et al. (2019) Xiaohua Zhai, Avital Oliver, Alexander Kolesnikov, and Lucas Beyer. 2019.  S4L: Self-Supervised Semi-Supervised Learning. In _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_.
  * Zhang et al. (2017) Jiani Zhang, Xingjian Shi, Irwin King, and Dit-Yan Yeung. 2017\.  Dynamic Key-Value Memory Networks for Knowledge Tracing. In _Proceedings of the 26th International Conference on World Wide Web, WWW, Perth, Australia_. 765–774.
  * Zhu (2015) Xiaojin Zhu. 2015\.  Machine Teaching: An Inverse Problem to Machine Learning and an Approach Toward Optimal Education.  _Proceedings of the AAAI Conference on Artificial Intelligence_ 29, 1 (2015).