# Topological-numerical analysis of a two-dimensional discrete neuron model

Paweł Pilarczyk  P. Pilarczyk is with Faculty of Applied Physics and
Mathematics & Digital Technologies Center, Gdańsk University of Technology,
ul. Narutowicza 11/12, 80-233 Gdańsk.
Orcid: 0000-0003-0597-697X<EMAIL_ADDRESS>,  Justyna Signerska-
Rynkowska∗ J. Signerska-Rynkowska is with Dioscuri Centre in Topological Data
Analysis, Institute of Mathematics of the Polish Academy of Sciences, ul.
Śniadeckich 8, 00-656 Warszawa;
and with Faculty of Applied Physics and Mathematics & BioTechMed Center,
Gdańsk University of Technology, ul. Narutowicza 11/12, 80-233 Gdańsk,
Orcid: 0000-0002-9704-0425<EMAIL_ADDRESS>and  Grzegorz Graff  G.
Graff is with Faculty of Applied Physics and Mathematics & BioTechMed Center,
Gdańsk University of Technology, ul. Narutowicza 11/12, 80-233 Gdańsk.
Orcid: 0000-0001-5670-5729<EMAIL_ADDRESS>
(Date: September 7, 2022)

###### Abstract.

We conduct computer-assisted analysis of the two-dimensional model of a neuron
introduced by Chialvo in 1995 (Chaos, Solitons & Fractals 5, 461–479). We
apply the method for rigorous analysis of global dynamics based on a set-
oriented topological approach, introduced by Arai et al. in 2009 (SIAM J.
Appl. Dyn. Syst. 8, 757–789) and improved and expanded afterwards.
Additionally, we introduce a new algorithm to analyze the return times inside
a chain recurrent set. Based on this analysis, together with the information
on the size of the chain recurrent set, we develop a new method that allows
one to determine subsets of parameters for which chaotic dynamics may appear.
This approach can be applied to a variety of dynamical systems, and we discuss
some of its practical aspects. The data and the software described in the
paper are available at http://www.pawelpilarczyk.com/neuron/.

###### Key words and phrases:

nonlinear neurodynamics, spiking-bursting oscillations, Conley index, Morse
decomposition, rigorous numerics, excitable systems, recurrence, computer-
assisted proof

###### 2020 Mathematics Subject Classification:

Primary: 37B35. Secondary: 37B30, 37M99, 37N25, 92-08.

∗Corresponding author.

Research supported by the National Science Centre, Poland, within the grants:
Sheng 1 2018/30/Q/ST1/00228 (for G.G.), SONATA 2019/35/D/ST1/02253 (for
J.S-R.), and OPUS 2021/41/B/ST1/00405 (for P.P.). J. S-R. also acknowledges
the support of Dioscuri program initiated by the Max Planck Society, jointly
managed with the National Science Centre (Poland), and mutually funded by the
Polish Ministry of Science and Higher Education and the German Federal
Ministry of Education and Research. Computations were carried out at the
Centre of Informatics Tricity Academic Supercomputer & Network.

##  1\. Introduction

With the increasing capabilities of contemporary computers, it is possible to
apply more and more computationally demanding methods to the analysis of
dynamical systems. Such methods may provide comprehensive overview of the
dynamics on the one hand, and thorough insight into specific features of the
system on the other hand.

In this paper, we discuss an application of a computationally advanced method
for the analysis of global dynamics of a system with many parameters. The
method was originally introduced in [1], and now we considerably enhance it by
introducing Finite Resolution Recurrence analysis, as explained below and in
Section 3.3. We apply this method to the two-dimensional discrete-time semi-
dynamical system introduced in [5] for modeling a single neuron. We describe
this model in Section 2.

Our approach uses rigorous numerical methods and a topological approach based
on the Conley index and Morse decompositions, and provides mathematically
validated results concerning the qualitative dynamics of the system. The main
idea is to cover the phase space (a subset of ${\mathbb{R}}^{n}$) by means of
a rectangular grid, and to use interval arithmetic to compute an outer
estimate of the map on the grid elements ($n$-dimensional rectangular boxes).
This construction gives rise to a directed graph, and fast graph algorithms
allow one to enclose all the recurrent dynamics in bounded subsets, further
called Morse sets, built of the grid elements, so that the dynamics outside
the collection of these subsets is gradient-like. The entire range of
parameters under consideration is split into classes in such a way that
parameters within one class yield equivalent dynamics. We outline this method
in Section 3.1.

Since existence of chaotic dynamics implies recurrence in large areas of the
phase space (existence of large “strange attractors”), construction of an
outer estimate for the chain recurrent set results in this case in just one
large isolating neighborhood, and therefore the approach based on constructing
a Morse decomposition provides very little information on the actual dynamics.
In order to address this problem, we introduce new algorithmic methods for the
analysis of the directed graph that represents the map in order to get insight
into the dynamics _inside_ this kind of a large Morse set. In particular, we
introduce the notion of Finite Resolution Recurrence (FRR for short) and we
analyze its variation within the set. We introduce this approach in Section
3.2.

After having introduced the computational method and having explained it from
the theoretical point of view, we show its practical application to an
effective computer-assisted comprehensive analysis of the various types of
dynamics that appear in the Chialvo model. We explain how to set up the
various parameters of the computations, and we discuss the results in Section
4.1. Finally, we develop certain indicators of chaotic dynamics that are based
on the FRR analysis and apply them in Sections 4.4, 4.5, and 5.

One could summarize the main ideas of our method for comprehensive analysis of
a dynamical system at finite resolution in the following way (see Figure 1):

  1. (a)

We split the dynamics into chain recurrent sets (Morse sets) between which the
dynamics is simple (gradient-like).

  2. (b)

We obtain some information about the Morse sets by looking at them from the
_outside_ (computing their Conley indices); in particular, a nontrivial Conley
index implies the existence of a non-empty isolated invariant subset of the
Morse set.

  3. (c)

Computation of (a) and (b) involves rigorous numerical methods (based on
interval arithmetic) and provides mathematically reliable results (computer-
assisted proof).

  4. (d)

We study the dynamics _inside_ the Morse sets using a new method: Finite
Resolution Recurrence analysis.

  5. (e)

We conduct (a)–(d) for a grid of parameters of the system varying in some
bounded ranges, and we gather the information in order to classify the types
of dynamics encountered, detect bifurcations, and search for chaotic dynamics.

In particular, by applying our approach to the classical Chialvo model (1) we
are able to obtain precise and comprehensive description of the dynamics for a
large range of parameters considered in [5].

As it will be made clear in the sequel, the set-oriented topological method
that we apply has several advantages over purely analytical methods, and over
plain numerical simulations as well. One could argue that analytical methods
typically focus on finding equilibria of the system and determining their
stability. On the other hand, numerical simulations are usually limited to
iterating individual trajectories and thus their ability is limited to finding
stable invariant sets only, not to mention their vulnerability to round-off or
approximation errors. In contrast to this, our approach detects all kinds of
recurrent dynamics (also unstable) in a given region of the phase space and
provides mathematically reliable results.

It is also worth pointing out that introducing advanced methods for
investigating invariant sets and their structure is especially desirable in
discrete systems already in dimension $2$, such as the Chialvo model. The
reason is that such systems are much more demanding than their ODE
“counterparts.” For example, in dimension $2$, the information on the
invariant (limit) sets can be concluded from the Poincaré-Bendixson Theorem,
the shape of stable and unstable manifolds of saddle fixed points, and other
elementary considerations. Indeed, as Chialvo noticed in [5], “Trajectories
associated with iterated maps are sets of discrete points, and not continuous
curves as in ODE. In two-dimensional ODE, orbits or stable and unstable
manifolds partition the phase space in distinct compact subsets with their
specific attractors. Structure of the stable sets might be more complicated
for 2D iterated maps.”

On the other hand, in the analysis of one-dimensional discrete models, one can
benefit from the theory of circle maps or the theory of interval maps; both
have undergone rapid development in recent decades. In particular, the theory
of $S$-unimodal maps can be successfully applied to obtain rigorous results
for the one-dimensional Chialvo model [21]. Unfortunately, in the discrete
setting, increasing the dimension from one to two makes the analysis
considerably harder, as no such powerful analytical tools exist for discrete
models in two dimensions. Hence there is a strong demand for reliable
computational techniques in studying (discrete) higher dimensional systems.
This demand has been one of our main motivations for conducting the research
described in this paper.

The structure of the paper is as follows. In Section 2, we describe a
$2$-dimensional discrete-time dynamical system introduced by Chialvo for
modeling an individual neuron. In Section 3, we explain the set-oriented
topological method for comprehensive analysis of a dynamical system, and we
introduce the Finite Resolution Recurrence (FRR) and its Variation (FRRV) as
new mathematical tools for deep analysis of recurrent dynamics at limited
resolution. In Section 4, we thoroughly discuss an application of the methods
introduced in Section 3 to the model introduced in Section 2. In particular,
we explain the various kinds of global dynamics that we encountered in the
phase space across the analyzed ranges of parameters, and we conduct the
Finite Resolution Recurrence analysis of the dynamics. Finally, in Section 5,
we use FRR as a tool in classification of dynamics, and we demonstrate its
usefulness as an indicator of the existence of chaotic dynamics.

##  2\. Model

The following model of a single neuron was proposed by D. Chialvo in [5]:

(1)
$\left\\{\begin{array}[]{rcl}x_{\\_}{n+1}&=&x_{\\_}n^{2}\exp(y_{\\_}n-x_{\\_}n)+k,\\\
y_{\\_}{n+1}&=&ay_{\\_}n-bx_{\\_}n+c.\end{array}\right.$

In this model, $x$ stands for the membrane (voltage) potential of a neuron. It
is the most important dynamical variable in all neuron models. However, in
order to model neuron kinetics in a more realistic way than by means of a
single variable, at least one other dynamical variable must be included in the
model. Therefore, the system (1) contains also $y$ that acts as a recovery-
like variable. There are four real parameters in this model: $k$ which can be
interpreted as an additive perturbation or a current input the neuron is
receiving, $a>0$ which is the time constant of neuron’s recovery, the
activation (voltage) dependence of the recovery process $b>0$, and the offset
$c$.

This discrete model, in which $x_{\\_}n$ and $y_{\\_}n$ are values of the
voltage and the recovery variable at the consecutive time units $n$, belongs
to the class of so-called map-based models. Such models have received a lot of
attention recently, and include the famous Rulkov models [30, 31, 32]) and
many others (see also review articles [9], [13] and references therein). For
completeness, we also mention the fact that neurons can be modeled by
continuous dynamical systems, i.e., ordinary differential equations, dating
back to the pioneering work of Hodgkin and Huxley [12], or by hybrid systems
(see e.g. [3, 14, 28, 29, 33, 36]). Models taking into account the propagation
of the voltage through synapses or models of neural networks often incorporate
PDEs and stochastic processes (see e.g. [4, 35] and references therein).
Although some might consider map-based models too simplified from the
biological point of view, their biological relevance in fact can be sometimes
satisfactory. Their important advantage is that they are often computationally
plausible as components of larger systems. Due to the tremendous complexity of
real neuronal systems, map-based models appeared as models that are simple
enough to be dealt with, yet able to capture the most relevant properties of
the cell. Map-based models can sometimes be seen as discretizations of ODE-
based models.

It seems that the Chialvo model (1) is not a direct discretization of any of
the popular ODE models. However, as noticed by Chialvo himself (see [5] and
references therein), the shape of nullclines reminds that of some two-
dimensional ODE excitable systems. Moreover, the equation (1) fulfills the
most common general form of map-based models of a neuron (compare with [13]):

