UNIVERSIDADE FEDERAL DO PARANÁ

Effects of the local dynamics in the synchronization of neural models

Bruno R. R. Boaretto

Curitiba
$24^{\mathrm{th}}$ February, 2022

UNIVERSIDADE FEDERAL DO PARANÁ

Effects of the local dynamics in the synchronization of neural models

Bruno R. R. Boaretto

Thesis presented as a requirement for the degree of PhD in Science in the
Graduate program in Physics, Exact Sciences Sector, of the Universidade
Federal do Paraná.

Advisor: Prof. Dr. Sergio Roberto Lopes

Curitiba
$24^{\mathrm{th}}$ February, 2022

Aos meus pais, Claudia Reichert e Jacir Boaretto.

AGRADECIMENTOS

Agradeço aos meus pais, Claudia Reichert e Jacir Boaretto, por proverem todo o
suporte e apoio necessários, permitindo que esta tese se concretizasse.
Agredeço também ao meu irmão, Marco, pelo convívio nos últimos anos.

Agradeço à Mariana, minha companheira, que acompanhou-me durante toda a essa
jornada, demonstrando compreensão, paciência e amor em todos os momentos.

Agradeço ao meu orientador, Prof. Dr. Sergio Roberto Lopes, pelo permanente
incentivo na execução desta tese e dos demais projetos produzidos desde a
iniciação científica. Homenageando-o, agradeço aos demais colaboradores que
com muita humildade e interesse acreditaram em mim e nas minhas capacidades.

Agradeço aos meus amigos e colegas de grupo, Roberto e Kalel, pela amizade e
cumplicidade nas discussões sobre a pesquisa e a todos os colegas do programa
pela amizade e apoio nos estágios desse processo.

Agradeço aos professores que leram revisaram e corrigiram atentamente esta
tese em suas diferentes versões: Profa. Dra. Sabrina Borges Lino Araujo
(UFPR), Prof. Dr. Ricardo Luiz Viana (UFPR), Prof. Dr. Marcus Werner Beims
(UFPR), Prof. Dr. Thiago de Lima Prado (UFPR), Prof. Dr. Giovani Vasconcelos
(UFPR), Prof. Dr. Iberê Luiz Caldas (USP), Prof. Dr. Antonio Marcos Batista
(UEPG), Prof. Dr. Elbert Einstein Nehrer Macau (UNIFESP).

Finalmente, agradeço à Universidade Federal do Paraná, ao Departamento de
Física e ao Programá de Pós-Graduaçao em Física por toda a estrutura cedida.
Registro meus agradecimentos também a CAPES pelo auxílio financeiro.

“O que aprendi, acima de tudo, é a seguir em frente, pois a grande ideia é a
de que, como o acaso efetivamente participa de nosso destino, um dos
importantes fatores que levam ao sucesso está sob o nosso controle: o número
de vezes que tentamos rebater a bola, o número de vezes que nos arriscamos, o
número de oportunidades que aproveitamos.”
Leonard Mlodinow - O Andar do Bêbado, 2009.

RESUMO

O comportamento cooperativo de neurônios e áreas neuronais associadas ao
comportamento de sincronização se apresenta como mecanismo fundamental para o
funcionamento cerebral. Além disso, níveis anormais de sincronização têm sido
relacionados a estados patológicos. Ao longo desta tese, abordam-se diferentes
fenômenos de sincronização que surgem por meio da dinâmica coletiva de modelos
de neurônios acoplados em uma rede. Primeiramente, mostra-se uma forte
correlação entre a dinâmica individual do neurônio com o comportamento global
da sincronização da rede, em que a periodicidade observada no neurônio isolado
é refletida em uma sincronização de fase ao considerar um acoplamento fraco.
Em segundo lugar, estuda-se o papel da biestabilidade na sincronização de uma
rede de neurônios idênticos, acoplados através de um esquema de campo médio.
Mostra-se que a simples existência de dois estados estáveis distintos pode
levar a rede a diferentes estados de sincronização, dependendo da
inicialização do sistema. Por fim, é investigado o mecanismo de sincronização
explosiva de uma rede neural complexa composta por neurônios não-idênticos. A
presença deste regime é acompanhada por um loop de histerese na dinâmica da
rede, à medida que o parâmetro de acoplamento é adiabaticamente aumentado e
reduzído. Demonstra-se que as transições de sincronização abruptas estão
associadas a rotas para o caos e que os mecanismos dinâmicos para a região de
biestabilidade são dados em termos de uma bifurcação de sela-nó e uma crise de
fronteira. Portanto, os resultados desta tese mostram uma riqueza de
comportamentos de sincronização associados a pequenas mudanças na dinâmica
neuronal, trazendo novos insights para o estudo teórico das redes neurais.

Palavras-chave: Redes neurais, Dinâmica local, Sincronização de Fase,
Transição de Sincronização, Sincronização Explosiva.

ABSTRACT

The cooperative behavior of neurons and neuronal areas associated with the
synchronization behavior proves to be a fundamental neural mechanism. In
addition, abnormal levels of synchronization have been related to unhealthy
neural behaviors. Throughout this thesis, it is explored different
synchronization phenomena which emerge through the collective dynamics of
models of neurons coupled in a network. Firstly, it is shown a strong
correlation between the individual dynamics of the neuron with the global
behavior of the synchronization, in which the periodicity seen in the isolated
neuron is reflected in a phase synchronization in the weak coupling region.
Secondly, it is studied the role of bistability in the synchronization of a
network of identical neurons coupled through a mean-field scheme. It is shown
that the simple existence of two distinct stable states can lead the network
to different states of synchronization, depending on the initialization of the
system. Lastly, it is investigated the mechanism for explosive synchronization
of a complex neural network composed of non-identical neurons. The presence of
this regime is accompanied by a hysteresis loop on the network dynamics as the
coupling parameter is adiabatically increased and decreased. It is 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. Therefore, the results of this
thesis show a richness of synchronization behaviors associated with small
changes of the neuronal dynamics bringing new insights to the theoretical
study of neural networks.

Keywords: Neural Network, Local Dynamics, Phase Synchronization,
Synchronization Transition, Explosive Synchronization.

###### List of Figures

  1. 2.1 Projection of the Lorenz attractor in the phase space and the sensibility to initial conditions.
  2. 2.2 Simple pendulum.
  3. 2.3 Generic phase space ($x_{1},x_{2}$).
  4. 2.4 Iterations of the logistic map for different values of $r$.
  5. 2.5 Temporal evolution of a bidimensional sphere of initial conditions.
  6. 2.6 Lyapunov exponent of the logistic map.
  7. 2.7 Period doubling route to chaos.
  8. 2.8 Intermittency route to chaos.
  9. 2.9 Three types of intermittency route to chaos.
  10. 2.10 Crisis route in the logistic map.
  11. 3.1 Basic anatomy of a neuron.
  12. 3.2 Cell membrane representation
  13. 3.3 Action potential recorded between inside and outside axon
  14. 3.4 Equivalent circuit model
  15. 3.5 Time-evolution of the variables of the Hodgkin-Huxley model.
  16. 3.6 Bifurcation diagram of the Hodgkin-Huxley model.
  17. 3.7 Time-evolution of the variables of the HH$\ell$ model.
  18. 3.8 Bifurcation diagram of the HH$\ell$ model.
  19. 3.9 Bifurcation diagram of the HR model.
  20. 3.10 Time evolution of the variables of the Chialvo model.
  21. 4.1 Illustration of the Königsberg bridges.
  22. 4.2 Examples of regular networks.
  23. 4.3 Examples of random networks.
  24. 4.4 Example of the Watts-Strogatz route.
  25. 4.5 Example of the Newman-Watts route.
  26. 4.6 Small-worldness coefficient.
  27. 5.1 Time evolution of the phases of identical Kuramoto oscillators and the relaxation time of the network.
  28. 5.2 Phase synchronization of the Kuramoto model.
  29. 6.1 Phase association of the HH$\ell$ model.
  30. 6.2 Neural dynamics of HH$\ell$ model as a function of ${g_{\mathrm{r}}}$.
  31. 6.3 Synchronization transition of individual regular neurons.
  32. 6.4 Synchronization transition of individual chaotic neurons.
  33. 6.5 Correlation between the individual dynamics and the phase synchronization varying the parameter space $g_{\mathrm{r}}\times g_{\mathrm{d}}$.
  34. 6.6 Correlation between the individual dynamics and the phase synchronization varying the parameter space $g_{\mathrm{sr}}\times g_{\mathrm{sd}}$.
  35. 6.7 Chaos induced by the application of external pulses.
  36. 6.8 Verification test with the Hindrmarsh-Rose model.
  37. 6.9 Synchronization transition for different networks’ topologies.
  38. 7.1 Bistability on the HH$\ell$ model.
  39. 7.2 Stability of each state.
  40. 7.3 The role of the bistability at the synchronization in a network of identical HH$\ell$ neurons.
  41. 7.4 Temporal evolution of the synchronization in the network at the situation (ii) with $\varepsilon=10^{-4}$.
  42. 7.5 Detailed analysis in the synchronization of situation (iii) with $\varepsilon=10^{-4}$.
  43. 7.6 Detailed analysis in the synchronization of situation (iii) with $\varepsilon=10^{-3}$.
  44. 7.7 A complete transition to state I depends on the initial conditions.
  45. 7.8 The dependence of synchronization features on the initial conditions analyzed for $\varepsilon=0.00016$.
  46. 7.9 Mean time $\langle\tau\rangle$ needed for the network to reach phase synchronization as a function of $\varepsilon$ for $100$ initial conditions.
  47. 8.1 Phase association in the Chialvo model.
  48. 8.2 The role of the topology in the synchronization behavior.
  49. 8.3 Numerical results for ${S_{\mathrm{h}}}$ as a function of ${p_{\mathrm{nw}}}$.
  50. 8.4 The explosive synchronization occurs due to a chaotic transition.
  51. 8.5 Illustration of the chaotic transient in a crisis.
  52. 8.6 Evidence of boundary crisis in the network.
  53. 8.7 Illustration of intermittent chaotic events.
  54. 8.8 Evidence of saddle node in the network.

###### List of Tables

  1. 3.1 Constants used for the simulation of the Hodgkin-Huxley model [33]
  2. 3.2 Values of the parameters of the HH$\ell$ model [76]
  3. 3.3 Constants used for the simulation of the Hindmarsh-Rose model [40].
  4. 7.1 Lyapunov spectra for IC-1 and IC-2 computed for two values of $g_{\mathrm{d}}$.

###### Contents

  1. 1 Introduction
  2. 2 Dynamical systems and chaos
    1. 2.1 Dynamical systems
    2. 2.2 Chaos in dynamical systems
    3. 2.3 Routes to chaos
      1. 2.3.1 The period doubling cascade
      2. 2.3.2 Intermittency route to chaos
      3. 2.3.3 Crisis
  3. 3 Neuronal models and synapses
    1. 3.1 The action potential
    2. 3.2 Neuron models
      1. 3.2.1 The Hodgkin-Huxley (HH) model
      2. 3.2.2 A Hodgkin-Huxley-like (HH$\ell$) model
      3. 3.2.3 The Hindmarsh-Rose model
      4. 3.2.4 The Chialvo model
    3. 3.3 Synaptic current
  4. 4 Complex networks
    1. 4.1 Graphs theory
    2. 4.2 Complex networks
  5. 5 Phase synchronization
  6. 6 The role of the individual dynamics in the synchronization processes of neural networks
    1. 6.1 Network properties
    2. 6.2 Results
      1. 6.2.1 External perturbation by electrical pulses
      2. 6.2.2 Verification test with the Hindmarsh-Rolse model
      3. 6.2.3 Topological effects
    3. 6.3 Discussions
  7. 7 Bistability in the synchronization of identical neurons
    1. 7.1 Bistability on the HH$\ell$ model
    2. 7.2 Network properties
    3. 7.3 Network results
      1. 7.3.1 Bistability with periodic-chaotic configuration (situation (ii))
      2. 7.3.2 Bistability with periodic-periodic configuration (situation (iii))
    4. 7.4 Discussions
  8. 8 Mechanism for explosive synchronization of neural networks
    1. 8.1 Network properties
    2. 8.2 Results
      1. 8.2.1 The dynamical mechanism of the loss of stability of the chaotic attractor
      2. 8.2.2 The dynamical mechanism of the loss of stability of the phase-synchronized attractor
    3. 8.3 Discussions
  9. 9 Conclusions and future perspectives
  10. A Kuramoto oscillators

##  Chapter 1 Introduction

The synchronization phenomenon is studied for centuries since the scientist
Christiaan Huygens reported his observations that two weakly coupled pendulum
clocks become synchronized in-phase [1] (English translation [2]). After that,
this phenomenon was detected in a wide range of biological systems [3, 4, 5],
as in the rhythmic flashing of fireflies [6], in the crickets which
synchronize their chirps by responding to the preceding chirp of their
neighbors [7], groups of women whose menstrual periods become mutually
synchronized [8], in the synchronous of rabbit sino-atrial pace-maker cells
[9], and also in the action potentials of the nervous system [10, 11].

A healthy human brain is composed of $\sim 10^{11}$ neuronal cells
interconnected by $\sim 10^{15}$ synapses creating groups of connected neurons
divided into brain regions, each with specific functions [12]. In particular,
the complex behaviors seen in the brain are directly related to the emergence
of spatial-temporal activation patterns that comes spontaneously as a result
of the cooperative interaction among neurons. The role of neuroscience is to
understand how does the behaviors produced by the brain like perception,
movement, language, thought, memory, etc; can be explained in terms of the
activity patterns of neurons since all the behavioral disorders that
characterize psychiatric illness are disturbances on the brain functioning
[12]. In this context, the cooperative behavior of neurons and neuronal areas
associated with the synchronization behavior proves to be a fundamental neural
mechanism [10]. It supports memory process [10, 13], information process [14],
and is relevant for many cognitive processes [15, 16]. In addition, abnormal
levels of synchronization have been related to unhealthy neural behaviors
[12]. While a high degree of synchronization is detected in epileptic seizures
[17], where the increase in the synchronization of some groups of neurons
generates seizure episodes, and in Parkinson’s disease, where there is an
excessive synchronization in basal ganglia [18, 19], reduced levels of
synchronization among cortical areas can be associated with brain disorders
such as autism [20], and Alzheimer’s disease [21].

From the theoretical point of view, the use of complex networks proves to be
useful for the study of an ensemble of coupled sites, in which the global
behavior of the system results from the interaction among the sites achieving
an ample possibility of phenomena even in the case of simple interactions [22,
23]. Particularly, for neuronal systems, each site of the network is composed
of a neuron and the edges of the network represent its synaptic connections.
This approach can be used to improve our knowledge and give theoretical
insights into the understating of brain functioning. In studies of neuronal
systems, there are always two critical choices: what model describes the
firing dynamics of each neuron and how the neurons are connected [24]. While
several works are focused on understanding the role of the connection
architecture, also called topology of connection, in the synchronization of
dynamical systems [25, 26, 27, 28, 29, 30, 31, 32], this thesis, it is studied
the sensibility of the synchronization features to the dynamics of the
isolated neurons.

In the last decades, after the success of the Hodgkin-Huxley model, being the
first quantitative description of the regenerative currents generating the
action potential [33, 34, 35], dozens of models were created to reproduce more
complicated firing patterns exhibited by neurons [24, 36]. From the simplest
models described by iterated maps [37, 38], to more sophisticated models of
several non-linear differential equations [35, 39, 40, 41, 42]. Often, by
changing parameters of the models it is possible to change the dynamical
behavior of the neuron, from regular activity to a chaotic one, where the
activation patterns of the neurons occur in a non-periodic way [33]. In
addition, some dynamical models also can exhibit multistable states, where a
neuron initialized with different initial conditions can present different
stable states, with different firing patterns, frequencies, regularity, and
chaoticity [43, 44, 45, 46].

Throughout this thesis, it is studied how the individual dynamics of the
neurons affect the synchronization of the neuronal network. For neurons that
individually exhibit chaotic behavior, synchronization occurs similarly to the
transition of chaotic oscillators known in the literature [47, 48]. However,
the knowledge about the synchronization of non-chaotic neurons that lose these
features when coupled, acquiring characteristics of chaotic neurons, is still
being studied [47]. It is shown that, in some cases, these neurons tend to
synchronize in phase to weak coupling regimes (a coupling parameter close to
zero) due to the influence of the regularity of the individual dynamics. This
synchronized state is lost with the increment of the coupling strength, due to
the influence of the collective behavior of the dynamics [49, 50, 51],
characterizing a non-monotonic evolution of the synchronization phenomenon as
a function of coupling.

