analyzed, which arises in the weak coupling regime when neurons that
individually exhibit non-chaotic behaviors are coupled. A network composed of
$1000$ neurons was simulated, firstly with the HH$\ell$ model, secondly with
the HR model. The bifurcation diagram of a coupled neuron shows that in the
lower coupling regime, the phase synchronization behavior occurs due to the
partial-regular behavior of the neurons. This phase synchronization signal is
lost with the increase of the coupling parameter since the synaptic current
gains relevance inducing chaoticity in the neurons. After that, for greater
values of coupling strength, the partial-regular behavior no longer exists,
and the neurons phase synchronize in a chaotic regime.

This phenomenon is explored with different approaches, by varying intrinsic
parameters of the model with the change of ionic conductances and disturbing
the neuron dynamics with external disturbances by applying an external pulsed
current. In both cases, it is possible to transition the neuronal dynamics of
one isolated neuron from periodic to chaotic, and consequently, change the
phase synchronization transition of a network composed of such neurons. In
this context, this bifurcation point can be achieved more easily with the
mutual change of more than one ion conductance, where small changes in one
specific ion-conductance may act as a catalyst to the second conductance
change [109]. On the other hand, the direct application of external currents
in the neuron can induce chaoticity in the regular dynamics, suppressing the
non-monotonic transition. It is important to mention that, for the parameters
used in this work, it was not possible to find a situation in which the pulsed
current applied to individually chaotic neurons induces a regular behavior in
neurons. Nevertheless, this possibility should not be discarded, as several
works use the application of external disturbances as a chaotic control
mechanism [111, 112, 113].

Also, it is important to notice that similar results can be obtained using the
other networks’ topologies of connections and other neuron models with the
same conditions, i.e., bursting behavior, and chemical coupling. Hence, this
correlation between the individual neuron dynamics and the collective behavior
of the network can be directly applied to the suppression (or optimization) of
phase synchronization of the network.

##  Chapter 7 Bistability in the synchronization of identical neurons

As shown in Chapter 2, nonlinear dynamical systems are known to exhibit a
complex behavior with the variation of a control parameter, going from a
stable periodic orbit to a chaotic state. Another interesting dynamical
behavior is called multistability, which is the coexistence of different
states for a given set of parameters [43, 60, 61]. This dynamical feature has
been studied for several years and it is observed in areas such as Physics
[114], Chemistry [115], climatology [116], and also neuroscience [117, 118].
In this sense, multistability can be associated with different states of the
brain and was proposed as a possible mechanism for memory storage and pattern
recognition [117, 118]. On the other hand, in simulated neuronal systems, this
coexistence can mean that a neuron may depict distinct stable states, with
different firing patterns, frequencies, regularity, and chaoticity [43, 44,
45, 46].

In this Chapter, it is studied how a bistable state can affect the
synchronization of a neural network. Firstly, it shows the existence of a
parameter region where the neuron exhibits bistable behavior, that is, a
neuron initialized with different initial conditions can present two different
stable states. Using the Lyapunov spectrum, it is detected that one state,
namely state I, is always periodic while the other state, namely state II,
depending on the parameters used, can exhibit periodic behavior or a chaotic
one. Furthermore, it is shown that state II is more sensitive to noise than
the state I, and the meantime that the system spends before escape from state
II to state I follow a Kramers law [119, 120]. After that, a network of
identical neurons is constructed using a generic mean-field electrical
coupling. If all the neurons are initialized in the periodic state, as shown
in Chapter 5 for identical oscillators, the network always reaches phase
synchronization, with more weakly-coupled networks needing more time to reach
the phase-synchronized state. But if there is bistability in the network, for
high coupling values the synchronization state is reached, but, it takes
longer compared to the case without bistability. On the other hand, in the
situations where both states are periodic, smaller values of coupling strength
are not able to make the network reach the phase synchronized state, and the
system depicts two groups with different synchronization features. In both
cases, interesting dynamic phenomena are observed such as chimera states and
anti-phase synchronization [121]. The results of this Chapter are published in
the article “Bistability in the synchronization of identical neurons.”
Physical Review E 104.2 (2021): 024204 [122].

###  7.1 Bistability on the HH$\ell$ model

The bistability in the HH$\ell$ model can be seen in Fig. 7.1. In panel (a) it
presents the bifurcation diagram of the IBI of one single neuron, using Eqs.
(3.18-3.22) for the same parameters of the Table 3.2, but varying the sodium
maximum conductance $g_{\mathrm{d}}$. The different colors refers to different
initial conditions
$\\{V(0),a_{\mathrm{d}}(0),a_{\mathrm{r}}(0),a_{\mathrm{sd}}(0),a_{\mathrm{sr}}(0)\\}$:
where the magenta dots represent $\\{-10,0,0,0,0.45\\}$, named IC-1, and the
cyan dots represent $\\{-70,0,0,0,0.45\\}$, named IC-2. It is noted that the
region $1.120\lesssim g_{\mathrm{d}}\lesssim 1.142$, the neuron is bistable,
since IC-1 (IC-2) leads to an upper (bottom) state called state I (state II).
The inner panel exhibits a magnification of the two states: while state I is
always periodic, the state II goes through a sequence of period doubling
bifurcation at $g_{\mathrm{d}}\gtrsim 1.1385$ and becomes chaotic due to
mechanisms discussed in Chapter 2.

Figure 7.1: Bistability on the HH$\ell$ model. (a) Bifurcation diagram of the
Inter-Burst-Intervals (IBI) as a function of $g_{\mathrm{d}}$, where two
different initial conditions are shown - cyan and magenta points. The
magnification in the inner panel exhibits the bistability region where state I
(magenta) and state II (cyan) are defined. The former is always periodic,
while the latter can depict periodic and chaotic behavior depending on the
$g_{\mathrm{d}}$ value (see inner panel). (b – c) Two-dimensional projection
$V(t)\times a_{\mathrm{sr}}(t)$ of the system’s phase portrait, at (b)
$g_{\mathrm{d}}=1.1350$ resulting in two periodic states, and (c)
$g_{\mathrm{d}}=1.1415$ a periodic (magenta) and chaotic (cyan) states are
detected. Each color corresponds to a distinct initial condition: magenta
(cyan) line with IC-1 (IC-2). The inner panel exhibits a magnification of the
chaotic attractor. (d) Maximum values of $a_{\mathrm{sr}}$ as a function of
$V(t=0)$ and $a_{\mathrm{sr}}(t=0)$ for
$a_{\mathrm{d}}(0)=a_{\mathrm{r}}(0)=a_{\mathrm{sd}}(0)=0$ using
$g_{\mathrm{d}}=1.1350$. The black (yellow) rectangle delimits each initial
condition regions to obtain the initialization at each state I (state II).

To proceed with the study of bistability, two values of conductances are
selected, $g_{\mathrm{d}}=1.1350$ (before the period doubling) and
$g_{\mathrm{d}}=1.1415$ (after the period doubling). Figure 7.1(b) depicts the
evolution of a two-dimensional projection of the phase portrait $V\times
a_{\mathrm{sr}}$ for $g_{\mathrm{d}}=1.1350$ where the magenta (cyan) line
represents the IC-1 (IC-2). The variations in the amplitude of
$a_{\mathrm{sr}}$ distinction the two states, while the state I
$\mathrm{max}(a_{\mathrm{sr}})\approx 0.47$, for state II
$\mathrm{max}(a_{\mathrm{sr}})\approx 0.43$. Figure 7.1(c)
$g_{\mathrm{d}}=1.1415$, it is observed that state II (cyan) depicts a greater
thickness, indicating chaotic dynamics. The inner panel exhibits a
magnification of this orbit, showing a projection of the chaotic attractor
[60]. Figure 7.1(d) displays the maximum values of $a_{\mathrm{sr}}$ in color
tones (from cyan $\approx 0.43$ to magenta $\approx 0.47$)
($g_{\mathrm{d}}=1.1350$) as a function of different initial conditions
$\\{V(0),0,0,0,a_{\mathrm{sr}}(0)\\}$, varying $V(0)$ and
$a_{\mathrm{sr}}(0)$. It is seen that the $\max(a_{\mathrm{sr}})$ can be used
to characterize each state, since the initial conditions leading to state I
(state II) are represented in magenta (cyan). The black and yellow rectangles
delimit the initial conditions used in this work to initialize the neurons in
each state: $(V(0),a_{\mathrm{sr}}(0))\in([-20,0],[0.40,0.48])$, for state I
and $(V(0),a_{\mathrm{sr}}(0))\in([-80,-60],[0.40,0.48])$ for II. It has to
mention that the results of panel (d) changes subtly with
$g_{\mathrm{d}}=1.1415$, but the rectangles still lead to their respective
states.

To verify the chaotic features of the states, Table 7.1 depicts the Lyapunov
spectrum for the two values of $g_{\mathrm{d}}$ and computed for both IC-1 and
IC-2. As mentioned in Chapter 2 the Lyapunov exponent computes the divergences
of orbits arbitrarily closed [59]. If at least one of the exponents
$\lambda_{i}>0$, it is an indication of chaos. For $g_{\mathrm{d}}=1.1350$
both initial conditions depict a non-chaotic behavior since the largest
Lyapunov exponent $\Lambda=\max\\{\lambda_{1},\cdots,\lambda_{5}\\}\approx 0$
while for $g_{\mathrm{d}}=1.1415$, the IC-2 produces a state with
$\Lambda\approx 10^{-4}>0$, evidencing its chaoticity. Besides this, the
smaller Lyapunov exponents indicate a difference in the stability of each
state, thereby showing that state I is more stable than state II, and state II
periodic is more stable than state II chaotic. In order to demonstrate the
convergence of the exponents, panels (a) and (b) of Fig. 7.2 present the time
evolution of the absolute value of the largest Lyapunov exponent ($\Lambda$)
for IC-1 (magenta) and IC-2 (cyan) after discarding $10^{5}$
$\text{\,}\mathrm{ms}$ to avoid transient effects. It is seen that, despite
the IC-2 with $g_{\mathrm{d}}=1.1415$ where $\Lambda$ is stable $\sim 10^{-4}$
(positive), for the other three cases $\Lambda$ converges asymptotically to
zero, exactly to zero if $t\rightarrow\infty$, denouncing the regular
dynamics.

Table 7.1: Lyapunov spectra for IC-1 and IC-2 computed for two values of $g_{\mathrm{d}}$. $g_{\mathrm{d}}=1.1350$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$
---|---|---|---|---|---
IC-1 | -0.000007 | -0.001657 | -0.102236 | -0.197086 | -5.466736
IC-2 | -0.000002 | -0.001059 | -0.122687 | -0.213510 | -5.420790
$g_{\mathrm{d}}=1.1415$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$
IC-1 | -0.000006 | -0.001818 | -0.099646 | -0.195736 | -5.471647
IC-2 | 0.000173 | -0.000017 | -0.122036 | -0.217536 | -5.418144
Figure 7.2: Stability of each state. Time evolution of the Largest Lyapunov
exponent for IC-1 and IC-2 (a) $g_{\mathrm{d}}=1.1350$ and (b)
$g_{\mathrm{d}}=1.1415$. Mean time that the system spends in state II before
escaping to the state I, as a function of the $1/\sigma^{2}$ ($\sigma$ is the
noise strength) considering $1000$ simulations and the threshold time of
$10^{8}$ $\text{\,}\mathrm{ms}$. The slope $U$ correspond to the height of the
potential barrier of a Kramers law
$\langle\tau(\sigma)\rangle\sim\exp(U/\sigma^{2})$. (c)
$g_{\mathrm{d}}=1.1350$ (periodic state) $U\approx 2.2\times 10^{-3}$. (d)
$g_{\mathrm{d}}=1.1415$ (chaotic state) $U\approx 2.0\times 10^{-6}$.

It is known that multistable systems are generally sensitive to noise [43,
120, 123, 124]. The mean time to perform a noise-induced transition from one
state to another gives an insight into the stability of the system. In this
sense, a current term is summed to Eq. (3.17), characterized by a noisy signal
which follows a random normal distribution with average $0$ and standard
deviation $\sigma$. It is noted that for the values of $\sigma$ and
considering a maximum simulation time of $10^{8}$ ms, state II switches to
state I only, and never the other way around. Panels (c) and (d) of Fig. 7.2
show the meantime that the system spends in state II before escaping to state
I as a function of the inverse of the noise strength $1/\sigma^{2}$
considering a $1000$ simulations. The time when a transition from state II to
state I occur is evaluated when $a_{\mathrm{sr}}$ crosses the threshold
$0.45$. The results show that the meantime to perform a transition follows a
Kramers law $\langle\tau(\sigma)\rangle\sim\exp(U/\sigma^{2})$ where $U$
corresponds to the height of the potential barrier [119, 120]. In panel (c)
$g_{\mathrm{d}}=1.1350$ (periodic state) $U\approx 2.2\times 10^{-3}$. In
panel (d) $g_{\mathrm{d}}=1.1415$ (chaotic state) $U\approx 2.0\times
10^{-6}$. Based on this test, the results indicate that the chaotic state II
switches more easily than the periodic state II.