(2) $\left\\{\begin{array}[]{rcl}x(t+1)&=&F\big{(}x(t),k\pm y(t)\big{)},\\\
y(t+1)&=&y(t)\mp\varepsilon\big{(}x(t)-qy(t)-\sigma\big{)}.\end{array}\right.$

Before we proceed with our analyses, let us briefly summarize main properties
of the dynamics of (1) that have already been described in the literature.
Note that, in general, the overall analysis of the phase plane dynamics for
the model (1) has not been conducted. There are, however, valuable
observations for some ranges of parameters supported by numerical simulations.

In the paper [5], in which this model was introduced, Chialvo discusses only
the case $k=0$, and then $k$ of small positive value with two prescribed
choices of the other parameters. For $k=0$, the point
$(x_{\\_}{\textrm{f}\,0},y_{\\_}{\textrm{f}\,0}):=(0,c/(1-a))$ is always a
stable (attracting) fixed point of the system (since the corresponding
eigenvalues are $0$ and $a<1$). For $k\neq 0$, Chialvo [5] treats only the
case when the phase portrait has exactly one equilibrium point (which happens,
e.g., for $a=0.89$, $b=0.6$ and $c=0.28$) and treats $k$ as the bifurcation
parameter, while $a$, $b$ and $c$ are usually kept constant. For this
particular choice of parameter values and small values of $k$, the unique
fixed point is globally attracting and this parameter regime is referred to as
_quiescent-excitable_ regime. For larger values of $k$ (e.g., for $k=0.1$),
the unique fixed point is no longer stable, and oscillatory solutions might
appear. This phenomenologically corresponds to the bifurcation from quiescent-
excitable to oscillatory solution. When the value of $k$ is increased a bit
more, chaotic-like behavior was observed in [5] for some values of the
parameters. For example, when $b$ is decreased from $0.6$ to $0.18$, and
$a=0.89$, $c=0.28$, $k=0.03$ then instead of periodic-like oscillations the
solution displays chaotic bursting oscillations with large spikes often
followed by a few oscillations of smaller amplitude, resembling so-called
mixed-mode oscillations (see Figure 10 in [5]).

The system (1) can have up to three fixed points and their existence and
stability as well as bifurcations were studied analytically in [15]. Numerical
simulations described in another work [38] suggest the existence of an
interesting structure in the $(a,b)$-parameter space (with $c=0.28$ and
$k=0.04$ fixed), including comb-shaped periodic regions (corresponding to
period-incrementing bifurcations), Arnold tongue structures (due to the
period-doubling bifurcations) and shrimp-shaped structures immersed in large
chaotic regions.

Let us also mention the fact that the recent work [21] studies in detail the
dynamics of the reduced Chialvo model, i.e., the evaluation of the membrane
voltage given by the first equation in (1), with $y_{\\_}n=\mathrm{const}$
treated as a parameter. These purely analytical studies take advantage of the
fact that the one-dimensional map $F_{\\_}y(x)=x^{2}\exp(y-x)+k$ of the
$x$-variable, restricted to the invariant interval of interest, is unimodal
with negative Schwarzian derivative, which makes it possible to use the well-
developed theory of S-unimodal maps.

Despite all these important observations mentioned above, it is clear that the
description of the dynamics of the model (1) in the existing literature is
very incomplete. In our research described in this paper, we aimed at
obtaining a better understanding of the two-dimensional model (1) in a
reasonable parameter range. We focused mainly on investigating the set of
parameters that covered most of the analyses conducted in [5]. Specifically,
we fixed $a=0.89$ and $c=0.28$, and we made the other two parameters vary. We
first studied the range $(b,k)\in\Lambda_{\\_}1:=[0,1]\times[0,0.2]$, and
later also its sub-region $\Lambda_{\\_}2:=[0,1]\times[0.015,0.030]$, see
Section 4 for the detailed results. In addition to typical behavior observed
in [5], including attracting points, attractor-repeller pairs consisting of a
point and a periodic orbit, and chaotic behavior as well, we detected many
regions with other types of interesting dynamics, especially for very small
values of the parameters $b$ and $k$.

##  3\. Method

In this section, we describe the method for computer-assisted analysis of
dynamics in a system with a few parameters, first introduced in [1] for
discrete-time dynamical systems, and further extended in [19] to flows. We
also introduce a new method based on the notion of Finite Resolution
Recurrence that provides insight into the dynamics inside chain recurrent
sets. This approach provides a considerable improvement, because – to the best
of our knowledge – in methods based on [1] introduced so far this kind of
analysis that would reveal the internal structure of chain recurrent sets was
never proposed.

The computations are conducted for entire intervals of parameters, and the
results are valid for each individual parameter in the interval. This allows
one to determine the dynamical features for entire parameter ranges if those
are subdivided into smaller subsets. By using interval arithmetic and
controlling the rounding of floating point numbers, the method provides
mathematically rigorous results (a.k.a.  _computer-assisted proof_).

The main steps of our method are summarized in Figure 1 and described in
Sections 3.1–3.3.

Figure 1. Schematic illustration of the main steps of our method. The
equations of the system, together with a chosen range of parameters
$\Lambda\subset\mathbb{R}^{m}$ and a specification of the phase space
($B\subset\mathbb{R}^{n}$), constitute input to a computer program: an
implementation of the method described in Sections 3.1–3.3. The computer
program can be treated as a “black box” that produces desired information
about the system as its output. Namely, the first thing to obtain is a
Continuation Diagram that partitions parameter space $\Lambda$ into classes
corresponding to equivalent dynamics. Next, for each particular value of the
parameters, a Conley-Morse graph and a corresponding phase-space portrait can
be obtained. This provides information on invariant sets. Finally, the
dynamics within individual Morse sets can be examined in more detail by means
of Recurrence Diagrams that allow one to deduce possible existence of periodic
attractors and chaotic dynamics.

###  3.1. Automatic analysis of global dynamics

While numerical simulations based on computing individual trajectories may
provide some insight into the dynamics, considerably better understanding may
be achieved by set-oriented methods in which entire sets are iterated by the
dynamical system. One of the first software packages that used this approach
was GAIO [10]. The first step is to partition the phase space into a
collection of bounded sets with simple structure (such as squares or cubes),
further called _grid elements_. By considering images of these sets, one can
represent a map that generates a discrete-time dynamical system as a directed
graph on grid elements. Numerical methods based on interval arithmetic provide
means for computing an outer enclosure of the map rigorously and effectively.
And here comes the key idea. Effective graph algorithms applied to such a
representation of the map make it possible to capture all the chain recurrent
dynamics contained within a given subset of the phase space. This dynamics is
captured in a collection of so-called _Morse sets_. In particular, this
construction proves that the dynamics in the remaining part of the phase space
is gradient-like (see also [2, 17]). By determining possible connections
between the chain recurrent components, one constructs so-called _Morse
decomposition_ , and the Conley index [8] provides additional information
about the invariant part of the Morse sets.

Formally, let $X$ be a topological space, and let $f\colon X\to X$ be a
continuous map. $S\subset X$ is called an invariant set with respect to $f$ if
$f(S)=S$. The invariant part of a set $N\subset X$ is an invariant set defined
as $\operatorname{Inv}N:=\bigcup\\{S\subset N:f(S)=S\\}$. An isolating
neighborhood is a compact set $N$ whose invariant part is contained in its
interior: $\operatorname{Inv}N\subset\operatorname{int}N$. A set $S$ is called
an isolated invariant set if $S=\operatorname{Inv}N$ for some isolating
neighborhood $N$.

A Morse decomposition of $X$ with respect to $f$ is a finite collection of
disjoint isolated invariant sets (called Morse sets)
$S_{\\_}1,\ldots,S_{\\_}p$ with a strict partial ordering $\prec$ on the index
set $\\{1,\ldots,p\\}$ such that for every $x\in
X\setminus(S_{\\_}1\cup\cdots\cup S_{\\_}p)$ and for every orbit
$\\{\gamma_{\\_}k\\}_{\\_}{k\in{\mathbb{Z}}}$ (that is, a bi-infinite sequence
for which $f(\gamma_{\\_}k)=\gamma_{\\_}{k+1}$) such that $\gamma_{\\_}0=x$
there exist indices $i\prec j$ such that $\gamma_{\\_}k\to S_{\\_}i$ as
$k\to\infty$ and $\gamma_{\\_}k\to S_{\\_}j$ as $k\to-\infty$.

Since it is not possible, in general, to construct numerically a valid Morse
decomposition of a compact set $B\subset{\mathbb{R}}^{n}$, we construct
isolating neighborhoods of the Morse sets instead. This is a family of
isolating neighborhoods $N_{\\_}1,\ldots,N_{\\_}p\subset B$ with a strict
partial ordering $\prec$ on the set of their indices such that the family
$\\{S_{\\_}i:=\operatorname{Inv}N_{\\_}i:i=1,\ldots,p\\}$ forms a Morse
decomposition of $\operatorname{Inv}B$ with the ordering $\prec$. The sets
$N_{\\_}i$, $i=1,\ldots,p$, will be called numerical Morse sets, and the
collection $N_{\\_}1,\ldots,N_{\\_}p$ is then a numerical Morse decomposition.

We visualize a numerical Morse decomposition by means of a directed graph that
corresponds to the transitive reduction of the relation $\prec$. Vertices in
this graph correspond to the numerical Morse sets, and a path from $N_{\\_}i$
to $N_{\\_}j$ indicates the possibility of existence of a connecting orbit
between them.

We construct numerical Morse sets as finite unions of small rectangles in
${\mathbb{R}}^{n}$ whose vertices form a regular mesh. Specifically, a
rectangular set in ${\mathbb{R}}^{n}$ is a product of compact intervals. Given
a rectangular set

$B=[a_{\\_}1,a_{\\_}1+\delta_{\\_}1]\times\cdots\times[a_{\\_}n,a_{\\_}n+\delta_{\\_}n]\subset{\mathbb{R}}^{n}$

and integer numbers $d_{\\_}1,\ldots,d_{\\_}n>0$, the set

${\mathcal{G}}_{\\_}{d_{\\_}1,\ldots,d_{\\_}n}(B):=\bigg{\\{}\prod_{\\_}{i=1}^{n}[a_{\\_}i+\frac{j_{\\_}i}{d_{\\_}{i}}\delta_{\\_}i,a_{\\_}i+\frac{j_{\\_}i+1}{d_{\\_}{i}}\delta_{\\_}i]:\\\
j_{\\_}i\in\\{0,\ldots,d_{\\_}i-1\\},i\in\\{1,\ldots,n\\}\bigg{\\}}$

is called the $d_{\\_}1\times\cdots\times d_{\\_}n$ uniform rectangular grid
in B. The grid elements are referred to by the $n$-tuples
$(j_{\\_}1,\ldots,j_{\\_}n)$. The $n$-tuple of integers
$(d_{\\_}1,\ldots,d_{\\_}n)$ is called the resolution in $B$. We shall often
write ${\mathcal{G}}(B)$ instead of
${\mathcal{G}}_{\\_}{d_{\\_}1,\ldots,d_{\\_}n}(B)$ for short.

Note that in the planar case, a numerical Morse decomposition in
$B\subset{\mathbb{R}}^{2}$ can be visualized as a digital raster image whose
pixels correspond to the individual boxes in each of the numerical Morse sets.
For convenience, each numerical Morse set can be plotted with a different
color. Obviously, this kind of visualization should be accompanied by a graph
that shows the ordering $\prec$.

A multivalued map
${\mathcal{F}}\colon{\mathcal{G}}(B)\multimap{\mathcal{G}}(R)$, where
$B\subset R\subset{\mathbb{R}}^{n}$ and ${\mathcal{G}}(R)$ is a uniform
rectangular grid containing ${\mathcal{G}}(B)$, is called a _representation_
of a continuous map $f\colon B\to R$ if the image $f(Q)$ of every grid element
$Q\in{\mathcal{G}}(B)$ is contained in the interior of the union of grid
elements in ${\mathcal{F}}(Q)$. If $|{\mathcal{A}}|$ denotes the union of all
the grid elements that belong to the set
${\mathcal{A}}\subset{\mathcal{G}}(R)$ then this condition can be written as
follows:

(3) $f(Q)\subset\operatorname{int}|{\mathcal{F}}(Q)|\text{ for every
}Q\in{\mathcal{G}}(B).$

A representation corresponds to a directed graph $G=(V,E)$ whose vertices are
grid elements and directed edges are defined by the mapping ${\mathcal{F}}$ as
follows: $(P,Q)\in E\iff Q\in{\mathcal{F}}(P)$. It is a remarkable fact that
the decomposition of $G$ into _strongly connected path components_ (maximal
collections of vertices connected in both directions by paths of nonzero
length) yields a numerical Morse decomposition in $B$, provided that
$N_{\\_}i\subset\operatorname{int}B$ for all the numerical Morse sets
$N_{\\_}i$; see [1, 2, 17] for justification.

Now consider a dynamical system that depends on $m$ parameters. Consider a
rectangular set $\Lambda\in{\mathbb{R}}^{m}$ of all the $m$ parameter values
of interest. Take a uniform rectangular grid
${\mathcal{G}}_{\\_}{s_{\\_}1,\ldots,s_{\\_}m}(\Lambda)$ for some positive
integers $s_{\\_}1,\ldots,s_{\\_}m$. Using interval arithmetic, one can
compute a representation ${\mathcal{F}}_{\\_}L$ valid for the maps
$f_{\\_}\lambda$ for all the parameters $\lambda\in L$. Then the numerical
Morse decomposition computed for $F_{\\_}L$ yields a collection of isolating
neighborhoods of a Morse decomposition for each $f_{\\_}\lambda$, where
$\lambda\in L$.

Given two parameter boxes $L_{\\_}1,L_{\\_}2\in{\mathcal{G}}(\Lambda)$, we use
the clutching graph introduced in [1, §3.2] to check if the numerical Morse
sets in the computed two numerical Morse decompositions are in one-to-one
correspondence. If this is the case then continuation of Morse decompositions
has been proved and we consider the dynamics found for the parameter boxes
$L_{\\_}1,L_{\\_}2$ equivalent. A visualization of the collection of
equivalence classes with respect to this relation is called a _continuation
diagram_.

###  3.2. Analysis of dynamics in individual Morse sets using the Conley index

In order to understand the dynamics around each numerical Morse set
${\mathcal{N}}_{\\_}i$ constructed by the method introduced in Section 3.1, we
check its stability by computing its forward image by ${\mathcal{F}}$ and
analyzing the part that “sticks out”:
${\mathcal{F}}({\mathcal{N}}_{\\_}i)\setminus{\mathcal{N}}_{\\_}i$. We say
that ${\mathcal{N}}_{\\_}i$ is attracting if
${\mathcal{F}}({\mathcal{N}}_{\\_}i)\subset{\mathcal{N}}_{\\_}i$. One can
prove that then $|{\mathcal{N}}|$ contains a local attractor, which justifies
this term. If
${\mathcal{F}}({\mathcal{N}}_{\\_}i)\not\subset{\mathcal{N}}_{\\_}i$ then we
say that $N$ is unstable. We qualify the kind of instability by computing the
Conley index using the approach introduced in [1, 23, 27].

The definition of the Conley index is based on the notion of an _index pair_.
This is a pair of sets $(P_{\\_}1,P_{\\_}2)$ such that $P_{\\_}1$ covers an
isolating neighborhood, and trajectories exit this neighborhood through
$P_{\\_}2$; see e.g. [1] for the precise definition. A few typical Conley
indices that might appear in the actual computations are shown in Figures 2
and 3. In the case of a flow, the homological Conley index is merely the
relative homology of the index pair. However, in the case of a map, one also
needs to consider the homomorhpism induced in homology by the map on the index
pair (denoted here by $H_{\\_}*(I_{\\_}P)$), with some reduction applied to
it; see e.g. [34] for the details. In order to simplify the representation of
the Conley index for a map, we compute the non-zero eigenvalues of the index
map, which is a weaker but easily computable invariant; see [1] for more
explanations on this approach.

$P_{\\_}1$(a)

An attracting fixed point
$P_{\\_}1$ is the disk shown in blue, $P_{\\_}2=\emptyset$
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong({\mathbb{Z}})$
index map at level $0$: $g_{\\_}1\mapsto g_{\\_}1$
eigenvalues: $\\{1\\}$ at level $0$
example map: $f(x,y)=(x/2,y/2)$

$f$$f$$P_{\\_}1$$P_{\\_}1$(b)

An attracting period-$2$ orbit
$P_{\\_}1$ is the union of the two disks shown in blue, $P_{\\_}2=\emptyset$
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong({\mathbb{Z}}^{2})$
index map at level $0$: $g_{\\_}1\mapsto g_{\\_}2$, $g_{\\_}2\mapsto g_{\\_}1$
eigenvalues: $\\{-1,1\\}$ at level $0$
example map: $f(x,y)=(-x(2-|x|),y/2)$
attracting periodic orbit for the example map: $(\pm 1,0)$

$f$$f$$f$$P_{\\_}1$$P_{\\_}1$$P_{\\_}1$(c)

An attracting period-$3$ orbit
$P_{\\_}1$ is the union of the three disks shown in blue, $P_{\\_}2=\emptyset$
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong({\mathbb{Z}}^{3})$
index map at level $0$: $g_{\\_}1\mapsto g_{\\_}2$, $g_{\\_}2\mapsto
g_{\\_}3$, $g_{\\_}3\mapsto g_{\\_}1$
eigenvalues: $\\{1,-\tfrac{1}{2}\pm\mathrm{i}\tfrac{\sqrt{3}}{{2}}\\}$ at
level $0$

$P_{\\_}1$(d)

An attracting invariant circle
$P_{\\_}1$ is the ring shown in blue, $P_{\\_}2=\emptyset$
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong({\mathbb{Z}},{\mathbb{Z}})$
index map at level $0$: $g_{\\_}1\mapsto g_{\\_}1$; at level $1$:
$g_{\\_}1\mapsto g_{\\_}1$
eigenvalues: $\\{1\\}$ at level $0$, $\\{1\\}$ at level $1$

Figure 2.  Typical Conley indices for stable isolated invariant sets that
appear in actual applied dynamical systems. The examples (a) and (d) may come
from a time-$t$ map for a flow (the index map is thus the identity).

$P_{\\_}1$$P_{\\_}2$(a)

A saddle fixed point (orientation preserving)
$P_{\\_}1$ is the blue rectangle together with its two sides $P_{\\_}2$
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong(0,{\mathbb{Z}})$
index map at level $1$: $g_{\\_}1\mapsto g_{\\_}1$
eigenvalues: $\\{1\\}$ at level $1$
example map: $f(x,y)=(2x,y/2)$

$f$$f$(b)

A saddle fixed point (orientation reversing)
$P_{\\_}1$ is the violet rectangle with its two sides $P_{\\_}2$
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong(0,{\mathbb{Z}})$
index map at level $1$: $g_{\\_}1\mapsto-g_{\\_}1$
eigenvalues: $\\{-1\\}$ at level $1$
example map: $f(x,y)=(-2x,y/2)$

$P_{\\_}1$$P_{\\_}2$(c)

A repelling fixed point (orientation preserving)
$P_{\\_}1$ is the blue disk together with its “boundary” $P_{\\_}2$
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong(0,0,{\mathbb{Z}})$
index map at level $2$: $g_{\\_}1\mapsto g_{\\_}1$
eigenvalues: $\\{1\\}$ at level $2$
example map: $f(x,y)=(2x,2y)$

$P_{\\_}1$$P_{\\_}2$(d)

A repelling invariant circle
$P_{\\_}1$ is the blue ring together with $P_{\\_}2$ shown in red
$H_{\\_}*(P_{\\_}1,P_{\\_}2)\cong(0,{\mathbb{Z}},{\mathbb{Z}})$
index map at level $1$: $g_{\\_}1\mapsto g_{\\_}1$; at level $2$:
$g_{\\_}1\mapsto g_{\\_}1$
eigenvalues: $\\{1\\}$ at level $1$, $\\{1\\}$ at level $2$

Figure 3.  Typical Conley indices for unstable isolated invariant sets that
appear in actual applied dynamical systems. The examples (a), (c), (d) may
come from a time-$t$ map for a flow (the index map is thus the identity). The
map in (b) flips the index pair horizontally and squeezes it vertically.

A selection of typical Conley indices is provided in Figures 2 and 3. In
particular, it is important to note that this index has a specific form for a
hyperbolic fixed point or a hyperbolic periodic orbit with a $d$-dimensional
unstable manifold. If we encounter one of these specific forms of the index
then we say that ${\mathcal{N}}_{\\_}i$ is of type of the corresponding
periodic point or orbit. Although in such a case $|{\mathcal{N}}_{\\_}i|$
indeed contains a periodic orbit of the expected period, it may turn out that
the stability of that orbit is different, and the dynamics inside the
numerical Morse set might be more complicated than it appears from the
outside. In particular, if $|{\mathcal{N}}_{\\_}i|\subset{\mathbb{R}}^{n}$ and
${\mathcal{N}}_{\\_}i$ is of type of a fixed point or a periodic orbit in
${\mathbb{R}}^{n}$ with $n$-dimensional unstable manifold then we say that
${\mathcal{N}}_{\\_}i$ is repelling. It is a crucial fact that if the Conley
index of ${\mathcal{N}}_{\\_}i$ is nontrivial then
$\operatorname{Inv}|{\mathcal{N}}_{\\_}i|\neq\emptyset$.

###  3.3. Finite resolution recurrence

While Conley index is a powerful topological tool that provides reliable
information about the isolating neighborhood, it sheds little light on the
actual dynamics inside the isolating neighborhood. This may be especially
disappointing if the neighborhood is large, consisting of hundreds of
thousands of grid elements. Since this set corresponds to a strongly connected
path component of the graph representation of the map, there exists a path in
the graph from every grid element to any other element, including a path back
to itself. In particular, a periodic orbit gives rise to a cycle in the graph
of the same length, and thus a fixed point (stable or not) yields a cycle of
length $1$. On the opposite, a path in the graph corresponds to a pseudo-orbit
for the underlying map: after each iteration we may need to switch to another
point within the image of the grid element in order to follow a path in the
graph by means of pieces of orbits.

Let us remark here that recurrence is one of fundamental properties of many
dynamical systems. In order to get insight into the recurrent dynamics inside
each numerical Morse set, in the remainder of this section, we introduce the
notion of Finite Resolution Recurrence (FRR for short) and an algorithm for
its analysis. This is a new method for the analysis of dynamics in a Morse set
by means of the distribution of minimum return times. For alternative
approaches to measure the recurrence in dynamical systems see, for example,
[22] and references therein.

In what follows, we define a finite resolution version of the notion of
recurrence, and we introduce an effective algorithm for its computation.

######  Definition 1.

Let ${\mathcal{F}}\colon{\mathcal{N}}\multimap{\mathcal{N}}$ be a multivalued
map on a set ${\mathcal{N}}\subset{\mathcal{G}}(B)$, and let
$Q\in{\mathcal{N}}$. The _recurrence time_ of $Q$ in ${\mathcal{N}}$ with
respect to ${\mathcal{F}}$ is defined as follows:
$\operatorname{rec}(Q):=\min\\{k>0:Q\in{\mathcal{F}}^{k}(Q)\\}$, with the
convention $\min\emptyset=\infty$.

If ${\mathcal{N}}$ is a numerical Morse set then one can restrict the codomain
of the multivalued map ${\mathcal{F}}$ on ${\mathcal{G}}(B)$ to
${\mathcal{N}}$ in order to obtain the multivalued map
${\mathcal{F}}\colon{\mathcal{N}}\multimap{\mathcal{N}}$ (with some images
possibly empty). Let $G=(V,E)$ be the directed graph that represents this
multivalued map. The recurrence time of $Q\in{\mathcal{N}}$ is the length of
the shortest non-trivial cycle passing through $Q$. If ${\mathcal{N}}$ is a
strongly connected path component of the graph $G$ then
$\operatorname{rec}(Q)$ is finite. The following algorithm allows one to
compute the recurrence times for all the elements of ${\mathcal{N}}$
effectively:

######  Algorithm 2 (computation of recurrence times).

function recurrence_times
---
input:
|  $G=(V,E)$: a directed graph
begin
|  for each $v\in V$:
|  |  for each $u\in V$:
|  |  |  $D_{\\_}{vu}:=$ length of the shortest path in $G$ from $v$ to $u$
|  for each $v\in V$:
|  |  if $(v,v)\in E$ then
|  |  |  $\operatorname{rec}(v):=1$
|  |  else
|  |  |  $\operatorname{rec}(v):=\min\\{D_{\\_}{vu}+D_{\\_}{uv}:u\in V\\}$
end.

In this algorithm, one first computes the shortest paths from every vertex
$v\in V$ to every other vertex $u\in V$. Then two possibilities are
considered. If there exists a self-loop from $v\in V$ to itself then $1$ is
recorded as the recurrence time for $v$. Otherwise, all the other loops go
through some other vertex $u\in V$, and the shortest possibilities are
explored to determine an optimal non-trivial loop from $v$ back to itself.

######  Theorem 3 (correctness and complexity of computation of recurrence
times).