On the other hand, for neurons that individually present periodic dynamics,
i.e., a well-defined frequency, the phase synchronization of a network with
identical neurons can be achieved for any non-zero coupling strength, with
more weakly-coupled networks needing more time to reach the phase-synchronized
state. In this sense, the existence of identical bistable neurons in the
network can delay the achievement of a complete-phase-synchronized state,
demanding transitions to a unique and identical state for all neurons in the
network. These transitions are induced by the coupling of the network and
occur according to the stability of each state, in which the less stable state
transitions to the more stable state.

Moreover, the synchronization of periodic neurons also can be disturbed with
the existence of heterogeneity in the network, in which non-identical-periodic
neurons are coupled. Neurons simulated with different parameters may depict
different activation patterns. This dissimilitude in the neuronal dynamics can
influence both frequency and amplitude of each neuron generating a chaotic
non-synchronized regime instead of a complete-periodic synchronization. It is
show that the increase of the coupling strength can transits the system
abruptly to the phase-synchronized state, characterizing an explosive
synchronization [52, 53, 54, 55, 56].

Given this context, this thesis aims to understand the effect of individual
dynamics in the synchronization of neural networks. This objective can be
divided into three independent topics, described as:

  * •

Understand the non-monotonic evolution of the synchronization as a function of
the coupling in neuronal networks where the chaoticity arises from the
synaptic currents.

  * •

Investigate the role of bistability in the synchronization of a network of
identical bursting neurons coupled through a generic electrical mean-field
scheme.

  * •

Study the mechanism for explosive synchronization of a complex neural network
composed of non-identical spiking neurons and coupled through a small-world
network.

This thesis is divided into two parts, the first part, from Chapters 2 to 5,
is devoted to introducing all the theoretical frameworks used in this work,
while the results and conclusions are presented in the second part from
Chapters 6 to 9.

Chapter 2 presents a general review of the main concepts related to dynamical
systems and nonlinear dynamics, like fixed points, stability of fixed points,
dynamical features of chaotic systems, Lyapunov exponent, and routes to chaos.
Chapter 3 is focused on the presentation of properties of physiological
neurons. Moreover, it is introduced the formulation ideas of dynamical systems
which display a qualitative behavior of real neurons. The models can be
described by differential equations, like the Hodgkin-Huxley model, and also
by iterated maps like the Chialvo model. Also, it is introduced the synaptic
structures that permit the communication between presynaptic and postsynaptic
neurons. After this, Chapter 4 shows some concepts of the graphs theory used
to construct the complex networks which rule the connections between the
neurons. Finally, Chapter 5 explores some properties of coupled oscillators,
and also it is presented a powerful tool used to quantify phase
synchronization, the Kuramoto order parameter.

In Chapter 6, it is studied the effect called non-monotonic synchronization,
where using a Hodgkin-Huxley-like (HH$\ell$) neuron, proposed by Braun et al.
[57], which in a periodic-bursting regime, it is possible to achieve phase
synchronization for small values of coupling. It shows a clear relation
between the individual dynamics of the isolated neuron with the type of phase
synchronization transition of the network. After that, it is shown that the
same correlation occurs for neurons under external perturbations and with
other similar neuronal models.

By investigating the parameter space of the HH$\ell$ model, it is found a
parameter region where the neuron exhibits bistability behavior. One of the
states named state I is always periodic, while the other state, named state
II, goes through a sequence of period-doubling bifurcations from periodic to
chaotic behavior. In this scenario, Chapter 7 investigates the stability of
each state, and how the existence of bistability can change the phase
synchronization of a network. To isolate this behavior, the neurons are
coupled through a generic mean-field scheme (all-to-all connection) where the
only difference between the neurons are the initial conditions. After all, the
simple existence of two distinct stable states in the network can produce
distinct synchronization states.

Chapter 8 analyzes the mechanism for an abrupt transition to phase
synchronization, called explosive synchronization, which is found in a complex
network composed of non-identical spiking neurons, simulated with the neuronal
model proposed by Chialvo [37]. It is shown that this regime is accompanied by
a hysteresis loop on the network dynamics as the coupling parameter is
adiabatically increased and decreased, characterizing a bistability regime.
The abrupt synchronization transitions are associated with routes to chaos,
and the dynamical mechanisms for the bistability region are given in terms of
a saddle-node bifurcation and a boundary crisis.

At last, the conclusions based on the results and open questions for further
investigation are presented in Chapter 9.

##  Chapter 2 Dynamical systems and chaos

This chapter presents some concepts about dynamical systems, which are
fundamental for the understanding of the time evolution of real systems, both
regular and chaotic systems. In this sense, Edward Lorenz has studied a
simplified mathematical model for atmospheric convection [58, 59, 60]. The
model is a system of three ordinary differential equations where $x$, $y$, and
$z$ denote variables proportional to convective intensity, horizontal, and
vertical temperature differences

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle\sigma(y-x),$ (2.1)
$\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle x(r-z)-y,$ (2.2)
$\displaystyle\dot{z}$ $\displaystyle=$ $\displaystyle xy-bz,$ (2.3)

where the $(\cdot)$ symbol represents the derivative of a function of time and
$\sigma$, $r$, and $b$ are parameters, called the Prandtl number, Rayleigh
number, and a geometric factor, respectively. This system becomes very popular
due to the “butterfly effect” which emphasizes the sensitive dependence on
initial conditions in which a small change in one state of a deterministic
nonlinear system can result in large differences in a later state. Figure 2.1
(a) presents the projection of the Lorenz attractor in the phase space, with
the shape that may also be seen to resemble a butterfly, and panels (b – d)
show the time evolution of two slightly different initial conditions. The
sensibility to initial conditions which evolve in such a way that their
trajectories diverge exponentially is one of the signatures of a chaotic
motion. This chapter shows how a dynamical system can achieve chaotic dynamics
and how to quantify it.

Figure 2.1: Projection of the Lorenz attractor in the phase space and the
sensibility to initial conditions. (a) Projection of the Lorenz attractor
$x\times z$. (b – d) Dynamical variables of the Lorenz system using
$x_{1}(0)=x_{2}(0)=0,\,y_{1}(0)=y_{2}(0)=1\,,z_{1}(0)=1.05,\,z_{2}(0)=1.051$
and $\sigma=16$, $r=45.92$ and $b=4$, with an integration timestep of $0.01$.
After some integration steps due to the chaotic behavior, even for close
initial conditions the orbits tend to diverge from each other.

###  2.1 Dynamical systems

A dynamical system is a concept in which a mathematical function describes the
behavior of particles under the action of a set of laws. In this sense, for a
given configuration, called initial condition, a dynamical system is
deterministic and presents only one possible solution. There are two ways to
describe the evolution of a dynamical system: with differential equations,
which describe the evolution of the system considering time as a continuous
variable, such systems are called flows; and through difference equations that
consider time as a discrete variable, such systems are called maps.

Considering a continuous system with $D$ dimensions characterized by the
variables $(x_{1},\,x_{2},$ $\,\cdots,\,x_{D})$, the temporal evolution of
these variables are described by the set of $D$ equations in which,

$\displaystyle\dot{x_{1}}$ $\displaystyle=$ $\displaystyle
f_{1}(x_{1},x_{2},\cdots,x_{D}),$ $\displaystyle\dot{x_{2}}$ $\displaystyle=$
$\displaystyle f_{2}(x_{1},x_{2},\cdots,x_{D}),$ $\displaystyle\vdots$
$\displaystyle\dot{x_{D}}$ $\displaystyle=$ $\displaystyle
f_{D}(x_{1},x_{2},\cdots,x_{D}),$

or rewriting

$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}).$ (2.4)

In this sense, if the elements of $\mathbf{f}$ presents products, power, and
functions of $\mathbf{x}$, such as $x_{i}x_{j}$, $x_{i}^{3}$, or even
$\cos(x_{i}x_{j})$ the system is characterized nonlinear, or linear, otherwise
[61]. The system is periodic if exists a time $\mathcal{T}$ where:

$\mathbf{x}(t+\mathcal{T})=\mathbf{x}(t).$ (2.5)

The simple pendulum is an example of a system where the time is a continuous
variable, represented in Fig. 2.2, where a particle of mass $m$ is attached to
a wire of length $\ell$ which under the action of the gravitational
acceleration oscillates with an angle $\theta(t)$ around a fixed point. Even
though it is two-dimensional system, when considering the constraint (wire),
the movement of the particle is restricted to an arc of angle $\theta$ of
radius $\ell$. Consequently

$\ddot{\theta}=-\frac{g}{\ell}\sin\theta,$ (2.6)

where $g$ is the absolute value of the gravitational acceleration, and $\ell$
is the length of the wire. The therm $\sin\theta$ characterizes the
nonlinearity of the system.

Figure 2.2: Simple pendulum. A particle with mass $m$ attached to an
inextensible wire with length $\ell$ that oscillates with an angle $\theta(t)$
around a fixed point due to the gravitational acceleration $\mathbf{g}$.

The nonlinearity in the equations makes it difficult to find an analytic
solution to the problem, that is, writing $\theta$ as a function of $t$. To
make the things simple, the small oscillation approximation is used, where
$\sin\theta\approx\theta$,

$\ddot{\theta}\approx-\frac{g}{\ell}\mathrm{\theta}.$ (2.7)

Considering the angular frequency $\omega=\sqrt{g/\ell}$, and as initial
conditions $\theta(0)=\theta_{0}$, and $\dot{\theta}(0)=0$, the analytical
solution is described by

$\theta(t)=\theta_{0}\cos(\omega t),$ (2.8)

and the period of oscillation is

$\mathcal{T}_{\mathrm{linear}}=2\pi\sqrt{\frac{\ell}{g}},$ (2.9)

such equation (Eq.(2.7)), now linear, is equivalent to the harmonic oscillator
equation, where the period of oscillation is independent of the amplitude of
the system.

The fact that a system does not present an analytical solution is not
necessarily a limiting factor. Using any integration method it is possible to
find a numerical solution (point-by-point) for the system. The nonlinear
pendulum is one of the special cases in which can be solved analytically, the
solution is given in terms of elliptic integrals [62], where the period of
oscillation is described by

$\mathcal{T}_{\mathrm{nonlinear}}=4\sqrt{\frac{\ell}{g}}\mathcal{F}\left(\sin\frac{\theta_{0}}{2},\frac{\pi}{2}\right),$
(2.10)

where $\mathcal{F}$ represents the elliptic function of Legendre of the first
kind [62] defined as

$\mathcal{F}(k,\varphi)=\int_{0}^{\varphi}\frac{1}{\sqrt{1-k^{2}\sin^{2}\theta}}d\theta.$
(2.11)

Considering $x_{1}=\theta$ and $x_{2}=\dot{\theta}$ the Eq. (2.7) that has a
second-order derivative, can be rewritten in terms of two first-order
derivatives.

$\displaystyle\dot{x_{1}}$ $\displaystyle=$ $\displaystyle x_{2},$ (2.12)
$\displaystyle\dot{x_{2}}$ $\displaystyle=$ $\displaystyle-\frac{g}{\ell}\sin
x_{1},$ (2.13)

that is, a two-variable system, where $x_{1}$ is the position of the particle,
and $x_{2}$ is the velocity. Assuming that the solution of this system for a
given initial condition is known, this solution will be a set of functions
$x_{1}(t)$ (position) and $x_{2}(t)$ (velocity). Considering an abstract space
with coordinates $(x_{1},\;x_{2})$, then the solution $(x_{1}(t),\;x_{2}(t))$
corresponds to a point moving along a curve in that space. Figure 2.3
represents this abstract space, known as phase space, in which the curve
represents the evolution of the system solution, called trajectory.

Figure 2.3: Generic phase space ($\mathbf{x_{1},x_{2}}$). For a set of initial
conditions ($x_{1}(0),\;x_{2}(0)$) the solution of the system
($x_{1}(t),\;x_{2}(t)$) corresponds to a curve in the phase-space called
trajectory.

The phase space is completely composed of trajectories and every point in
phase space can be considered an initial condition. When analyzing the system
equations, several restrictions or transitions can be found in its dynamics in
phase space. For example, there are solutions known as equilibrium solutions
where for $\mathbf{x}(t)=\mathbf{x}^{*}$ the temporal evolution of the system
is

$\frac{d\mathbf{x}}{dt}\bigg{|}_{\mathbf{x}=\mathbf{x}^{*}}=\mathbf{f}(\mathbf{x}^{*})=0$
(2.14)

and the point $\mathbf{x}^{*}$ is denominated equilibrium point. When the
trajectory reaches such a point it is confined to it. In this sense, the
equilibrium points of a system can be stable or unstable. When disturbing the
system in the vicinity of the point, if the trajectory is attracted (repelled)
the equilibrium point is considered stable (unstable). For the case of the
pendulum, the equilibrium points are
$\mathbf{x}^{*}=(x_{1}^{*}=j\pi,x_{2}^{*}=0)$ where $j\in\mathbb{Z}$. However,
physically there is no difference between $0$ and $2\pi$, so only two
equilibrium points are considered $\mathbf{x}^{(*,1)}=(0,0)$ and
$\mathbf{x}^{(*,2)}=(\pi,0)$, that is, the lower and upper extremes of the
pendulum, respectively. As this is an idealized system (without friction, or
damping), intuitively, it is possible to notice that there is a difference
between the two equilibrium points. When disturbing the system at the
equilibrium point $\mathbf{x}^{(*,2)}$, that is, the pendulum at the top end,
the system is automatically repelled. On the other hand, when perturbing the
point $\mathbf{x}^{(*,1)}$, the trajectory will be kept very close to it. Note
that for this case without damping, the system will never stabilize at the
equilibrium point again, because even if it returns to the point $x_{1}=0$,
the velocity will be $x_{2}\neq 0$.

In general, the stability of the equilibrium points can be determined by
analyzing the vicinity of the points. Considering a equilibrium point
$\mathbf{x}^{*}$ of a $D$-dimensional system, and a perturbation in the
vicinity of the point $\delta\mathbf{x}=\mathbf{x}(t)-\mathbf{x}^{*}$, the
temporal evolution of the disturbance can be described in terms of

$\dot{\delta\mathbf{x}}=\frac{d}{dt}(\mathbf{x}(t)-\mathbf{x}^{*}),$

and since $\mathbf{x}^{*}$ is a constant
$\dot{\delta\mathbf{x}}=\dot{\mathbf{x}}$, and
$\dot{\delta\mathbf{x}}=\mathbf{f}(\mathbf{x})=\mathbf{f}(\mathbf{x}^{*}+\delta\mathbf{x})$.
The Taylor series expansion is

$\mathbf{f}(\mathbf{x}^{*}+\delta\mathbf{x})=\mathbf{f}(\mathbf{x}^{*})+\frac{\partial\mathbf{f}(\mathbf{x})}{\partial\mathbf{x}}\bigg{|}_{\mathbf{x}^{*}}\delta\mathbf{x}+\mathcal{O}(\delta\mathbf{x}^{2}),$

where $\mathbf{f}(\mathbf{x}^{*})=0$ due to the equilibrium definition and
$\mathcal{O}(\delta\mathbf{x}^{2})$ representing the higher-order terms of
$\delta\mathbf{x}$. Disregarding the terms of order greater than or equal to
2, as long as the disturbance is sufficiently small,

$\dot{\delta\mathbf{x}}\approx\frac{\partial\mathbf{f}(\mathbf{x})}{\partial\mathbf{x}}\bigg{|}_{\mathbf{x}^{*}}\delta\mathbf{x},$
(2.15)

the partial derivative of the equation is known as Jacobian matrix and is
defined by

$\displaystyle\mathcal{J}(\mathbf{x}^{*})\equiv\frac{\partial\mathbf{f}(\mathbf{x})}{\partial\mathbf{x}}\bigg{|}_{\mathbf{x}^{*}}=\begin{bmatrix}\frac{\partial
f_{1}}{\partial x_{1}}&\frac{\partial f_{1}}{\partial
x_{2}}&\cdots&\frac{\partial f_{1}}{\partial x_{D}}\\\ \frac{\partial
f_{2}}{\partial x_{1}}&\frac{\partial f_{2}}{\partial
x_{2}}&\cdots&\frac{\partial f_{2}}{\partial x_{D}}\\\
\vdots&\vdots&\ddots&\vdots\\\ \frac{\partial f_{D}}{\partial
x_{1}}&\frac{\partial f_{D}}{\partial x_{2}}&\cdots&\frac{\partial
f_{D}}{\partial x_{D}}\\\
\end{bmatrix}_{x_{1}^{*},x_{2}^{*},\cdots,x_{D}^{*}},$ (2.16)

so

$\dot{\delta\mathbf{x}}=\mathcal{J}(\mathbf{x}^{*})\delta\mathbf{x},$ (2.17)