###  7.2 Network properties

To explore the role of bistability in synchronization, is considered a neural
network composed of $N=100$ HH$\ell$ identical neurons. The evolution of the
membrane potential for the $i$-th neuron is given by

$C_{\mathrm{M}}\frac{dV_{i}}{dt}=-I_{i,\mathrm{d}}-I_{i,\mathrm{r}}-I_{i,\mathrm{sd}}-I_{i,\mathrm{sr}}-I_{i,\mathrm{l}}+I_{i,\mathrm{syn}},$
(7.1)

where the ionic currents are given by the Eqs. (3.18-3.22) and the synaptic
current is characterized by a mean-field coupling (all-to-all)

$I_{i,\mathrm{syn}}=\frac{\varepsilon}{N}\sum_{j=1}^{N}(V_{j}-V_{i})=\varepsilon(\langle
V\rangle-V_{i}),$ (7.2)

where $\varepsilon$ is the coupling (synaptic) strength measured in
$\text{\,}\mathrm{m}\mathrm{S}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ which, for
simplicity, is subsequently omitted. And

$\langle V\rangle=\frac{1}{N}\sum_{i=1}^{N}V_{i}$ (7.3)

is the mean-field of the network. This configuration is chosen because allows
the isolation of the bistability effect since $\varepsilon$ and $\langle
V\rangle$ are the same for all neurons, the only difference between them is
the initial conditions.

Moreover, the network is artificially subdivided in two groups
$\Omega_{1}=\\{1,\cdots,N_{1}\\}$ and $\Omega_{2}=\\{N_{1}+1,\cdots,N\\}$ with
$N_{1}$ and $N_{2}$ neurons, respectively. With this definition, it is
possible to adapt the Kuramoto order parameter (defined to the whole network)
to analyze the phase synchronization level of each group separately. To do
this, Eq. (5.5), for the $\ell$-th group is

$R_{\ell}=\left|\frac{1}{N_{\ell}}\sum_{j\in\Omega_{\ell}}e^{i\theta_{j}(t)}\right|,\;\;\;\;\;\ell=1,2.$
(7.4)

Finally, to measure if there is a difference in the synchronization of each
group is defined the absolute difference between $R_{1}$ and $R_{2}$:

$\Delta R=|R_{2}-R_{1}|.$ (7.5)

The initial conditions for the network simulations were selected to be in a
random position on the attractor of each state. To do so, an uncoupled neuron
is simulated for $10^{6}$ $\text{\,}\mathrm{ms}$ with random initial
conditions according to the rectangles of Fig. 7.1(d) (for each state). To
select the initial conditions for each neuron it is selected a random time
instant of the last $5\times 10^{5}$ $\text{\,}\mathrm{ms}$ of the simulation.
This approach avoids any initial synchronization bias ($R(t=0)\approx 1$) and
allows the initialization of each desired state. All simulations are performed
using a time limit of $10^{8}$ $\text{\,}\mathrm{ms}$.

###  7.3 Network results

The investigation of the role of the bistability in the network is made by
considering three different situations:

  * (i) - All neurons are in the state I (periodic) with $g_{\mathrm{d}}=1.1350$ (absence of bistability);

  * (ii) - Half of neurons of the network in each group ($N_{1}=N_{2}=N/2$), where neurons in $\Omega_{1}$ are initialized in state I (periodic) and $\Omega_{2}$ in the state II (chaotic) with $g_{\mathrm{d}}=1.1415$;

  * (iii) - Half of neurons of the network in each group ($N_{1}=N_{2}=N/2$), where neurons $\Omega_{1}$ are initialized in the state I (periodic) and $\Omega_{2}$ in the state II (periodic) with $g_{\mathrm{d}}=1.1350$.

In Fig. 7.3, it is evaluated the time-evolution of the Kuramoto order
parameter $R(t)$ for the three different values of coupling $\varepsilon$,
$10^{-3}$ (black line), $10^{-4}$ (red line), and $10^{-5}$ (blue line). Panel
(a) shows the result for the situation (i) (absence of bistability), for the
three values of $\varepsilon$ the system monotonically evolves from a non-
synchronized state $R\approx 0$ to a complete phase synchronized state $R=1$,
with the time needed to reach phase synchronization decreasing as the coupling
strength increases. A similar scenario is observed in panel (b), where
situation (ii) is studied, but the transition is not monotonic as in panel
(a). Panel (c) depicts the result for situation (iii) where the stable phase
synchronized state is achieved only for $\varepsilon=10^{-3}$, otherwise, for
$\varepsilon=10^{-5},\,10^{-4}$ the system stays in an oscillatory state where
$R$ vary between $0$ and $1$. This behavior is maintained for a long period as
observed in panel (d), which shows the last $10^{5}$ of the simulation. The
difference between the panels highlights that the bistability can influence
the final state of the system. While situation (i) the network depicts the
same behavior of periodic coupled oscillators, the existence of two different
states in the network makes things more complicated. To reach the synchronized
state, the coupling first needs to induce a transition in the neurons leading
them to a final unique state, with the same frequency, and only then, it can
align their phases. The details of these processes are analyzed in the next
sections.

Figure 7.3: The role of the bistability at the synchronization in a network of
identical HH$\ell$ neurons. Kuramoto order parameter $R(t)$ for three values
of coupling $\varepsilon=10^{-3}$ (black lines), $\varepsilon=10^{-4}$ (red
lines), and $\varepsilon=10^{-5}$ (blue lines). Three network configurations
are used: (a) - a network with all neurons in state $\mathrm{I}$; (b) - a
mixed network with half of neurons in state $\mathrm{I}$ (periodic) and half
in state $\mathrm{II}$ (chaotic); (c) - a mixed network with half of neurons
in state $\mathrm{I}$ (periodic) and half in state $\mathrm{II}$ (periodic).
(d) Magnification of the last $10^{5}$ times of panel (c).

####  7.3.1 Bistability with periodic-chaotic configuration (situation (ii))

In this section, it is considered a network with $g_{\mathrm{d}}=1.1415$ with
neurons $\Omega_{1}$ initialized in state I (periodic), and neurons of
$\Omega_{2}$ in state II (chaotic). As shown in Fig. 7.3 (b) the network
presents the same behavior for the three $\varepsilon$ values. In order to
study this situation it is used $\varepsilon=10^{-4}$ (red line). The results
are presented in Fig. 7.4 wherein panel (a), the black line in panel
corresponds to $R$ (whole network), the magenta and cyan lines represent the
local order parameters $R_{1}$ and $R_{2}$, respectively, and the orange line
presents $\Delta R$. As expected the network starts from the non-syncrhronized
state and, as the system evolves, $R_{1}$ increases slowly while $R_{2}$
remains close to $0$, leading to a local maximum of $\Delta R$. After the
first group reach the synchronized state $R_{1}=1$, $R_{2}$ quickly rises to
$1$. The colored arrows above panel (a) are three instants to be analyzed
separately in raster plots in panels (b – d). In this sense, each dot
corresponds to the beginning of a spike for each neuron, evaluated when
$V_{i}(t)$ reaches $-20$ $\text{\,}\mathrm{mV}$ (with a positive derivative).
Panel (b) shows the raster plot for a non-synchronized state of the network.
Panel (c) presents an interesting behavior where the first group is
synchronized ($R_{1}\approx 1$) while the second one does not ($R_{2}\approx
0$). This behavior is called chimera, where the system displays the
coexistence of one coherent-phase-locked group with an incoherent-non-
synchronized one [121]. In this case, the chimera is transient [125, 126] and
disappears when the second group gains synchronization, shown in panel (d),
where the network as a whole is synchronized.

Figure 7.4: Temporal evolution of the synchronization in the network at the
situation (ii) with $\mathbf{\bm{\varepsilon}=10^{-4}}$. (a) $R$, $R_{1}$,
$R_{2}$, and $\Delta R$ (black, magenta, cyan, and orange lines,
respectively). (b - d) Raster plots of the network where each dot corresponds
to the beginning of a spike. The colored arrows indicate the times when the
raster plots are obtained, which match with the dots’ colors.

####  7.3.2 Bistability with periodic-periodic configuration (situation (iii))

In this section, it is consider a network with $g_{\mathrm{d}}=1.1350$ with
neurons $\Omega_{1}$ initialized in the state I (periodic), and neurons of
$\Omega_{2}$ in the state II (periodic). As shown in panels (c) and (d) of
Fig. 7.3, the final state of the network depends on the $\varepsilon$. Firstly
it is studied the non-synchronizaed state where $\varepsilon=10^{-4}$ (red
line of Fig. 7.3 (c - d)). Figure 7.5 (a) depicts the synchronization features
of the network, where $R$, $R_{1}$, $R_{2}$, and $\Delta R$ are represented by
the black, magenta, cyan, and orange lines, respectively. The network goes
from the non-synchronized state ($R\approx 0$) and, as time evolves, the $R$
oscillates with increasing amplitude. On the other hand, $R_{1}$ and $R_{2}$
grows, but $R_{1}$ increases more quickly than $R_{2}$, leading to a local
maximum of $\Delta R$. After that, $R_{2}$ also increases, leading both
$R_{1}$ and $R_{2}$ close to $1$, but $R_{1}>R_{2}$. Panel (b) depicts the
last $10^{5}$ $\text{\,}\mathrm{ms}$ of the simulation, suggesting a beating
process of $R$.

Figure 7.5: Detailed analysis in the synchronization of situation (iii) with
$\bm{\varepsilon=10^{-4}}$. (a) $R$, $R_{1}$, $R_{2}$ and $\Delta R$ (black,
magenta, cyan, and orange lines) as a function of time. (b) Last $10^{5}$
$\text{\,}\mathrm{ms}$ of the simulation, where an oscillatory behavior of $R$
is evidenced while each group assumes phase synchronization separately. (c -
f) Raster plots obtained for the network in the times indicated by the colored
arrows.

In order to understand this oscillatory behavior of $R$ in Fig. 7.5, $4$ time
instants of the simulation are selected, represented by the colored arrows
above panels (a) and (b). Panels (c - f) represent the raster plots
corresponding to the colors of the arrows. First, at the panel (c) an
incoherent behavior is observed in the red dots. As time evolves, in panel (d)
the network depicts a transient chimera state since $\Delta R$ increases.
Eventually, each group depicts phase synchronization separately, since $R_{1}$
and $R_{2}$ show values close to one. However, each group evolves following
its frequency, which leads to momentary non-synchronization (e) and phase-
synchronization (f) (the oscillatory behavior of $R$).

The detailed analysis for the situation (iii) with $\varepsilon=10^{-3}$
(black line of Fig. 7.3 (c - d)) is presented in Fig. 7.6. The color scheme
follows the same one of the previous figure. The network starts in a non-
synchronized case and $R$ gains amplitude in an oscillatory way, reaching the
phase synchronized asymptotic state ($R=1$) as time evolves. While $R_{1}$
approaches monotonically the state of phase synchronization, $R_{2}$
approaches in an oscillatory manner. In addition, for $t\sim 2\times 10^{5}$
$\text{\,}\mathrm{ms}$, the $R_{2}$ momentarily loses synchronization and a
peak can be observed in $\Delta R\approx 0.9$. Three-time instants (colored
arrows) are selected to evaluate the raster plots. Panel (b) depicts the non-
synchronized case. Panel (c) depicts the raster plot for the network in the
time instant where the $R_{2}$ loses synchronization and $\Delta R$ assumes a
maximum characterizing a transient chimera state [125]. This behavior occurs
for a short period. Panel (d) shows the raster plot for the network when phase
synchronization is reached and maintained until the end of the simulation.
Different from the previous case, with $\varepsilon=10^{-4}$, here, the
coupling is strong enough to induce a transition in the network which yields a
phase synchronized state, the oscillatory behavior of the neurons in
$\Omega_{2}$ indicates that it undergoes a series of transitions to the state
of the $\Omega_{1}$.

Figure 7.6:  Detailed analysis in the synchronization of situation (iii) with
$\bm{\varepsilon=10^{-3}}$. (a) $R$, $R_{1}$, $R_{2}$, and $\Delta R$ (black,
magenta, cyan, and orange lines) as a function of time. (b - d) Raster plots
for the network where the dots’ colors match the arrows’ colors, which
indicate the analyzed times.

To understand the influence of the coupling strength $\varepsilon$ in the
final state of the network in Fig. 7.7 it is simulated $100$ different
initials conditions for the situation (iii). At the panel (a) is computed the
number of simulations $\mathcal{N}$ where there is a complete transition from
state II to state I and the entire network reaches phase synchronization. The
transition from state II to state I is recorded when $a_{\mathrm{sr}}$ crosses
the threshold 0.45 for the first time. For $\varepsilon>0.00030$ (dark-gray
area) the network depicts phase synchronization for all simulations
$\mathcal{N}=100$, for $0.00014<\varepsilon<0.00030$ (light-gray area), a
fraction of simulations induce the network to present phase synchronization
$0<\mathcal{N}<100$, and, at last, for $\varepsilon<0.00014$, no transitions
are observed $\mathcal{N}=0$. At panel (b) it is computed the time when the
last neuron initialized in state II transitions to state I, named
$\tau_{\mathrm{last}}$. In this case, the higher the coupling, the smaller the
time at which this happens.