Algorithm 2 called with a directed graph $G=(V,E)$ computes the recurrence
times for all the vertices $v\in V$. It can reach this goal within as little
as $O(|V||E|+|V|^{2}\log|V|)$ time and using $O(|V|^{2})$ memory.

###### Proof.

Let us first notice that Algorithm 2 computes the recurrence times for all the
vertices. Indeed, a cycle through a vertex $v\in V$ is either a self-loop (of
length 1), or a cycle that runs through another vertex $u\in V$. Then the
algorithm clearly computes the minimum length of any nontrivial cycle that
runs through $v$.

Let us now analyze the computational complexity of Algorithm 2. Dijkstra’s
algorithm can be used to compute simultaneously all the lengths of the
shortest paths that start at a fixed vertex $v\in V$. Its time complexity is
$O(|E|+|V|\log|V|)$ when implemented using the Fibonacci heap. Therefore, the
time complexity of the first loop is $O(|V||E|+|V|^{2}\log|V|)$, because
Dijkstra’s algorithm must be run for every vertex $v\in V$ as a source. The
time complexity of the second loop is $O(|V|^{2})$, because its body is run
for every $v\in V$, and the minimum is computed over all $u\in V$.

The memory complexity is $O(|V|^{2})$, or even $\Theta(|V|^{2})$, because the
entire matrix $D$ must be stored after having been computed in the first loop
for the purpose of being used in the second loop. ∎