where $\mathcal{J}(\mathbf{x}^{*})$ is the Jacobian matrix evaluated at the
equilibrium points $\mathbf{x}^{*}$. According to the theory of differential
equations, the solution can be described in terms of the complex eigenvalues
of the Jacobian matrix. If the real part of all eigenvalues of $\mathcal{J}$
is negative, the point is stable. However, if at least the real part of one of
the eigenvalues is positive, the point is unstable [60].

For discrete systems, instead of differential equations, it is used the
iterative maps

$\mathbf{x}_{t+1}=\mathbf{f}(\mathbf{x}_{t}),$ (2.18)

in which the index $t$ represents the $t$-th iteration of the map. That is,
the state of the variable $\mathbf{x}_{t+1}$ is described in terms of
$\mathbf{x}_{t}$, where $t\in\mathbb{N}$. Compared to a continuous system

$\frac{d\mathbf{x}}{dt}=\lim_{\Delta t\rightarrow 0}\frac{\mathbf{x}(t+\Delta
t)-\mathbf{x}(t)}{\Delta t}\approx\frac{\mathbf{x}(t+\Delta
t)-\mathbf{x}(t)}{\Delta t},$

that is, ignoring the limit of $\Delta t\rightarrow 0$, Eq. (2.4)

$\frac{d\mathbf{x}}{dt}=\mathbf{f}(\mathbf{x})\approx\frac{\mathbf{x}(t+\Delta
t)-\mathbf{x}(t)}{\Delta t},$

isolating the term $\mathbf{x}(t+\Delta t)$

$\mathbf{x}(t+\Delta t)\approx\mathbf{x}(t)+\mathbf{f}(\mathbf{x})\Delta t,$

that is, on the left, there is a state $\mathbf{x}(t+\Delta t)$, and on the
right-side terms that depend only on $\mathbf{x}(t)$, something analogous to
Eq. (2.18). Therefore, maps can be understood as discrete approximations of
continuous systems. In addition, the fact that an integrator is not required
to perform the numerical calculations the maps are known to exhibit high
performance in the computation of the system.

The classic example of a map is the logistic map [63], which, despite being a
discrete and one-dimensional model, by varying the single parameter of the
map, different dynamical regimes can be reached, having applications in
several areas such as physics [60, 61], biology [64], economy [65],
electronics [66], etc. The map is described by the following equation

$x_{t+1}=rx_{t}(1-x_{t}),$ (2.19)

where $r$ is a positive parameter, and $x_{t}$ is the ratio of existing
population to the maximum possible population which varies from $0\leq
x_{t}\leq 1$ if $0\leq r\leq 4$. In this sense, the logistic equation
describes population growth, where the growth rate is controlled by the
parameter $r$. In Fig. 2.4 it is studied the iterations for the logistic map
considering different values of $r$. In panel (a) $r=0.9$, the map starts from
the initial condition and evolves to a $x_{t}=0$, that is, the extinction of
the population. In panel (b) $r=1.5$, the map evolves to an equilibrium point
$x_{t}=1/3$. In panel (c) $r=3.5$, the system evolves to an orbit of period
$2$. In panel (d) $r=4$ the dynamics of the map presents a great variability
of values of $x_{t}$, indicating an absence of period.

Figure 2.4: Iterations of the logistic map for different values of
$\mathbf{r}$. (a) For $r=0.9$ after a transient time the system reaches
$x_{t}=0$. (b) For $r=1.5$ the system reaches a fixed point. (c) For $r=3.5$
the system presents a $2$-periodic dynamics. (d) For $r=4$, there is no
specific period. (e) Bifurcation diagram of the logistic map, after discarding
the first $100$ iterations to avoid the transient effects, is plotted for all
iterations $x_{t}$ as a function of $r$.

A more robust way to understand the role of a control parameter is through a
bifurcation diagram, which is a graphical representation of the qualitative
behavior of a dynamical variable ($x_{t}$) as a function of a control
parameter ($r$). Figure 2.4 (e) presents the bifurcation diagram as a function
of $r$, in which for each value of $r$, the system evolves by $500$
iterations. Discarding the first $100$ iterations, more than the necessary
time for the system to reach the equilibrium point, known as transient time,
from this point all map iterations are plotted as a function of $r$. For $r<1$
the system evolves to $x_{t}=0$. For $1<r<3$ the system reaches an equilibrium
point. For $3<r<1+\sqrt{6}$ the orbit that was periodic is doubled, generating
an orbit of period 2, where $x_{t+2}=f(f(x_{t}))$. For $1+\sqrt{6}<r\lesssim
3.54409$ each orbit splits into two, generating orbits of period $4$. As
$r>3.54409$ is increased orbits of period $8$ appear, and in sequence
$16,\,32,\,\dots$ This cascade of period doublings occurs repeatedly until
$r\approx 3.56995$ where an infinite number of values of $x_{t}$ appear. This
aperiodicity is one of the characteristics of chaotic behavior. For
$r=1+\sqrt{8}$ the system loses chaoticity, and presents a $3$-period
dynamics, these orbits become $6$-periodic orbit, until for $3.57\lesssim
r\leq 4$ the system presents, again, chaotic behavior [60].

For discrete systems, a fixed point is defined as the point that maps itself,
from Eq. (2.18)

$\mathbf{x}^{*}=\mathbf{f}(\mathbf{x}^{*}).$ (2.20)

Considering a point in the vicinity of the fixed point
$\delta\mathbf{x}_{t}=\mathbf{x}_{t}-\mathbf{x}^{*}$, and expanding
$\mathbf{f}(\mathbf{x})$ in Taylor series

$\mathbf{f}(\mathbf{x}_{t})=\mathbf{f}(\mathbf{x}^{*}+\delta\mathbf{x}_{t})=\mathbf{f}(\mathbf{x}^{*})+\frac{\partial\mathbf{f}(\mathbf{x}_{t})}{\partial\mathbf{x}_{t}}\bigg{|}_{\mathbf{x^{*}}}\delta\mathbf{x}_{t}+\mathcal{O}(\delta\mathbf{x}_{t}^{2}),$

disregarding the second-order (and higher) terms, and isolating the derivative

$\mathbf{f}(\mathbf{x}^{*}+\delta\mathbf{x}_{t})-\mathbf{f}(\mathbf{x}^{*})\approx\frac{\partial\mathbf{f}(\mathbf{x}_{t})}{\partial\mathbf{x}_{t}}\bigg{|}_{\mathbf{x^{*}}}\delta\mathbf{x}_{t},$

and

$\mathbf{x}_{t+1}=\delta\mathbf{x}_{t+1}+\mathbf{x}^{*},$ (2.21)

isolating the term $\delta\mathbf{x}_{t+1}$ and substituting the Eqs. (2.18)
and (2.20) at Eq. (2.21),

$\delta\mathbf{x}_{t+1}=\frac{\partial\mathbf{f}(\mathbf{x}_{t})}{\partial\mathbf{x}_{t}}\bigg{|}_{\mathbf{x^{*}}}\delta\mathbf{x}_{t},$
(2.22)

and taking into account that the term of the partial derivatives is the
Jacobian matrix of the map, evaluated at the fixed point $\mathbf{x}^{*}$

$\delta\mathbf{x}_{t+1}=\mathcal{J}(\mathbf{x}^{*})\delta\mathbf{x}_{t},$
(2.23)

and, analogously to a flow, the stability of the fixed points can be studied
from the eigenvalues of the Jacobian matrix, in which, if the modules of the
eigenvalues are $<1$, it implies that the trajectory converges to the fixed
point for $t\rightarrow\infty$ (stable). If at least one of the eigenvalues is
$>1$, $\delta\mathbf{x}_{t}$ diverges. In this sense, considering the case of
the logistic map Eq. (2.19), the fixed points satisfy the following equation,

$x^{*}=rx^{*}(1-x^{*}),$

so, the logistic map presents two fixed points $x^{(*,1)}=0$ and
$x^{(*,2)}=1-1/r$. Due to the restriction on the map domain ($0\leq x_{t}\leq
1$), $x^{(*,2)}$ is fixed point for $r\geq 1$.

Differentiating the Eq. (2.19)

$f^{\prime}(x_{t})=r(1-2x_{t}),$ (2.24)

and, the first fixed point,

$f^{\prime}(x^{(*,1)}=0)=r,$ (2.25)

and $x^{(*,1)}$ is stable for $r<1$, and unstable $r>1$. The second fixed
point,

$f^{\prime}(x^{(*,2)}=1-1/r)=2-r,$ (2.26)

that is, $x^{(*,2)}$ is stable for $1<r<3$, and unstable for $r>3$. It is
noteworthy that for the study of orbits with period $>1$, it is necessary to
analyze the stability of the fixed points of subsequent iterations, e.g.,
$x_{t+2}=f(f(x_{t}))$ for the case of period orbits $2$, where a fixed point
is the point that maps itself every two iterations $x^{*}=f(f(x^{*}))$.

###  2.2 Chaos in dynamical systems

As shown at the beginning of this chapter with the Lorenz system, one of the
main characteristics of a chaotic motion is the sensitivity to initial
conditions, being two conditions arbitrarily close, after a sufficiently long
time the trajectories tend to diverge exponentially from each other [60, 61,
67]. In nonlinear dynamics, systems that exhibit such behavior are often
studied.

Considering a $D$-dimensional system
$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})$. Given an arbitrary initial
condition $\mathbf{x}(0)$, exists a hypersphere of radius $\Gamma_{0}$ of
infinitesimal volume centered on $\mathbf{x}(0)$ that involves infinitely
close initial conditions to $\mathbf{x}(0)$. The temporal evolution of this
system comes from both $\mathbf{x}(t)$ and the rate of expansion (or
retraction) of the axis of the hypersphere, transforming it into an ellipsoid
whose directions are given by the set of $D$ vectors
$\mathbf{\Gamma}_{1}(t),\mathbf{\Gamma}_{2}(t),\cdots,\mathbf{\Gamma}_{D}(t)$,
as illustrated in the two-dimensional example in Fig. 2.5.

Figure 2.5: Temporal evolution of a bidimensional sphere of initial
conditions. For a set of initial conditions in a sphere of the radios
$\mathbf{\Gamma}_{0}$, the expansion and the retraction of the directions
$\mathbf{\Gamma}_{1}(t)$, $\mathbf{\Gamma}_{2}(t)$ indicates the chaotic
divergence in the initial conditions.

To distinguish regular and chaotic dynamics for a multidimensional system, the
Lyapunov spectrum is computed using Benettin’s algorithm [59, 68] and
including a Gram–Schmidt re-orthonormalization procedure [59]. The method
consists of the computation in parallel of the evolution of the system
$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})$ and the set of vectors
$\mathbf{\Gamma}(t)$. It is defined a matrix $\mathcal{M}$ in which the
elements $\mathcal{M}_{ij}$ describe the $j$-th component of the $i$-th vector
$\mathbf{\Gamma}_{i}(t)$,

$\displaystyle\mathcal{M}(t)=\begin{bmatrix}\mathcal{M}_{1,1}(t)&\mathcal{M}_{1,2}(t)&\cdots&\mathcal{M}_{1,D}(t)\\\
\mathcal{M}_{2,1}(t)&\mathcal{M}_{2,2}(t)&\cdots&\mathcal{M}_{2,D}(t)\\\
\vdots&\vdots&\ddots&\vdots\\\
{\mathcal{M}_{D,1}(t)}&{\mathcal{M}_{D,2}(t)}&\cdots&{\mathcal{M}_{D,D}(t)}\\\
\end{bmatrix},$ (2.27)

where the temporal evolution of the matrix is based on the Jacobian matrix

$\dot{\mathcal{M}}=\mathcal{J}(\mathbf{x})\mathcal{M},$ (2.28)

and the initialization of $\mathcal{M}(0)$ is equal to a identity matrix
$\mathbf{I}$.

In order to kept the vectors $\mathbf{\Gamma}(t)$ in an orthogonal direction,
in each integration step it is used the Gram-Schmidt reorthonormalization
procedure on the vector frame [59]. After discard the transient effects, the
Lyapunov spectrum are evaluated using the norm of the vectors

$\lambda_{i}=\frac{1}{t_{\mathrm{f}}-t_{\mathrm{i}}}\sum_{t=t_{\mathrm{i}}}^{t_{\mathrm{f}}}\ln|\mathbf{\Gamma}_{i}(t)|,\;\;\;\;\;\;\;\;\;\;\;\;\;i=1,2,\cdots,D,$
(2.29)

where $t_{\mathrm{i}}$ and $t_{\mathrm{f}}$ are the initial and final times of
computation, respectively, to obtain the stationary solutions of the dynamical
system. If at least one of the exponents $\lambda_{i}>0$, it usually taken as
an indication that the system is chaotic, since there is an exponential
divergence at the $i$-th direction, even if $\lambda_{j}\leq
0,\,\forall\,j\neq i$. In this sense, it is useful to define the largest
Lyapunov exponent ($\Lambda$) where

$\Lambda=\mathrm{max}\\{\lambda_{1},\lambda_{2},\cdots,\lambda_{D}\\},$ (2.30)

which means that if $\Lambda>0$ the system depicts a chaotic behavior.

In the case of a unidimensional discrete system described by Eq. (2.18), the
($n+1$)-th iteration of the map can be written in terms of the initial
condition $x_{0}$

$x_{n+1}=f(x_{n})=f^{n}(x_{0}).$ (2.31)

Consider a nearby point $x_{0}+\delta_{0}$, where the initial separation
$\delta_{0}$ is infinitesimal and $\delta_{n}$ the separation after $n$
iterations which can be described by

$\delta_{n}=f^{n}(x_{0}+\delta_{0})-f^{n}(x_{0}),$ (2.32)

and assuming that the separation evolves exponentially

$|\delta_{n}|\approx|\delta_{0}|e^{n\lambda},$ (2.33)

isolating $\lambda$

$\lambda\approx\frac{1}{n}\ln\left|\frac{\delta_{n}}{\delta_{0}}\right|,$

and by Eq. (2.32)

$\lambda\approx\frac{1}{n}\ln\left|\frac{f^{n}(x_{0}+\delta_{0})-f^{n}(x_{0})}{\delta_{0}}\right|,$
(2.34)

taking the limit $\delta_{0}\rightarrow 0$ Eq. (2.34) yields

$\lambda\approx\frac{1}{n}\ln\left|\lim_{\delta_{0}\rightarrow
0}\frac{f^{n}(x_{0}+\delta_{0})-f^{n}(x_{0})}{\delta_{0}}\right|=\frac{1}{n}\ln|(f^{n}(x_{0}))^{\prime}|,$
(2.35)

expanding the logarithm term using the chain rule

$\lambda\approx\frac{1}{n}\ln\left|\prod_{i=0}^{n-1}f^{\prime}(x_{i})\right|=\frac{1}{n}\sum_{i=0}^{n-1}\ln|f^{\prime}(x_{i})|,$

and, if the limit of $n\rightarrow\infty$ exists, this limit is defined to be
the Lyapunov exponent of the map

$\lambda=\lim_{n\rightarrow\infty}\left\\{\frac{1}{n}\sum_{i=0}^{n-1}\ln|f^{\prime}(x_{i})|\right\\}.$
(2.36)

In the same way that for the Lyapunov spectrum, a positive exponent
$\lambda>0$ is a signature of chaos. In the case of the logistic map, explored
in the last section,

$\lambda=\lim_{n\rightarrow\infty}\left\\{\frac{1}{n}\sum_{i=0}^{n-1}\ln|r(1-2x_{i})|\right\\},$
(2.37)

the lyapunov exponent depends on the control parameter $r$, which is expected
since as shown in Fig. 2.4, different values of $r$ depict different dynamical
behaviors. Figure 2.6 depicts the bifurcation diagram (left scale) and the
lyapunov exponent $\lambda$ (red line, right scale) as a function of $r$. The
chaotic dynamics takes place at $r\approx 3.57$ due to a known route to chaos
called period doubling cascade [60] which will be presented in more details in
the next section.

Figure 2.6: Lyapunov exponent of the logistic map. Bifurcation diagram of the
logistic map (left scale) and the Lyapunov exponent $\lambda$ (right scale).

###  2.3 Routes to chaos

As shown in the last section, the chaotic behavior can emerge in the dynamical
system with the simplest change of a parameter, which causes a change in the
stability of fixed points of the system. In this section, it is explored three
of the most famous route to chaos: period doubling, intermittency, and crisis
[60].

####  2.3.1 The period doubling cascade