Figure 7.7: A complete transition to state I depends on the initial
conditions. (a) Number of simulations where the network as a whole depicts
phase synchronization, meaning all neurons transition to state I, as a
function of $\varepsilon$. For $\varepsilon<0.00014$, no transitions are
observed. For $0.00014<\varepsilon<0.00030$ (light-gray area), a fraction of
simulations induce the network to present phase synchronization, and, at last,
for $\varepsilon>0.00030$ (dark-gray area), the network depicts phase
synchronization for all simulations. (b) The time when the last neuron
transitions to state I occur ($\tau_{\mathrm{last}}$) as a function of
$\varepsilon$. Here, $100$ simulations with different initial conditions are
performed.

An interesting behavior can be seen at the light-gray-area
$0.00014<\varepsilon<0.00030$, where the final state of the network depends on
the initial condition even for the same value of $\varepsilon$. In this sense,
in Fig. 7.8 it is studied the sensibility of initial conditions for
$\varepsilon=0.00016$. Panels (a – c) depict the time evolution of $R(t)$
(left scale) and the number of neurons in state I ($n_{1}(t)$, magenta) and
state II ($n_{2}(t)$, cyan) (right scale). Panel (a) presents an initial
condition where $R$ evolves to an oscillatory behavior and $n_{1}=n_{2}=50$
during all the simulation. Panel (b) depicts an initial condition where $R$
evolves to a phase synchronized state and at the end of simulation $n_{1}=100$
and $n_{2}=0$, a complete transition. An interesting behavior is described in
panel (c), for this initial condition, a partial transition to state I is
observed. In this case, $R$ oscillates between $0$ and $1$ as $t$ evolves, but
after a certain time, the oscillation becomes restricted between $0.5$ and
$1$. To understand this behavior it is selected a time instant (red arrow) to
analyze the phase synchronization and the spatiotemporal of the network. Panel
(d) shows $R$, $R_{1}$, $R_{2}$ and $\Delta R$, where it is observed that
neurons in $\Omega_{1}$ are phase synchronization while $R_{2}$ oscillates
between $0$ and $1$. The raster plot of the network (bursts only), depicted in
panel (e), shows the coexistence of two frequencies in the $\Omega_{2}$,
explaining the oscillations in $R_{2}$ and $R$. In this situation, chimera
states occur: they appear when the second group is non-synchronized, then
disappear when it synchronizes and reappears later again.

Figure 7.8: The dependence of synchronization features on the initial
conditions analyzed for $\bm{\varepsilon=0.00016}$. (a), (b), and (c) present
the Kuramoto order parameter for the entire network (left scale) and the
number of neurons in each state as a function of $t$ (right scale). Here,
different initial conditions are considered, in which no transition (a),
complete transition (b), and a partial transition (c) from state II to state I
are observed. (d) Synchronization features of the network for the case of
partial transition considering the time region indicated by the red arrow. (e)
Raster plot of the network (bursts only) for the same condition.

In Fig. 7.9 it is studied the meantime to reach a synchronized state ($R\geq
0.99$) $\langle\tau\rangle$ as a function of $\varepsilon$ considering $100$
initial conditions where different network configurations are considered,
varying the number of neurons initialized in each state $[N_{1},N_{2}]$. The
black dots correspond to the situation (i) (all neurons in the state I with
$g_{\mathrm{d}}=1.1350$) where the mean time decays with a power-law
$f(\varepsilon)\sim\varepsilon^{-1}$ represented by the black solid line. It
has to be mentioned that this is the same situation of periodic oscillators
studied in Chapter 5, the relaxation time is inversely proportional to the
coupling parameter. The squares represent configurations where state I is
periodic and state II is chaotic ($g_{\mathrm{d}}=1.1415$) in red [50:50]
(situation (ii)), purple [60:40], green [75,25], and cyan [90:10], and the
triangles both states are periodic ($g_{\mathrm{d}}=1.1350$) in in blue
[50:50] (situation (iii)), brown [60:40], gray [75:25], and yellow [90:10]. In
a general point of view the existence of bistability in the network increases
the time necessary to the system reach phase synchronization, greater is the
proportion of neurons in the state I closest is the result to the black line.
In particular, the absence of triangles on the left side of the figure can be
explained since for coupling values lower than a critical value no transition
is observed, consequently there is no phase synchronization to compute.

Figure 7.9:  Mean characteristic time $\bm{\langle\tau\rangle}$ needed for the
network to reach phase synchronization as a function of $\bm{\varepsilon}$ for
$\bm{100}$ initial conditions. The initialization of the network is made in
different configurations of states $[N_{1}:N_{2}]$. The black dots correspond
to all neurons in state I (situation (i)), following a power-law curve (solid
black line) $f(\varepsilon)\sim\varepsilon^{-1}$. The squares represent
configurations where state II is chaotic, in red [50:50] (situation (ii)),
purple [60:40], green [75:25], and cyan [90:10]. For the triangles, state II
is periodic, in blue [50:50] (situation (iii)), brown [60:40], gray [75:25],
and yellow [90:10].

###  7.4 Discussions

Throughout this chapter, the effects of bistability on the synchronization
process of neural networks were analyzed. The HH$\ell$ neuron presents a
bistability region when the sodium conductance is varied. It is shown that one
of the states, namely state I, is always periodic, while the other state,
namely state II depending on $g_{\mathrm{d}}$ can be periodic or chaotic. By
evaluating the Lyapunov spectrum and with noise application it was found that
state I is always more stable than state II, and state II periodic is more
stable than state II chaotic.

To understand the impact of bistability in the synchronization, it was built a
network of identical neurons with a global coupling, where the only difference
between neurons is their initial conditions. In this sense, three network
configurations have been considered: (i) all neurons in the state I; (ii) half
of the network in the state I and half in state II (chaotic); (iii) half of
the network in the state I and half in state II (periodic). The existence of
bistability in the network increases, on average, the time to reach a
synchronized state. This can be explained due to the fact that the phase
synchronization is only possible with a unique and final state. Hence, a
transition from the initial states to a final and unique state is mandatory.
The transitions occur from state II (periodic or chaotic) to state I
(periodic). If the state II is periodic small coupling values may not induce
the transition, resulting in a local synchronization in each of the groups
(not global synchronization), or even in a partial transition of neurons. At
last, using different fractions of neurons in each state, the time the network
takes to reach phase synchronization decreases, on average, as the coupling
strength increases and the behavior gets closer to configuration (i) as the
fraction of neurons initialized in the state I increase.

In summary, bistability plays an important role in the synchronization of
neural networks. The simple existence of two distinct stable states can lead
the network to different states of synchronization, depending on the
initialization of the system: from a non-synchronized state to a complete
phase synchronized state. Bistability also leads to the existence of a variety
of chimera states where the network displays the coexistence of one coherent-
phase-locked group with an incoherent-non-synchronized one [121]. These
synchronization and chimera states occur due to the difference in the
stability of the states of the uncoupled neurons, thus highlighting the
importance of the individual neuronal dynamics.

##  Chapter 8 Mechanism for explosive synchronization of neural networks

The phenomenon of explosive synchronization in oscillators was first studied
in chaotic networks by Gomez et al. [52], in which, for the particular case of
Kuramoto oscillators coupled by a network with scale-free topology [22], when
the natural frequency $\omega_{i}$ of each oscillator is given by the number
of connections $\omega_{i}=n_{i}$, the transition from the non-synchronized
state to the synchronized state occurs abruptly. This abrupt phase-variation
regarding the synchronization is named explosive synchronization. This
behavior occurs due to the existence of hubs (sites with a high degree of
connectivity), a characteristic of the topology of the scale-free network
[127]. From the point of view of neural networks, this effect is not
necessarily reproduced. This chapter explores how explosive synchronization
appears in a network of spiking neurons. It is considered an ensemble of
Chialvo neurons coupled in a network generated with the Newman-Watts route
[91]. For different values of the connection probability $p_{\mathrm{nw}}$
(which controls the number of shortcuts added in the network) it is possible
to found different synchronization transitions, in particular, the behavior of
explosive synchronization [52, 128, 129]. In addition to this abrupt phase
transition, a range of coupling strength values depicts a bistability
behavior, where for the same parameters there is a coexistence of both
synchronized and non-synchronized states. In the end, it is shown that the
dynamical mechanisms of this bistability are described by a saddle-node
bifurcation and a boundary crisis [60, 71, 72]. These results are published in
the article “Mechanism for explosive synchronization of neural networks.”
Physical Review E 100.5 (2019): 052301, [130].

###  8.1 Network properties

A network with $N=10000$ Chialvo neurons is considered. The neurons are
coupled using a complex network generated with the Newman-Watts route (with
$n_{0}=40000$ connections for $p_{\mathrm{nw}}=0$). The dynamics of the $i$-th
neuron is described as

$\displaystyle x_{i,t+1}$ $\displaystyle=$ $\displaystyle
x_{i,t}^{2}\exp(y_{i,t}-x_{i,t})+k_{i}+\frac{\varepsilon}{\bar{n}}\sum_{j=1}^{N}e_{i,j}x_{j,t},$
(8.1) $\displaystyle y_{i,t+1}$ $\displaystyle=$ $\displaystyle
ay_{i,t}-bx_{i,t}+c,$ (8.2)

where $x_{i,t}$ and $y_{i,t}$ are the activation and recovery variables,
respectively as defined in Chapter 3. $a=0.89$, $b=0.6$, and $c=0.28$ are
constant parameters and $k_{i}$ acts like an addictive disturbance in the
neuron that affects both the amplitude of oscillation and the frequency of
neurons. $\bar{n}$ is equivalent to the average number of connections on the
network, $e_{i,j}$ the element of the connection matrix, and $\varepsilon$ is
the coupling parameter. As shown in Chapter 3, for this set of parameters this
model depicts periodic spikes, different from chaotic oscillators, the
synchronization state of periodic oscillators is achieved for any
$\varepsilon>0$. This synchronization can be avoided with the introduction of
a dissimilitude parameter $\sigma=0.001$. In this sense, for each neuron, the
parameter $k_{i}$ is randomly selected between $[0.03,\,0.03+\sigma]$ imposing
a network with non-identical neurons.

To compute the occurrence of the spikes of each neuron, and further define a
phase $\theta_{i}$ in the Chialvo model, a threshold is defined with a value
of $x_{0}=0.5$, which means that a spike occurs when $x_{i,t}$ crosses $x_{0}$
(with positive derivative). After computing all the spikes of all the neurons,
with Eq. (5.11), it is possible to evaluate $\theta_{i}$, and hence, the phase
synchronization with Eq. (5.12) using $t_{\mathrm{i}}=150000$ and
$t_{\mathrm{f}}=200000$. Figure 8.1 depicts at panel (a) the time-evolution of
$x_{i,t}$ for $\varepsilon=0$, where the red dashed line delimits $x_{0}$. At
panel (b) is computed the $\cos(\theta_{i})$, at each spike the phase is
$2\pi$ multiple and $\cos(\theta_{i})=1$.

Figure 8.1: Phase association in the Chialvo model. (a) Dynamics of the
Chialvo model for $a=0.89$, $b=0.6$, $c=0.28$, and $k_{i}=0.03$. The red
dashed line delimits $x_{0}=0.5$. (b) The cosine of the phase of the neuron,
where a spike begins $\theta_{i}$ is $2\pi$ multiple and $\cos(\theta_{i})=1$.

###  8.2 Results

The topological effects of the probability of connection $p_{\mathrm{nw}}$ are
presented in Fig. 8.2, by evaluating the Kuramoto order parameter $\langle
R\rangle$ as a function of $\varepsilon$ for different values of
$p_{\mathrm{nw}}$. To produce this results, it is considered a continuation
process as follows: for $\varepsilon=0$ a random initial conditions is
selected for the neurons in the network. After evaluating $\langle R\rangle$
for this $\varepsilon$, the coupling is increased adiabatically from
$\varepsilon$ to $\varepsilon+\delta\varepsilon$ (with
$\delta\varepsilon=0.001$). For this new $\varepsilon$, the initial conditions
are the final conditions of the previous $\varepsilon$, without rebooting the
system. The increment is made until $\varepsilon=0.5$. After that, it is made
the adiabatic decrease ($\delta\varepsilon=-0.001$) until $\varepsilon=0$.
Therefore, each $\varepsilon$ value is studied two times, in the forward
direction (positive variations) and backward direction (negative variations).
If at each $\varepsilon$ value the system was restarted, i.e. without the
continuation method previously described, there would be a probability that
the chosen initial condition would fall into one of two possible states. In
this situation, it has been found that the backward direction is more likely
than the forward direction.