At this point we emphasize the fact that knowledge of the recurrence time of a
grid element provides certain rigorous information about the dynamics of the
points that belong to the grid element. Specifically, if
$\operatorname{rec}(Q)=r$ then for every $x\in Q$, we know that if the
trajectory $x=x_{\\_}0,x_{\\_}1,\ldots,x_{\\_}{r-1}$ stays in
$|{\mathcal{N}}|$ then $x_{\\_}i\notin Q$ for $i=1,\ldots,r-1$. In particular,
there is no periodic orbit contained in $|{\mathcal{N}}|$ that goes through
$Q$ whose period is below $r$. On the other hand, if $\operatorname{rec}(Q)=r$
then we know that there exists a _$\delta$ -pseudo-orbit_
$x_{\\_}0,\ldots,x_{\\_}r$, with
$\delta=\max_{\\_}{Q\in{\mathcal{N}}}\operatorname{diam}(|{\mathcal{F}}(Q)|)$,
that begins and ends in $Q$; this is a sequence of points
$x_{\\_}i\in|{\mathcal{N}}|$ such that $x_{\\_}0,x_{\\_}r\in Q$, and
$\operatorname{dist}(f(x_{\\_}{i-1}),x_{\\_}{i})<\delta$ for all
$i\in\\{1,\ldots,r\\}$.

In order to illustrate the usefulness of recurrence times, we show three
examples. Figure 4 shows a numerical Morse set constructed for the well-known
$2$-dimensional Hénon map [11] that exhibits chaotic dynamics. The recurrence
times are low, and different values are scattered unevenly throughout the
entire set. It seems that some orbits with low periods were captured
correctly, especially the fixed point and the period-two orbit. Figure 5 shows
an example of a large numerical Morse set computed for the non-linear Leslie
population model discussed in [1]. The recurrence times reveal its internal
structure. Indeed, separate isolating neighborhoods for a fixed point in the
middle and period-$3$ orbits can be found for nearby parameters or when
conducting the computation at a much finer resolution. Figure 6, on the other
hand, shows a numerical Morse set that is an isolating neighborhood for a
time-$t$ discretization of the Van der Pol oscillator flow on
${\mathbb{R}}^{2}$ for the parameters for which an attracting periodic
trajectory is observed. In this example, one can clearly see the recurrence
time $1$ that corresponds to the fixed point in the middle, and high
recurrence times around $50$ that capture the stable periodic trajectory.

Figure 4. Recurrence diagram for a numerical Morse set constructed for the
Hénon map that encloses a chaotic attractor, and a histogram that shows the
amounts of the various recurrence times encountered. The fixed points and
period-$2$ orbits are clearly visible, as well as the complicated inner
structure of the set.

Figure 5. Recurrence diagram for a numerical Morse set constructed for the
Leslie population model map discussed in [1], and a histogram that shows the
amounts of the various recurrence times encountered. Inner structure of the
large isolating neighborhood is revealed, with a fixed point and period-$3$
orbits.

Figure 6. Recurrence diagram for a numerical Morse set constructed for a
time-$t$ map for the Van der Pol oscillator with a periodic attracting
trajectory, and a histogram that shows the amounts of the various recurrence
times encountered. The recurrence times for the map vary from $1$ around the
fixed point at the origin up to $59$ at the border of the set.

Based on the illustrations, one may conjecture that high local variation in
the recurrence time is an indicator of complicated dynamics, such as chaos. In
the next subsection, we introduce notion of variation of the Finite Resolution
Recurrence that we further use to effectively quantify this local variation in
recurrence time.

###  3.4. Finite Resolution Recurrence Variation

Let us recall the notion of variation of a function in multiple variables
introduced by Vitali [37]. It is a generalization of the well known notion of
variation of a real-valued function in one variable.

######  Definition 4 (see [37]).

Let $f\colon B\to{\mathbb{R}}$ be a function defined on a rectangular set
$B=[a_{\\_}1,b_{\\_}1]\times\cdots\times[a_{\\_}n,b_{\\_}n]\subset{\mathbb{R}}^{n}$.
For $k\in\\{1,\ldots,n\\}$ and $h_{\\_}k>0$, define

$\Delta^{k}_{\\_}{h_{\\_}k}(f,x):=f(x_{\\_}1,\ldots,x_{\\_}k+h_{\\_}k,\ldots,x_{\\_}n)-f(x_{\\_}1,\ldots,x_{\\_}k,\ldots,x_{\\_}n),$

and recursively:

$\displaystyle\Delta_{\\_}{h_{\\_}1}(f,x)$ $\displaystyle:=$
$\displaystyle\Delta^{1}_{\\_}{h_{\\_}1}(f,x),$
$\displaystyle\Delta_{\\_}{h_{\\_}1,\ldots,h_{\\_}k}(f,x)$ $\displaystyle:=$
$\displaystyle\Delta^{k}_{\\_}{h_{\\_}k}\big{(}\Delta_{\\_}{h_{\\_}1,\ldots,h_{\\_}{k-1}},x\big{)}.$

_Vitali variation_ of $f$ on $B$ is defined as the supremum of the sums

(4)
$\sum_{\\_}{i_{\\_}1=1}^{N_{\\_}1}\cdots\sum_{\\_}{i_{\\_}n=1}^{N_{\\_}n}\left|\Delta_{\\_}{h_{\\_}1^{i_{\\_}1},\ldots,h_{\\_}n^{i_{\\_}n}}\big{(}f,(x_{\\_}1^{i_{\\_}1},\ldots,x_{\\_}n^{i_{\\_}n})\big{)}\right|$

over all the possible finite subdivisions of $[a_{\\_}1,b_{\\_}1]$, …,
$[a_{\\_}n,b_{\\_}n]$, where $h^{i_{\\_}j}_{\\_}k$ is the difference between
the subdivision points $x_{\\_}{k}^{1},\ldots,x_{\\_}{k}^{N_{\\_}k+1}$ of
$[a_{\\_}k,b_{\\_}k]$.

Note that the recurrence time function is constant on the interior of each
grid element, and can be set to the minimum of the values of the intersecting
grid elements at each boundary point. Therefore, for the practical computation
of this variation, we only need to check the differences in the values of the
function on adjacent grid elements. In two dimensions, the formula on a grid
of $(N+1)\times(M+1)$ rectangular grid elements reduces to the following:

(5)
$\sum_{\\_}{i=1}^{N}\sum_{\\_}{j=1}^{M}\left|f(x_{\\_}{i+1},x_{\\_}{j+1})-f(x_{\\_}{i+1},x_{\\_}j)-f(x_{\\_}i,x_{\\_}{j+1})+f(x_{\\_}i,x_{\\_}j)\right|,$

and we consider only those values for which the four grid elements are all in
the set that we analyze. We make this more precise in the following.

######  Definition 5 (Finite Resolution Recurrence Variation (FRRV)).

Let ${\mathcal{N}}\subset{\mathcal{G}}(B)$ and let
$f\colon{\mathcal{N}}\to{\mathbb{R}}$. Given an $n$-tuple of integers
$(j_{\\_}1,\ldots,j_{\\_}n)$, denote by $Q(j_{\\_}1,\ldots,j_{\\_}n)$ the grid
element in ${\mathcal{G}}(B)$ referred to by this $n$-tuple of integers. Let
$I({\mathcal{N}})$ denote the set of those grid elements
$Q(j_{\\_}1,\ldots,j_{\\_}n)\in{\mathcal{N}}$ for which all the grid elements
of the form $Q(j_{\\_}1+\delta_{\\_}1,\ldots,j_{\\_}n+\delta_{\\_}n)$, where
$\delta_{\\_}i\in\\{0,1\\}$, are in ${\mathcal{N}}$. Following the idea of
Definition 4, define

$\Delta_{\\_}{k}\big{(}f,Q(j_{\\_}1,\ldots,j_{\\_}n)\big{)}:=f\big{(}Q(j_{\\_}1,\ldots,j_{\\_}k+1,\ldots,j_{\\_}n)\big{)}-f\big{(}Q(j_{\\_}1,\ldots,j_{\\_}k,\ldots,j_{\\_}n)\big{)},$

and recursively:

$\Delta_{\\_}{1,\ldots,k}\big{(}f,Q(j_{\\_}1,\ldots,j_{\\_}n)\big{)}:=\Delta_{\\_}{k}\big{(}\Delta_{\\_}{1,\ldots,{k-1}},Q(j_{\\_}1,\ldots,j_{\\_}n)\big{)}.$

_Finite Resolution Recurrence Variation_ of $f$ on ${\mathcal{N}}$ is defined
as the sum

(6) $\operatorname{FRRV}(f,{\mathcal{N}})=\sum_{\\_}{K\in
I({\mathcal{N}})}\left|\Delta_{\\_}{1,\ldots,n}(f,K)\right|$

###  3.5. Normalized Finite Resolution Recurrence Variation

It is clear that the variation grows with the increase in the size of the set,
because there are more and more adjacent grid elements along the border with
the same difference. Therefore, we propose to normalize the resulting
variation by dividing it by a quantity that corresponds to the diameter of the
set on which the variation is computed, so that the normalized variation is
independent of the resolution at which we investigate the dynamics. In an
$n$-dimensional system, we propose to divide by
$\sqrt[n]{\operatorname{card}{\mathcal{N}}}$.