In a doubling period bifurcation, a stable fixed point loses stability to an
attracting 2-period orbit as the parameter reaches a critical value [60]. An
illustration of this route can be seen in panel (a) of Fig. 2.7. In the
logistic map, the fixed point $x^{(*,2)}=1-\frac{1}{r}$ is stable for $1<r<3$
since $|f^{\prime}(x^{*,2})|=|2-r|$ is $<1$ at this range. For $r>3$ the fixed
point loses stability since $|f^{\prime}(x^{*,2})|>1$, and, simultaneously,
$f(f(x))$ leads the creation of a 2-periodic stable fixed point (which is not
a fixed points for $f(x)$). As $r$ is increased to $1+\sqrt{(}6)\approx
3.449$, the 2-period fixed point loses stability leading the creation to a
4-period fixed point (due to the stability of $f^{4}(x)$), etc. This process
continues indefinitely; at each bifurcation, the periodic orbit is replaced by
a new attracting periodic orbit of twice the period, producing an infinite
cascade of period doublings with ranges $r_{m-1}<r<r_{m}$, in which a $2^{m}$
orbit is stable. Panel (b) of Fig. 2.7 presents the evolution of successive
doubling of periods on a logarithmic scale. Since the length of $r$ of the
range of stability for an orbit of period $2^{m}$ decreases as $m$ increases,
there is a saturation value $r_{\infty}$ where there is an accumulation point
of an infinite number of period-doubling bifurcations

$r_{\infty}=\lim_{m\rightarrow\infty}r_{m}\approx 3.5699,$ (2.38)

with the property of

$\lim_{m\rightarrow\infty}\frac{r_{m}-r_{m-1}}{r_{m+1}-r_{m}}\equiv\hat{\delta}\approx
4.6692,$ (2.39)

where $\hat{\delta}$ is called Feigenbaum constant [60, 69]. It is remarkable
that $\hat{\delta}$ is a universal constant for functions approaching chaos
via period doubling.

Figure 2.7: Period doubling route to chaos. (a) A stable point loses stability
when a control parameter reaches a critical value, creating a 2-periodic
stable fixed point. (b) The evolution of successive doubling of periods. The
length of $r$ of the range of stability for an orbit of period $2^{m}$
decreases as $m$ increases, collapsing an accumulation of an infinite number
of period doublings when $r_{\infty}$ making the system converges to a chaotic
attractor.

####  2.3.2 Intermittency route to chaos

Considering a generic system, the intermittency route to chaos refers to how a
periodic orbit is replaced by a chaotic attractor when the control parameter,
namely $r$ reaches a critical value $r^{*}$ [60, 70]. Supposing that the
periodic orbit exists for $r>r^{*}$, for values $r\lesssim r^{*}$ the periodic
orbit no longer exists and it is possible to see “nearly-period” orbits which
are intermittently interrupted by a finite chaotic behavior. An example of
this behavior can be seen in the logistic map where for
$r^{*}=1+2\sqrt{2}\approx 3.8284$ a 3-periodic orbit is created, and the
chaotic behavior occurs for $r<r^{*}$. Figure 2.8 (a) exhibits a magnification
of the bifurcation diagram of the logistic map in the critical value $r^{*}$
when the 3-period orbit appears. Panels (b – d) show an illustration of the
intermittent behavior, where the system presents a 3-periodic orbit for
$r=2.83$ (b). At panel (c) for $r=2.8282$, it is possible to see a remnant
nearly-3-period (approximately regular), this orbit is interrupted by a
chaotic behavior that occurs intermittently. And panel (d) or $r=2.82$ there
is no apparent evidence of the 3-period orbit anymore.

Figure 2.8: Intermittency route to chaos. (a) Bifurcation diagram for the
logistic map as a function of $r$. (b) 3-periodic orbit for $r=2.83$. (c) For
$r=2.8282\approx r^{*}$ a remnant nearly-3-period orbit is seen (the
3-periodic orbit no longer exists), this orbit is interrupted by chaos that
occurs intermittently. (d) For $r=2.82$ there is no apparent evidence of the
3-period orbit anymore.

Different from the doubling period cascade, where the stable fixed point loses
stability and creates a 2-periodic stable fixed point, and the chaotic
attractor is achieved due to the successive creation of an infinite number of
period-doubling bifurcations. In this route, the stable fixed point either
becomes unstable or is destroyed as the control parameter reaches a critical
parameter. In this sense, three types of intermittency are distinguished
corresponding to three types of bifurcations [70]. Figure 2.9 depicts an
illustration of the three types of intermittency route to chaos. Type I:
saddle-node bifurcation, where stable and unstable fixed points of a dynamical
system collide and annihilate each other; Type II: Hopf bifurcation, where a
stable fixed point of a dynamical system loses stability; Type III: inverse
period-doubling bifurcation, where there is a stable fixed point with
2-periodic unstable fixed point, the stable fixed point collides with the
unstable fixed points and lost its stability.

Figure 2.9: Three types of intermittency route to chaos. In this route to
chaos, the chaotic state is reached in three different ways. Figure inspired
by Fig. 8.6 of Ref [60].

It is possible to define $\langle\tau_{\mathrm{int}}(r)\rangle$ as the mean
characteristic time of occurrence of an intermittency behavior with

$\lim_{r\rightarrow r^{*}}\langle\tau_{\mathrm{int}}(r)\rangle=\infty$ (2.40)

where $r$ is studied in the chaotic regime. In this sense, for each type of
intermittency case, this mean time decays with

$\langle\tau_{\mathrm{int}}(r)\rangle\sim\begin{dcases}|r-r^{*}|^{-1/2}\;\;&\mathrm{for\;Type\;I},\\\
|r-r^{*}|^{-1}\;\;&\mathrm{for\;Type\;II\;and\;III}.\end{dcases}$ (2.41)

####  2.3.3 Crisis

In the crisis route, the chaotic attractor is changed with the variations in
the control parameter [60]. In particular, in the boundary crisis, the chaotic
attractor is annihilated with the collision of an unstable periodic orbit on
its basin boundary [60, 71, 72]. For example, for a dynamical system that is
chaotic for $r<r^{*}$ and periodic, otherwise. For values $r$ sightly greater
than the critical value the attractor no longer exists but is replaced by a
chaotic saddle that does not attract trajectories but allows a chaotic
transient. This means that an initialization of the system in some region of
the phase space inside the chaotic saddle spend some time behaving as chaotic
until leaving the attractor and never return. The mean transient lifetime that
the system spent near the chaotic transient
$\langle\tau_{\mathrm{crisis}}\rangle$ becomes longer as $r\rightarrow r^{*}$,
the decay with the distance of the critical point following

$\langle\tau_{\mathrm{crisis}}(r)\rangle\propto|r-r^{*}|^{-\kappa}$ (2.42)

where $\kappa$ is called critical exponent of the crisis [60]. A simple
example of a crisis is given by the logistic map, where the chaotic attractor
exists for $r=4$, but for $r>4$ the chaotic attractor is replaced by a chaotic
transient and for a initial condition inside the basin of attraction, i.e.
$x_{0}$ $\in[0,1]$, after $\langle\tau_{\mathrm{crisis}}\rangle$ the system
diverges to $-\infty$. Figure 2.10 presents the numerical result of
$\langle\tau_{\mathrm{crisis}}(r)\rangle$ for the logistic map, in this case
as well as for other one dimensional maps the critical exponent of the crisis
is $\kappa=1/2$ (solid magenta line) [60]. In the case of multidimensional
system is expected $\kappa>1/2$ [60].

Figure 2.10: Crisis route in the logistic map. The mean transient life time of
that the system spent near the chaotic transient
$\langle\tau_{\mathrm{crisis}}\rangle$. The theoretical result
$\propto|r-r^{*}|^{-1/2}$ is plotted as a line.

##  Chapter 3 Neuronal models and synapses

This chapter presents some properties of the neuron, the main cell of the
nervous system [12, 33]. Then, it is presented the Hodgkin-Huxley (HH) model
[35], where the equations to modeling the time-evolution of the membrane
potential are created based on an electrical circuit. After that, it is
presented one of the adaptations of the HH model proposed by Braun et al. [41,
57], here called the Hodgkin-Huxley-like (HH$\ell$) model, where the
adaptation consists of the addition of two ionic currents and some
temperature-dependence parameters, that makes it possible the neuron to depict
a different dynamical feature, called burst dynamics. Then, it is presented
the Hindmarsh-Rose (HR) model [40], a model of neuron of three coupled first-
order differential equations. Therefore, the model proposed by Dante R.
Chialvo [37], in which, with two discrete equations, the dynamics of the
action potential seen in the neuron can be reproduced. Lastly, it is briefly
presented the equations that ruled the synaptic interactions among neurons,
which enables the transfer of information between the presynaptic neuron to
the postsynaptic neuron [12].

###  3.1 The action potential

The neuron is one of the main cells of the nervous system, being responsible
for the conduction of electrical impulses. In a general point of view, the
neuron can be divided into three regions: soma, which contains the nuclei and
cytoplasm; dendrites, whose main function is the reception of signals which
comes from other neurons; axons, which carry the signals from the soma to
other neurons [12] (Fig. 3.1).

Figure 3.1: Basic anatomy of a neuron. A typical neuron can be separated into
three basic regions: the soma, dendrites, and axons [12]. This figure is
adapted from [73].

The main function of the neuron is the transmission of action potentials,
which are electrical signals which propagate information at the nerve system.
The action potential consists of a depolarization followed by repolarization,
depicting a spike shape. These variations at the membrane potential occur due
to the variations of the ion concentrations between the intracellular and
extracellular media. In this sense, the membrane potential can be defined as

$V(t)=V_{\mathrm{int}}-V_{\mathrm{ext}},$ (3.1)

where $V_{\mathrm{int}}$ and $V_{\mathrm{ext}}$ are the potentials of the
intracellular and extracellular media, respectively.

The ions permeate the neural membrane by proteins which work as ion channels.
There are gated channels and non-gated channels. While non-gated channels are
always open, allowing the entry and exit of ions from the intracellular side,
gated channels have potential-dependence gates, that is, the permissiveness of
the channels depends directly on the potential of the membrane. The
predominant ions found on either side of the cell membrane are potassium ions
($\mathrm{K}^{+}$), sodium ($\mathrm{Na}^{+}$), and chlorine
($\mathrm{Cl}^{-}$). One illustration of the ionic channels is presented in
Fig. 3.2, where one channel is open allowing the exit of a K+ ion, the other
is closed, and a non-gated channel allows the entrance of a Cl- ion.

Figure 3.2: Cell membrane representation. The cell membrane contains proteins,
that function as channels that allow the ions to move through it. The ion flux
changes the concentration of the intracellular (and extracellular) side,
inducing variations in the membrane potential.

At the equilibrium point also called the resting state, there are no ionic
changes between the media. In this case, the intracellular side contains in
majority K+ ions in comparison with Na+ and Cl- that are more abundant in the
extracellular side [33]. Hence, the membrane potential is determined primarily
by the K+ resting potential,

$V(t)=V_{\mathrm{int}}-V_{\mathrm{ext}}=V_{\mathrm{eq}}\approx-$70\text{\,}\mathrm{mV}$.$

When the cell is stimulated above a threshold, the Na+ channels open allowing
the entrance of ions inside the cell, this influx of Na+ tends to depolarizes
the cell, resulting in a positive variation in $V(t)$. This abrupt increase of
$V$ inverts the polarity of the cell, closing the Na+ channels and opening the
K+ channels, allowing the efflux of K+ ions to repolarizes the cell, which
takes the membrane potential to a level below the $V_{\mathrm{eq}}$
(hyperpolarization). Until the K+ channels close up again, the membrane is in
a refractory stage. During this time, the cell pumps the exchange excess Na+
ions inside the cell with excess K+ ions outside the cell, and a new cycle can
be started if the stimulus is kept [33]. The Fig. 3.3 represents the first
intracellular record of an action potential [12] of the giant axon of squids
Loligo forbesi, published by Hodgkin and Huxley at 1939 [34] where temporal
pulses was applied each $2\text{\,}\mathrm{ms}$.

Figure 3.3: Action potential of the giant axon of squids. Action potential
between inside and outside of axon was recorded by Hodgkin and Huxley at
$1939$ using a micro-electrode consisted of a glass tube filled with seawater.
Figure taken from [34]. Adapted.

###  3.2 Neuron models

####  3.2.1 The Hodgkin-Huxley (HH) model

Alan Lloyd Hodgkin and Andrew Fielding Huxley were the firsts to describe
mathematically a regenerative current that generates an action potential. They
were awarded the Nobel Prize in Physiology or Medicine in 1963, together with
Sir John Carew Eccles for their discoveries about the ionic mechanisms
involved in the excitation and inhibition of the membrane of nerve cells [74].

To describe the action potential of the cell, Hodgkin and Huxley have used a
circuit model, considering the contribution of two ionic currents (K+ and Na+)
plus a passive current which takes the contribution of the non-gated channels
and the less abundant ions in the media. The circuit, illustrated in Fig. 3.4,
is composed of three components: resistors, which represent the ion channels;
sources, which represent the ion concentration gradient; and capacitors, which
represent the charge that can be stored in the membrane [33].

Figure 3.4: Equivalent circuit model. The circuit consists of three
components: resistors, which represent the ion channels; sources, which
represent the ion concentration gradient; and capacitors, which represent the
charge that can be stored in the membrane [33].

The charge stored in the capacitor is proportional to the membrane potential,
where the proportionality constant is the capacitance

$Q=C_{\mathrm{M}}V,$ (3.2)

where $C_{\mathrm{M}}$ is the specific capacitance of the membrane, measured
in
$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{F}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$.
$V$ is the membrane potential measured in $\text{\,}\mathrm{mV}$.
Differentiating both sides in relation to time

$\frac{dQ}{dt}=i_{\mathrm{cap}}=C_{\mathrm{M}}\frac{dV}{dt}$ (3.3)

where $i_{\mathrm{cap}}$ is the specific current referring to the capacitance
divided by the area.

Applying Kirchhoff’s law to the circuit at Fig. 3.4 implies that

$0=i_{\mathrm{cap}}+I_{\mathrm{K}}+I_{\mathrm{Na}}+I_{\mathrm{leak}}-I(t),$
(3.4)

where $I_{\mathrm{K}}$ and $I_{\mathrm{Na}}$ are the ionic currents,
$I_{\mathrm{leak}}$ is the passive current, and $I(t)$ is an external applied
current. Since the currents are Ohmic it can be rewritten as
$I_{\mu}=g_{\mu}(V-E_{\mu})$. Substituting $i_{\mathrm{cap}}$ with Eq.(3.3)

$C_{\mathrm{M}}\frac{dV}{dt}=-g_{\mathrm{K}}(V-E_{\mathrm{K}})-g_{\mathrm{Na}}(V-E_{\mathrm{Na}})-g_{\mathrm{leak}}(V-E_{\mathrm{leak}})+I(t),$
(3.5)

in which $E_{\mathrm{K}}$, $E_{\mathrm{Na}}$ and $E_{\mathrm{leak}}$ are the
resting potential of each channel, and $g_{\mathrm{K}}$, $g_{\mathrm{Na}}$ and
$g_{\mathrm{leak}}$ are the conductances of each channel.

The conductances of the gated channels are voltage-dependent and are directly
related to the permissively of the ion channels, which means that it rules the
probability of opening/close the channels. Using voltage-clamp techniques
Hodgkin and Huxley have fitted the equations for each conductance as

$\displaystyle g_{\mathrm{K}}$ $\displaystyle=$
$\displaystyle\overline{g}_{\mathrm{K}}n^{4},$ (3.6) $\displaystyle
g_{\mathrm{Na}}$ $\displaystyle=$
$\displaystyle\overline{g}_{\mathrm{Na}}m^{3}h,$ (3.7) $\displaystyle
g_{\mathrm{leak}}$ $\displaystyle=$
$\displaystyle\overline{g}_{\mathrm{leak}},$ (3.8)

where $\overline{g}_{\mathrm{K}}$, $\overline{g}_{\mathrm{Na}}$ and
$\overline{g}_{\mathrm{leak}}$ are the maximum values of conductance of each
current, $n$, $m$ are activation functions and $h$ an inactivation function,
these functions are related to the probabilities of opening and closing the
channels.

Considering $\alpha_{\mu}(V)$ as the probability of the channel being open,
and $\beta_{\mu}(V)$ closed (with $\mu=m,\,n,\,h$) these functions satisfy the
following equations

$\frac{d\mu}{dt}=\alpha_{\mu}(V)(1-\mu)-\beta_{\mu}(V)\mu,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mu=n,{m},{h},$
(3.9)

in which $\alpha_{\mu}$ and $\beta_{\mu}$ are described