Figure 8.2: The role of the topology in the synchronization behavior. $\langle
R\rangle$ as a function of the coupling $\varepsilon$ for a fixed
$\sigma=0.001$ and $5$ $p_{\mathrm{nw}}$ values. The coupling strength is
evolved adiabatically ($\delta\varepsilon=0.001$) in two different directions:
forward (increment) and backward (decrement). (a) For $p_{\mathrm{nw}}=0.10$,
the network does not depict phase synchronization. (b) For a larger value of
$p_{\mathrm{nw}}=0.15$, a large interval is observed where explosive
synchronization is likely to occur, the vertical dotted (dashed) line
represents the transition points of $\varepsilon^{*}$
($\varepsilon^{\dagger}$). (c) and (d) for $p_{\mathrm{nw}}=0.25$ and
$p_{\mathrm{nw}}=0.35$, explosive synchronization still may exist but for a
smaller range of $\varepsilon$. (e) For larger values of
$p_{\mathrm{nw}}=0.45$ the transition from unsynchronized to phase-
synchronized state is smooth. The inner panel displays a magnification of the
transition region.

At panel (a) $p_{\mathrm{nw}}=0.1$ the network is kept in the non-synchronized
state in both directions due to the small number of non-local connections
which forbids the phase synchronization. Panel (b) $p_{\mathrm{nw}}=0.15$, for
the forward direction (cyan triangles), as the coupling is increased there is
a non-synchronized region ($\langle R\rangle\approx 0$) which is kept until
the critical coupling $\varepsilon=\varepsilon^{\dagger}$ when the network
transition abruptly to the synchronized state ($\langle R\rangle\approx 1$).
On the other hand, for the backward direction (magenta triangles), the network
starts at the synchronized state, which is maintained until
$\varepsilon=\varepsilon^{*}\neq\varepsilon^{\dagger}$. Hence, the forward and
backward directions are different, creating a hysteretic behavior that
characterizes a bistability region delimited
$\varepsilon^{*}<\varepsilon<\varepsilon^{{\dagger}}$, where at this range it
is possible to achieve both synchronized and non-synchronized states depending
on the initial conditions. Panels (c) and (d) where $p_{\mathrm{nw}}=0.25$ and
$0.35$, respectively, present the same behavior of panel (b), but for distinct
values of critical couplings $\varepsilon^{*},\,\varepsilon^{\dagger}$. At
panel (e) $p_{\mathrm{nw}}=0.45$ the transition is not explosive, for this
case, the phase transition occurs smoothly in both directions, which can be
easily seen in the magnification at the inner panel.

The synchronization features caused by the topological changes can be
summarized with the definition of the hysteretic area

$S_{\mathrm{h}}(p_{\mathrm{nw}})=\left|\int_{\mathrm{forward}}\langle
R\rangle(\varepsilon,p_{\mathrm{nw}})\,d\varepsilon-\int_{\mathrm{backward}}\langle
R\rangle(\varepsilon,p_{\mathrm{nw}})\,d\varepsilon\right|.$ (8.3)

Figure 8.3 presents $S_{\mathrm{h}}$ as a function of $p_{\mathrm{nw}}$
considering $0\leq\varepsilon\leq 0.5$, the error bars are the standard
deviation over $50$ different simulations. As discussed in Fig. 8.2, it can be
identified three different regions. The region (I) $0\leq p_{\mathrm{nw}}\leq
0.1$ with $S_{\mathrm{h}}\approx 0$ since the network does not synchronize.
The region (II) $0.11\leq p_{\mathrm{nw}}\lesssim 0.4$ with $S_{\mathrm{h}}>0$
where the network presents its bistability behavior, and $S_{\mathrm{h}}$
decreases as $p_{\mathrm{nw}}$ increases. Region (III) $\varepsilon>0.4$ with
$S_{\mathrm{h}}=0$ different from region (I), $S_{\mathrm{h}}\approx 0$ due
the smoothly transition in both directions.

Figure 8.3: Numerical results for $\mathbf{S_{\mathrm{h}}}$ as a function of
$\mathbf{p_{\mathrm{nw}}}$. The vertical lines delimit three different
scenarios: (I) the network does not show synchronization and only a chaotic-
asymptotic state is observed; (II) the network depicts explosive
synchronization as the result of attraction basin changes. At the same time,
the network displays a hysteretic loop; (III) the network shows a smooth
transition from the unsynchronized state to the synchronized state. The error
bars are the standard deviation over $50$ samples.

To investigate the bistability region, Fig. 8.4 explores the neuron dynamics
and the synchronization behavior of the network for a fixed probability of
$p_{\mathrm{nw}}=0.15$, where it is computed the maxima of $y_{\mathrm{max}}$
variable of one arbitrary neuron (left scale), comparing with the $\langle
R\rangle$ (right scale) of the network. In the forward direction (a), the
dissimilitude among neurons allows a chaotic-non-synchronized state, until
$\varepsilon=\varepsilon^{\dagger}$, where the chaotic attractor is subtly
replaced by a stable periodic orbit for $\varepsilon>\varepsilon^{\dagger}$.
In the backward direction (b), the synchronized state is maintained from
$\varepsilon=0.5$ until $\varepsilon=\varepsilon^{*}$, after that, for
$\varepsilon<\varepsilon^{*}$ (inner panel) the coupling is not sufficiently
strong enough to keep the non-identical neurons in the same period,
consequently, the synchronized state is lost. Hence, abrupt transitions in the
phase synchronization occurs due to the chaotic transitions when the coupling
strength $\varepsilon$ reaches a critical value, $\varepsilon^{\dagger}$ in
the forward direction, and $\varepsilon^{*}$ in the backward direction. In the
next section, it is investigated the dynamical mechanisms of the chaotic
transitions of both directions.

Figure 8.4: The explosive synchronization occurs due to a chaotic transition.
Bifurcation diagram of the $y_{\mathrm{max}}$ of an arbitrary neuron (left
scale), $\langle R\rangle$ of the network as a function of $\varepsilon$
(right scale). (a) The coupling is adiabatically increased, the chaotic-non-
synchronized state is kept until $\varepsilon=\varepsilon^{\dagger}\approx
0.4$ where the chaotic attractor loses stability to a periodic orbit. (b) The
coupling is adiabatically decreased (from $\varepsilon=0.5$ to $0$), the
periodic-synchronized state exists for
$\varepsilon\approx\varepsilon^{*}\approx 0.01$ (inner panel) in which the
non-synchronized state is reached.

####  8.2.1 The dynamical mechanism of the loss of stability of the chaotic
attractor

The physical mechanism responsible for the abrupt transition in the forward
direction, where the network changes abruptally from the non-synchronized to
the synchronized state which occurs for
$\varepsilon=\varepsilon^{\dagger}(p_{\mathrm{nw}})$, can be described in
terms of a route to chaos called boundary crisis [60, 71, 72]. As described in
Chapter 2, a boundary crisis is characterized by the collision of a chaotic
attractor with an unstable periodic orbit [60]. In this sense, the attractor
loses stability when the control parameter reaches a critical point
$\varepsilon=\varepsilon^{{\dagger}}$. For
$\varepsilon>\varepsilon^{{\dagger}}$ the chaotic attractor no longer exists
but is replaced by a chaotic transient. If it is considered an initial
condition in the basin of attraction (which exists only for
$\varepsilon<\varepsilon^{{\dagger}}$) and increase the
$\varepsilon\gtrsim\varepsilon^{{\dagger}}$, for a period, the orbit is
similar to the chaotic attractor, however, after a chaotic transient, called
$\tau_{\mathrm{crisis}}$, the system loses this behavior [60]. The Fig. 8.5
illustrates this chaotic transient where it is calculated the mean-field of
the networks, defined as

$X_{t}=\frac{1}{N}\sum_{i=1}^{N}x_{i,t}.$ (8.4)

With the initial condition of the network being chosen within the chaotic
saddle ($\varepsilon\lesssim\varepsilon^{{\dagger}}$). In panel (a)
$\varepsilon\lesssim\varepsilon^{{\dagger}}$ the mean-field shows an
incoherent low amplitude behavior, related to the out-of-sync state of the
network. At panel (b) $\varepsilon\gtrsim\varepsilon^{{\dagger}}$ the chaotic
attractor no longer exists, and after a transient time,
$\tau_{\mathrm{crisis}}\approx 3\times 10^{4}$ the amplitude of the mean-field
increases, showing periodic behavior related to the synchronization state. In
panel (c) $\varepsilon>\varepsilon^{{\dagger}}$ (greater than in panel (b)),
the coupling value is already far from the critical point, making the
transient time extremely short $\tau_{\mathrm{crisis}}\approx 0.5\times
10^{4}$.

Figure 8.5: Illustration of the chaotic transient in a crisis. Time evolution
of the mean-field of the network $X_{t}=(1/N)\sum_{i=1}^{N}x_{i,t}$ with a
initial condition inside the chaotic saddle, where (a)
$\varepsilon\lesssim\varepsilon^{{\dagger}}$, (b)
$\varepsilon\gtrsim\varepsilon^{{\dagger}}$ and (c)
$\varepsilon>\varepsilon^{{\dagger}}$.

In order to evidence the post-crisis behavior, Fig. 8.6 presents how this
chaotic transient $\langle\tau_{\mathrm{crisis}}\rangle$ changes as a function
of $|\varepsilon-\varepsilon^{\dagger}|$. To do so, firstly an initial
condition is evolved from $\varepsilon=0$ until
$\varepsilon<\varepsilon^{\dagger}$, to ensure that the initial condition
belongs to the chaotic saddle. After that, for this initial condition
($\varepsilon<\varepsilon^{\dagger}$), it is computed how much time the system
spends until reach the periodic synchronized state ($R=0.9$) for values
$\varepsilon>\varepsilon^{\dagger}$. The $\langle\cdot\rangle$ represents the
average over $50$ distinct orbits initialized in the chaotic saddle, using a
representative value of $p_{\mathrm{nw}}=0.15$. In a crisis, it is expected a
power-law decay
$\langle\tau_{\mathrm{crisis}}\rangle\sim|\varepsilon-\varepsilon^{\dagger}|^{-\kappa}$,
where $\kappa$ is namely the crisis exponent, and $\kappa>0.5$ for
multidimensional systems [60, 71, 72, 131]. For sufficient close values
$|\varepsilon-\varepsilon^{\dagger}|<2\times 10^{-2}$ the power-law decay
occurs with $\kappa=0.85$. For $|\varepsilon-\varepsilon^{\dagger}|>2\times
10^{-2}$ a power-law decay is still expected, but due the numerical
impossibility of start the system exactly in the saddle point
$\varepsilon^{\dagger}$ the power-law decay is replaced by an exponential
decay.

Figure 8.6: Evidence of boundary crisis in the network. Average chaotic
transient $\langle\tau_{\mathrm{crisis}}\rangle$ as a function of
$|\varepsilon-\varepsilon^{\dagger}|$ for the system initialized at the
chaotic saddle $\varepsilon\lesssim\varepsilon^{\dagger}$ (but, iterated with
$\varepsilon>\varepsilon^{\dagger}$). The average is evaluated over $50$
distinct initialization. For values sufficient close to the critical point
$\varepsilon^{\dagger}$ the chaotic transient decays with a power-law
$\langle\tau_{\mathrm{crisis}}\rangle\sim|\varepsilon-\varepsilon^{\dagger}|^{-\kappa}$
with $\kappa=0.85$.

####  8.2.2 The dynamical mechanism of the loss of stability of the phase-
synchronized attractor

For the backward direction, the network transition from the synchronized
periodic state to a non-synchronized chaotic one, that is, the network
presents a periodic dynamics for
$\varepsilon<\varepsilon^{*}(p_{\mathrm{nw}})$, the periodic orbit loses
stability to a chaotic one due a saddle-node bifurcation [60]. In such
transition, for $\varepsilon\lesssim\varepsilon^{*}$ intermittency can be
detected in the traces of periodic orbit due to the quasi-stable character of
the synchronized state before the saddle-node bifurcation of the periodic-
synchronized state [60]. Figure 8.7 presents an illustration of the
intermittent chaotic interruptions for $\varepsilon<\varepsilon^{*}$. It is
calculated the maxima values of the $y_{i,t}=y_{\mathrm{max}}$ variable for
three random neurons in the network. At panel (a)
$\varepsilon>\varepsilon^{*}$, as expected, the three neurons present periodic
dynamics. In the other panels (b), (c), and (d), $\varepsilon<\varepsilon^{*}$
(in descending order), the neurons exhibit some escapes of the periodic
orbits. The greater the distance between $\varepsilon$ and $\varepsilon^{*}$,
the smaller is the time between the escapes. It is possible to define a
characteristic time of intermittent events $\tau_{\mathrm{int}}$ where,

$\lim_{\varepsilon\rightarrow\varepsilon^{*}}\tau_{\mathrm{int}}=+\infty.$

As discussed in Chapter 2, if the transition is characterized by a saddle-node
bifurcation, the occurrence of the escapes decays with a power-law exponent
[60] in which

$\langle\tau_{\mathrm{int}}\rangle\sim|\varepsilon-\varepsilon^{*}|^{-1/2}.$
(8.5)