Moreover, instead of the absolute value of the variation of the recurrence, we
are more interested in the relative variation in comparison to the average
recurrence times. For example, the recurrence times in the Van der Pol system
(see Figure 6) are considerably larger than in the Hénon system (see Figure
4), and thus small relative fluctuations in the recurrence times might yield
considerably higher values of the variation in the former map than in the
latter. In order to overcome this problem, it seems reasonable to further
normalize the variation by dividing it by the mean recurrence time
encountered. We thus propose the following.

######  Definition 6 (NFRRV).

Let ${\mathcal{N}}\subset{\mathcal{G}}(B)$. Let
$f\colon{\mathcal{N}}\to{\mathbb{R}}$. Let
$\overline{\operatorname{rec}}({\mathcal{N}})$ denote mean recurrence in
${\mathcal{N}}$, that is,
$\sum_{\\_}{Q\in{\mathcal{N}}}\operatorname{rec}(Q)/\operatorname{card}{\mathcal{N}}$.
_Normalized Finite Resolution Recurrence Variation (NFRRV)_ of $f$ in
${\mathcal{N}}$ is the following quantity:

(7)
$\operatorname{NFRRV}(f,{\mathcal{N}}):=\frac{\operatorname{FRRV}(f,{\mathcal{N}})}{\overline{\operatorname{rec}}({\mathcal{N}})\sqrt[n]{\operatorname{card}{\mathcal{N}}}},$

where $n$ is the dimension of the phase space.

We show the results of computation of Finite Resolution Recurrence (FRR) in
Section 4.4 and its variation (FRRV), also in the normalized form of NFRRV, in
Section 4.5.

##  4\. Results and discussion

In this section, we describe an application of the methods introduced in
Section 3 to the Chialvo model of a neuron, explained in Section 2, and we
discuss the results.