$\displaystyle\alpha_{\mathrm{n}}(V)$ $\displaystyle=$
$\displaystyle\frac{0.01(V+55)}{(1-\exp[-(V+55)/10])},$ (3.10)
$\displaystyle\alpha_{\mathrm{m}}(V)$ $\displaystyle=$
$\displaystyle\frac{0.1(V+40)}{(1-\exp[-(V+40)/10])},$ (3.11)
$\displaystyle\alpha_{\mathrm{h}}(V)$ $\displaystyle=$ $\displaystyle
0.07\exp[-(V+65)/20],$ (3.12) $\displaystyle\beta_{\mathrm{n}}(V)$
$\displaystyle=$ $\displaystyle 0.125\exp[-(V+65)/80],$ (3.13)
$\displaystyle\beta_{\mathrm{m}}(V)$ $\displaystyle=$ $\displaystyle
4\exp[-(V+65)/18],$ (3.14) $\displaystyle\beta_{\mathrm{h}}(V)$
$\displaystyle=$ $\displaystyle\frac{1}{(1+\exp[-(V+35)/10])},$ (3.15)

taking into account that $V$ is measured in $\text{\,}\mathrm{mV}$, and, for
simplicity, the dimensions required in the other values are omitted.

Substituting the Eqs. (3.6–3.8) at Eq. (3.5), it is achieved the Hodgkin-
Huxley equations [33]

$C_{\mathrm{M}}\frac{dV}{dt}=-\overline{g}_{\mathrm{K}}n^{4}(V-E_{\mathrm{K}})-\overline{g}_{\mathrm{Na}}m^{3}h(V-E_{\mathrm{Na}})-\overline{g}_{\mathrm{leak}}(V-E_{\mathrm{leak}})+I(t).$
(3.16)

Figure 3.5: Time-evolution of the variables of the Hodgkin-Huxley model. For
$I=$10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$
the membrane potential presents a sequence of action potentials.

The Fig. 3.5 demonstrates the dynamics of the HH model using the constants of
Table 3.1 and
$I=$10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$.
Each panel depicts the time evolution of a variable of the model. $V$ is
presented in panel (a) and exhibits a sequence of action potentials (spikes).
The other panels (b - d) depict the activation variables $n$ and $m$, and
inactivation $h$, respectively. It is noted that when $V$ reaches a threshold,
$n$ and $m$ increases quickly reaching their maximum values. In contrast, $h$
has the opposite effect, characterizing its inactivation effect, repolarizing
the cell starting a new cycle. The neuron continues with this dynamics of
periodic spikes while $I$ is kept active.

To understand the role of the applied current $I$, Figure 3.6 presents the
bifurcation diagram of the inter-spike-interval (ISI) of the neuron. The
occurrence of a spike is defined when $V$ reaches a threshold value of $-20$
$\text{\,}\mathrm{mV}$ (with positive derivative). After that, it is evaluated
the time between two consecutive spikes as a function of $I$ disregarding the
first second to avoid the transient effect. It is noted that for $I<8.39$ the
current is not strong enough to keep the neuron spiking. For $I\approx 8.39$
the model pass through a Hopf bifurcation [60] and the ISI jumps from $\infty$
to $15.55\text{\,}\mathrm{ms}$. Panel (b) and (c) depict the time evolution of
$V$ for $I=5$ and $I=10$, respectively. (b) The neuron depicts one spike, but
after this, the equilibrium state is recovered since the current is not strong
enough to keep the neuron spiking, as seen in panel (c).

Figure 3.6: Bifurcation diagram of the Hodgkin-Huxley model. (a) Inter-spike-interval (ISI) as a function of the applied current $I$, for $I\approx 8.39$ there is a Hopf bifurcation point where the ISI jumps from $\infty\rightarrow$15.55\text{\,}\mathrm{ms}$$. (b) For $I=5$ after a transient time, the membrane potential reaches an equilibrium point. (c) $I=10$ the membrane potential depicts periodic spikes. Table 3.1: Constants used for the simulation of the Hodgkin-Huxley model [33] Membrane capacitance ($\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{F}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$) | $C_{\mathrm{M}}$ | $1$
---|---|---
Maximum conductances ($\text{\,}\mathrm{m}\mathrm{S}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$) | $\overline{g}_{\mathrm{Na}}$ | $120$
$\overline{g}_{\mathrm{K}}$ | $36$
$\overline{g}_{\mathrm{leak}}$ | $0.3$
Resting potentials ($\text{\,}\mathrm{mV}$) | $E_{\mathrm{Na}}$ | $50$
$E_{\mathrm{K}}$ | $-77$
$E_{\mathrm{leak}}$ | $-54.4$

####  3.2.2 A Hodgkin-Huxley-like (HH$\ell$) model

The model proposed by Braun et al. [41] consists in an adaptation of the HH
model in order to reproduce similar patterns observed in thermally sensitive
electroreceptors of the catfish (Ictaluris nebulosus) [57]. This dynamical
behavior is called burst, is characterized by a silent phase of near-steady-
state resting behavior alternating with an active phase of rapid spike
oscillations. Bursting occurs in the activity of some thalamic cells, e.g.,
can implicate in the generation of sleep rhythms, whereas patients with
parkinsonian tremors exhibit increased bursting activity in neurons within the
basal ganglia [33].

The adaptation consists of the addition of two sub-threshold currents with
slower frequency, and some temperature-dependence parameters. The main
equation of the model is defined as

$C_{\mathrm{M}}\frac{dV}{dt}=-I_{\mathrm{d}}-I_{\mathrm{r}}-I_{\mathrm{sd}}-I_{\mathrm{sr}}-I_{\mathrm{l}},$
(3.17)

where, again, $C_{\mathrm{M}}$ is the membrane capacitance; $I_{\mathrm{d}}$
and $I_{\mathrm{r}}$ represent the classical HH ionic currents related to Na+
and K+, respectively. $I_{\mathrm{sd}}$ and $I_{\mathrm{sr}}$ are sub-
threshold currents of depolarization and repolarization, respectively, which
can be related to the Ca2+ ions [75]. $I_{\mathrm{l}}$ represents the passive
current which takes the contribution of the non-gated channels. The currents
are described by

$\displaystyle I_{\mathrm{d}}$ $\displaystyle=$ $\displaystyle\varrho
g_{\mathrm{d}}a_{\mathrm{d}}(V-E_{\mathrm{d}}),$ (3.18) $\displaystyle
I_{\mathrm{r}}$ $\displaystyle=$ $\displaystyle\varrho
g_{\mathrm{r}}a_{\mathrm{r}}(V-E_{\mathrm{r}}),$ (3.19) $\displaystyle
I_{\mathrm{sd}}$ $\displaystyle=$ $\displaystyle\varrho
g_{\mathrm{sd}}a_{\mathrm{sd}}(V-E_{\mathrm{sd}}),$ (3.20) $\displaystyle
I_{\mathrm{sr}}$ $\displaystyle=$ $\displaystyle\varrho
g_{\mathrm{sr}}a_{\mathrm{sr}}(V-E_{\mathrm{sr}}),$ (3.21) $\displaystyle
I_{\mathrm{l}}$ $\displaystyle=$ $\displaystyle
g_{\mathrm{l}}(V-E_{\mathrm{l}}),$ (3.22)

in which $g_{\mathrm{d}}$, $g_{\mathrm{r}}$, $g_{\mathrm{sd}}$,
$g_{\mathrm{sr}}$, and $g_{\mathrm{l}}$ are the maximum values of the
respective conductances, and $a_{\mathrm{d}}$, $a_{\mathrm{r}}$,
$a_{\mathrm{sd}}$, and $a_{\mathrm{sr}}$ are the activation functions of the
ion channels. $\varrho$ is a temperature-dependent parameter defined by

$\varrho=1.3^{\frac{T-T_{0}}{\tau_{0}}},$ (3.23)

where $T$ is the temperature of the system, $T_{0}$ and $\tau_{0}$ are
constants.

The temporal evolution of each activation term follows a differential equation