Figure 8.7: Illustration of intermittent chaotic events. Maxima values of the
variable $y_{i,t}$ for three random neurons in the network. (a)
$\varepsilon>\varepsilon^{*}$ where $y_{\mathrm{max}}$ is constant for the
three neurons (periodic dynamics). (b), (c) and (d)
$\varepsilon<\varepsilon^{*}$, (in descending order), the neurons present some
escape events from the periodic orbit which occur intermittently.

To study the mean time of the intermittent escapes
$\langle\tau_{\mathrm{int}}\rangle$, it is computed the times in which
$y_{i,\mathrm{max}}$ of the $i$-th neuron deviates from two standard
deviations below the mean. In Fig. 8.8 it is presented the
$\langle\tau_{\mathrm{int}}\rangle$ value is an average for $50$ different
initial conditions, as function of $|\varepsilon-\varepsilon^{*}|$. The
magenta line represents the theoretic curve
$\propto|\varepsilon-\varepsilon^{*}|^{-1/2}$, characterizing a saddle node
bifurcation due to the intermittency transition of type I [60].

Figure 8.8: Evidence of saddle node in the network. Average intermittent time
of the saddle-node bifurcation $\langle\tau_{\mathrm{int}}\rangle$ as a
function of the distance from the critical coupling
$|\varepsilon-\varepsilon^{*}|$ for different $50$ initial conditions. The
magenta line represents the theoretic curve
$\propto|\varepsilon-\varepsilon^{*}|^{-1/2}$.

###  8.3 Discussions

Throughout this chapter, it was explored different features of synchronization
using a Newman-Watts complex network composed of non-identical spiking
neurons. The dissimilitude is necessary to avoid the phase synchronization
behavior detected in periodic neurons. While a network with lower connectivity
no longer reaches a synchronized state $p_{\mathrm{nw}}\leq 0.1$, as the
number of connections increases the synchronized state is achieved in
different ways. For $0.1<p_{\mathrm{nw}}\lesssim 0.4$ the transition is
characterized by an explosive synchronization, and $p_{\mathrm{nw}}>0.4$ there
is a continuous phase-transition scenario. In particular, the presence of the
explosive synchronization is accompanied by a hysteresis loop on the network
dynamics as the coupling parameter is adiabatically increased and decreased.
It was shown that the abrupt synchronization transitions are associated with
routes to chaos. The dynamical mechanisms for the bistability region, are
given in terms of a saddle-node bifurcation and a boundary crisis [60].

At last, different from the original findings of Gomez et al. [52] where the
explosive synchronization occurs due to the existence of hubs (sites with a
high degree of connectivity), a characteristic of the topology of the scale-
free network, the explosive synchronization can be achieved with a set of non-
identical periodic neurons coupled in a complex network (in a small-world
regime). It is concluded that this scenario is sufficient for the presence of
explosive synchronization and opens a new approach with application in neural
networks.

##  Chapter 9 Conclusions and future perspectives

Throughout this thesis, the cooperative behavior of dynamical systems which
mimics the membrane potential of neuronal cells has been studied, in
particular, it is explored how the individual dynamics of neurons affect the
global behavior of phase synchronization of the network. In this sense,
different approaches have been used with distinct neuronal models and coupling
schemes.

The main results of this thesis have been presented in Chapter 6 where the
emergence of phase synchronization in weakly-coupled neurons is explored. This
phenomenon occurs when HH$\ell$ bursting neurons that individually exhibit
periodic behavior are coupled. Despite the mechanisms of the non-monotonic
transitions are still unknown, one of the possible explanations for this
observation is due to the interplay of the neural activity patterns with the
synaptic current because the collective dynamics induced by the synaptic
currents are not strong enough to dismiss the periodicity of the neurons, but
is sufficient to lead the network to a phase synchronized state (local maximum
of phase synchronization). The phase synchronization is lost with the
increment of the coupling parameter since the chaoticity induced by the
synaptic current desynchronizes the network but is recovered for higher values
of coupling.

Particularly interesting, the results have shown a possibility to control the
phase synchronization in weakly-coupled neurons which occurs by varying
parameters of the model allowing the individual neuron dynamics to migrate
from a regular to chaotic behavior and, consequently, changing the entire
transition scenario of the phase synchronization even when just tiny changes
in one ion-conductance is performed. In this sense, the ion conductance’s
variations can be understood as the blocking or activation of ion channels,
and small changes in one specific ion conductance may act as a catalyst to the
second conductance change. This is one of the working mechanisms of drugs used
to normalize neuronal functioning to unhealthy neural behaviors. Such results
may be of great interest in researches on drugs to control individual and
collective behavior on the brain. This phenomenon has been explored with
different approaches, by varying different parameters of the model and
disturbing the neuron dynamics with the application of an external pulsed
current. To corroborate these findings, the same non-monotonic transition has
been found in another neuronal model under the same dynamical conditions, and
in complex networks with different topologies. The hypothesis that such an
effect can be found in networks of nonlinear oscillators models which not
reproduces the neuronal dynamics is not discarded.

Moreover, the same HH$\ell$ model presents a range of parameters where a
bistable state has been found. Hence, two identical neurons initialized with
distinct initial conditions depict different dynamical behavior. The role of
bistability in the phase synchronization of a network has been explored in
Chapter 7. The simplest coexistence of neurons in different states can lead
the network to different states of phase synchronization, which depends only
on the initialization of the system. In general, the existence of bistability
at the network delays the occurrence of phase synchronization, since phase
synchronization is only possible if the synaptic current is sufficiently
strong to induce a transition from the less-stable state to a final and unique
state. Otherwise, a local synchronization is achieved where neurons
initialized in their respective states may synchronize with neurons in the
same state preventing a global synchronized stable state. Furthermore, during
the synaptic-induced transitions, interesting dynamical phenomena are observed
such as chimera states and anti-phase synchronization.

Lastly, in Chapter 8, a discrete model of spiking neurons coupled in a small-
world network has been studied. The dissimilitude among neurons allows a
chaotic-non-synchronized regime which is kept with the adiabatic increment of
the coupling until reaches a critical coupling and abruptly synchronizes. On
the other hand, the adiabatic decrements of the coupling maintain the network
in a periodic-synchronized state which is lost for a different critical
coupling value, characterizing a hysteretic loop. This phenomenon, called
explosive synchronization, has been explained with Kuramoto oscillators
coupled in a scale-free topology of connections with the existence of hubs
(one of the characteristics of this topology). However, the heterogeneity
produced by the dissimilitude of the neurons is sufficient to achieve such a
transition. In conclusion, the dynamical mechanisms for the bistability region
are given in terms of a boundary crisis where the chaotic-non-synchronized
loses stability when the chaotic attractor collides with a non-stable fixed
point and the periodic-synchronized state loses stability to saddle-node
bifurcation.

The results of this thesis offer various directions for further research. It
will be worthwhile to describe the mechanisms of the non-monotonic
synchronization reported in neurons with individual non-chaotic behaviors but
sensitive to transition to chaotic behaviors when coupled. In this sense, the
use of the master stability functions [132] and other tools used in data
analysis [133, 134, 135] allow the achievement of a general description to
understand the dynamical properties of the synaptic current, which allows the
non-monotonic synchronization. Also, this may offer a way of approaching
similar observations in more general oscillators’ systems. On the other hand,
from the point of view of the bistability found in the neuron model, the
analysis of the network with noise application can present an interesting
result and will be explored in a future opportunity, since it is known that
the noise may disrupt the synchronized state. However, the noise can induce
transitions among neuronal states facilitating a possible state of phase
synchronization. Lastly, the explosive synchronization found in the non-
identical spiking neurons of Chialvo, future works will be devoted to the
generalization of this behavior to other spiking neuronal models and
oscillators to generalize this phenomenon.