The computations were carried out at the Centre of Informatics Tricity
Academic Supercomputer & Network (https://task.gda.pl/en/) on a cluster
equipped with Intel® Xeon® CPU E5-2670 v3 at 2.30GHz running 64-bit GNU/Linux
with kernel 3.10.0. In the paper we report CPU time in terms of single CPU-
hours, and we also indicate the maximum amount of memory used by a single
process.

All the real numbers that we mention in the paper should be understood as
approximations of the actual binary floating-point numbers stored in the
computer. For example, if we say that we set a parameter value to $0.1$ then
this means that the actual value used in the computations was a binary
representation of $0.1$, which is close but not exactly equal $0.1$. Also, if
we say that we obtained a value of $0.30$ then this means that the actual
number obtained was a binary floating-point number whose approximation within
the provided precision was $0.30$, that is, the number was in $[0.295,0.305]$.

###  4.1. Application of the automatic analysis method to the Chialvo model of
a neuron

An implementation of the method for the analysis of global dynamics introduced
in Section 3 that we applied in this paper is available at [25]. The software
can run on a single computer or on a cluster, using the parallelization scheme
introduced in [24]; in particular, the software is capable of using multiple
CPU cores if run on a single machine in an appropriate way. The software uses
the Computer Assisted Proofs in Dynamics (CAPD) library [7, 18] for rigorous
numerics. Conley indices are computed using the Computational Homology Project
(CHomP) software [6] based on [16] with further improvements [23, 26, 27].

We applied this method to the system (1) in the following setup. We fixed the
parameters $a:=0.89$ and $c:=0.28$, and we made the two other parameters vary:
$(b,k)\in\Lambda_{\\_}1:=[0,1]\times[0,0.2]$. Our choice of the parameters and
their ranges was based on what Chialvo indicated in his paper [5]. We took the
$100\times 100$ uniform rectangular grid in $\Lambda_{\\_}1$ to compute the
continuation diagram.

Determining a bounded rectangular subset $B$ of the phase space that contained
the constructed numerical Morse decomposition in its interior consisted of two
steps. We first did a series of numerical simulations to make sure that at
least most prominent attractors would be captured. For that purpose, we took
$30\times 30=900$ uniformly spaced initial conditions from
$[0,100]\times[-100,100]$, and we considered $11\times 11=121$ uniformly
spaced parameter values from the ranges of interest. For each of the $121$
parameter combinations, we iterated each of the $900$ initial conditions
$10{,}000$ times by the map, and then we analyzed the ranges of the variables
for the next $100$ iterations. We found out that $x\in[0.000,7.102]$ and
$y\in[-0.575,2.545]$. In the second step, we took
$B_{\\_}0:=[-0.1,10]\times[-5,5]$, and we conducted low-resolution test
computation with the set-oriented analysis method, using the $256\times 256$
uniform grid in the phase space (without computation of Conley indices to save
time). This computation took about $10$ minutes on a laptop PC. The
constructed numerical Morse sets for all the parameter boxes considered were
contained in $[-0.022,8.55]\times[-4.54,2.78]$. We thus took
$B:=[-0.1,9]\times[-5,3]$ for the final computation. We remark that it would
be ideal to obtain some bounds for the recurrent dynamics analytically, but
the only obvious bound in our case was $x\geq 0$, and the absorbing region
$D^{+}$ determined in [9, p. 1643] turned out not to contain all the recurrent
dynamics we were interested in and thus could not be used for our purpose.

Once the set $B$ has been fixed, we conducted the computations with the
$1024\times 1024$ uniform rectangular grid in $B$. We explain that due to the
gradual refinements scheme [1, §4.2], the phase space must be subdivided into
a $d\times\cdots\times d$ uniform rectangular grid, where $d$ is a power of
$2$. The computation completed within $29$ CPU-hours, and each process used no
more than $1.8$ GB of memory. All the numerical Morse sets fit in
$[-0.003,7.24]\times[-1.35,2.58]$. Between $1$ and $4$ numerical Morse sets
were found in each numerical Morse decomposition, typically $1.24\pm 0.5$
(average $\pm$ standard deviation). Either one or two of them were attracting,
typically $1.015\pm 0.122$. The sizes of the total of $12{,}378$ Morse sets
were between $1$ and $106{,}978$ grid elements, typically $10{,}863\pm
19{,}065$, with the median of $1{,}193$. Ten continuation classes were found,
including one consisting of a single parameter box. The resulting continuation
diagram is shown in Figure 7.

Figure 7. Continuation diagram for the Chialvo model with $a=0.89$, $c=0.28$,
and $(b,k)\in\Lambda_{\\_}1=[0,1]\times[0,0.2]$ split into the $100\times 100$
uniform rectangular grid.

Since most changes in global dynamics appear to take place for a limited range
of $k$, we conducted an analogous computation with the range of parameters
limited to $(b,k)\in\Lambda_{\\_}2:=[0,1]\times[0.015,0.030]$. The set
$\Lambda_{\\_}2$ was split into the $200\times 75$ uniform rectangular grid.
The computation took a little over $34$ CPU-hours and memory usage did not
exceed $396$ MB per process. All the numerical Morse sets were contained in
$[0.015,7.07]\times[-1.33,2.56]$. The corresponding continuation diagram is
shown in Figure 8; note, however, that the colors chosen for the continuation
classes may differ from those shown in Figure 7.

Figure 8. Continuation diagram for the Chialvo model with $a=0.89$, $c=0.28$,
and $(b,k)\in\Lambda_{\\_}2=[0,1]\times[0.015,0.030]$ split into the
$200\times 75$ uniform rectangular grid. The areas (a)–(g) correspond to the
areas marked with the same labels in Figure 7, but note that the colors are
different in most cases.

The results of the computations are available in [25] for interactive
browsing. Clicking a point in the continuation diagram launches a page with
the phase space portrait of the numerical Morse decomposition computed for the
specific rectangle of parameters, as well as a visualization of the
corresponding Conley-Morse graph. The details shown in the visualization are
briefly explained in Figures 9 and 10.

$N_{\\_}0$and $N_{\\_}1$$N_{\\_}2$$N_{\\_}3$Colors assigned toconsecutive
numerical Morse sets Figure 9. Sample phase space portrait computed for the
Chialvo model, as shown in the interactive visualization available at [25],
with the gray bounding box added for clarity. The color boxes in the lower
left corner indicate the color coding of consecutive numerical Morse sets. The
numbering starts with $0$. Note the barely visible tiny sets $N_{\\_}0$
(black) and $N_{\\_}1$ (blue), both located in the left top corner of the
figure. consecutive number of the numerical Morse setand the number of grid
elements in the seta set without incoming connecting orbitsa set with incoming
and outgoing connecting orbitshomology of the index paireigenvalues for each
nontrivial mapan attractor with proven
${\mathcal{F}}({\mathcal{N}}_{\\_}i)\subset{\mathcal{N}}_{\\_}i$map at each
non-trivial homology level(each homology generator is denotedby an integer,
starting with 1) Figure 10. Sample Conley-Morse graph computed for the Chialvo
model, as shown in the interactive visualization available at [25], with the
gray background added for clarity. Attracting numerical Morse sets are
indicated by yellow rectangles (it has been proved that
${\mathcal{F}}({\mathcal{N}}_{\\_}i)\subset{\mathcal{N}}_{\\_}i$). Numerical
Morse sets that have no trajectories coming from other Morse sets in the
decomposition are indicated by green rectangles. Pass-through sets (with
orbits coming in from other numerical Morse sets and orbits leaving towards
other numerical Morse sets) are indicated by ovals.

###  4.2. Types of dynamics and bifurcations found in the system

A comprehensive overview of the types of dynamics that were found in the model
can be seen in Figure 11. For each continuation class, a simplified Conley-
Morse graph is shown. An overview of bifurcations that were detected in the
system, on the other hand, is better visible in Figure 12, where the Conley-
Morse graphs were joined by edges whenever the corresponding parameter regions
were adjacent (this adjacency can be seen in Figures 8 and 11). Let us now
discuss the Morse decompositions found in the different continuation classes
shown in Figure 8, and also the bifurcations that were observed.

Figure 11. Conley-Morse graphs for the continuation classes shown in Figure 8
for the Chialvo model with $a=0.89$, $c=0.28$, and
$(b,k)\in\Lambda_{\\_}2=[0,1]\times[0.015,0.030]$. The background color of
each frame is the same as that of the corresponding continuation class in
Figure 8.

Figure 12. Continuation graph between the continuation classes shown with the
corresponding Conley-Morse graphs for the Chialvo model with $a=0.89$,
$c=0.28$, and $(b,k)\in\Lambda_{\\_}2=[0,1]\times[0.015,0.030]$. The
background color of each frame is the same as that of the corresponding
continuation class in Figure 8. Connections corresponding to classes
intersecting by a single vertex, such as between (a) and (f), are neither
shown here nor discussed in the text for the sake of clarity.

In Region (a), there is exactly one attracting neighborhood, so the detected
dynamics is very simple. However, the constructed numerical Morse set is of
different size and shape, depending on the specific part of the region: it is
small for lower values of $b$, and suddenly increases in size above Region
(e), as shown in Figure 13.

Figure 13. Morse decompositions for selected parameter boxes in Region (a)
shown in Figure 8. The parameter boxes have integer coordinates $(21,56)$ and
$(36,56)$, and correspond to $(b,k)\in[0.105,0.110]\times[0.0262,0.0264]$, and
$(b,k)\in[0.180,0.185]\times[0.0262,0.0264]$, respectively. There is exactly
one attracting numerical Morse set in each case. Note the different scale in
both plots.

When the parameter $b$ is increased to move from Region (a) to Region (b), the
internal structure of the large isolating neighborhood is revealed, and in
Region (b) one can see it split into an attractor–repeller pair: a small
repeller ($308$ boxes) surrounded by a circle-shaped attractor (almost
$31{,}000$ boxes); see Figure 14. The Morse graph shows the Conley indices
computed for the numerical Morse sets. The exit set of the small set surrounds
it: the relative homology is like for the pointed sphere, with the identity
index map. The exit set of the large set is empty, and the index map shows
that the orientation is preserved. This situation corresponds to what is shown
in Figure 3 (c).

Figure 14. The Morse graph (left) and the Morse decomposition (right) computed
for a sample parameter box taken from Region (b) shown in Figure 8. The
parameter box with integer coordinates $(56,56)$ was chosen, which corresponds
to $(b,k)\in[0.280,0.285]\times[0.0262,0.0264]$. There is an attractor in the
shape of a circle, and a small repeller inside.

On the other hand, when we move from Region (a) through Region (c) down to
Region (d) by decreasing the parameter $k$, we observe a numerical version of
the saddle-node bifurcation. A new numerical Morse set appears in Region (c)
with trivial index, which then splits into two Morse sets: one with one
unstable direction (a saddle) and one attractor. These features can be derived
from the Conley index; see Figure 15 and compare the indices to the ones shown
in Figures 3 (a) and 2 (a).

Figure 15. The Morse graph (left) and the Morse decomposition (right) computed
for a sample parameter box taken from Region (d) shown in Figure 8. The
parameter box with integer coordinates $(24,16)$ was chosen, which corresponds
to $(b,k)\in[0.120,0.125]\times[0.0182,0.0184]$. There is a small
attractor–repeller pair and another attractor; this is a case of bi-stability.

When the parameter $b$ is increased to move from Region (d) to Region (e), the
newly created saddle joins the large attractor, and a large numerical Morse
set appears; see Figure 16. The Conley index of this large numerical Morse
set, however, is trivial, which suggests that it might contain no non-empty
invariant set. Its existence is most likely due to the dynamics slowing down
in preparation for another bifurcation. Such a bifurcation indeed appears if
we increase $b$ further to enter Region (f). The large numerical Morse set
splits into a saddle and a repeller, which is another version of the saddle-
node bifurcation; see Figure 17 and compare the indices to the ones shown in
Figures 3 (a) and 3 (c). Increasing the parameter $b$ further makes these two
sets collapse in Region (g) and disappear in Region (h).

Figure 16. The Morse graph (left) and the Morse decomposition (right) computed
for a sample parameter box taken from Region (e) shown in Figure 8. The
parameter box with integer coordinates $(35,23)$ was chosen, which corresponds
to $(b,k)\in[0.175,0.180]\times[0.0196,0.0198]$. The Conley index of the large
numerical Morse set is trivial. The tiny attractor is barely visible at the
top left corner of the figure.

Figure 17. The Morse graph (left) and the Morse decomposition (right) computed
for a sample parameter box taken from Region (f) shown in Figure 8. The
parameter box with integer coordinates $(55,23)$ was chosen, which corresponds
to $(b,k)\in[0.275,0.280]\times[0.0196,0.0198]$. There is a repeller (shown in
black), a saddle (shown in blue) and a small attractor (shown in red).

An interesting and somewhat unusual bifurcation occurs when one decreases the
parameter $k$ to move from Region (b) to Region (f). The circle-shaped
attractor observed in Region (b) splits into a node-type attractor and a
saddle, while the node-type repeller inside persists. Apparently, this might
be a saddle–node bifurcation. Right after the transition, a small neighborhood
of the attractor appears close to the large circular isolating neighborhood of
the saddle, and the latter one suddenly shrinks with further decrease in $b$;
see Figure 18.

Figure 18. Numerical Morse decompositions in transition from Region (b) to
Region (f) shown in Figure 8. The leftmost plot corresponds to the parameter
box with integer coordinates $(62,35)$, locatd in Region (b), the other two –
parameter boxes $(62,30)$ and $(62,29)$, respectively, both located in Region
(f). The actual parameters are: $b\in[0.310,0.315]$ in all the three cases,
and $k\in[0.0208,0.0210]$, $k\in[0.0220,0.0212]$, or $k\in[0.0210,0.0212]$,
respectively. In the middle plot, the attractor shown in red is barely visible
at the top left corner, surrounded by the blue set. Note that the blue set is
topologically a circle, but its leftmost edge is very thin and might not be
well visible in the plots.

There are also a few additional regions in the continuation diagram for very
small values of the parameter $b$ that can be seen in Figure 8 and can be
investigated with the interactive continuation diagram available at [25]. When
decreasing $b$ from Regions (a), (c) and (d) to Regions (a′), (c′) and (d′),
respectively, that is to $b\in[0.010,0.015]$, an attracting isolating
neighborhood splits into a period-two attracting orbit and a saddle in the
middle, with the map reversing the orientation, like in a typical period-
doubling bifurcation; compare the indices to the ones shown in Figures 2 (b)
and 3 (b). When $b$ is decreased even further, at $b\in[0.005,0.010]$, the
period-doubling bifurcation is undone, and the two numerical Morse sets again
become one. For the lowest values of $b$, that is, $b\in[0,0.005]$, an
additional numerical Morse set appears that looks like a layer on top of the
attracting numerical Morse set. Its Conley index is trivial, and thus its
appearance is most likely due to slow-down in the dynamics; see Figure 19.

Figure 19. Numerical Morse decompositions for very small values of $b$ shown
in Figure 8 as vertical stripes along the left-hand side edge of the diagram.
The parameter boxes chosen for the plots are (from left to right): $(0,51)$,
$(1,51)$, and $(2,51)$. The actual parameters are: $k\in[0.0252,0.0254]$ in
all the three cases, and $b\in[0,0.005]$, $b\in[0.005,0.01]$, or
$b\in[0.01,0.015]$, respectively. The thin numerical Morse set (drawn in
black) in the leftmost plot has trivial Conley index. The isolating
neighborhood in the middle plot is an attractor. The Morse decomposition in
the rightmost plot looks like a period-two attracting orbit with a saddle in
the middle.

We would like to point out the fact that our discussion of the dynamics and
bifurcations was only based on isolating neighborhoods and their Conley
indices. The actual dynamics might be much more subtle and complicated, and
therefore, any statements about possible hyperbolic fixed points are merely
speculations. Moreover, the actual bifurcations may take place for some nearby
parameters, at locations that are somewhat shifted from the lines shown in
Figures 7 and 8. Nevertheless, if the isolating neighborhoods are small then,
from the point of view of applications in which the accuracy is limited and
there is some noise or other disturbances, the numerical results shed light
onto the global dynamics and our discussion explains it in terms of simple
models that can be built with hyperbolic fixed points. This is not true,
however, in the cases in which the computations yield large isolating
neighborhoods. Indeed, the computational method introduced in [1] does not
provide any means for understanding the dynamics inside such sets, apart from
what can be deduced from the knowledge of their Conley indices. We proposed
some methods for this purpose in Section 3.2 and we show their application in
Section 4.4.

###  4.3. Sizes of invariant sets in the Chialvo model

The diagrams in Figures 20 and 21 show the total size of all the numerical
Morse sets constructed for all the parameter combinations considered. These
diagrams complement the corresponding continuation diagrams (Figures 7 and 8)
in providing the information about the global overview of the dynamics.

A larger numerical Morse set allows more room for fluctuations or even
appearance of complicated dynamics in a real system that is approximated by
the mathematical model. In particular, one can notice that very large
numerical Morse sets, consisting of some $80{,}000$ grid elements or more,
appear especially in two regions of the parameters: $(b,k)\approx(0.5,0.03)$
and $(b,k)\approx(0.2,0.03)$. Such parameters may make the model resistant to
purely analytical investigation due to its complexity, making our methods a
better fit for the purpose of understanding the dynamics. One may also
speculate that in the case of larger sets the dynamics is less predictable and
thus chaotic dynamics might emerge.

Although it is not entirely obvious, it may be possible that some numerical
Morse decompositions fall in the same continuation class even if the sizes of
the numerical Morse sets being matched differ considerably. Indeed, this
happens in our case. For example, sudden change in the size of the numerical
Morse sets sometimes can be found for $(b,k)\approx(0.8,0.025)$ or
$(b,k)\approx(0.17,0.026)$. The reason in all the observed cases seems to be
the emergence of cyclic behavior that corresponds to the spiking-bursting
oscillations in the Chialvo model.

An important observation is that sets of parameters that yield very large
numerical Morse sets sometimes span across adjacent continuation classes. This
can be observed, for example, for $(b,k)\approx(0.22,0.025)$, where the
Conley-Morse diagram changes only because the unstable fixed point can be
isolated from the large attracting isolating neighborhood. This results in a
qualitative change in the perception of the dynamics at the prescribed
resolution, switching from Region (a) to Region (b), even though the actual
behavior of the vast majority of trajectories might still be similar.

Figure 20. The size of the union of all the numerical Morse sets found in the
phase space for the corresponding parameters in $\Lambda_{\\_}1$ in the
Chialvo model. The black lines indicate borders between different continuation
classes; see Figure 7 for the corresponding continuation diagram. Figure 21.
The size of the union of all the numerical Morse sets found in the phase space
for the corresponding parameters in $\Lambda_{\\_}2$ in the Chialvo model. The
black lines indicate borders between different continuation classes; see
Figure 8 for the corresponding continuation diagram.

###  4.4. Recurrence analysis of large invariant sets in the Chialvo model

Through numerical simulations, we identified three kinds of dynamics that
yield large numerical Morse sets in the Chialvo model. They are shown in
Figure 22. The initial condition $(2,1.8)$ was taken in all the cases.
Plotting the trajectories for $100{,}000$ iterations was started after $1$
million of initial iterations to allow the trajectories settle down on the
attractors; we remark that $10{,}000$ initial iterations were not enough for
the winding orbit. In addition to them, an attracting fixed point is shown
that appears for another combination of the parameters.

The fixed point around $(0,2.5)$ is shown in Figure 22 in black, the prominent
periodic orbit in orange, an orbit that seems to follow a chaotic attractor is
shown in partly transparent blue, and a winding periodic orbit with weak
attraction is shown in red. This last orbit does not look like part of a
chaotic trajectory, because the few dozens points in the figure actually
correspond to $100{,}000$ iterates, so this is most likely a periodic
attractor. Note the small differences in the bifurcation parameters that yield
the qualitatively different asymptotic behavior of the orbits.

Figure 22. Four specific types of attracting orbits observed in the Chialvo
model, discussed in Section 4.4.

Recurrence diagrams for the numerical Morse sets constructed for the three
combinations of parameters shown in Figure 22 that yield large sets are shown
in Figures 23–25. Since the time complexity of the algorithm is worse than
$O(|V|^{2})$, we conducted the computations at a relatively low resolution in
the phase space in order to quickly obtain the results (in less than $2$
minutes each) and to clearly illustrate the results. The constructed numerical
Morse sets consist of $8{,}265$, $7{,}740$, and $6{,}740$ grid elements,
respectively. Instead of the exact values of the parameters $a$ and $k$,
intervals of width $0.0001$ were taken in order to make the computations more
realistic.

Recurrence histograms are shown along with each recurrence diagram in Figures
23–25. There are some differences between these histograms that reflect subtle
differences in recurrence diagrams. They are discussed in the next sections.

Figure 23. Recurrence diagram with histogram computed for the large numerical
Morse set found for the Chialvo map with $a\approx 0.89$, $b=0.18$, $c=0.28$,
and $k\approx 0.03$, for which $\operatorname{NFRRV}=1.447$. This is the range
of parameters that corresponds to the chaotic orbit observed in Figure 22.

Figure 24. Recurrence diagram with histogram computed for the large numerical
Morse set found for the Chialvo map with $a\approx 0.9$, $b=0.18$, $c=0.28$,
and $k\approx 0.022$, for which $\operatorname{NFRRV}=1.661$. This is the
range of parameters that corresponds to the clear periodic orbit observed in
Figure 22.

Figure 25. Recurrence diagram with histogram computed for the large numerical
Morse set found for the Chialvo map with $a\approx 0.9$, $b=0.18$, $c=0.28$,
and $k\approx 0.03$, for which $\operatorname{NFRRV}=1.598$. This is the range
of parameters that corresponds to the periodic orbit winding multiple times
observed in Figure 22.

###  4.5. Computation of the Normalized Finite Resolution Recurrence Variation
(NFRRV)

Let us begin by considering the Finite Resolution Recurrence Variation (FRRV)
introduced in Section 3.4. In order to check how this quantity changes with
the change of the resolution at which a numerical Morse set is computed, we
conducted the following experiment. We chose the six specific dynamical
systems whose FRR was discussed so far: the three sample systems described in
Section 3.3 and the three cases of the Chialvo model shown in Section 4.4. For
each of these systems, we constructed a numerical Morse set at a few different
resolutions. The computed FRRV is shown in Figure 26 as a function of the size
of the set, counted in terms of the number of grid elements. One can
immediately notice the increasing trend in all the cases, which justifies the
need for normalization introduced in Section 3.5.

Figure 26. Finite Resolution Recurrence Variation ($\operatorname{FRRV}$)
computed for six sample systems at several different resolutions. Figure 27.
Normalized finite resolution recurrence variation ($\operatorname{NFRRV}$)
computed for six sample systems at several different resolutions.

Figure 27 shows the Normalized Finite Resolution Recurrence Variation (NFRRV)
computed for the six sets. One can see that the values computed for the
systems that experience complicated dynamics (Hénon, Leslie) still have an
increasing trend, which is clearly due to the fact that with the increase in
the resolution, more and more details of the dynamics are revealed. However,
the values of $\operatorname{NFRRV}$ computed for the Van der Pol oscillator,
as well as for the Chialvo model, are essentially constant, with the latter
somewhat above the former, which indicates a slightly more complicated
dynamics. The fact that the Chialvo model in the parameter regimes
corresponding to chaos and long periodic orbits shows smaller values of
$\operatorname{NFRRV}$ than the Hénon attractor or the chaotic Leslie model
might also be due to the fact that the Chialvo model behaves like a fast-slow
system, in which trajectories spend long time in some parts of the phase space
and very quickly pass through others, and this behavior is common for almost
all trajectories.

A probably counter-intuitive observation is that out of the three sets
constructed for the Chialvo model, the lowest values of $\operatorname{NFRRV}$
are encountered by the set with chaotic dynamics. Since the chaotic trajectory
wanders round the set in rotational motion, its effect is the actual averaging
of recurrence times, which can be seen in Figure 23 as compared to Figures 24
and 25. This is in contrast to the Hénon attractor in which the trajectories
tend to run apart each other due to the positive Lyapunov exponent.

Figure 28 shows the results of comprehensive computation of NFRRV for the
Chialvo model for a wide range of parameters. Since the computation of FRR is
considerably more demanding (in terms of CPU time and memory usage) than the
set-oriented analysis of dynamics discussed in Section 4.1, we conducted the
computations for a somewhat narrower range of the parameters
$(b,k)\in\Lambda_{\\_}3:=[0,0.5]\times[0.017,0.027]$ subdivided into
$200\times 50$ boxes, and with the resolution reduced to $256\times 256$. The
phase space was taken as $B_{\\_}3:=[-0.1,7.5]\times[-1.3,2.7]$. Like
previously, we fixed $a:=0.89$ and $c:=0.28$. The computation time was a
little over $372$ CPU-hours and the memory usage did not exceed $491$ MB per
process. The continuation diagram that we obtained was similar to the relevant
part of the diagram shown in Figure 8, and the Conley-Morse graphs with phase
space images can be browsed at [25], additionally with recurrence diagrams and
the corresponding histograms created for the largest numerical Morse set
observed in each computation. Since the value of NFRRV computed for small sets
(say, with $\operatorname{card}{\mathcal{N}}<1000$) does not correspond to the
intuitive interpretation, the values computed for all the numerical Morse sets
of less than $1000$ elements were shown in gray in Figure 28. We use the
results of these computations in Section 5.

Figure 28. Normalized Finite Resolution Recurrence Variation
($\operatorname{NFRRV}$) computed for the Chialvo model (1) with $a=0.89$,
$b\in[0,0.5]$ (horizontal axis), $c=0.28$, and $k\in[0.017,0.027]$ (vertical
axis). The black lines indicate borders between different continuation
classes. Values computed for numerical Morse sets of less than $1{,}000$
elements were grayed out.

###  4.6. Attractor sizes in numerical simulations

In order to find regions of parameters for which chaotic attractors might
exist, we conducted the following numerical simulations. We took a set of
evenly spaced $400\times 400$ values of
$(b,k)\in\Lambda_{\\_}3=[0,0.5]\times[0.017,0.027]$. For each of the
parameters, we computed four trajectories, starting at $(0.01,0.01)$,
$(1.0,1.0)$, $(10.0,2.0)$, and $(2.0,10.0)$, respectively. We discarded the
first $10{,}000$ iterations of each trajectory in order to make it settle on a
hypothetical attractor. In order to obtain an approximation of this attractor,
we then plotted the next $10{,}000$ points at the resolution $4\times 4$ times
higher than in the computation upon which Figure 28 is based. The largest
numbers of boxes obtained in these attractor plots are shown in Figure 29 for
each parameter selection separately. Whenever a large chaotic attractor
exists, the numbers that we obtained should be large. Such regions of
parameters are clearly seen in Figure 29. Note, however, that this condition
is not sufficient, because a very long periodic trajectory might also yield a
large set. Nevertheless, we use the results of these simulations in the next
section as a reference to the regions of parameters where chaotic dynamics
might be present.

Figure 29. Maximal size of a cover of an attractor found in numerical
simulations for the Chialvo model (1) with $a=0.89$, $b\in[0,0.5]$ (horizontal
axis, $400$ individual values), $c=0.28$, and $k\in[0.017,0.027]$ (vertical
axis, $400$ individual values). The black lines indicate borders between
different continuation classes.

##  5\. Recurrence as a Tool in Classification of Dynamics

In order to demonstrate the usefulness of the quantities that were introduced
in Sections 3.3–3.5 (specifically, FRR and NFRRV), we conducted the following
two analyses. We used the data computed in Section 4.5 to group parameters by
the shape of the corresponding recurrence histogram (discussed in Section
5.1), and by the value of NFRRV together with the median of FRR (discussed in
Section 5.2). We used the unsupervised machine learning density-based spatial
clustering algorithm DBSCAN that finds prominent clusters in the data and
leaves the remaining points as “noise”. The results of this kind of clustering
are shown in Figures 30 and 33, respectively.

We would like to point out the fact that our analysis of variation of
recurrence times within large invariant sets is aimed at providing some
quantification of chaos in the Chialvo model and developing guidelines for how
to analyze other systems. So far, to the best of our knowledge, chaos in the
Chialvo model has not been studied analytically with the exception of the work
[15] where the authors studied chaos in the sense of Marotto (existence of the
so-called snap-back repeller) which is closer to the notion of Li-Yorke chaos
than chaos in the sense of Devaney (see e.g. [39]). However, the existence of
Marotto chaos in their work is mostly relevant to the case $k<0$, thus not
much connected with our results. The existence of the snap-back repeller in
biological systems was examined in [20] with the use of computations conducted
in interval arithmetic.

###  5.1. Recurrence histograms

In our first analysis, we took the recurrence histograms, such as those shown
in Figures 23, 24, 25. We removed the bar corresponding to recurrence $1$ from
each of the histograms, and we scaled the horizontal axis to make each
histogram consist of precisely $10$ bars. We then normalized the histograms so
that the sum of the heights of the bars was equal $1$. Both the original and
reduced histograms can be browsed at [25]. We equipped the $10$-dimensional
space of all such histograms with the $l^{1}$ metric (also known as the
Manhattan metric, or the taxicab metric). We restricted our attention to all
the histograms obtained for the largest numerical Morse set at each of the
$200\times 50$ parameter boxes, provided that this set consisted of at least
$1{,}000$ boxes in the phase space. Then we ran the DBSCAN clustering
algorithm on this collection of histograms. This algorithm finds core samples
of high density and expands clusters from them. It takes two parameters: the
maximum distance $\varepsilon$ between two samples for one to be considered in
the neighborhood of the other, and the number $p$ of samples in a neighborhood
of a point to be considered as a core point. We tried all the combinations of
$\varepsilon\in\\{0.1,0.2,\ldots,1.3,1.4\\}$ and
$p\in\\{50,100,150,200,250,300\\}$. The most useful clustering was obtained
for $\varepsilon=0.2$ and $p=150$, and is shown in Figure 30.

Figure 30. Clusters of similar recurrence histograms found by the DBSCAN
algorithm in the Chialvo model (1) with $a=0.89$, $b\in[0,0.5]$ (horizontal
axis), $c=0.28$, and $k\in[0.017,0.027]$ (vertical axis). The black lines
indicate borders between different continuation classes. Histograms for
numerical Morse sets of less than $1{,}000$ elements were not taken into
consideration.

It is interesting to see the histograms that correspond to the consecutive
classes; they are shown in Figure 31. The first histogram (Class 1) shows the
presence of many boxes in the numerical Morse set with a wide range of
recurrence times. The shape of the histogram resembles normal distribution
skewed to the right and thus is most similar to the histogram obtained for the
Hénon map (see Figure 4). This shape of a histogram might thus be an indicator
of chaotic dynamics. Indeed, this region of parameters coincides with the
neighborhood of part of the region in which large attractors were found in
numerical simulations (see Figure 29) that also indicate likely chaotic
dynamics.

The second histogram (Class 2) has three major peaks: the highest one in the
middle, another at high recurrence, and a considerably less prominent one for
some lower recurrence values. This result reflects the internal structure of
the numerical Morse set that can be seen in Figure 32, and both the histogram
and the recurrence diagram resemble the situation observed in the Leslie
system shown in Figure 5.

The next four histograms are similar to each other, with the high peak around
$34$, but with higher and higher recurrence values appearing and thus pushing
the peak to the left (it gets shifted from the 3rd position from the right to
the 4th, 5th, and 6th positions, respectively). All the cases are somewhat
similar to the Van der Pol oscillator case shown in Figure 6, and might
indicate simple, non-chaotic dynamics, at least as perceived at the finite
resolution.

Figure 31.  Reduced recurrence histograms computed for representative
parameter boxes from the classes shown in Figure 30.

Figure 32. Recurrence diagram with histogram computed for the large numerical
Morse set found for the Chialvo map with $a=0.89$, $b\in[0.225,0.2275]$,
$c=0.28$, and $k\in[0.0232,0.0234]$. These parameters correspond to the box
no. $(90,31)$ of parameters in the computation described in Section 4.5.

###  5.2. Clustering based on FRR and NFRRV

In our second analysis, we took the pairs of two numbers: NFRRV and the median
of FRRs computed for the largest numerical Morse set at each of the $200\times
50$ parameter boxes. We only considered numerical Morse sets of at least
$1{,}000$ boxes in the phase space. We standardized each of the two variables
(subtracted the mean and divided by the standard deviation). We tried DBSCAN
with the same values of $\varepsilon$ and $p$ as in the first analysis, but we
first multiplied the variables by $7$ in order to compensate for the lower
dimension of the space. We used the $l^{1}$ metric. The most reasonable
clustering was obtained for $\varepsilon=0.8$ and $p=100$, and is shown in
Figure 33.

Figure 33. Clusters of similar values of NFRRV together with the median of FRR
found by the DBSCAN algorithm in the Chialvo model (1) with $a=0.89$,
$b\in[0,0.5]$ (horizontal axis), $c=0.28$, and $k\in[0.017,0.027]$ (vertical
axis). The black lines indicate borders between different continuation
classes. Histograms for numerical Morse sets of less than $1{,}000$ elements
were not taken into consideration.

The characteristics of each of the six clusters found in this computation are
shown in Table 1. One can notice that the features that distinguish Cluster 1
from the others is its lowest FRR together with relatively high NFRRV, which
might be an indicator of chaotic dynamics. Indeed, low values of FRR were
found in the Hénon system (see Figure 4) and in the Leslie model (see Figure
5), with very high values of NFRRV (see Figure 27). Like in the first
analysis, this region of parameters coincides with the neighborhood of part of
the region in which large attractors were found in numerical simulations (see
Figure 29) that also indicate likely chaotic dynamics.

Cluster 3 has the second lowest median FRR; however, its NFRRV values are
relatively low. Even though numerical simulations found some relatively large
attractors in some of the parameters in Cluster 3 (see Figure 29), the FRR and
NFRRV values do not suggest the existence of chaotic dynamics and comply with
findings of our first analysis.

An interesting observation is that the considerably higher values of NFRRV in
Clusters 5 and 6 provide means to distinguish the different types of dynamics
found in the different parts of Clusters 4 and 5 found in the first analysis
(see Figure 30).

Clusters 2 and 4 have considerably lower values of NFRRV and higher median FRR
than were found in the other clusters. Their location is just below the
possible chaotic dynamics found in the numerical simulations (see Figure 29)
and thus the dynamics might be more similar to what was found in the Van der
Pol oscillator (see Figure 6).

cluster | size | NFRRV | median FRR
---|---|---|---
1 | 1133 | $1.378\pm 0.215$ | $22.675\pm 0.648$
2 | 123 | $0.753\pm 0.046$ | $35.203\pm 0.584$
3 | 752 | $0.927\pm 0.080$ | $31.134\pm 1.479$
4 | 327 | $0.635\pm 0.046$ | $38.780\pm 1.384$
5 | 315 | $1.287\pm 0.067$ | $34.032\pm 1.060$
6 | 116 | $1.084\pm 0.045$ | $40.000\pm 0.643$
Table 1.  Characteristics of the clusters found using the following two
variables: NFRRV and median FRR. The size is given in the number of parameter
boxes, the other values are provided as mean $\pm$ standard deviation.

##  6\. Final remarks

To sum up, the method that we propose in this paper allows one to obtain an
overview of global dynamics in a given system, to arrange the set of
parameters depending on the type and complexity of dynamics they yield, and to
compare the observed dynamics to other known systems. In particular, it is
possible to identify ranges of parameters in which complicated dynamics is
likely to occur. All this can be done through automated computation that
requires little human effort.

We would like to emphasize the fact that the mathematical framework and the
software introduced in this paper are not limited to dimension two. However,
we showed an application of this method to a two-dimensional system with two
varying parameters for the sake of the ease of visualization of the results.
An important objective for further work might thus be to understand the
results of the computations without the need for visualizing all the details.
The idea of Conley-Morse graphs meets this goal as far as understanding the
global dynamics in the phase space is considered. However, getting hold of the
interplay between changes in the various parameters of the system and the
corresponding changes in the dynamics might be a more demanding task.

By applying this method to the analysis of the Chialvo model, we obtained a
comprehensive overview of its global dynamics across a wide range of
parameters. We also used a heuristically justified method based on the
computation of recurrence to distinguish between the occurrence of periodic
and chaotic dynamics. The challenge for the future research would be to
provide a complete mosaic of bifurcation patterns. As it was illustrated in
Figure 22, this would be a demanding task, and thus we may repeat after
Chialvo [5] that more “work is still needed to fully understand the
bifurcation structure of Equation (1).”

## References

  * [1] Z. Arai, W. Kalies, H. Kokubu, K. Mischaikow, H. Oka, P. Pilarczyk. A database schema for the analysis of global dynamics of multiparameter systems. SIAM J. Appl. Dyn. Syst. 8 (2009), 757–789.
doi:10.1137/080734935

  * [2] H. Ban, W. Kalies. A computational approach to Conley’s decomposition theorem. Journal of Computational and Nonlinear Dynamics 1 (2006), 312–319.
doi:10.1115/1.2338651

  * [3] R. Brette, W. Gerstner. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94 (2005), 3637–3642.
doi:10.1152/jn.00686.2005

  * [4] A.N. Burkitt. A Review of the Integrate-and-fire Neuron Model: I. Homogeneous Synaptic Input. Biol Cybern 95 (2006), 1–19.
doi:10.1007/s00422-006-0068-6

  * [5] D. Chialvo. Generic excitable dynamics on a two-dimensional map. Chaos, Solitons & Fractals 5 (1995), 461–479.
doi:10.1016/0960-0779(93)E0056-H

  * [6] Computational Homology Project, http://www.pawelpilarczyk.com/chomp/ (accessed on September 7, 2022).
  * [7] Computer Assisted Proofs in Dynamics group, http://capd.ii.uj.edu.pl/ (accessed on September 7, 2022).
  * [8] C. Conley. Isolated invariant sets and the Morse index. CBMS Regional Conference Series in Mathematics Vol. 38, American Mathematical Society, Providence, R.I., 1978.
  * [9] M. Courbage, V.I. Nekorkin. Map based models in neurodynamics. Internat. J. Bifur. Chaos Appl. Sci. Engrg. 20 (2010), 1631–1651.
doi:10.1142/S0218127410026733

  * [10] M. Dellnitz, G. Froyland, O. Junge. The algorithms behind GAIO—set oriented numerical methods for dynamical systems. In Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems, B. Fiedler, ed., Springer-Verlag, New York, 2001, pp. 145–174.
doi:10.1007/978-3-642-56589-27

  * [11] M. Hénon. A two-dimensional mapping with a strange attractor. Communications of Mathematical Physics 50 (1976) 69–77.
doi:10.1007/BF01608556

  * [12] A.L. Hodgkin, A.F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117 (1952) 500–544.
doi:10.1113/jphysiol.1952.sp004764

  * [13] B. Ibarz, J.M. Casado, M.A.F. Sanjuán. Map-based models in neuronal dynamics. Physics Reports 501 (2011), 1–74
doi:10.1016/j.physrep.2010.12.003

  * [14] E. Izhikevich. Simple model of spiking neurons. IEEE Trans. Neural Netw. 14 (2003), 1569–1572.
doi:10.1109/TNN.2003.820440

  * [15] Z. Jing, J. Yang, W. Feng. Bifurcation and chaos in neural excitable system. Chaos, Solitons & Fractals 27 (2006), 197–215.
doi:10.1016/j.chaos.2005.04.060

  * [16] T. Kaczynski, K. Mischaikow, and M. Mrozek. Computational homology. Volume 157 of Applied Mathematical Sciences, Springer-Verlag, New York, 2004.
  * [17] W.D. Kalies, K. Mischaikow, R.C.A.M. VanderVorst. An algorithmic approach to chain recurrence. Found. Comput. Math. 5 (2005), 409–449.
doi:10.1007/s10208-004-0163-9

  * [18] T. Kapela, M. Mrozek, D. Wilczak, P. Zgliczyński. CAPD::DynSys: a flexible C++ toolbox for rigorous numerical analysis of dynamical systems. Communications in Nonlinear Science and Numerical Simulation 101 (2021), 105578.
doi:10.1016/j.cnsns.2020.105578

  * [19] D.H. Knipl, P. Pilarczyk, G. Röst. Rich bifurcation structure in a two-patch vaccination model. SIAM J. Appl. Dyn. Syst. 14 (2015), 980–1017.
doi:10.1137/140993934

  * [20] K.-L. Liao, C.-W. Shih, C.-J. Yu. The snapback repellers for chaos in multi-dimensional maps. J. Comput. Dyn. 5 (2018), no. 1–2, 81–92.
doi:10.3934/jcd.2018004

  * [21] F. Llovera Trujillo, J. Signerska-Rynkowska, P. Bartłomiejczyk. Periodic and chaotic dynamics in a map-based neuron model. submitted
doi:10.48550/arXiv.2111.14499

  * [22] N. Marwan, C.L. Webber. Mathematical and Computational Foundations of Recurrence Quantifications. In Webber, Jr. C., Marwan N. (eds) Recurrence Quantification Analysis. Understanding Complex Systems. (2015) Springer, Cham.
doi:10.1007/978-3-319-07155-81

  * [23] K. Mischaikow, M. Mrozek, P. Pilarczyk. Graph approach to the computation of the homology of continuous maps. Found. Comput. Math. 5 (2005), 199–229.
doi:10.1007/s10208-004-0125-2

  * [24] P. Pilarczyk. Parallelization method for a continuous property. Found. Comput. Math. 10 (2010), 93–114.
doi:10.1007/s10208-009-9050-8

  * [25] P. Pilarczyk. Topological-numerical analysis of a two-dimensional discrete neuron model. Data and software (2022). http://www.pawelpilarczyk.com/neuron/ (accessed on September 7, 2022).
  * [26] P. Pilarczyk, P. Real. Computation of cubical homology, cohomology, and (co)homological operations via chain contraction. _Adv. Comput. Math._ 41 (2015), 253–275.
doi:10.1007/s10444-014-9356-1

  * [27] P. Pilarczyk, K. Stolot. Excision-preserving cubical approach to the algorithmic computation of the discrete Conley index. Topology Appl. 155 (2008), 1149–1162.
doi:10.1016/j.topol.2008.02.003

  * [28] J.E. Rubin, J. Signerska-Rynkowska, J. Touboul, A. Vidal. Wild oscillations in a nonlinear neuron model with resets: (I) Bursting, spike-adding, and chaos. Discrete Contin. Dyn. Syst. Ser. B. 22 (2017), 3967–4002.
doi:10.3934/dcdsb.2017204

  * [29] J.E. Rubin, J. Signerska-Rynkowska, J. Touboul. Type III responses to transient inputs in hybrid nonlinear neuron models. SIAM J. Appl. Dyn. Syst. 20 (2021), 953–980.
doi:10.1137/20M1354970

  * [30] N.F. Rulkov. Regularization of synchronized chaotic bursts. Phys. Rev. Lett. 86 (2001), 183–186.
doi:10.1103/PhysRevLett.86.183

  * [31] N.F. Rulkov. Modeling of spiking-bursting neural behaviour using two-dimensional map. Physical Review E 65 (2002), 041922.
doi:10.1103/PhysRevE.65.041922

  * [32] A.L. Shilnikov, N.F. Rulkov. Subthreshold oscillations in a map-based neuron model. Phys. Lett. A 328 (2004), 177–184.
doi:10.1016/j.physleta.2004.05.062

  * [33] J. Signerska-Rynkowska. Analysis of interspike-intervals for the general class of integrate-and-fire models with periodic drive. Math. Model. Anal. 20 (2015), 529–551.
doi:10.3846/13926292.2015.1085459

  * [34] A. Szymczak. The Conley index for discrete semidynamical systems. Topology Appl. 66 (1995), 215–240.
doi:10.1016/0166-8641(95)0003J-S

  * [35] N. Torres, M.J. Cáceres, B. Perthame, D. Salort. An elapsed time model for strongly coupled inhibitory and excitatory neural networks. Physica D: Nonlinear Phenomena 425 (2021), 132977.
doi:10.1016/j.physd.2021.132977

  * [36] J. Touboul, R. Brette. Spiking dynamics of bidimensional integrate-and-fire neurons. SIAM J. Appl. Dyn. Syst. 8 (2009), 1462–1506.
doi:10.1137/080742762

  * [37] Vitali variation. Encyclopedia of Mathematics. http://encyclopediaofmath.org/index.php?title=Vitali_variation&oldid=44408 (accessed on September 7, 2022)
  * [38] F. Wang, H. Cao. Mode locking and quasiperiodicity in a discrete-time Chialvo neuron model. Commun. Nonlinear Sci. Numer. Simul. 56 (2018), 481–489.
doi:10.1016/j.cnsns.2017.08.027

  * [39] Y. Zhao, L. Xie, K.F.C. Lingli. An improvement on Marotto’s theorem and its applications to chaotification of switching systems. Chaos Solitons Fractals 39 (2009), no. 5, 2225–2232.
doi:10.1016/j.chaos.2007.06.109