$\displaystyle\frac{da_{\mathrm{d}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\phi}{\tau_{\mathrm{d}}}(a_{\mathrm{d},\infty}-a_{\mathrm{d}}),$
(3.24) $\displaystyle\frac{da_{\mathrm{r}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\phi}{\tau_{\mathrm{r}}}(a_{\mathrm{r},\infty}-a_{\mathrm{r}}),$
(3.25) $\displaystyle\frac{da_{\mathrm{sd}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\phi}{\tau_{\mathrm{sd}}}(a_{\mathrm{sd},\infty}-a_{\mathrm{sd}},)$
(3.26) $\displaystyle\frac{da_{\mathrm{sr}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\phi}{\tau_{\mathrm{sr}}}(-\eta I_{\mathrm{sd}}-\gamma
a_{\mathrm{sr}}),$ (3.27)

in which $\phi$ represents another temperature-dependent parameter given by
$\phi=3^{(T-T_{0})/\tau_{0}}$; $\tau_{\mathrm{d}}$, $\tau_{\mathrm{r}}$,
$\tau_{\mathrm{sd}}$, $\tau_{\mathrm{sr}}$ are characteristic times which
correspond to each activation function [41]. The term $\eta$ is a factor which
relates the mixed Na/Ca current to the increment of the intracellular Ca2+,
and $\gamma$ is the decrease rate of Ca2+ [75].

And finally the functions

$\displaystyle a_{\mathrm{d},\infty}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\exp[-s_{\mathrm{d}}(V-V_{0\mathrm{d}})]},$ (3.28)
$\displaystyle a_{\mathrm{r},\infty}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\exp[-s_{\mathrm{r}}(V-V_{0\mathrm{r}})]},$ (3.29)
$\displaystyle a_{\mathrm{sd},\infty}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\exp[-s_{\mathrm{sd}}(V-V_{0\mathrm{sd}})]},$ (3.30)

where $s_{\mathrm{d}}$, $s_{\mathrm{r}}$, $s_{\mathrm{sd}}$,
$V_{0\mathrm{d}}$, $V_{0\mathrm{r}}$ and $V_{0\mathrm{sd}}$ are parameters.
The Fig. 3.7 presents the time-evolution of the $5$ variables of the model,
using the parameters of the Table 3.2. In panel (a) the membrane potential $V$
depicts a bursting oscillation, a sequence of $4$ spikes followed by a resting
time. At panels (b) and (c) the activation variables $a_{\mathrm{d}}$ and
$a_{\mathrm{r}}$ present a similar dynamics. In panels (d) and (e) the
activation variables $a_{\mathrm{sd}}$ and $a_{\mathrm{sr}}$ show a slower
oscillation which follow the bursting oscillation.

Figure 3.7: Time-evolution of the variables of the HH$\ell$ model. For the values of the Table 3.2 the membrane potential $V(t)$ (a) depicts the burst dynamics: a sequence of spikes followed by a resting time. The activation variables $a_{\mathrm{d}}$ (b) and $a_{\mathrm{r}}$ (c) follow the $V$ fast dynamics, while $a_{\mathrm{sd}}$ (d) and $a_{\mathrm{sr}}$ (e) present a slower dynamics. In particular, $a_{\mathrm{sr}}$ oscillates with the bursting frequency slightly influenced by the spike occurrences. Table 3.2: Values of the parameters of the HH$\ell$ model [76] Membrane Capacitance ($\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{F}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$) | $C_{\mathrm{M}}$ | $1$
---|---|---
Maximum conductances ($\text{\,}\mathrm{m}\mathrm{S}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$) | ${g}_{\mathrm{d}}$ | $1.5$
${g}_{\mathrm{r}}$ | $2$
${g}_{\mathrm{sd}}$ | $0.25$
${g}_{\mathrm{sr}}$ | $0.4$
${g}_{\mathrm{l}}$ | $0.1$
Characteristic times ($\text{\,}\mathrm{ms}$) | $\tau_{\mathrm{d}}$ | $0.05$
$\tau_{\mathrm{r}}$ | $2$
$\tau_{\mathrm{sd}}$ | $10$
$\tau_{\mathrm{sr}}$ | $20$
Reversal potentials $($\text{\,}\mathrm{mV}$)$ | $E_{\mathrm{d}}$ | $50$
$E_{\mathrm{r}}$ | $-90$
$E_{\mathrm{sd}}$ | $50$
$E_{\mathrm{sr}}$ | $-90$
$E_{\mathrm{l}}$ | $-60$
$V_{0\mathrm{d}}$ | $-25$
$V_{0\mathrm{r}}$ | $-25$
$V_{0\mathrm{sd}}$ | $-40$
Temperature parameters ($\text{\,}\mathrm{\SIUnitSymbolCelsius}$) | $T_{0}$ | $25$
$T$ | $13$
$\tau_{0}$ | $10$
Other parameters | $s_{\mathrm{d}}$ | $0.25\text{\,}\mathrm{m}\mathrm{V}^{-1}$
$s_{\mathrm{r}}$ | $0.25\text{\,}\mathrm{m}\mathrm{V}^{-1}$
$s_{\mathrm{sd}}$ | $0.09\text{\,}\mathrm{m}\mathrm{V}^{-1}$
$\eta$ | $0.012\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$
$\gamma$ | $0.17$

Considering different parameters it is possible to change the dynamics of the
model. Figure 3.8 displays the bifurcation diagram of the ISI as a function of
the temperature $T$ measured in $\text{\,}\mathrm{\SIUnitSymbolCelsius}$. The
spike is computed in the same way as the HH model when $V$ crosses $-20$
$\text{\,}\mathrm{mV}$ (with positive derivative). The bifurcation diagram
represents a classical route to chaos called doubling period [60]. The first
doubling occurs at $T\approx 6.765$ and the second $T\approx 7.195$, and for
$T>7.31$ the system reach the chaoticity [77]. For $T\approx 10.66$ the ISI
depicts extremely higher values, in Ref. [77], the authors show evidences that
this behavior is associated with a homoclinic bifurcation. For $T\gtrsim 13$
the system loses its chaoticity but presents two branches of ISI associate
with the fast and slow temporal scales, associated with the spikes and bursts
respectively. Panels (b – d) depict the dynamics of $V$ for $T=5,\,10,\,15$,
respectively. For $T=5$ the neuron depicts periodic spikes in a similar way
that the HH model. For $T=10$ the neuron presents chaotic bursts, and $T=15$
periodic bursts. It should be noted that the transitions to chaos can also be
achieved with the change of other parameters, as the maximum conductances for
example.

Figure 3.8: Bifurcation diagram of the HH$\ell$ model. (a) Inter-spike-
interval (ISI) as a function of the temperature $T$ measured
$\text{\,}\mathrm{\SIUnitSymbolCelsius}$. The membrane potential dynamics
consists of (b) periodic spikes $T=5$, (c) chaotic bursts $T=10$, and (d)
periodic bursts $T=15$.

####  3.2.3 The Hindmarsh-Rose model

In addition to the model described in the previous section, in the literature,
it is possible to find several models that present the bursting dynamics [33,
38, 40, 42]. Most of the models are purely dynamical, which means that,
besides the dynamics, there is no relation with the real neuron. In this
sense, the Hindmarsh-Rose (HR) model [40] is composed of three dimensionless
variables ($x,y,z$) in which the $x$ variable depicts a similar-bursting
dynamics observed in the membrane potential of real neurons [33]. The model is
described by three differential equations

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle
y-ax^{3}+bx^{2}-z+\mathcal{I},$ (3.31) $\displaystyle\frac{dy}{dt}$
$\displaystyle=$ $\displaystyle c-dx^{2}-y,$ (3.32)
$\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle
r[s(x-x_{\mathrm{r}})-z],$ (3.33)

where $a$, $b$, $c$, $d$, $r$, $x_{\mathrm{r}}$, $s$, and $\mathcal{I}$ are
parameters of the model. The $\mathcal{I}$ parameter acts similarly to an
external current being applied to the neuron, often used as a control
parameter [40]. In this sense, Fig. 3.9 depicts a bifurcation diagram of the
ISI as a function of $\mathcal{I}$, using the other parameters of Table 3.3.
The spike is evaluated when $x$ crosses $0$ (with positive derivative). This
model depicts periodic bursts for $\mathcal{I}<2.92$, chaotic bursts for
$2.92<\mathcal{I}\lesssim 3.4$, and periodic spikes for $\mathcal{I}\gtrsim
3.4$. Panels (b – d) present the dynamics of the three variables of the model
for $\mathcal{I}=2.92$ to clarify the chaoticity of the bursting dynamics of
this model.

Table 3.3: Constants used for the simulation of the Hindmarsh-Rose model [40]. $a=1$ | $b=3$ | $c=1$ | $d=5$ | $x_{\mathrm{r}}=-8/5$ | $s=4$ | $r=0.006$
---|---|---|---|---|---|---
Figure 3.9: Bifurcation diagram of the HR model. (a) Inter-spike-interval
(ISI) as a function of the applied current $\mathcal{I}$. The time of each
spike is calculated when $x$ crosses $0.0$ with a positive derivative. (b – d)
Dynamics of each variable of the model, $x$, $y$ and $z$ respectively. Using
$\mathcal{I}=2.92$ and parameters of Table 3.3.

####  3.2.4 The Chialvo model

The model proposed by Dante R. Chialvo [37] contains two dimensionless map
equations described by

$\displaystyle x_{t+1}$ $\displaystyle=$ $\displaystyle
x_{t}^{2}\exp(y_{t}-x_{t})+k,$ (3.34) $\displaystyle y_{t+1}$ $\displaystyle=$
$\displaystyle ay_{t}-bx_{t}+c,$ (3.35)

where $x_{t}$ acts like the potential of a membrane, and $y_{t}$ like an
recovery variable. The model depends of four parameters, $a$, $b$, $c$ e $k$.
In particular $k$ acts like an external current. In Fig. 3.10, the time
evolution of both variables $x_{t}$ and $y_{t}$ is shown using $a=0.89$,
$b=0.6$, $c=0.28$, for different values of $k$. For $k<0.03$ (blue line) the
system is at an equilibrium state. For $k\geq 0.03$ the dynamics of the model
consists of periodic spikes, but for higher values of $k$, smaller is the
amplitude of the spikes and higher the frequency [37].

Figure 3.10: Time evolution of the variables of the Chialvo model. Evolution
of the variables $x_{t}$ (a) and $y_{t}$ (b). For $a=0.89$, $b=0.6$, $c=0.28$,
and different values of $k$. For $k<0.03$ the systems stays at a fixed point
while $k>0.03$ the model depicts a periodic motion.

###  3.3 Synaptic current

The synapse is the region responsible for carrying electrical signals
initiated in the presynaptic neuron that propagate within a postsynaptic
neuron [12]. This connection is often represented by a current term called
synaptic current. The synapses generally consist of three components: the axon
terminals of the presynaptic neuron, a target on the postsynaptic neuron, and
a zone of apposition. The structure of the apposition rules the type of the
synapse, which can be electrical or chemical synapses [12]. While electrical
synapses provide almost instantaneous signal transmission, chemical synapses
can amplify the signal.

The electrical synapses are characterized by a direct interaction among ions
of the two cells, where the presynaptic terminal and the postsynaptic cell are
in very close apposition at regions called gap junctions. The signal produced
by the action potential of the presynaptic neuron reaches the postsynaptic
neuron at a gap junction which connects both neurons [12]. In this sense, the
postsynaptic neuron receives a synaptic current described as

$I_{\mathrm{syn}}=g_{\mathrm{syn}}(V_{\mathrm{pre}}-V_{\mathrm{post}}),$
(3.36)

where $g_{\mathrm{syn}}$ is the conductance of the channels,
$V_{\mathrm{pre}}$ and $V_{\mathrm{post}}$ are the membrane potentials of the
presynaptic and postsynaptic neuron, respectively.

In a chemical synapse, the connection does not occur directly. The action
potential of the presynaptic neuron leaves neurotransmitters, the most common
in cortical neurons are glutamate and $\gamma$-Aminobutyric acid (GABA) [33].
These neurotransmitters diffuse through the synaptic cleft to reach the
postsynaptic neuron, inducing the opening or the closing of channels. For
these cases, the synaptic current is given by

$I_{\mathrm{syn}}=g_{\mathrm{syn}}r(t)(E_{\mathrm{syn}}-V_{\mathrm{post}}),$
(3.37)

where $E_{\mathrm{syn}}$ is the reversal synaptic potential which
characterizes if the synapse is excitatory or inhibitory, and $r(t)$ is a
function that simulates the neurotransmitter kinetics [78].

A way of modeling $r(t)$ is considered a sigmoidal function [33]

$r(t)=\frac{1}{1+\exp(-\lambda_{\mathrm{s}}(V_{\mathrm{pre}}(t)-\beta_{\mathrm{s}}))},$
(3.38)

where $\lambda_{\mathrm{s}}$ and $\beta_{\mathrm{s}}$ are constants related to
chemical synapses, adjusted for the respective neuronal model.

Another way to model $r(t)$ is with a kinetic function, which takes into
account the fraction of open channels that allow the transmission of
neurotransmitters,

$\frac{dr}{dt}=\left(\frac{1}{\tau_{\mathrm{r}}}-\frac{1}{\tau_{\mathrm{d}}}\right)\frac{1-r}{1+\exp[-s_{0}(V_{\mathrm{pre}}(t)-V_{0})]}-\frac{r}{\tau_{\mathrm{d}}},$
(3.39)

in which $\tau_{\mathrm{r}}$ and $\tau_{\mathrm{d}}$ are characteristic times,
$s_{0}=1/$\text{\,}\mathrm{mV}$$ is a unitary constant, and $V_{0}$ is a
reversal potential [78].

##  Chapter 4 Complex networks

The study of the global behavior of a group of connected entities is made
using a complex network. A dynamical system is called complex if its final
state does not depend on the existence of a central controller, i.e., the
cooperative behavior depends only on the interactions among entities [22]. In
this chapter, it is introduced some concepts of the graphs theory, which is
used in the network study. The way that the connections are distributed
characterizes the network connection topology. In this thesis, all the
networks were generated using a Python language package called “NetworkX” [79]
which was developed for exploration and analysis of networks and network
algorithms.

###  4.1 Graphs theory

The graph theory is a section of discrete mathematics that was developed in
part by Leonard Euler to solve the challenge of the Königsberg bridges. The
Pregel river, when cross the Königsberg city is ramified creating two islands,
which were connected to the city by seven bridges. The challenge consisted at
to cross all the bridges without repetition, independently where the journey
was started or finished. Figure 4.1 depicts one adaptation of the original
illustration of the challenge made by Euler [80, 81].

Figure 4.1: Illustration of the Königsberg bridges. The Pregel River
ramifications and form two islands $A$ and $D$, which are connected to the
city ($B$ and $C$) by the $7$ bridges $a$, $b$, $c$, $d$, $e$, $f$ and $g$.
Figure taken from [81], an adaptation of the original work of Euler [80].

The idea used to solve this challenge was to represent each earth region and
the bridges as sites connected by edges. It was found that, since all the
earth regions have an odd number of bridges, there is no solution for this
challenge. The challenge could be solved if any of the bridges were removed.
The way in that Euler obtained his results gives a start to the study of graph
theory [82].

A graph $\mathcal{G}$ is an ordered pair of distinct sets $\mathcal{V}$ and
$\mathcal{E}$, where $\mathcal{V}$ is a non-empty discrete set of elements
called sites or nodes, and $\mathcal{E}$ a subset constituted of ordered pair
of $\mathcal{V}$ elements, called edges or connections. For two given elements
of $\mathcal{V}$, the pair $(v_{i},\,v_{j})\in\mathcal{E}$ is a line which
connects $v_{i}$ and $v_{j}$. The connections among sites can be binary or
weighted: for binary connections, $(v_{i},\,v_{j})=1$ if the sites are
connected or $0$, otherwise; when weighted, some connections are more relevant
than others $(v_{i},\,v_{j})\in\mathbb{R}$. On the other hand, the connections
can be directed and non-directed: when non-directed the connections are
reciprocal which means that, if $v_{i}$ is connected to $v_{j}$, implies that
$v_{j}$ is connected to $v_{i}$; when directed, the connections are not
reciprocal. The notation $\mathcal{G}(N,n)$ represents a graph with $N$ sites
and $n$ connections. The graph $\mathcal{G}(N,n)$ can be represented as a
squared matrix $G$, called connection matrix, where $N$ is the number of rows
and columns, and $n$ is the number of non-zero elements. In this sense, the
element $e_{ij}$ of the matrix $G$ represents the connection between the site
$i$ and $j$.

###  4.2 Complex networks

The graph theory can be used in the study of networks, in particular, neural
networks, where each node of the network is a neuron and their connections are
the synapses. The way of the connections are distributed in the network
determines its topology which plays a role in the global dynamics of the
network [25]. In this sense, besides the size of the network $N$ and the
number of connections $n$, it is possible to classify the topology of the
network using some properties.

The average shortest path length $\mathcal{L}$ is a quantity that refers to
the shortest path between the $i$-th site and $j$-th site. This path relates
to the number of sites that the information needs to pass to reach the target
[22]. If $\mathcal{L}_{ij}$ is the shortest path between $i$ and $j$, the
average shortest path length of the network is

$\mathcal{L}=\dfrac{1}{N(N-1)}\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\mathcal{L}_{ij}.$ (4.1)

The clustering coefficient $\mathcal{C}$ is a quantity of how the sites tend
to cluster. The evaluation is based on the number of sites which groups in
trios [22, 83]. Considering three sites, a closed trio is when all sites are
connected while an open trio one of the connections is missing. The clustering
coefficient is defined by the ratio between the number of closed trios divided
by the total number of trios (open and closed),

$\mathcal{C}=\frac{\mathcal{T}_{\mathrm{\triangle}}}{\mathcal{T}_{\mathrm{\triangle}}+\mathcal{T}_{\land}},$
(4.2)

where $\mathcal{T}_{\mathrm{\triangle}}$ and $\mathcal{T}_{\land}$ are the
number of closed trios and open trios, respectively. If
$\mathcal{T}_{\triangle}\gg\mathcal{T}_{\land}$, $\mathcal{C}\approx 1$, which
means in a clustered network [22].

The simplest network topology is called a regular network, where all the sites
depict the same number of connections, as is the case of a first-neighborhood
network, where the $i$-th site is connected to the $i+1$ and $i-1$ sites
(closest neighbors), or the second-neighborhood network, and global network
where all the possible connections exist. Figure 4.2 depicts an example of
these three regular networks with $N=20$ sites.

Figure 4.2: Examples of regular networks. Using $20$ sites. (a) First-
neighborhood network. (b) Second-neighborhood network. (c) Global network
(complete network).

Another topology often used in the literature is the random topology, where
their connections are randomly distributed [22]. One of the possibilities to
build such network is using the Erdos-Rényi algorithm [84], which consider a
probability of connection $0\leq p_{\mathrm{rand}}\leq 1$, where the number of
connections in the network is given by
$n_{\mathrm{rand}}=p_{\mathrm{rand}}N(N-1)$. Figure 4.3 depicts examples of
random networks with $N=20$ and (a) $p_{\mathrm{rand}}=0.05$ and $n=20$
connections, in this case the small number of connections makes possible to
find isolated sites from the network, (b) $p_{\mathrm{rand}}=0.5$ and $n=190$
connections, and (c) $p_{\mathrm{rand}}=1$ all possible connections.

Figure 4.3: Examples of random graphs. Using $20$ sites. (a) For
$p_{\mathrm{rand}}=0.05$ there are $20$ connections. (b) For
$p_{\mathrm{rand}}=0.5$ there are $190$ connections. (c) For
$p_{\mathrm{rand}}=1$ there are $n(n-1)$ connections.

However, regular and random networks are idealization cases, real networks are
believed to be between these extremes of order and randomness [25]. Watts and
Strogatz have found that with the substitution of local connections for random
connections in a regular network, it is possible to create a network with
shortcuts that presents the following properties: low average shortest path
length and a high clustering coefficient (in comparison with random networks).
These networks were called small-world networks [85]. In addition to being
networks optimized from the point of view of information, propagation [85], it
was found that some real networks have the same topological characteristics of
a small-world network, like social networks such as the world-wide-web [86,
87, 88], electric power grids [87], and even neural networks, as in the case
of the nematode nervous system C. Elegans [85, 89], and in the anatomy of the
human brain and other mammals [90].

One of the ways to generate a small-world network is using the Watts-Strogatz
algorithm [85]. Starting with a regular network with $N$ sites and $n_{0}$
connections, the idea is to replace $n_{\mathrm{ws}}$ connections randomly.
Controlling $n_{\mathrm{ws}}$ with a probability $p_{\mathrm{ws}}$, where
$n_{\mathrm{ws}}=p_{\mathrm{ws}}n$. For lower values of
$p_{\mathrm{ws}}\approx 0$ the network depicts the same features of a regular
one, and higher values $p_{\mathrm{ws}}\approx 1$ of a random one. But for
intermediate values of $p_{\mathrm{ws}}$ the number of shortcuts can put the
network in the small-world regime, presenting a high cluster coefficient and
low average shortest path length. Figure 4.4 presents an example of the Watts-
Strogatz route, at panel (a) a second-neighborhood regular network
($p_{\mathrm{ws}}=0$), (b) a small-world network ($p_{\mathrm{ws}}=0.05$), and
(c) a random network ($p_{\mathrm{ws}}=1$).

Figure 4.4: Example of the Watts-Strogatz route. (a) Regular network
$p_{\mathrm{ws}}=0$, (b) Small-world network $p_{\mathrm{ws}}=0.05$, (c)
Random network $p_{\mathrm{ws}}=1$. The small-world regime is achieved between
the extremes of order and randomness.

However, in the Watts-Strogatz route, since the connections are randomly
replaced, it is possible to disconnect a neuron from the network. To avoid
this, Newmann and Watts proposed an alternative algorithm to achieve the
small-world regime [91], the idea is to add random connections instead of
replacing [22]. For this route, the number of connections depends on a
probability $p_{\mathrm{nw}}$, if $n_{0}$ is the number of connections of the
regular network, the new number of connections is given by
$n_{\mathrm{nw}}=(1+p_{\mathrm{nw}})n_{0}$. This means that for each existing
connection ($n_{0}$), a new random connection can be added with a probability
$p_{\mathrm{nw}}$. The Fig. 4.5 represents an example of the Newman-Watts
route, at the panel (a) a regular network ($p_{\mathrm{nw}}=0$), (b) a small-
world network ($p_{\mathrm{nw}}=0.1$) and (c) a network which is no longer in
the small-world regime.

Figure 4.5: Example of the Newman-Watts route. (a) Regular network
$p_{\mathrm{nw}}=0$, (b) Small-world network $p_{\mathrm{nw}}=0.2$, (c)
Network with $p_{\mathrm{nw}}=1$ in this case the network is no longer at the
small-world regime. Figure 4.6: Small-worldness coefficient. $\mathcal{S}$ as
a function of the probability $p_{\mathrm{ws}}$ (first row), $p_{\mathrm{nw}}$
(second row), for networks with $N=100$ (a) and (d), $N=1000$ (b) and (e),
$N=10000$ (c) and (f). Each color refers to an initial number of connections
on a regular network. The error bars are the standard deviation over $10$
different initialization.

Despite the fact that there is no ideal way to characterize the small-world
regime [86, 92], it is possible to define a quantity to quantify the small-
worldness of a generic network using

$\mathcal{S}(p)=\frac{(\mathcal{C}(p)/\mathcal{C}_{\mathrm{r}})}{(\mathcal{L}(p)/\mathcal{L}_{\mathrm{r}})}$
(4.3)

where $\mathcal{L}$ and $\mathcal{C}$ are the average shortest path length and
the clustering coefficient of a generic network, and
$\mathcal{L}_{\mathrm{r}}$ and $\mathcal{C}_{\mathrm{r}}$ for a random
network. A network is commonly classified as a small-world if $\mathcal{S}\gg
1$, which happens when $\mathcal{C}>\mathcal{C}_{\mathrm{r}}$ and
$\mathcal{L}<\mathcal{L}_{\mathrm{r}}$. Figure 4.6 presents $\mathcal{S}$ as a
function of the probability $p_{\mathrm{ws}}$ (at the first row) and
$p_{\mathrm{nw}}$ (at the second row) for networks with different sizes
$N=100$ (a) and (d), $N=1000$ (b) and (e), $N=10000$ (c) and (f). The Watts-
Strogatz route has a maximum of $\mathcal{S}$ at $p_{\mathrm{ws}}\approx 0.1$,
and for the Newman-Watts the maximum occurs $p_{\mathrm{nw}}\approx 0.3$.
Therefore, the network is in a small-world regime for values close to these
maxima.

##  Chapter 5 Phase synchronization

Synchronization, this basic nonlinear phenomena, detected in the $17$th
century when the scientist Christiaan Huygens reported his observations that
two very weakly coupled pendulum clocks become synchronized in phase [1]
(English translation [2]). In the context of dynamical systems, many different
synchronization states have been studied: complete synchronization, also
called identical synchronization, is the simplest form of synchronization, it
consists of the exact convergence of all the trajectories of the system to a
unique synchronization manifold; frequency synchronization consists in the
frequency locking of the system, where all elements of the system evolve with
the same periodicity; and phase synchronization were beyond the frequencies
the phases are also locked [93]. In this chapter, it is introduced some
concepts of the synchronization of dynamical systems, based on the Kuramoto
oscillator model. In particular, it shows a powerful tool to quantify phase
synchronization, called the Kuramoto order parameter, which can be generalized
to measure phase synchronization of distinct oscillator models [5], and in
special, to neuronal models.

The Kuramoto model consists of an ensemble of $N$ oscillators with a nonlinear
coupling where the phase of the $i$-th oscillator evolves as

$\dot{\theta}_{i}=\omega_{i}+\frac{\varepsilon}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})\,\,\,\,\,\,\,\,\,\,i=1,\cdots,N,$
(5.1)