## Bibliography

  * [1] C. Huygens  “Horologium oscillatorium”  Parisiis, France, 1673
  * [2] C. Huygens  “The pendulum clock”  In _Trans RJ Blackwell, The Iowa State University Press, Ames_ Springer, 1986
  * [3] L. M. Pecora and T. L. Carroll  “Synchronization in chaotic systems”  In _Physical Review Letters_ 64.8 APS, 1990, pp. 821
  * [4] A. Goldbeter  “Biochemical oscillations and cellular rhythms: the molecular bases of periodic and chaotic behaviour”  Cambridge University Press, 1997
  * [5] A. Pikovsky, M. Rosenblum and J. Kurths  “Synchronization: a universal concept in nonlinear sciences”  Cambridge University Press, 2003
  * [6] J. Buck and E. Buck  “Synchronous fireflies”  In _Scientific American_ 234.5 JSTOR, 1976, pp. 74–85
  * [7] T. J. Walker  “Acoustic synchrony: two mechanisms in the snowy tree cricket”  In _Science_ 166.3907 American Association for the Advancement of Science, 1969, pp. 891–894
  * [8] M. K. McClintock  “Menstrual synchrony and suppression”  In _Nature_ 229.5282 Springer, 1971, pp. 244–245
  * [9] J. Jalife  “Mutual entrainment and electrical coupling as mechanisms for synchronous firing of rabbit sino-atrial pace-maker cells.”  In _The Journal of Physiology_ 356.1 Wiley Online Library, 1984, pp. 221–243
  * [10] J. Fell and N. Axmacher  “The role of phase synchronization in memory processes”  In _Nature Reviews Neuroscience_ 12.2 Nature Publishing Group, 2011, pp. 105–118
  * [11] L. Glass  “Synchronization and rhythmic processes in physiology”  In _Nature_ 410.6825 Nature Publishing Group, 2001, pp. 277–284
  * [12] E. R. Kandel, J. H. Schwartz, T. M. Jessell, S. A. Siegelbaum and A. J. Hudspeth  “Principles of Neural Science”  McGraw-hill New York, 2013
  * [13] W. Klimesch  “Memory processes, brain oscillations and EEG synchronization”  In _International Journal of Psychophysiology_ 24.1-2 Elsevier, 1996, pp. 61–100
  * [14] P. R. Roelfsema, A. K. Engel, P. König and W. Singer  “Visuomotor integration is associated with zero time-lag synchronization among cortical areas”  In _Nature_ 385.6612 Nature Publishing Group, 1997, pp. 157–161
  * [15] E. Rodriguez, N. George, JP. Lachaux, J. Martinerie, B. Renault and F. J. Varela  “Perception’s shadow: long-distance synchronization of human brain activity”  In _Nature_ 397.6718 Nature Publishing Group, 1999, pp. 430–433
  * [16] J. F. Cavanagh, M. X. Cohen and J. J. B. Allen  “Prelude to and resolution of an error: EEG phase synchrony reveals cognitive control dynamics during action monitoring”  In _Journal of Neuroscience_ 29.1 Soc Neuroscience, 2009, pp. 98–105
  * [17] F. Mormann, K. Lehnertz, P. David and C. E. Elger  “Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients”  In _Physica D: Nonlinear Phenomena_ 144.3 Elsevier, 2000, pp. 358–369
  * [18] C. Hammond, H. Bergman and P. Brown  “Pathological synchronization in Parkinson’s disease: networks, models and treatments”  In _Trends in Neurosciences_ 30.7 Elsevier, 2007, pp. 357–364
  * [19] O. V. Popovych and P. A Tass  “Control of abnormal synchronization in neurological disorders”  In _Frontiers in Neurology_ 5 Frontiers, 2014, pp. 268
  * [20] I. Dinstein, K. Pierce, L. Eyler, S. Solso, R. Malach, M. Behrmann and E. Courchesne  “Disrupted neural synchronization in toddlers with autism”  In _Neuron_ 70.6 Elsevier, 2011, pp. 1218–1225
  * [21] M. D. Greicius, G. Srivastava, A. L. Reiss and V. Menon  “Default-mode network activity distinguishes Alzheimer’s disease from healthy aging: evidence from functional MRI”  In _Proceedings of the National Academy of Sciences_ 101.13 National Acad Sciences, 2004, pp. 4637–4642
  * [22] N. Boccara  “Modeling complex systems”  Springer Science & Business Media, 2010
  * [23] D. G. Green  “Emergent behavior in biological systems”  In _Complex systems: from biology to computation_ Citeseer, 1993, pp. 24–35
  * [24] E. M. Izhikevich  “Which model to use for cortical spiking neurons?”  In _IEEE Transactions on Neural Networks_ 15.5 IEEE, 2004, pp. 1063–1070
  * [25] S. H. Strogatz  “Exploring complex networks”  In _Nature_ 410.6825 Nature Publishing Group, 2001, pp. 268–276
  * [26] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno and C. Zhou  “Synchronization in complex networks”  In _Physics Reports_ 469.3 Elsevier, 2008, pp. 93–153
  * [27] F. A. Rodrigues, T. K. D. M. Peron, P. Ji and J. Kurths  “The Kuramoto model in complex networks”  In _Physics Reports_ 610 Elsevier, 2016, pp. 1–98
  * [28] A. Townsend, M. Stillman and S. H. Strogatz  “Dense networks that do not synchronize and sparse ones that do”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 30.8 AIP Publishing LLC, 2020, pp. 083142
  * [29] R. C. Budzinski, B. R. R. Boaretto, K. L. Rossi, T. L. Prado, J Kurths and S. R. Lopes  “Nonstationary transition to phase synchronization of neural networks induced by the coupling architecture”  In _Physica A: Statistical Mechanics and its Applications_ 507 Elsevier, 2018, pp. 321–334
  * [30] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado and S. R. Lopes  “Investigation of details in the transition to synchronization in complex networks by using recurrence analysis”  In _Mathematical and Computational Applications_ 24.2 Multidisciplinary Digital Publishing Institute, 2019, pp. 42
  * [31] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado and S. R. Lopes  “Phase synchronization and intermittent behavior in healthy and Alzheimer-affected human-brain-based neural network”  In _Physical Review E_ 99.2 APS, 2019, pp. 022402
  * [32] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado, R. L. Viana and S. R. Lopes  “Synchronous patterns and intermittency in a network induced by the rewiring of connections and coupling”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 29.12 AIP Publishing LLC, 2019, pp. 123132
  * [33] G. B. Ermentrout and D. H. Terman  “Mathematical foundations of neuroscience”  Springer Science & Business Media, 2010
  * [34] A. L. Hodgkin and A. F. Huxley  “Action potentials recorded from inside a nerve fibre”  In _Nature_ 144.3651 Nature Publishing Group, 1939, pp. 710–711
  * [35] A. L. Hodgkin and A. F. Huxley  “A quantitative description of membrane current and its application to conduction and excitation in nerve”  In _The Journal of Physiology_ 117.4 Wiley Online Library, 1952, pp. 500–544
  * [36] E. M. Izhikevich  “Dynamical systems in neuroscience”  MIT press, 2007
  * [37] D. R. Chialvo  “Generic excitable dynamics on a two-dimensional map”  In _Chaos, Solitons & Fractals_ 5.3-4 Elsevier, 1995, pp. 461–479
  * [38] N. F. Rulkov, M. M. Sushchik, L. S. Tsimring and H. D. I. Abarbanel  “Generalized synchronization of chaos in directionally coupled chaotic systems”  In _Physical Review E_ 51.2 APS, 1995, pp. 980
  * [39] R. FitzHugh  “Impulses and physiological states in theoretical models of nerve membrane”  In _Biophysical Journal_ 1.6 Elsevier, 1961, pp. 445–466
  * [40] J. L. Hindmarsh and R. M. Rose  “A model of neuronal bursting using three coupled first order differential equations”  In _Proceedings of the Royal Society of London B: Biological Sciences_ 221.1222 The Royal Society, 1984, pp. 87–102
  * [41] H. A. Braun, M. T. Huber, M. Dewald, K. Schäfer and K. Voigt  “Computer simulations of neuronal signal transduction: the role of nonlinear dynamics and noise”  In _International Journal of Bifurcation and Chaos_ 8.05 World Scientific, 1998, pp. 881–889
  * [42] E. M. Izhikevich  “Simple model of spiking neurons”  In _IEEE Transactions on neural networks_ 14.6 IEEE, 2003, pp. 1569–1572
  * [43] U. Feudel  “Complex dynamics in multistable systems”  In _International Journal of Bifurcation and Chaos_ 18.06 World Scientific, 2008, pp. 1607–1626
  * [44] J. Foss, A. Longtin, B. Mensour and J. Milton  “Multistability and delayed recurrent loops”  In _Physical Review Letters_ 76.4 APS, 1996, pp. 708
  * [45] M. Sainz-Trapága, C. Masoller, H. A. Braun and M. T. Huber  “Influence of time-delayed feedback in the firing pattern of thermally sensitive neurons”  In _Physical Review E_ 70.3 APS, 2004, pp. 031904
  * [46] J. Ma and J. Wu  “Multistability in spiking neuron models of delayed recurrent inhibitory loops”  In _Neural Computation_ 19.8 MIT Press, 2007, pp. 2124–2148
  * [47] K. Xu, J. P. Maidana, S. Castro and P. Orio  “Synchronization transition in neuronal networks composed of chaotic or non-chaotic oscillators”  In _Scientific Reports_ 8.1 Nature Publishing Group, 2018, pp. 1–12
  * [48] Y. Kuramoto  “Self-entrainment of a population of coupled non-linear oscillators”  In _International symposium on mathematical problems in theoretical physics_ , 1975, pp. 420–422  Springer
  * [49] T. L. Prado, S. R. Lopes, C. A. S. Batista, J. Kurths and R. L. Viana  “Synchronization of bursting Hodgkin-Huxley-type neurons in clustered networks”  In _Physical Review E_ 90 American Physical Society, 2014, pp. 032818
  * [50] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado and S. R. Lopes  “Detection of nonstationary transition to synchronized states of a neural network using recurrence analyses”  In _Physical Review E_ 96.1 APS, 2017, pp. 012320
  * [51] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado, J Kurths and S. R. Lopes  “Suppression of anomalous synchronization and nonstationary behavior of neural network under small-world topology”  In _Physica A: Statistical Mechanics and its Applications_ 497 Elsevier, 2018, pp. 126–138
  * [52] J. Gómez-Gardenes, S. Gómez, A. Arenas and Y. Moreno  “Explosive synchronization transitions in scale-free networks”  In _Physical Review Letters_ 106.12 APS, 2011, pp. 128701
  * [53] I. Leyva, R. Sevilla-Escoboza, J. M. Buldú, I. Sendina-Nadal, J. Gómez-Gardeñes, A. Arenas, Y. Moreno, S. Gómez, R. Jaimes-Reátegui and S. Boccaletti  “Explosive first-order transition to synchrony in networked chaotic oscillators”  In _Physical Review Letters_ 108.16 APS, 2012, pp. 168702
  * [54] P. Ji, T. K. DM. Peron, P. J. Menck, F. A. Rodrigues and J. Kurths  “Cluster explosive synchronization in complex networks”  In _Physical Review Letters_ 110.21 APS, 2013, pp. 218701
  * [55] Y. Zou, T. Pereira, M. Small, Z. Liu and J. Kurths  “Basin of attraction determines hysteresis in explosive synchronization”  In _Physical Review Letters_ 112.11 APS, 2014, pp. 114102
  * [56] V. Avalos-Gaytán, J. A. Almendral, I. Leyva, F. Battiston, V. Nicosia, V. Latora and S. Boccaletti  “Emergent explosive synchronization in adaptive complex networks”  In _Physical Review E_ 97.4 APS, 2018, pp. 042301
  * [57] H. A Braun, K. Schäfer, K. Voigt, R. Peters, F. Bretschneider, X. Pei, L. Wilkens and F. Moss  “Low-dimensional dynamics in sensory biology 1: thermally sensitive electroreceptors of the catfish”  In _Journal of Computational Neuroscience_ 4.4 Springer, 1997, pp. 335–347
  * [58] E. N. Lorenz  “Deterministic nonperiodic flow”  In _Journal of Atmospheric Sciences_ 20.2, 1963, pp. 130–141
  * [59] A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano  “Determining Lyapunov exponents from a time series”  In _Physica D: Nonlinear Phenomena_ 16.3 Elsevier, 1985, pp. 285–317
  * [60] E. Ott  “Chaos in dynamical systems”  Cambridge University Press, 2002
  * [61] S. H. Strogatz  “Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering (Studies in Nonlinearity)”  Westview Press, 2001
  * [62] J. B. Marion  “Classical dynamics of particles and systems”  Academic Press, 2013
  * [63] R. M. May  “Simple mathematical models with very complicated dynamics”  In _The Theory of Chaotic Attractors_ Springer, 2004, pp. 85–93
  * [64] B. E. Kendall and G. A. Fox  “Spatial structure, environmental heterogeneity, and population dynamics: Analysis of the coupled logistic map”  In _Theoretical Population Biology_ 54.1 Elsevier, 1998, pp. 11–37
  * [65] J. Miśkiewicz and M. Ausloos  “A logistic map approach to economic cycles.(I). The best adapted companies”  In _Physica A: Statistical Mechanics and its Applications_ 336.1-2 Elsevier, 2004, pp. 206–214
  * [66] M. Suneel  “Electronic circuit realization of the logistic map”  In _Sadhana_ 31.1 Springer, 2006, pp. 69–78
  * [67] A. J. Lichtenberg and M. A. Lieberman  “Regular and chaotic dynamics”  Springer Science & Business Media, 2013
  * [68] G. Benettin, L. Galgani, A. Giorgilli and J. M. Strelcyn  “Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory”  In _Meccanica_ 15.1 Springer, 1980, pp. 9–20
  * [69] M. J. Feigenbaum  “Quantitative universality for a class of nonlinear transformations”  In _Journal of Statistical Physics_ 19.1 Springer, 1978, pp. 25–52
  * [70] Y. Pomeau and P. Manneville  “Intermittent transition to turbulence in dissipative dynamical systems”  In _Communications in Mathematical Physics_ 74 Springer, 1980, pp. 189––197
  * [71] C. Grebogi, E. Ott and J. A. Yorke  “Critical exponent of chaotic transients in nonlinear dynamical systems”  In _Physical Review Letters_ 57.11 APS, 1986, pp. 1284
  * [72] C. Grebogi, E. Ott, F. Romeiras and J. A. Yorke  “Critical exponents for crisis-induced intermittency”  In _Physical Review A_ 36.11 APS, 1987, pp. 5365
  * [73] “Free pik repository, Designed by brgfx”, https://www.freepik.com/free-vector/stem-cell-diagram-white-background_2480958.htm
  * [74] “The Nobel Prize in Physiology or Medicine 1963” Accessed: 2020-05-16, https://www.nobelprize.org/prizes/medicine/1963/summary/
  * [75] P. R. Shorten and D. J. N. Wall  “A Hodgkin–Huxley model exhibiting bursting oscillations”  In _Bulletin of Mathematical Biology_ 62.4 Elsevier, 2000, pp. 695–715
  * [76] S. Postnova, C. Finke, W. Jin, H. Schneider and H. A. Braun  “A computational study of the interdependencies between neuronal impulse pattern, noise effects and synchronization”  In _Journal of Physiology-Paris_ 104.3-4 Elsevier, 2010, pp. 176–189
  * [77] U. Feudel, A. Neiman, X. Pei, W. Wojtenek, H. Braun, M. Huber and F. Moss  “Homoclinic bifurcation in a Hodgkin–Huxley model of thermally sensitive neurons”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 10.1 American Institute of Physics, 2000, pp. 231–239
  * [78] A. Destexhe, Z. F. Mainen and T. J. Sejnowski  “An efficient method for computing synaptic conductances based on a kinetic model of receptor binding”  In _Neural Computation_ 6.1 MIT Press, 1994, pp. 14–18
  * [79] A. A. Hagberg, D. A. Schult and P. J. Swart  “Exploring Network Structure, Dynamics, and Function using NetworkX”  In _Proceedings of the 7th Python in Science Conference_ , 2008, pp. 11 –15
  * [80] L. Euler  “Solutio problematis ad geometriam situs pertinentis”  In _Commentarii academiae scientiarum Petropolitanae_ , 1741, pp. 128–140
  * [81] G. Alexanderson  “About the cover: Euler and Königsberg’s Bridges: A historical view”  In _Bulletin of the american mathematical society_ 43.4, 2006, pp. 567–573
  * [82] D. B. West  “Introduction to graph theory”  Prentice hall Upper Saddle River, 2001
  * [83] J. Saramäki, M. Kivelä, J. Onnela, K. Kaski and J. Kertesz  “Generalizations of the clustering coefficient to weighted complex networks”  In _Physical Review E_ 75.2 APS, 2007, pp. 027105
  * [84] P. Erdös and A. Rényi  “On random graphs”  In _Publicationes mathematicae_ 6.26, 1959, pp. 290–297
  * [85] D. J. Watts and S. H. Strogatz  “Collective dynamics of ’small-world’ networks”  In _Nature_ 393.6684 Nature Publishing Group, 1998, pp. 440–442
  * [86] V. Latora and M. Marchiori  “Efficient behavior of small-world networks”  In _Physical Review Letters_ 87.19 APS, 2001, pp. 198701
  * [87] D. J. Watts  “Small worlds: the dynamics of networks between order and randomness”  Princeton University Press, 2004
  * [88] H. Ebel, J. Davidsen and S. Bornholdt  “Dynamics of social networks”  In _Complexity_ 8.2 Wiley Online Library, 2002, pp. 24–27
  * [89] T. B. Achacoso and W. S. Yamamoto  “AY’s Neuroanatomy of C. elegans for Computation”  CRC Press, 1991
  * [90] D. S. Bassett and E. D. Bullmore  “Small-world brain networks”  In _The Neuroscientist_ 12.6 Sage Publications Sage CA: Thousand Oaks, CA, 2006, pp. 512–523
  * [91] M. E. J. Newman and D. J. Watts  “Renormalization group analysis of the small-world network model”  In _Physics Letters A_ 263.4-6 Elsevier, 1999, pp. 341–346
  * [92] Q. K. Telesford, K. E. Joyce, S. Hayasaka, J. H. Burdette and P. J. Laurienti  “The ubiquity of small-world networks”  In _Brain connectivity_ 1.5, 2011, pp. 367–375
  * [93] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares and C. S. Zhou  “The synchronization of chaotic systems”  In _Physics Reports_ 366.1-2 Elsevier, 2002, pp. 1–101
  * [94] D. A. Wiley, S. H. Strogatz and M. Girvan  “The size of the sync basin”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 16.1 American Institute of Physics, 2006, pp. 015103
  * [95] S. W. Son, H. Jeong and H. Hong  “Relaxation of synchronization on complex networks”  In _Physical Review E_ 78.1 APS, 2008, pp. 016106
  * [96] M. V. Ivanchenko, G. V. Osipov, V. D. Shalfeev and J. Kurths  “Phase synchronization in ensembles of bursting oscillators”  In _Physical Review Letters_ 93.13 APS, 2004, pp. 134101
  * [97] O. V. Popovych and P. A. Tass  “Desynchronizing electrical and sensory coordinated reset neuromodulation”  In _Frontiers in Human Neuroscience_ 6 Frontiers, 2012, pp. 58
  * [98] W. Löscher and D. Schmidt  “New horizons in the development of antiepileptic drugs”  In _Epilepsy Research_ 50.1-2 Elsevier, 2002, pp. 3–16
  * [99] W. Löscher, H. Klitgaard, R. E. Twyman and D. Schmidt  “New avenues for anti-epileptic drug discovery and development”  In _Nature Reviews Drug Discovery_ 12.10 Nature Publishing Group, 2013, pp. 757–776
  * [100] J. S. Perlmutter and J. W. Mink  “Deep brain stimulation”  In _Annu. Rev. Neurosci._ 29 Annual Reviews, 2006, pp. 229–257
  * [101] C. A. S. Batista, S. R. Lopes, R. L. Viana and A. M. Batista  “Delayed feedback control of bursting synchronization in a scale-free neuronal network”  In _Neural Networks_ 23.1 Elsevier, 2010, pp. 114–124
  * [102] V. H. P. Louzada, N. A. M. Araújo, J. S. Andrade and H. J. Herrmann  “How to suppress undesired synchronization”  In _Scientific Reports_ 2.1 Nature Publishing Group, 2012, pp. 1–5
  * [103] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado and S. R. Lopes  “Synchronization domains in two coupled neural networks”  In _Communications in Nonlinear Science and Numerical Simulation_ 75 Elsevier, 2019, pp. 140–151
  * [104] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado, J. Kurths and S. R. Lopes  “Neuron dynamics variability and anomalous phase synchronization of neural networks”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 28.10 AIP Publishing, 2018, pp. 106304
  * [105] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado, J. Kurths and S. R. Lopes  “Protocol for suppression of phase synchronization in Hodgkin–Huxley-type networks”  In _Physica A: Statistical Mechanics and its Applications_ 528 Elsevier, 2019, pp. 121388
  * [106] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado and S. R. Lopes  “Suppression of phase synchronization in scale-free neural networks using external pulsed current protocols”  In _Mathematical and Computational Applications_ 24.2 Multidisciplinary Digital Publishing Institute, 2019, pp. 46
  * [107] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado and S. R. Lopes  “Temperature dependence of phase and spike synchronization of neural networks”  In _Chaos, Solitons & Fractals_ 123 Elsevier, 2019, pp. 35–42
  * [108] K. L. Rossi, R. C. Budzinski, B. R. R. Boaretto, T. L. Prado, U. Feudel and S. R. Lopes  “Phase-locking intermittency induced by dynamical heterogeneity in networks of thermosensitive neurons”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 31.8 AIP Publishing LLC, 2021, pp. 083121
  * [109] B. R. R. Boaretto, C. Manchein, T. L. Prado and S. R. Lopes  “The role of individual neuron ion conductances in the synchronization processes of neuron networks”  In _Neural Networks_ 137 Elsevier, 2021, pp. 97–105
  * [110] Zack B. and Helen B.-S.  “High frequency deep brain stimulation and neural rhythms in Parkinson’s disease”  In _Neuropsychology Review_ 25.4 Springer, 2015, pp. 384–397
  * [111] E. Ott, C. Grebogi and J. A. Yorke  “Controlling chaos”  In _Physical Review Letters_ 64.11 APS, 1990, pp. 1196
  * [112] T. Shinbrot, C. Grebogi, J. A. Yorke and E. Ott  “Using small perturbations to control chaos”  In _Nature_ 363.6428 Nature Publishing Group, 1993, pp. 411–417
  * [113] S. J. Schiff, K. Jerger, D. H. Duong, T. Chang, M. L. Spano and W. L. Ditto  “Controlling chaos in the brain”  In _Nature_ 370.6491 Nature Publishing Group, 1994, pp. 615–620
  * [114] F. T. Arecchi, R. Meucci, G. Puccioni and J. Tredicce  “Experimental evidence of subharmonic bifurcations, multistability, and turbulence in a Q-switched gas laser”  In _Physical Review Letters_ 49.17 American Physical Society, 1982, pp. 1217
  * [115] N. Ganapathisubramanian and K. Showalter  “Bistability, mushrooms, and isolas”  In _The Journal of Chemical Physics_ 80.9 American Institute of Physics, 1984, pp. 4177–4184
  * [116] D. Paillard  “The timing of Pleistocene glaciations from a simple multiple-state climate model”  In _Nature_ 391.6665 Nature Publishing Group, 1998, pp. 378–381
  * [117] John A. Hertz  “Introduction to the theory of neural computation”  CRC Press, 2018
  * [118] C. C. Canavier, D. A. Baxter, J. W. Clark and J. H. Byrne  “Nonlinear dynamics in a model neuron provide a novel mechanism for transient synaptic inputs to produce long-term alterations of postsynaptic activity”  In _Journal of Neurophysiology_ 69.6 American Physiological Society Bethesda, M. D., 1993, pp. 2252–2257
  * [119] P. Hänggi, P. Talkner and M. Borkovec  “Reaction-rate theory: fifty years after Kramers”  In _Reviews of Modern Physics_ 62.2 APS, 1990, pp. 251
  * [120] S. Kraut, U. Feudel and C. Grebogi  “Preference of attractors in noisy multistable systems”  In _Physical Review E_ 59.5 APS, 1999, pp. 5253
  * [121] D. M. Abrams and S. H. Strogatz  “Chimera states for coupled oscillators”  In _Physical Review Letters_ 93.17 APS, 2004, pp. 174102
  * [122] B. R. R. Boaretto, R. C. Budzinski, K. L. Rossi, C. Manchein, T. L. Prado, U. Feudel and S. R. Lopes  “Bistability in the synchronization of identical neurons”  In _Physical Review E_ 104.2 APS, 2021, pp. 024204
  * [123] P. R. Sharma, M. D. Shrimali, A. Prasad and U. Feudel  “Controlling bistability by linear augmentation”  In _Physics Letters A_ 377.37 Elsevier, 2013, pp. 2329–2332
  * [124] A. N. Pisarchik and U. Feudel  “Control of multistability”  In _Physics Reports_ 540.4 Elsevier, 2014, pp. 167–218
  * [125] A. Zakharova, M. Kapeller and E. Schöll  “Chimera death: Symmetry breaking in dynamical networks”  In _Physical Review Letters_ 112.15 APS, 2014, pp. 154101
  * [126] F. P. Kemeth, S. W. Haugland, L. Schmidt, I. G. Kevrekidis and K. Krischer  “A classification scheme for chimera states”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 26.9 AIP Publishing LLC, 2016, pp. 094815
  * [127] R. Albert and A.-L. Barabási  “Statistical mechanics of complex networks”  In _Reviews of Modern Physics_ 74.1 APS, 2002, pp. 47
  * [128] X. Zhang, S. Boccaletti, S. Guan and Z. Liu  “Explosive synchronization in adaptive and multilayer networks”  In _Physical Review Letters_ 114.3 APS, 2015, pp. 038701
  * [129] W. Zhou, J. Yang, L. Zhou and D. Tong  “Stability and Synchronization Control of Stochastic Neural Networks”  Berlin: Springer, 2015
  * [130] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado and S. R. Lopes  “Mechanism for explosive synchronization of neural networks”  In _Physical Review E_ 100.5 APS, 2019, pp. 052301
  * [131] G. T. Kubo, R. L. Viana, S. R. Lopes and C. Grebogi  “Crisis-induced unstable dimension variability in a dynamical system”  In _Physics Letters A_ 372.34 Elsevier, 2008, pp. 5569–5574
  * [132] L. M. Pecora and T. L. Carroll  “Master stability functions for synchronized coupled systems”  In _Physical Review Letters_ 80.10 APS, 1998, pp. 2109
  * [133] F. Chollet  “Deep Learning with Python”  Manning Publications Co., 2017
  * [134] C. Bandt and B. Pompe  “Permutation entropy: a natural complexity measure for time series”  In _Physical Review Letters_ 88.17 APS, 2002, pp. 174102
  * [135] N. Marwan, M. C. Romano, M. Thiel and J. Kurths  “Recurrence plots for the analysis of complex systems”  In _Physics Reports_ 438.5-6 Elsevier, 2007, pp. 237–329
  * [136] Y. Kuramoto  “Cooperative dynamics of oscillator communitya study based on lattice of rings”  In _Progress of Theoretical Physics Supplement_ 79 Narnia, 1984, pp. 223–240
  * [137] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort and R. Spigler  “The Kuramoto model: A simple paradigm for synchronization phenomena”  In _Reviews of Modern Physics_ 77.1 APS, 2005, pp. 137

##  Appendix A Kuramoto oscillators

The Kuramoto model consist of the ensemble of $N$ oscillators with phase
$\theta_{i}$, where the dynamics of each oscillator depends on a natural
frequency and a nonlinear coupling [136, 137]. The evolution of the phase of
each oscillator is described by

$\frac{d\theta_{i}}{dt}=\omega_{i}+\frac{1}{N}\sum_{j=1}^{N}\varepsilon_{ij}\sin(\theta_{j}-\theta_{i}),\;\;\;i=1,2,\cdots,N,$
(A.1)

in which $\omega_{i}$ is the natural frequency of each oscillator, and
$\varepsilon_{ij}$ is the parameter which controls the magnitude of the
coupling between the $i$-th and the $j$-th oscillator. In the particular case
of global coupling where $\varepsilon_{ij}=\varepsilon\;\forall\;i,j$,

$\frac{d\theta_{i}}{dt}=\omega_{i}+\frac{\varepsilon}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}),\;\;\;i=1,2,\cdots,N,$
(A.2)

where $\varepsilon$ is the intensity of the coupling. The natural frequencies
{$\omega_{i}$} are distributed following a probability density function
$\varkappa(\omega)$, which is assumed to be symmetric for a given frequency
$\Omega$, that is,

$\varkappa(\Omega-\omega)=\varkappa(\Omega+\omega).$ (A.3)

Kuramoto have defined an order parameter to quantify the phase synchronization
of the oscillators

$R(t)e^{i\psi(t)}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}(t)},$ (A.4)

where $R(t)$ is the modulus of the order parameters and $\psi(t)$ is the
average circular frequency of the oscillators. Multiplying both sides of Eq.
(A.4) by $e^{-i\theta_{i}}$

$Re^{i(\psi-\theta_{i})}=\frac{1}{N}\sum_{j=1}^{N}e^{i(\theta_{j}-\theta_{i})},$

and taking the imaginary part

$R\sin(\psi-\theta_{i})=\frac{1}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}).$
(A.5)

Replacing Eq.(A.5) at Eq. (A.2),

$\frac{d\theta_{i}}{dt}=\omega_{i}+\varepsilon
R\sin(\psi-\theta_{i}),\;\;\;i=1,2,\cdots,N,$ (A.6)

that is, the phase of each oscillator is attracted by the frequency $\psi$ and
the intensity of the coupling is proportional to the order parameter $R$.

The Eq. (A.4) can be rewritten as [137]

$Re^{i\psi}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}=\int_{-\pi}^{\pi}e^{i\theta}\left(\frac{1}{N}\sum_{j=1}^{N}\delta(\theta-\theta_{j})\right)d\theta.$
(A.7)

Considering infinity oscillators $N\rightarrow\infty$,

$Re^{i\psi}=\int_{-\pi}^{\pi}\int_{-\infty}^{+\infty}e^{i\theta}\rho(\theta,\omega,t)\varkappa(\omega)d\omega
d\theta.$ (A.8)

in which the system can be described in terms of the probability density
$\rho(\theta,\omega,t)$. Where $\rho(\theta,\omega,t)d\theta$ represents the
fraction of oscillators with frequency $\omega$ which are found between
$\theta$ and $\theta+d\theta$ in a given time instant $t$.

Since the oscillators are moving in a unitary circle $[-\pi,\pi]$ with angular
velocity $\dot{\theta_{i}}$ described by the Eq. (A.6), the probability
density $\rho$ obeys the continuity equation