where $\omega_{i}$ is the natural frequency of the oscillator and
$\varepsilon$ the coupling strength. For $\varepsilon=0$ each oscillator
evolves linearly which his own frequency
$\theta_{i}(t)=\omega_{i}t+\theta_{i}(0)$, where $\theta_{i}(0)$ is the
initial condition of the $i$-th oscillator.

For the particular case of identical oscillators
$\omega_{i}=\omega_{0},\;\forall\;i$, and, without losing generality
$\omega_{0}=0$ (since for $\omega_{0}\neq 0$, it is always possible to choose
a frame that rotates at the same frequency $\omega_{0}$), therefore

$\dot{\theta}_{i}=\frac{\varepsilon}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})\,\,\,\,\,\,\,\,\,\,i=1,\cdots,N.$
(5.2)

With the definition of a vector
${\bm{\theta}}=(\theta_{1},\theta_{2},\cdots,\theta_{N})$, it is shown that
Eq. (5.2) is a gradient equation $\dot{\bm{\theta}}=-\nabla\mathcal{U}$ [94],
where the potential function is given by

$\mathcal{U}=-\frac{\varepsilon}{2N}\sum_{i=1}^{N}\sum_{j=1}^{N}\cos(\theta_{j}-\theta_{i}),$
(5.3)

which means that the trajectories of such systems flow monotonically at the
potential surface and asymptotically approaching the equilibrium point
$\dot{\bm{\theta}}=0$. For $\varepsilon>0$, any symmetric initial condition
around the circle $[0,2\pi]$ will lead the system to a steady-state solution,
called twisted states [94], since the sinusoidal sum among all oscillators
mutually cancel. For non-symmetric initial conditions, the coupling drives the
system to the phase-locking stable state, where
$\theta_{1}=\theta_{2}=\cdots=\theta_{N}$, also called complete phase
synchronized state.

Figure 5.1 depicts the time evolution of the phases of $N=100$ identical
oscillators, at panel (a) a symmetric initial condition given by a uniform
distribution $\theta_{i}(0)=2\pi i/N$, where even for $\varepsilon\neq 0$ the
solution is an unstable equilibrium state. At panels (b – d) a random initial
condition is studied ($\theta_{i}(0)\in[0,2\pi]$), the phases asymptotically
approach to the phase-synchronized state, for
$\varepsilon=10^{-4},\;10^{-3},\;10^{-2}$, respectively. The time required for
the system to reach the synchronized state is called the relaxation time
$\tau$ [27, 95]. $\tau$ is presented in panel (e) using different values of
$N$, it is observed that $\tau$ decreases with the increase of the coupling
parameter in such a way that $\tau(\varepsilon)\propto\varepsilon^{-1}$.
Hence, despite of the particular cases of $\varepsilon=0$ and the twisted
states, a network of identical coupled oscillators always reach the
synchronized state.

Figure 5.1: Time evolution of the phases of identical Kuramoto oscillators and
the relaxation time of the network. (a) For a initial condition which follows
a uniform distribution, i.e., $\theta_{i}(0)=2\pi i/N$, the symmetry of the
system cancels the network dynamics, which is kept in a steady state called
twisted state. (b – d) For a random initial condition
($\theta_{i}(0)\in[0,2\pi]$), the phases evolve asymptotically to the phase-
locking state (synchronized state), where
$\theta_{1}=\theta_{2}=\cdots=\theta_{N}$, with
$\varepsilon=10^{-4},\;10^{-3},\;10^{-2}$, respectively. Greater the coupling
strength $\varepsilon$, fastest is the time to reach the synchronized state
$\tau$ (relaxation time). (e) $\langle\tau\rangle$ as a function of
$\varepsilon$ for different $N$ values (number of oscillators), where a clear
inverse relation is observed $\tau(\varepsilon)\propto\varepsilon^{-1}$ . The
$\langle\cdot\rangle$ denotes the average over 1000 realizations.

On the other hand, for non-identical oscillators, the natural frequency of the
oscillators can be described according to a probability density
$\varkappa(\omega)$ which is defined to be symmetric as a function of a mean
frequency $\Omega$,

$\varkappa(\Omega-\omega)=\varkappa(\Omega+\omega).$ (5.4)

With this exact configuration, Kuramoto had defined a complex mean-field
$\mathcal{Z}$, defined as

$\mathcal{Z}=R(t)e^{i\psi(t)}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}(t)},$
(5.5)

where $\psi$ is the circular average frequency. The absolute value of
$\mathcal{Z}$, called Kuramoto order parameter $R$, quantifies the degree of
phase coherence in the oscillators and can be used as an indicator of phase
synchronization.

Considering $N\rightarrow\infty$ oscillators, the phase transition from the
non-synchronized to the phase-synchronized state occurs in a critical value
$\varepsilon^{*}$ given by

$\varepsilon^{*}=\frac{2}{\pi\varkappa(0)},$ (5.6)

and the modulus of the Kuramoto order parameter is equal to

$R=\begin{dcases}0,\;\;\;&\mathrm{for}\;\;\varepsilon\leq\varepsilon^{*},\\\
\sqrt{\frac{-16}{\pi(\varepsilon^{*})^{3}\varkappa^{\prime\prime}(0)}\left(1-\frac{\varepsilon^{*}}{\varepsilon}\right)},\;&\mathrm{for}\;\;\varepsilon>\varepsilon^{*},\end{dcases}$
(5.7)

which means that, if $\varepsilon$ is lower than a critical value
$\varepsilon^{*}$, $R\approx 0$ (desynchronization) and for
$\varepsilon>\varepsilon^{*}$, $R(\varepsilon)$ approaches asymptotically to
$1$, and is equal to $1$ at the limit of $\varepsilon\rightarrow\infty$. The
results of the Eqs. (5.6) and (5.7) are deduced in the appendix A.

Considering the case where natural frequency of the ensemble of oscillators
that follows a zero-centered Cauchy–Lorentz distribution,

$\varkappa(\omega)=\frac{\zeta}{\pi(\zeta^{2}+\omega^{2})},$ (5.8)

where $\zeta$ is related to the width of the distribution. In this sense,

$\varkappa(0)=\frac{1}{\pi\zeta},\;\;\varkappa^{\prime\prime}(\omega)=-\frac{2\zeta(\zeta^{2}-3\omega^{2})}{\pi(\zeta^{2}+\omega^{2})^{3}}\Rightarrow\varkappa^{\prime\prime}(0)=-\frac{2}{\pi\zeta^{3}},$

therefore the critical coupling for the Cauchy-Lorentz distribution

$\varepsilon^{*}=2\zeta,$ (5.9)

and the order parameter, for $\varepsilon>\varepsilon^{*}$,

$R=\sqrt{\frac{-16}{\pi(8\zeta^{3})(-2/\pi\zeta^{3})}\left(1-\frac{\varepsilon^{*}}{\varepsilon}\right)}=\sqrt{1-\frac{\varepsilon^{*}}{\varepsilon}}.$
(5.10)

Figure 5.2 presents the Kuramoto order parameter $R$ as a function of the
coupling $\varepsilon$ considering a Cauchy-Lorenz distribution (Eq. (5.8))
with $\zeta=0.5$ which gives a critical coupling $\varepsilon^{*}=1$. The
dashed line represents the exact solution of $R$ (Eq. 5.10) and the solid
lines show the numerical results for different values of $N$. The greater the
value of $N$, closer is the numerical result to the exact solution.

Figure 5.2: Phase synchronization of the Kuramoto model. Kuramoto order
parameter $R$ as a function of the coupling $\varepsilon$ using a Cauchy-
Lorenz distribution with $\zeta=0.5$, which gives $\varepsilon^{*}=1$.

Despite the order parameter being created to calculate the phase
synchronization of Kuramoto oscillators, this method can be effectively
applied to distinct oscillators and also neurons regardless of periodicity
presence. To do this, it is necessary to associate a phase to the dynamical
model. The phase can be obtained with a projection of the dynamical system
$\theta(t)=h(\mathbf{x})$ and also with the definition of an appropriate
Poincare section where the orbit crosses once for each $2\pi$ rotation. In
this sense, the phase can be defined

$\theta_{i}(t)=2\pi k_{i}+2\pi\frac{t-t_{{k,i}}}{t_{k+1,i}-t_{{k,i}}},\hskip
14.22636ptt_{k,i}\leq t<t_{k+1,i},$ (5.11)

where $t_{k,i}$ represents the $k$th time where the $i$-th system crosses the
Poincare surface. Therefore, for every $t=t_{k}$, the second term of the
equation vanishes and the phase is equal to $2\pi k_{i}$, for the other
instants of time, the second term interpolates these times until the next
cross occurs in $t=t_{k+1}$. The phase is increased by a factor of $2\pi$
every $t_{k+1,i}-t_{k,i}$ [96]. After this association, it is possible to
quantify the phase synchronization of the system using the Eq. (5.5) for the
order parameter. Different from the Kuramoto model, a dynamical system might
not present an stationary state of $R(t)=R$. To measure phase synchronization
in distinct models, the average order parameter is defined, which consists of
the temporal average of the order parameter module

$\langle
R\rangle=\frac{1}{t_{\mathrm{f}}-t_{\mathrm{i}}}\sum_{t=t_{\mathrm{i}}}^{t_{\mathrm{f}}}R(t),$
(5.12)

where $t_{\mathrm{i}}$ and $t_{\mathrm{f}}$ are the start and end times of the
order parameter computation, respectively. If $\langle R\rangle=1$ represents
a completely phase-synchronized state, in which all elements start to cross
the surface at the same time. On the other hand, $\langle R\rangle=0$ means
that each element in the network has a corresponding pair that is completely
out-of-phase. This can correspond to a completely incoherent state (completely
unsynchronized) or a state with clusters of in-phase neurons that are anti-
phase between themselves. If the $N$ phases were to be randomly distributed,
the result would be $\langle R\rangle\sim\sqrt{1/N}$ [26].

##  Chapter 6 The role of the individual dynamics in the synchronization
processes of neural networks

Most of the neuronal behaviors produced by the brain are directly related to
the collective patterns of activation of groups of neurons, generating
(partially) synchronized dynamics of these groups [12]. However, the brain
activity can be disturbed by pathological states caused by some neuropathies,
as in the case of epilepsy, where the increase in neuron synchronization
generates seizure episodes [17], and in Parkinson’s disease, where an
excessive synchronization in the basal ganglia is evidenced [18, 19, 97]. In
this sense, drugs and other treatments are used to normalize brain activity,
allowing the neurons to depolarize and repolarize in a healthy way [98, 99,
100]. In this scenario, the use of simulated neural networks have proven to be
very useful to understand the synchronization mechanism to optimize or
suppress the synchronization of neurons [51, 101, 102].

This chapter is devoted to studying the influence of the dynamics of the
neurons in their synchronization processes of a neural network with HH$\ell$
neurons. It is shown that exists a strong correlation between the individual
dynamics of the neuron with the type of synchronization of the network. For
neurons that individually exhibit chaotic behavior, phase synchronization
occurs analogously to the transition of chaotic oscillators known in the
literature [48]. However, when non-chaotic neurons are coupled, the network
may present phase synchronization to weak coupling regimes (a coupling
parameter close to zero) due to the influence of the regularity of the
individual dynamics, and the increase of the coupling parameter can decrease
the phase synchronization of the network [49, 50, 51], characterizing a non-
monotonic evolution of synchronization as a function of coupling.

This phenomenon has already been explored in other works using different
approaches: network of networks [49, 103], detection of nonstationarity [50],
suppression of phase synchronization [51, 104, 105, 106] and by temperature
changes [107, 108]. Here, the dynamics of the neurons are changed by
variations in the ion conductances of the model, and with the application of
an external pulsed current to the neurons, both situations show similar
results. To confirm the generality of this behavior, a verification test is
successfully performed with another neuron model in the same conditions. In
conclusion, the occurrence of this non-monotonic phase transition occurs due
to an interplay between the individual-regular behavior of the neurons and the
influence of the synaptic current. Hence, it is possible to change the phase
synchronization in the weak coupling regime just with small variations to
achieve a chaotic transition in the neuron dynamics. Most of the results of
this chapter are published in the article “The role of individual neuron ion
conductances in the synchronization processes of neuron networks” Neural
Networks 137 (2021) 97–105 [109].

###  6.1 Network properties

It is considered a network with $N=1000$ HH$\ell$ neurons coupled in a small-
world network, generated with the Watts-Strogatz route with $n_{0}=6000$
connections and $p_{\mathrm{ws}}=0.1$. The membrane potential of the $i$-th
neuron is described 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}},$
(6.1)

where the ionic currents are given by the Eqs. (3.18-3.22). For this specific
case, the neurons are coupled with a chemical synapse

$I_{i,\mathrm{syn}}=\frac{\varepsilon}{\bar{n}}(E_{\mathrm{syn}}-V_{i})\sum\limits_{j=1}^{N}e_{i,j}r_{j}(t),$
(6.2)

with $\varepsilon$ being the coupling strength parameter, $\bar{n}$ is the
average of the number of connections in the network, $E_{\mathrm{syn}}$ is the
synaptic reversal potential. It is used $E_{\mathrm{syn}}=20$
$\text{\,}\mathrm{mV}$, since $E_{\mathrm{syn}}>V_{i}$, the synaptic term
$I_{\mathrm{syn}}$ is always positive featuring an excitatory synapse.
$e_{i,j}$ gives the element of the connection matrix, so that if $i$ and $j$
are connected $e_{i,j}=1$, otherwise $e_{i,j}=0$. The variable $r_{i}$
represents the fraction of bound receptors in the synapse where the kinetics
model is described by [78]

$\frac{dr_{i}}{dt}=\left(\frac{1}{\tau_{r}}-\frac{1}{\tau_{d}}\right)\frac{1-r_{i}}{1+\exp[-s_{0}(V_{i}-V_{0})]}-\frac{r_{i}}{\tau_{d}},$
(6.3)

where $s_{0}$ is an unitary constant $1\text{\,}\mathrm{m}\mathrm{V}^{-1}$,
$V_{0}=-20$ $\text{\,}\mathrm{mV}$, $\tau_{r}=0.5$ $\text{\,}\mathrm{ms}$, and
$\tau_{d}=8$ $\text{\,}\mathrm{ms}$ are constants.

Figure 6.1 presents the dynamics of the membrane potential (a) and the
variable $a_{\mathrm{sd}}$ (b) for one neuron simulated using the values of
Table 3.2. The temporal dynamic is composed of a set of (four) spikes followed
by a resting time, characterizing bursting dynamics. The bursting dynamic is
maintained for all parameters used in this Chapter. The $a_{\mathrm{sr}}$
variable, described in Eq. (3.27), is used as an auxiliary variable to compute
the bursting time since the local minimum of $a_{\mathrm{sr}}$ coincides with
the beginning of a burst. After compute all the bursts of all the neurons,
using Eq. (5.11) it is possible to evaluate $\theta_{i}$ and the phase
synchronization $\langle R\rangle$ with Eq. (5.12), here, considering
$t_{\mathrm{i}}=150$ $\text{\,}\mathrm{s}$ (avoiding transient effects) and
$t_{\mathrm{f}}=250$ $\text{\,}\mathrm{s}$ as a final time of computation. The
initial conditions for each neuron are randomly distributed to avoid any
initial synchronization bias. Each synchronization result ($\langle R\rangle$)
is a mean value over $10$ distinct initial conditions and $10$ different
networks preserving the same connection properties.

Figure 6.1: Phase association of the HH$\ell$ model. Time evolution of the:
(a) membrane potential of the neuron, (b) $a_{\mathrm{sr}}$ variable, using
the parameters of Table 3.2. The local minimums of $a_{\mathrm{sr}}$ (dashed
line) can be used as an auxiliary to compute the beginning of a burst.

###  6.2 Results