$\frac{\partial\rho}{\partial
t}+\frac{\partial}{\partial\theta}(\upsilon\rho)=0,$ (A.9)

where $\upsilon$ is the angular velocity
$\upsilon=\dot{\theta}=\omega+\varepsilon R\sin(\psi-\theta)$. The probability
density presents the normalization condition described by

$\int_{-\pi}^{\pi}\rho(\theta,\omega,t)d\theta=1.$ (A.10)

One of the possible solutions for this system is to consider $R=0$. This case
corresponds to the incoherent solution where the infinite oscillators are
uniformly distributed around the circle $\rho=1/(2\pi)$. Thus, we say that for
$R=0$ we have complete desynchronization, in which for each $i$-th oscillator
with a phase $\theta_{i}=\theta$ there will be a $j$-th oscillator in anti-
phase $\theta_{j}=\theta+\pi$. For the case of complete synchronization, where
all oscillators have the same phase, that is, $\theta_{i}=\psi$, and,
consequently, the module $R=1$ due

$|R(t)|\cdot|e^{i\psi(t)}|=\left|\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}(t)}\right|=\left|\frac{1}{N}(e^{i\theta_{1}(t)}+e^{i\theta_{2}(t)}+\cdots+e^{i\theta_{N}(t)})\right|=\frac{1}{N}\cdot
N|e^{i\psi(t)}|=1.$

To study the cases where $0<R<1$, the stationary solution of the Eq. (A.6),

$\frac{d\theta_{i}}{dt}=\omega_{i}+\varepsilon
R\sin(\psi-\theta_{i})=0,\Longrightarrow\varepsilon
R\sin(\theta_{i}-\psi)=\omega_{i}.$

For oscillators that have a natural frequency $|\omega_{i}|<\varepsilon R$,
these oscillators are locked at
$\theta_{i}-\psi=\arcsin({\omega_{i}/\varepsilon R})$. If
$|\omega_{i}|>\varepsilon R$ the oscillators rotate incoherently around the
circle with a probability density $\rho$ constant over time, that is
$\frac{\partial\rho}{\partial t}=0$. Therefore, from Eq. (A.9)
$\upsilon\rho=C$ (constant), that is

$\rho=\frac{C}{|\omega-\varepsilon R\sin(\theta_{i}-\psi)|},$

where $C$ can be evaluated using the normalization condition of Eq. (A.10),

$\int_{-\pi}^{\pi}\rho(\theta,\omega,t)d\theta=1,\Rightarrow\int_{-\pi}^{\pi}\frac{d\theta}{|\omega-\varepsilon
R\sin(\theta-\psi)|}=\frac{1}{C},$

being $\omega>\varepsilon R$, it is possible to remove the absolute value

$\int_{-\pi}^{\pi}\frac{d\theta}{\omega-\varepsilon R\sin(\theta-\psi)},$

and is used the solution

$\int\frac{dx}{a-b\sin
x}=\frac{-2\mathrm{arctg}\left(\frac{b-a\mathrm{tg}(\frac{x}{2})}{\sqrt{a^{2}-b^{2}}}\right)}{\sqrt{a^{2}-b^{2}}}+C_{0},$

where $a=\omega$, $b=\varepsilon R$, $x=\theta-\psi$ and $dx=d\theta$, so,

$\displaystyle\mathrm{if}\;x$ $\displaystyle=$
$\displaystyle+\pi\rightarrow\frac{-2\cdot\left(\frac{-\pi}{2}\right)}{\sqrt{\omega^{2}-(\varepsilon
R)^{2}}}=\frac{\pi}{\sqrt{\omega^{2}-(\varepsilon R)^{2}}},$
$\displaystyle\mathrm{if}\;x$ $\displaystyle=$
$\displaystyle-\pi\rightarrow\frac{-2\cdot\left(\frac{\pi}{2}\right)}{\sqrt{\omega^{2}-(\varepsilon
R)^{2}}}=\frac{-\pi}{\sqrt{\omega^{2}-(\varepsilon R)^{2}}},$

hence,

$C=\frac{\sqrt{\omega^{2}-(\varepsilon R)^{2}}}{2\pi},$ (A.11)

therefore, the probability density can be written as

$\displaystyle\rho(\theta,\omega)=\begin{dcases}\delta\left[\theta-\psi-\arcsin\left(\frac{\omega}{\varepsilon
R}\right)\right],\;\;\;|\omega|<\varepsilon R\\\ \frac{C}{|\omega-\varepsilon
R\sin(\theta-\psi)|},\;\;\;\;\;\;\;\;\;\;\;\;\mathrm{otherwise.}\end{dcases}$
(A.12)

$\rho$ can be replaced using Eq. (A.8) for $R$,

$\begin{split}R&=\int_{-\pi}^{\pi}\int_{|\omega|<\varepsilon
R}e^{i\theta}\delta\left[\theta-\psi-\arcsin\left(\frac{\omega}{\varepsilon
R}\right)\right]\varkappa(\omega)d\omega d\theta\;+\\\
&+\int_{-\pi}^{\pi}\int_{|\omega|>\varepsilon
R}e^{i\theta}\frac{C}{|\omega-\varepsilon
R\sin(\theta-\psi)|}\varkappa(\omega)d\omega d\theta.\end{split}$ (A.13)

And with the symmetry of the probability density of $\rho$

$\begin{split}\rho(\theta+\pi,-\omega)&=\frac{C}{|-\omega-\varepsilon
R\sin(\theta+\pi-\psi)|}=\frac{C}{|-\omega+\varepsilon
R\sin(\theta-\psi)|}=\\\ &=\frac{C}{|-(\omega-\varepsilon
R\sin(\theta-\psi))|}=\frac{C}{|\omega-\varepsilon
R\sin(\theta-\psi)|}=\rho(\theta,\omega).\end{split}$

Assuming that $\varkappa(\omega)=\varkappa(-\omega)$, and with
$\rho(\theta+\pi,-\omega)=\rho(\theta,\omega)$ the second term of the integral
(for $|\omega|>\varepsilon R$)

$\begin{split}\int_{-\pi}^{\pi}\int_{|\omega|>\varepsilon
R}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta&=\underbrace{\int_{-\pi}^{\pi}\int_{-\infty}^{-\varepsilon
R}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta}_{\omega=-\omega,\;\theta=\theta+\pi}+\\\
&+\int_{-\pi}^{\pi}\int_{\varepsilon
R}^{+\infty}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta,\end{split}$ $\begin{split}\int_{-\pi}^{\pi}\int_{|\omega|>\varepsilon
R}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta&=\int_{-\pi}^{\pi}\int_{+\infty}^{\varepsilon
R}e^{i\theta}e^{i\pi}\rho(\theta+\pi,-\omega)\varkappa(-\omega)(-d\omega)d\theta+\\\
&+\int_{-\pi}^{\pi}\int_{\varepsilon
R}^{+\infty}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta,\end{split}$

after the substitutions, the first term of the integral gains three signals,
one due the $e^{i\pi}=-1$, other from $-d\omega$, and the last from the
inversion of the integral limits at $\omega$. Finally,

$\begin{split}\int_{-\pi}^{\pi}\int_{|\omega|>\varepsilon
R}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta=&-\int_{-\pi}^{\pi}\int_{\varepsilon
R}^{+\infty}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta\;+\\\ &+\int_{-\pi}^{\pi}\int_{\varepsilon
R}^{+\infty}e^{i\theta}\rho(\theta,\omega)\varkappa(\omega)d\omega
d\theta=0.\end{split}$

Consequently, the oscillators with $|\omega|>\varepsilon R$ do not contribute
to the synchronization of the system. Remaining

$R=\int_{-\pi}^{\pi}\int_{|\omega|<\varepsilon
R}e^{i\theta}\delta\left[\theta-\psi-\arcsin\left(\frac{\omega}{\varepsilon
R}\right)\right]\varkappa(\omega)d\omega d\theta,$

and rewriting,

$R=\int_{-\varepsilon R}^{\varepsilon
R}\varkappa(\omega)d\omega\int_{-\pi}^{\pi}(\cos\theta+i\sin\theta)\delta\left[\theta-\psi-\arcsin\left(\frac{\omega}{\varepsilon
R}\right)\right]d\theta,$

being $\delta(x)$ and $\cos(x)$ functions of even parity, and $\sin(x)$ with
odd parity. By symmetry, the $\sin$ therm vanishes, hence,

$R=\int_{-\varepsilon R}^{\varepsilon
R}\varkappa(\omega)d\omega\int_{-\pi}^{\pi}\cos\theta\delta\left[\theta-\psi-\arcsin\left(\frac{\omega}{\varepsilon
R}\right)\right]d\theta$

solving the integral in $\theta$, using the filtering property of the
$\delta$, and translating the system to $\theta=\theta-\psi$

$R=\int_{-\varepsilon R}^{\varepsilon
R}\varkappa(\omega)d\omega\cos{\left(\arcsin\left(\frac{\omega}{\varepsilon
R}\right)\right)},$

and replacing

$\begin{split}\omega&=\varepsilon R\sin\theta,\;\;\Rightarrow
d\omega=\varepsilon R\cos\theta d\theta,\\\
\mathrm{If}&\;\;\omega=-\varepsilon R\rightarrow\theta=-\frac{\pi}{2},\\\
\mathrm{If}&\;\;\omega=+\varepsilon R\rightarrow\theta=+\frac{\pi}{2},\\\
\end{split}$

the integral

$R=\varepsilon
R\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos{{}^{2}\theta}\varkappa(\varepsilon
R\sin\theta)d\theta,\;\;\;\;0<R<1.$ (A.14)

Considering the critical coupling $\varepsilon^{*}$ in which

$\lim_{R\rightarrow 0^{+}}\varkappa(\varepsilon R\sin\theta)=\varkappa(0),$

hence,

$1=\varepsilon^{*}\varkappa(0)\underbrace{\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos{{}^{2}\theta}d\theta}_{=\;\pi/2},$

that is, there is an expression for the critical coupling that limits the
region of desynchronization $0\leq\varepsilon\leq\varepsilon^{*}$. Region
which depends only on the probability density of natural frequencies, which it
is assumed to be unimodal and symmetrical

$\varepsilon^{*}=\frac{2}{\pi\varkappa(0)}.$ (A.15)

For $\varepsilon\gtrsim\varepsilon^{*}$, close to the synchronization
transition point, using Eq. (A.14)

$1=\varepsilon\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos{{}^{2}\theta}\varkappa(\varepsilon
R\sin\theta)d\theta,$

expanding the Taylor function $\varkappa(\omega)$ centered in $0$

$\varkappa(\omega)=\varkappa(0)+\varkappa^{\prime}(0)(\omega)+\frac{\varkappa^{\prime\prime}(0)}{2}(\omega^{2})+\cdots,$

it is know that $\varkappa(\omega)=\varkappa(-\omega)$, hence, the first
derivative $\varkappa^{\prime}(0)=0$, and disregarding the cubic therms, the
integral

$\begin{split}1=\varepsilon\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos{{}^{2}\theta}\left[\varkappa(0)+\frac{\varkappa^{\prime\prime}(0)}{2}(\varepsilon
R\sin\theta)^{2}\right]d\theta&=\varepsilon\varkappa(0)\overbrace{\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos{{}^{2}\theta}d\theta}^{=\;\pi/2}\;+\\\
&+\frac{\varepsilon^{3}R^{2}\varkappa^{\prime\prime}(0)}{2}\overbrace{\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos{{}^{2}\theta}\sin{{}^{2}\theta}d\theta}^{=\;\pi/8},\end{split}$

soon,

$1=\varepsilon\frac{\pi\varkappa(0)}{2}+\frac{\pi\varepsilon^{3}R^{2}\varkappa^{\prime\prime}(0)}{16},$

and replacing Eq. (A.15)

$1=\frac{\varepsilon}{\varepsilon^{*}}+\frac{\pi\varepsilon^{3}R^{2}\varkappa^{\prime\prime}(0)}{16}\Rightarrow\varepsilon^{*}=\varepsilon+\frac{\pi\varepsilon^{3}\varepsilon^{*}R^{2}\varkappa^{\prime\prime}(0)}{16}\Rightarrow\varepsilon-\varepsilon^{*}=-\frac{\pi\varepsilon^{3}\varepsilon^{*}R^{2}\varkappa^{\prime\prime}(0)}{16}.$

Isolating $R$, finally the expression for the modulus of the order parameter

$R=\sqrt{\frac{-16}{\pi\varepsilon^{2}\varepsilon^{*}\varkappa^{\prime\prime}(0)}\left(1-\frac{\varepsilon^{*}}{\varepsilon}\right)},$
(A.16)

and since it is an approximation for $\varepsilon\gtrsim\varepsilon^{*}$, it
is written the first therm as a $\varepsilon^{*}$ dependence only [137]

$R=\sqrt{\frac{-16}{\pi(\varepsilon^{*})^{3}\varkappa^{\prime\prime}(0)}\left(1-\frac{\varepsilon^{*}}{\varepsilon}\right)}.$
(A.17)