The effects of variations in neuron dynamics due to changes in the ion
conductance can be illustrated by a bifurcation diagram of the Inter-Burst-
Interval (time between two successive bursts, IBI). Figure 6.2 (a) depicts the
IBI (x-axis) for one single isolated neuron ($\varepsilon=0$) as a function of
$g_{\mathrm{r}}$ (y-axis). It is observed for $1.90<g_{\mathrm{r}}\lesssim
2.03$ that the neuron dynamics has a periodic behavior. When $g_{\mathrm{r}}$
is increased above $1.95$ a cascade of period-doubling bifurcations occurs
producing infinite orbits and finally, for $g_{\mathrm{r}}\gtrsim 2.03$, the
large variability of IBI illustrates that the neuron dynamics is chaotic [60].
At panel (b), it is computed the Kuramoto order parameter $\langle R\rangle$
in color codes from blue tones ($\langle R\rangle=0$, non-synchronized) to red
tones ($\langle R\rangle=1$, complete phase synchronized), as function of
$\varepsilon$ and $g_{\mathrm{r}}$. It is noted that for values of
$g_{\mathrm{r}}>2.03$, the transition for the non-synchronized state to the
synchronized state is very similar to the transition seen in chaotic
oscillators [48, 101], which means that for $\varepsilon\lesssim 0.01$,
$R\approx 0$, and for $\varepsilon>0.01$ the network gains phase
synchronization, $R\approx 1$. On the other hand, for lower values of
$g_{\mathrm{r}}$, the network starts to synchronize for very low values of
$\varepsilon$, this behavior is associated with the fact that the neurons,
which individually present regular dynamics, are susceptible to synchronize
since the synaptic current is not strong enough to disrupt the periodic
features of the neurons. Hence, due to the frequency-locked phenomenon, the
synaptic current acts to aligning the neurons in phase generating the partial
phase synchronization behavior. As the coupling strength increases, the
synaptic current gains relevance, inducing chaoticity in the neurons,
resulting in a sharp decrease in the phase synchronization level,
characterizing the non-monotonic evolution of the synchronization.

Figure 6.2: Neural dynamics of HH$\ell$ model as a function of
$\bm{g_{\mathrm{r}}}$. (a) Bifurcation diagram of the inter-burst-interval
(IBI) for the individual neuron as a function of the maximum conductance
$g_{\mathrm{r}}$. (b) The Kuramoto order parameter of the network $\langle
R\rangle$ as a function of the coupling parameter ($\varepsilon$) and the
conductance $g_{\mathrm{r}}$.

Figure 6.3 shows how this transition occurs for $g_{\mathrm{r}}=1.92$ (regular
individual dynamics), where at the panel (a), the IBI for one arbitrary neuron
of the network is plotted as a function of $\varepsilon$ (left scale) with
$\langle R\rangle$ which is presented in red (right scale). In this case, a
weak coupling initially leads the network to a partial phase-synchronized
state $\langle R\rangle\approx 0.8$, but as the coupling is increased the
phase synchronization is decreased and for even stronger coupling a chaotic-
phase-synchronized behavior is achieved. It is important to mention that a
weak chaotic synchronized dynamic is still preserved for a moderate dispersion
of the IBI and some intervals of the coupling parameter. Panels (b – e)
present raster plots of the network where each dot corresponds to the
beginning of a burst. In panel (b), $\varepsilon=0.002$, an incoherent
behavior is observed. For panel (c) $\varepsilon=0.004$, some vertical
structures are noticed in the raster plot, this characterizes a partial phase
synchronization behavior, where bursts of neurons occur at close time
instants. In panel (d), $\varepsilon=0.008$ it shows the decrease of the phase
synchronization, and (e) $\varepsilon=0.020$ the chaotic-phase-synchronized
behavior. On the other hand, Fig. 6.4 depicts the same configuration, but with
$g_{\mathrm{r}}=2.05$, which confirms that once the individual dynamics is
already chaotic only the final transition from desynchronized state to the
chaotic-phase-synchronized state is observed. This can be seen at the
monotonic evolution of $R$ (panel (a)), and at the raster plots (panels (b –
e)).

Figure 6.3: Synchronization transition of individual regular neurons. (a)
Bifurcation diagram of the IBI (left scale) for an $i$-th coupled neuron in
the network, as a function of the coupling parameter ($\varepsilon$), the red
curve, represents the Kuramoto order parameter of the network (right scale)
using $g_{\mathrm{r}}=1.92$. (b – e) Raster plots of the network where each
dot corresponds to the beginning of a burst ($t_{k,i}$). The colored arrows
indicate the coupling values when the raster plots are obtained, which match
with the colors of the dots. Figure 6.4: Synchronization transition of
individual chaotic neurons. (a) Bifurcation diagram of the IBI (left scale)
for an $i$-th coupled neuron in the network, as a function of the coupling
parameter ($\varepsilon$), and the red curve, represents the Kuramoto order
parameter of the network (right scale) using $g_{\mathrm{r}}=2.05$. (b – e)
Raster plots of the network where each dot corresponds to the beginning of a
burst ($t_{k,i}$). The colored arrows in panel (a) indicate the coupling
values when the raster plots are obtained, which match with the colors of the
dots.

Another option to study the individual dynamics of the neuron, and correlates
to the type of transition of the phase synchronization, is comparing the
largest Lyapunov exponent ($\Lambda$) of the individual neuron with the
synchronization features of the network. The panel (a) of Fig. 6.5, presents
the largest Lyapunov exponent ($\Lambda$) presented in Chapter 2 of one
isolated neuron for a parameter space $g_{\mathrm{r}}\times g_{\mathrm{d}}$
using a grid of 800 $\times$ 800 parameter-pairs. The color bar codifies the
value of $\Lambda$ going from black $\Lambda\leq 0$ to hot colors $\Lambda>0$.
It is selected two sets of 4 points each in the parameter space where the
colored triangles represent points where $\Lambda\leq 0$,
$g_{\mathrm{r}}\times g_{\mathrm{d}}=\\{(1.90,1.30)$, $(2.00,1.50)$,
$(2.00,1.68)$, $(2.24,1.60)\\}$ and the colored dots represent points where
$\Lambda>0$, $g_{\mathrm{r}}\times g_{\mathrm{d}}=\\{(1.80,1.40)$,
$(1.90,1.60)$, $(2.10,1.50)$, $(2.20,1.70)\\}$. In panels (b) and (c) it is
measured the phase synchronization transition for each selected point,
plotting $\langle R\rangle$ as a function of $\varepsilon$ for each triangle
(b) and dot (c) point (with the respective colors). It is noted that at the
panel (b), the $4$ curves depict a non-monotonic transition, which means that
for lower coupling values there is a local maximum of phase synchronization,
for the magenta line it happens at $\varepsilon=0.006$, brown
$\varepsilon=0.004$, purple and yellow $\varepsilon=0.002$. On the other hand,
at the panel (c), all $4$ curves the intrinsic chaoticity of the neurons is
sufficient to avoid phase synchronization for the weak coupling regime, and
the network presents phase synchronization only for couplings greater than a
critical coupling $\varepsilon^{*}$ which its value depends on the parameter
space. It is noted that a synchronization transition can be easily achieved
for weaker values of $\varepsilon$ with the change of more than one parameter
since it is necessary to achieve a bifurcation point in the parameter space to
change the individual dynamics. In these cases, small changes in one specific
ion-conductance may act as a catalyst to the second conductance change [109].

Figure 6.5: Correlation between the individual dynamics and the phase
synchronization varying the parameter space $\bm{g_{\mathrm{r}}\times
g_{\mathrm{d}}}$. (a) Largest Lyapunov Exponent ($\Lambda$) obtained for one
isolated neuron in the parameter space $g_{\mathrm{r}}\times g_{\mathrm{d}}$.
Two sets of $4$ points are selected inside chaotic $\Lambda>0$ (colored dots)
and regular $\Lambda\approx 0$ (colored triangles) parameter space regions.
(b, c) $\langle R\rangle$ as a function of $\varepsilon$ for the two sets of
$4$ points of the parameter space, (b) regular dynamics (triangles in (a)) and
(c) chaotic dynamics (dots in (a)).

To investigate further the relation between the individual dynamics with the
phase synchronization, assuming that other conductances also can play a
similar role in the network, it is fixed $g_{\mathrm{r}}$ and $g_{\mathrm{d}}$
in $2.00$ and $1.50$ (original values of Table 3.2 [76]), but with the
variations of the slow currents of the model $g_{\mathrm{sr}}$ and
$g_{\mathrm{sd}}$. The results are depicted in Fig. 6.6 (the same methodology
of Fig. 6.5 is used). Panel (a) $\Lambda$ is plotted as a function of
$g_{\mathrm{sr}}$ and $g_{\mathrm{sd}}$, and, once again panels (b) and (c)
depict the $\langle R\rangle$ as a function of $\varepsilon$ for two sets of
$4$ points selected in the parameters space where (b) corresponds to the
colored triangles (regular individual dynamics) $g_{\mathrm{sr}}\times
g_{\mathrm{sd}}=\\{(0.395,0.210)$, $(0.395,0.235)$, $(0.400,0.250)$,
$(0.410,0.275)\\}$ and (c) the colored dots (chaotic individual dynamics)
$g_{\mathrm{sr}}\times g_{\mathrm{sd}}=\\{(0.385,0.200)$, $(0.395,0.220)$,
$(0.396,0.248)$, $(0.400,0.260)\\}$. Even for very similar parameters, as the
purple triangle ($0.400,0.250$) and the the blue dot ($0.396,0.248$) the phase
synchronization evolves in each case in completely different ways when the
lower coupling regime is considered.

Figure 6.6: Correlation between the individual dynamics and the phase
synchronization varying the parameter space $\bm{g_{\mathrm{sr}}\times
g_{\mathrm{sd}}}$. (a) $\Lambda$ obtained for one isolated neuron in the
parameter space $g_{\mathrm{sr}}\times g_{\mathrm{sd}}$. Two sets of $4$
points are selected inside chaotic $\Lambda>0$ (colored dots) and regular
$\Lambda\approx 0$ (colored triangles) parameter space regions. (b, c)
$\langle R\rangle$ as a function of $\varepsilon$ for the two sets of $4$
points of the parameter space, (b) regular dynamics (triangles in (a)) and (c)
chaotic dynamics (dots in (a)).

####  6.2.1 External perturbation by electrical pulses

In the last section, it was shown a correlation between the individual
dynamics of the neuron with the transition of the phase synchronization of the
network. Precisely, the synchronization transition of periodic neurons is
different from chaotic ones. In this sense, instead of change the parameters
of the neuron to achieve a bifurcation transition, it is possible to induce
chaoticity with the application of external perturbations. To do this, it is
made the addition of an external current $\xi(t)$ in Eq. (6.1), which
describes a sequence of electrical pulses

$\xi(t)=\frac{\xi_{0}}{2}+\sum_{m=1,3,5,\cdots}^{\infty}\frac{2\xi_{0}}{m\pi}\sin(2m\pi\nu
t),$ (6.4)

where $\xi_{0}$ is the amplitude of the pulse measured in
$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$,
$\nu$ is the frequency measured in $\text{\,}\mathrm{Hz}$, and $m$ is an odd
integer, and for the simulations it is used the $m<1000$ terms. The choice for
a pulsed current is inspired in the Deep Brain Stimulation procedure, where
electrodes are placed deep in the brain and are connected to a stimulator
device. When turned on, the stimulator emits electrical pulses, which in high
frequencies ($>100$ $\text{\,}\mathrm{Hz}$ [110] in comparison with the
frequency of the inter-bursts $\sim 1$ $\text{\,}\mathrm{Hz}$) suppresses
tremor symptoms associated with essential tremor or Parkinson disease [100].

The results are depicted in Fig. 6.7, where panel (a) presents the IBI of the
individual neuron as a function of the amplitude of the applied pulse
$\xi_{0}$ with a fixed frequency of $\nu=100$ $\text{\,}\mathrm{Hz}$, for
fixed parameters as in Table 3.2. It is seen that the regular dynamic
exhibited by the neuron is kept until $\xi_{0}\approx 0.02$, where a great
variability of IBI appears in the dynamics. At panel (b), it is measured the
phase synchronization ($\langle R\rangle$) as a function of the coupling
parameter $\varepsilon$ and $\xi_{0}$. Similarly, the individual dynamics of
the neuron rule the synchronization transition, since for $\xi_{0}<0.02$ there
is a local maximum of synchronization for $\varepsilon\approx 0.006$, which is
associated with the non-monotonic evolution of $\langle R\rangle$. On the
other hand, for $\xi_{0}>0.02$, the individual neuron presents chaotic
dynamics, and the phase synchronization is reached only for coupling values
greater than a critical value $\varepsilon^{*}\approx 0.01$, corroborating
with the results of the previous section.

Figure 6.7: Chaos induced by the application of external pulses. (a)
Bifurcation diagram of IBI for an individual neuron as a function of the
amplitude of the current $\xi_{0}$. (b) The Kuramoto order parameter of the
network $\langle R\rangle$ as a function of $\varepsilon$ and $\xi_{0}$.

####  6.2.2 Verification test with the Hindmarsh-Rolse model

To generalize the non-monotonic of the synchronization phenomenon, the same
situation is studied as a verification test using a second network in the same
conditions, i.e. same size, connections, and topology, but simulating the
Hindmarsh-Rose neuron model described by Eqs. (3.31 – 3.33), which now reads

$\displaystyle\frac{dx_{i}}{dt}$ $\displaystyle=$ $\displaystyle
y_{i}-ax_{i}^{3}+bx_{i}^{2}-z_{i}+\mathcal{I}+\mathcal{I}_{i,\mathrm{syn}},$
(6.5) $\displaystyle\frac{dy_{i}}{dt}$ $\displaystyle=$ $\displaystyle
c-dx_{i}^{2}-y_{i},$ (6.6) $\displaystyle\frac{dz_{i}}{dt}$ $\displaystyle=$
$\displaystyle r[s(x_{i}-x_{\mathrm{r}})-z_{i}],$ (6.7)

where $a$, $b$, $c$, $e$, $r$, $s$, and $x_{\mathrm{r}}$ are constants whose
values are depicted in Table 3.3, $\mathcal{I}$ is a free parameter, and
$\mathcal{I}_{i,\mathrm{syn}}$ is the dimensionless “synaptic” parameter
defined as

$\mathcal{I}_{i,\mathrm{syn}}=\frac{\varepsilon}{\bar{n}}(x_{\mathrm{s}}-x_{i})\sum_{j=1}^{N}e_{i,j}\frac{1}{1+\exp[-\lambda_{\mathrm{s}}(x_{j}-\beta_{\mathrm{s}})]}.$
(6.8)

where $\varepsilon$ is the coupling parameter, $\bar{n}$ is the average of
connections, and $e_{i,j}$ is the element of the connection matrix. $x_{s}$,
$\lambda_{s}$ and $\beta_{s}$ are constants related to the chemical synapses,
as described in Chapter 4 where it is used $x_{s}=2$, $\lambda_{s}=30$, and
$\beta_{s}=0$. To evaluate the bursting times $t_{k,i}$ to compute the phase
of the $i$-th neuron ($\theta_{i}$), it is used the local minimum values of
$z_{i}(t)$ which correlates with the beginning of a burst, as shown in Fig.
3.9 (analogously to the $a_{\mathrm{sr}}$ variable of the HH$\ell$ model). The
dynamic of the individual neuron is depicted in panel (a) of Fig. 6.8 where is
computed the IBI as a function of $\mathcal{I}$. The neuron presents a
periodic behavior for $\mathcal{I}<2.92$, and vast variability of IBI (chaotic
dynamics) $\mathcal{>}2.92$. At panel (b) it is evaluated the $\langle
R\rangle$ in color codes as a function of $\varepsilon$ and $\mathcal{I}$. In
a similar way to the results of the previous sections, the transition from
non-synchronized state to the synchronized state varies with the individual
neuron dynamics: non-monotonic transition for $\mathcal{I}<2.92$ and monotonic
for $\mathcal{I}>2.92$.

Figure 6.8: Verification test with the Hindmarsh-Rose model. (a) Bifurcation
diagram of IBI for an individual Hindrmarsh-Rose neuron as a function of the
parameter $\mathcal{I}$. The Kuramoto order parameter of the network $\langle
R\rangle$ as a function of $\varepsilon$ and $\mathcal{I}$.

####  6.2.3 Topological effects

In this project, it is opted to couple the neurons in a complex network with a
small-world topology. This choice is justified since, as discussed in Chapter
4, the small-world features are detected in real neural networks. However, one
can believe that the non-monotonic evolution of the network can arise due to
the small-world topology. To disprove this hypothesis, in Fig. 6.9, it is
studied $\langle R\rangle$ as a function of the coupling $\varepsilon$ for
different networks’ topologies generated with approximately the same number of
connections ($n_{\mathrm{sw}}=n_{\mathrm{random}}=6000$ and
$n_{\mathrm{sf}}=5982$). In panel (a) $g_{\mathrm{r}}=1.92$ (individually
periodic) is used and the phase synchronization transitions are non-monotonic
and panel (b) $g_{\mathrm{r}}=2.05$ (individually chaotic) is used and the
transitions are monotonic. These results confirm that the same effects of the
individual dynamics can be detected in different complex networks’ topologies.

Figure 6.9: Synchronization transition for different networks’ topologies.
$\langle R\rangle$ as a function of $\varepsilon$ (a) for
$g_{\mathrm{r}}=1.92$, and (b) $g_{\mathrm{r}}=2.05$ comparing different
topologies the same number of connections. The error bars are the standard
deviation over $10$ different realizations.

###  6.3 Discussions

Throughout this chapter, the phenomenon of non-monotonic synchronization was