# The generation of diverse traveling pulses and its solution scheme in an
excitable slow-fast dynamics

Arnab Mondal1, Argha Mondal2,3,∗, M. A. Aziz-Alaoui4, Ranjit Kumar Upadhyay1,
Sanjeev Kumar Sharma1, Chris G. Antonopoulos3 1 Department of Mathematics &
Computing, Indian Institute of Technology (Indian School of Mines), Dhanbad
826004, India
2 Department of Mathematics, Sidho-Kanho-Birsha University, Purulia 723104,
WB, India
3 Department of Mathematical Sciences, University of Essex, Wivenhoe Park CO4
3SQ, UK
4 Normandie Univ, UNIHAVRE, LMAH, FR-CNRS-3335, ISCN, 76600 Le Havre, France

###### Abstract

In this paper, we report on the generation and propagation of traveling pulses
in a homogeneous network of diffusively coupled, excitable, slow-fast
dynamical neurons. The spatially extended system is modelled using the nearest
neighbor coupling theory, in which the diffusion part measures the spatial
distribution of the coupling topology. We derive analytically the conditions
for traveling wave profiles that allow the construction of the shape of
traveling nerve impulses. The analytical and numerical results are used to
explore the nature of the propagating pulses. The symmetric or asymmetric
nature of the traveling pulses is characterized and the wave velocity is
derived as a function of system parameters. Moreover, we present our results
for an extended excitable medium by considering a slow-fast biophysical model
with a homogeneous, diffusive coupling that can exhibit various traveling
pulses. The appearance of series of pulses is an interesting phenomenon from
biophysical and dynamical perspective. Varying the perturbation and coupling
parameters, we observe the propagation of activities with various amplitude
modulations and transition phases of different wave profiles that affect the
speed of the pulses in certain parameter regimes. We observe different types
of traveling pulses, such as envelope solitons and multi-bump solutions and
show how system parameters and the coupling play a major role in the formation
of different traveling pulses. Finally, we obtain the conditions for stable
and unstable plane waves.

###### keywords:

FitzHugh-Rinzel model, traveling waves, complex Ginzburg-Landau equation,
perturbation theory, solution scheme, stability, instability

††journal: Journal of LaTeX Templates††footnotetext: Corresponding author
Email address<EMAIL_ADDRESS>
> The emergence of diverse traveling wave profiles might reflect the rhythmic
> properties in signal processing and functional connectivity in different
> brain areas. In [1], we have shown the existence of solitary wave profiles
> in the Hindmarsh-Rose model using the tanh method. Here, we use a different
> analytical approach to find the traveling wave profiles for a slow-fast
> FitzHugh-Rinzel model. We observe different types of traveling pulses, such
> as envelope solitons and multi-bump solutions and show how system parameters
> and the coupling play a major role in the formation of various traveling
> pulses. Finally, we obtain the conditions for which plane waves become
> stable or unstable.

##  1 Introduction

Excitable biophysical mechanisms show various dynamical characteristics and
can be studied from a dynamical systems perspective. Neurons receive incoming
sensory inputs, encode them into various biophysical variables and produce
relevant outputs [2, 3]. Theoretical modeling and numerical simulations are
prominent tools in understanding various functional mechanisms in neural
computations [4, 5, 6]. Hindmarsh-Rose (HR) [7], Izhikevich [2, 3], FitzHugh-
Nagumo (FHN) [8] and FitzHugh-Rinzel (FHR) [9, 10] systems are such types of
biophysical models that produce diverse firing properties. In extended spatial
systems of coupled neuronal populations, a transmembrane potential difference
can travel across neurons by means of traveling wave propagation [11, 12, 13,
14, 15]. This may cause the appearance and disappearance of neuronal responses
with traveling pulses. Such types of dynamics are responsible for playing a
major role in the biodynamics of signal processing [16].

Traveling pulses present a characteristic feature of certain neurological
disorders in humans, including epilepsy [17] and migraines [18]. Therefore,
investigating the mechanisms underlying wave propagation in neuronal tissue is
important in understanding both normal and pathological brain states. The
spatial propagation of wavy profiles in coupled neuronal networks can emerge
due to variations in the electrical activities at the single cell level. One
of the key challenges is to analyze the intrinsic dynamics of neuronal
activities and its characteristics [19, 20, 21]. The emergence of diverse
traveling wave profiles is often associated with self-organization and
interactions within the network and might reflect the rhythmic properties in
signal processing, and functional connectivity in different brain areas [22,
23, 24].

The main goal here is to study the emergence of diverse traveling pulses using
a mathematical approach for spatially excitable media in the context of
excitable biophysical systems. We combine the existence and propagation of
traveling pulses with speed equations and the stability of the solutions.
Traveling pulses are special types of solutions to a reaction diffusion
system. Equilibrium solutions can lead to two types of traveling waves known
as traveling fronts and pulse solutions. A traveling front connects distinct
steady states at both ends, however, a traveling pulse connects the same
steady state [25, 26].

Our study is motivated by theoretical and experimental works that assume the
existence of localized excitations in certain neuronal populations [27, 28,
29, 30, 1]. Many research articles have investigated the existence of envelope
solitons of the nerve impulse, not only in diffusively connected networks but
also in memristive neuron models [29, 30]. In [1], we have shown the existence
of solitary wave profiles in the HR model using the tanh method [31, 32].
Further, a special type of traveling pulse, i.e., the multi-bump or multi-
pulse solution has been studied in [26]. Specific mechanisms for traveling
pulses have been examined mathematically. Excitable neuron models can produce
two types of traveling pulses: envelope solitons and envelope of multi-bump
solutions. Envelope solitons are observed in nonlinear and weakly dispersive
systems. The wave equation in the small amplitude limit is very important to
describe envelope solitons. The disappearance of bumps changes the propagation
activity. The individual bump solution in multi-pulse wave profiles are non-
propagating and transient [26]. However, there arises important questions
about the existence of analytical, exact or approximate solutions to the
excitable models as it is interesting to test the numerical solutions with
theoretical approaches.

Here, we investigate different forms of nerve impulses in terms of traveling
wave profiles in a diffusively coupled network of slow-fast spiking-bursting
FHR type 3 neuron models. The main objectives are analytical and numerical
approaches to find the dynamical behavior of modulated traveling pulses in an
excitable diffusive network of slow-fast dynamics. We use the multiple-scale
expansion method in a semi-discrete approximation to obtain the modified
complex Ginzburg-Landau (CGL) equation using perturbation theory. Next, we
derive the envelope of the traveling pulses to observe the collective
dynamics, i.e., nerve impulses. Finally, we derive the role of this type of
connected oscillators with identical coupling and the effects of small
perturbations on the solution of the wave profiles to observe both the
envelope of multipulses and single solitary pulses.

The study of the propagation of various traveling pulses in such networks may
be relevant, e.g., in brain pathologies [33] and different functional
mechanisms especially in cortical areas [29]. Our work describes an explicit
analytical scheme and numerical approaches to analyze the diffusively coupled
FHR models for a range of coupling strengths. First the network can be
considered in the Lienard form and then it can be transformed into an extended
complex Ginzburg-Landau equation (CGLE) using perturbation theory. We apply
the multiple-scale expansion theory with a semi-discrete approximation. The
reduced CGLE follows the time evolution of the modulated wave propagation in
the network. To obtain the traveling wave solution, we derive the expression
of the nerve impulse using the membrane voltage variable.

The paper is organized as follows. In Sec. 2, we describe preliminaries about
the slow-fast subsystems and present the FHR type 3 neuron model. We
demonstrate the theory to derive the approximate solution using semi-discrete
approximations. In Sec. 3, we reduce the FHR model to the CGLE and find in
Sec. 4, the solution to the latter. We also present different types of
traveling pulses, such as envelope solitons and multi-bump solutions. Further,
we discuss how these depend on system parameters. In Sec. 5, we obtain the
conditions for which the plane wave becomes stable or unstable. Finally, we
conclude our results in Sec. 6 following a brief discussion.

##  2 The biophysical excitable model and the approach to finding its solution

###  2.1 Formulation of the coupled network model

The dynamics of the minimal mathematical model to explore different time
scales of slow and fast membrane voltage processes can be defined as
$dx/dt=F(x,z)$, $dz/dt=\mu G(x,z)$, where $x\in R^{2}$ and $z\in R$ with
$0<\mu\ll 1$. The fast and slow subsystems are 2D and 1D abstractions of the
dynamics while facilitating the geometrical features in the phase space. The
dynamics of the system evolves slowly around the slow variable when bursting
interacts dynamically with the fast subsystem. The transition between
quiescent and firing states can be studied by performing a bifurcation
analysis where the bifurcation parameter is $\mu$. The slow-fast model can be
presented as a generic model of excitabilities with diverse neuronal
responses. It exhibits different nonlinear properties that are governed by the
set of three coupled ODEs in Eq. (1). Our primary goal is to study the
dynamics on a spatially extended medium with system (1) exhibiting various
oscillatory firing patterns.

The FHR model was introduced by FitzHugh and Rinzel [9, 10, 34, 35, 36], and
is an extended modification of the 2D classical FHN model [8]. A $z$-shaped
slow manifold can be observed in the case of the traditional bursting model,
while the FHR model of type 3 bursting does not have slow manifolds [37]. The
most important variable is the membrane voltage $v$. The third variable of the
model enables us to understand the control of the resting period between the
generation of two action potentials. The equations of an excitable network of
$N$ neurons with gap junction coupling are given by

$\displaystyle\dot{v}_{k}(t)$ $\displaystyle=w_{k}-4(v_{k}^{3}-v_{k})-z_{k}+I$
$\displaystyle+D(v_{k+1}-2v_{k}+v_{k-1}),$ $\displaystyle\dot{w}_{k}(t)$
$\displaystyle=-(4v_{k}+1+w_{k}),$ (1) $\displaystyle\dot{z}_{k}(t)$
$\displaystyle=\mu(1.25v_{k}-(z_{k}-z_{0})/4),$

where $k=1,\ldots,N$, $v$ is the membrane voltage and $w$ the recovery
variable. $I$ measures the magnitude of the external input stimulus current.
$z$ is a slow variable, and $z_{0}$ the slow subsystem parameter that controls
the nature of bursting. The small parameter $0<\mu\ll 1$ determines the pace
of the slow variable $z$ and $D$ is the coupling strength of the gap
junctions.

We have already studied the dynamics of system (1) in [38]. In this work, we
consider the case of two nearest-neighbors coupled with weak coupling. Neurons
can connect with their nearest neighbors electrically. The weak coupling
reflects the situation that appears in the study of bursting-like activity in
$\beta$-cell islets of the pancreas that secrete insulin [29, 39, 40].
Biophysically, the neuronal variability may reflect the characteristics at
different levels of certain receptors or differences in regulatory effects
that can be induced by internal or external neuron-modulatory processes.

In the following, we will use the semi-discrete approximation [30, 29] to
derive analytically the localized traveling pulses that propagate in a
diffusively connected excitable network of FHR neurons (1). In order to
investigate the resulting interactions, the system is considered in such a way
that each neuronal node follows the same coupling topology. The dynamics of
spike generation can be described by the dynamics of the FHR model (1).

###  2.2 The solution method

Large scale coupled neuronal networks can generate various spatial states with
different structures. One of these states give rise to traveling pulses. Here
we derive analytically the emergence of traveling pulses using the semi-
discrete approach [30, 29]. In particular, we transform system (1) in the wave
form by reducing its first and second equations into the second order
differential equation

$\displaystyle\ddot{v}_{k}+f(v_{k},\dot{v}_{k},z_{k})$
$\displaystyle=D(v_{k+1}-2v_{k}+v_{k-1})$
$\displaystyle+D(\dot{v}_{k+1}-2\dot{v}_{k}+\dot{v}_{k-1}),$ (2)
$\displaystyle\dot{z}_{k}$ $\displaystyle=g(v_{k},z_{k}),$

where the first equation is a second-order ordinary differential equation.
Equation (2) is also known as the Lienard form. The associated constant
parameters are to be derived during the study of the semi-discrete
approximation. We note that this transformation does not affect the behavior
of the system. We can easily decompose the system and obtain the transformed
slow-fast system.

To study different traveling wave profiles, we treat the slow-fast system as a
discrete coupled system. The evolution of modulated waves in the network is
described by a modified CGLE. To find the traveling wave profiles, we need to
reduce the system into a CGLE and find its solution using the semi-discrete
approximation [30, 29]. We consider the new variables $m_{k}$ and $n_{k}$ as

$\displaystyle v_{k}$ $\displaystyle=\epsilon m_{k}$ (3) $\displaystyle z_{k}$
$\displaystyle=\epsilon n_{k},$ (4)

where $0<\epsilon\ll 1$. Then, system (2) becomes

$\displaystyle\ddot{m}_{k}+f_{1}(m_{k},\dot{m}_{k},n_{k})$
$\displaystyle=D(m_{k+1}-2m_{k}+m_{k-1})$
$\displaystyle+D(\dot{m}_{k+1}-2\dot{m}_{k}+\dot{m}_{k-1}),$ (5)
$\displaystyle\dot{n}_{k}$ $\displaystyle=g_{1}(m_{k},n_{k}).$

According to the method, which is a perturbation technique, the carrier waves
are kept discrete whereas the amplitude is considered in the continuum limit
[29, 30]. The semi-discrete approximation approach is used to study plane wave
modulation due to the nonlinear terms in the system. The main objective to use
the multiple-scale expansion method is to find the solution $v(x,t)$ in terms
of the new independent space and time variables $X_{k}$ and $T_{k}$, that we
discuss next. Applying the approximation, the method allows one to study the
modulation of a plane wave due to nonlinear effects [29].

In this context, we introduce the new space and time variables
$X_{k}=\epsilon^{k}x$ and $T_{k}=\epsilon^{k}t$, where $x$ and $t$ are the
space and time variables. Then, we consider

$\displaystyle v(x,t)$
$\displaystyle=\sum_{k=1}^{\infty}\epsilon^{k}m_{k}(X_{0},X_{1},X_{2},\ldots
T_{0},T_{1},T_{2},\ldots),$ $\displaystyle z(x,t)$
$\displaystyle=\sum_{k=1}^{\infty}\epsilon^{k}n_{k}(X_{0},X_{1},X_{2},\ldots
T_{0},T_{1},T_{2},\ldots).$

The partial derivatives with respect to the new space and time variables are
given by

$\displaystyle\frac{\partial}{\partial t}$
$\displaystyle=\frac{\partial}{\partial
T_{0}}+\epsilon\frac{\partial}{\partial
T_{1}}+\epsilon^{2}\frac{\partial}{\partial T_{2}}+\ldots$
$\displaystyle\frac{\partial}{\partial x}$
$\displaystyle=\frac{\partial}{\partial
X_{0}}+\epsilon\frac{\partial}{\partial
X_{1}}+\epsilon^{2}\frac{\partial}{\partial X_{2}}+\ldots.$

Finally, to obtain the solution, we use the new form of the voltage and slow
variables, the time and space derivatives and compare various terms at
different orders of $\epsilon$.

##  3 Equation of motion of the amplitude (CGLE)

To study the envelope solitons of wave equations in the small amplitude limit,
we find the solution to this equation analytically and compute it numerically.
To proceed, it is important to transform system (1) into the wave form. To do
so, we introduce a second-order differential equation in $v_{k}$ for the first
and second equations in model (1) and obtain the system of ODEs

$\displaystyle\ddot{v}_{k}+f(v_{k},\dot{v}_{k},z_{k})$
$\displaystyle=D(v_{k+1}-2v_{k}+v_{k-1})$
$\displaystyle+D_{1}(\dot{v}_{k+1}-2\dot{v}_{k}+\dot{v}_{k-1}),$
$\displaystyle\dot{z}_{k}$ $\displaystyle=g(v_{k},z_{k}),$ (6)

where
$f(v_{k},\dot{v}_{k},z_{k})=E_{0}v_{k}+(E_{1}+E_{2}v_{k}^{2})\dot{v}_{k}+\frac{{{E_{2}}}}{3}v_{k}^{3}+E_{3}z_{k}+I_{0}$,
$g(v_{k},z_{k})=\mu(1.25v_{k}-(z_{k}-z_{0})/4)$, $E_{0}=1.25\mu$, $E_{1}=-3$,
$E_{2}=12,$ $E_{3}=1-\frac{\mu}{4}$, $I_{0}=1+\frac{\mu z_{0}}{4}-I$ and
$D_{1}=D$. This transformation does not change the system (1), it only reduces
its dimension from 3 to 2. The transformed system (6) is equivalent to system
(1) and we can obtain approximate solutions of the nonlinear system (6) using
a perturbation technique. As we need to derive the solution in an extended
diffusive medium, we will use the perturbation $v_{k}=\epsilon m_{k}$,
$z_{k}=\epsilon n_{k}$ and perturb the parameters $E_{1}$, $E_{3}$ and $D_{1}$
up to order $\epsilon^{2}$. Then, system (6) becomes

$\displaystyle\ddot{m}_{k}+f_{1}(m_{k},\dot{m}_{k},n_{k})$ $\displaystyle=$
$\displaystyle D(m_{k+1}-2m_{k}+m_{k-1})$ $\displaystyle+$
$\displaystyle\epsilon^{2}D_{1}(\dot{m}_{k+1}-2\dot{m}_{k}+\dot{m}_{k-1}),$
$\displaystyle\dot{n}_{k}$ $\displaystyle=$ $\displaystyle
g_{1}(m_{k},n_{k}),$ (8)

where
$f_{1}(m_{k},\dot{m}_{k},n_{k})=E_{0}m_{k}+\epsilon^{2}(E_{1}+E_{2}m_{k}^{2})\dot{m}_{k}+\epsilon^{2}\frac{E_{2}}{3}m_{k}^{3}+\epsilon^{2}E_{3}n_{k}$,
$g_{1}(m_{k},n_{k})=-sn_{k}+E_{0}m_{k}$ and $s=\frac{\mu}{4}$. Equation (3)
generates the dynamical behavior of the membrane voltages in coupled FHR
neurons. Equation (8) reflects the associated coupling term that describes the
dynamics of the bursting variable. Since we are dealing with a weakly coupled
diffusively connected network with nonlinear excitations, it is possible to
use the semi-discrete approximation. In the following section, we will explore
analytically and numerically the traveling pulses as envelope solitons. We
will demonstrate the method to obtain the modified CGLE. Next, we reduce it to
a CGLE by using multiple-scale expansions to find its solution.

We consider the solution of Eqs. (3)-(8) in the following forms

$\displaystyle m_{k}$
$\displaystyle=A_{k}e^{i\theta_{k}}+\overline{A}_{k}e^{-i\theta_{k}}+\epsilon\Big{(}B_{k}+C_{k}e^{2i\theta_{k}}+\overline{C}_{k}e^{-2i\theta_{k}}\Big{)},$
(9) $\displaystyle n_{k}$
$\displaystyle=F_{k}e^{i\theta_{k}}+\overline{F}_{k}e^{-i\theta_{k}}+\epsilon(G_{k}+H_{k}e^{2i\theta_{k}}+\overline{H}_{k}e^{-2i\theta_{k}}),$

where $\theta_{k}=qk-\omega t$. $q$ and $\omega$ indicate the normal mode wave
vector and angular velocity, respectively. Here, we use the continuum limit
approximation on the amplitudes $A_{k}(t)$, $B_{k}(t)$ and $C_{k}(t)$
($F_{k}(t)$, $G_{k}(t)$ and $H_{k}(t)$, respectively) as they change slowly
with respect to space and time. Using the continuum limit approximation [29]
on the amplitudes $A_{k}(t)$, $B_{k}(t)$ and $C_{k}(t)$, the amplitudes become
$A(X_{1},X_{2},T_{1},T_{2})$, $B(X_{1},X_{2},T_{1},T_{2})$ and
$C(X_{1},X_{2},T_{1},T_{2})$ (similarly for $F$, $G$ and $H$). Using Taylor
expansions, $A_{k\pm 1}$ is given by

$A_{k\pm 1}=A\pm\epsilon\frac{\partial A}{\partial
X_{1}}\pm\epsilon^{2}\frac{\partial A}{\partial
X_{2}}+\frac{\epsilon^{2}}{2}\frac{\partial^{2}A}{\partial
X_{1}^{2}}+O(\epsilon^{3}).$ (10)

The first and second-order temporal derivatives are given by

$\frac{\partial A_{k}}{\partial t}=\epsilon\frac{\partial A}{\partial
T_{1}}+\epsilon^{2}\frac{\partial A}{\partial T_{2}}+o(\epsilon^{3})$ (11)

and

$\frac{\partial^{2}A_{k}}{\partial
t^{2}}=\epsilon^{2}\frac{\partial^{2}A}{\partial T_{1}^{2}}+o(\epsilon^{3}),$
(12)

respectively. Similar expressions hold for the amplitudes $B$, $C$, $F$, $G$
and $H$. The first and second-order temporal derivatives of $m_{k}$ with
respect to the new space and time variables and amplitudes $A$, $B$, $C$ are
given by

$\displaystyle\dot{m}_{k}$ $\displaystyle=\left(\epsilon\frac{\partial
A}{\partial T_{1}}+\epsilon^{2}\frac{\partial A}{\partial T_{2}}-i\omega
A\right)e^{i\theta_{k}}$
$\displaystyle+\left(\epsilon\frac{\partial\bar{A}}{\partial
T_{1}}+\epsilon^{2}\frac{\partial\bar{A}}{\partial
T_{2}}+i\omega\bar{A}\right)e^{-i\theta_{k}}+\epsilon^{2}\frac{\partial
B}{\partial T_{1}}$ $\displaystyle+\left(\epsilon^{2}\frac{\partial
C}{\partial T_{1}}-\epsilon 2i\omega C\right)e^{2i\theta_{k}}$
$\displaystyle+\left(\epsilon^{2}\frac{\partial\bar{C}}{\partial
T_{1}}+\epsilon 2i\omega\bar{C}\right)e^{-2i\theta_{k}}+o(\epsilon^{3}),$
$\displaystyle\ddot{m}_{k}$
$\displaystyle=\left(\epsilon^{2}\frac{\partial^{2}A}{\partial
T_{1}^{2}}-\epsilon 2i\omega\frac{\partial A}{\partial
T_{1}}-\epsilon^{2}2i\omega\frac{\partial A}{\partial
T_{2}}-\omega^{2}A\right)e^{i\theta_{k}}$
$\displaystyle+\left(\epsilon^{2}\frac{\partial^{2}\bar{A}}{\partial
T_{1}^{2}}+\epsilon 2i\omega\frac{\partial\bar{A}}{\partial
T_{1}}+\epsilon^{2}2i\omega\frac{\partial\bar{A}}{\partial
T_{2}}-\omega^{2}\bar{A}\right)e^{-i\theta_{k}}$
$\displaystyle+\left(-\epsilon^{2}4i\omega\frac{\partial C}{\partial
T_{1}}-\epsilon 4\omega^{2}C\right)e^{2i\theta_{k}}$
$\displaystyle+\left(\epsilon^{2}4i\omega\frac{\partial\bar{C}}{\partial
T_{1}}-\epsilon 4\omega^{2}\bar{C}\right)e^{-2i\theta_{k}}+o(\epsilon^{3}).$

Similarly, the first order temporal derivatives of $n_{k}$ with respect to the
new space and time variables and amplitudes $F$, $G$, $H$ are given by

$\displaystyle\dot{n}_{k}$ $\displaystyle=\left(\epsilon\frac{\partial
F}{\partial T_{1}}+\epsilon^{2}\frac{\partial F}{\partial T_{2}}-i\omega
F\right)e^{i\theta_{k}}$
$\displaystyle+\left(\epsilon\frac{\partial\bar{F}}{\partial
T_{1}}+\epsilon^{2}\frac{\partial\bar{F}}{\partial
T_{2}}+i\omega\bar{F}\right)e^{-i\theta_{k}}$
$\displaystyle+\epsilon^{2}\frac{\partial G}{\partial
T_{1}}+\left(\epsilon^{2}\frac{\partial H}{\partial T_{1}}-\epsilon 2i\omega
H\right)e^{2i\theta_{k}}$
$\displaystyle+\left(\epsilon^{2}\frac{\partial\bar{H}}{\partial
T_{1}}+\epsilon 2i\omega\bar{H}\right)e^{-2i\theta_{k}}+o(\epsilon^{3}).$

Substituting the above equations in Eqs. (3) - (8) (see Appendix A) and
equating the coefficients of different orders in $\epsilon$, we obtain the
following important relations:

The dispersion relation of linear pulses by comparing the coefficient of
$e^{\pm i\theta_{k}}$ in Eq. (3) (see Eq. (A1)) is given by

$\omega^{2}=E_{0}+4D\sin^{2}\left(\frac{q}{2}\right).$ (13)

Figure 1: The effect of $q$ and $D$ on $\omega$ and $v_{g}$. (a) Dispersion
relation of linear pulses for $\mu=0.003$ and $D=0.01$. (b) The effect of the
diffusion coefficient $D$ on the velocity of the pulses for $\mu=0.003$,
$q=1.5$. At smaller diffusion coefficient, velocity $v_{g}$ changes very
rapidly, whereas for larger diffusion coefficient, the rate of change of
velocity is very low. The insets show a zoom-in of the behavior in the black
rectangles.

We have also shown the dispersion relation numerically in Fig. (1)(a). The
dispersion relation of the traveling pulses in the network depends on the
associated parameters of the coupled system.

Similarly, comparing the coefficient of $\epsilon e^{i\theta_{k}}$ in Eq. (3)
(see Eq. (A1)), we obtain

$\frac{\partial A}{\partial T_{1}}+v_{g}\frac{\partial A}{\partial X_{1}}=0,$
(14)

where $v_{g}$ is known as the group velocity, given by

$v_{g}=\frac{D\sin q}{\omega}.$ (15)

Clearly, the diffusion coefficient $D$ plays a major role in controlling the
velocity of the pulses. We have shown the effects of diffusive coupling on the
velocity of the pulses (Fig. 1(b)). Interestingly, for smaller diffusion
coefficient, the velocity $v_{g}$ changes very rapidly, whereas for higher
diffusion coefficient, the rate of change of the velocity is very low. We have
shown a zoomed version in Fig. 1(b). The dispersion equation is related to the
parameters associated with the excitable network. The wave velocity depends on
the diffusive property of the membrane dynamics, as the nerve impulses move
faster according to the ionic movements in the cell membrane.

Equating the coefficients without exponential terms and $\epsilon
e^{2i\theta_{k}}$ in Eq. (3), we obtain $B=0$ and $C=0$, respectively. As a
result of the coupling between the membrane voltage and slow variables, we get

$F=\frac{E_{0}(s+i\omega)}{s^{2}+\omega^{2}}A$

by equating the coefficients of $e^{i\theta_{k}}$ in Eq. (8) (see Eq. (A2)).

The terms depending on $\epsilon^{2}e^{i\theta_{k}}$ have the following
relation in Eq. (3) (see Eq. (A1))

$\displaystyle\frac{\partial^{2}A}{\partial T_{1}^{2}}-2i\omega\frac{\partial
A}{\partial T_{2}}$ $\displaystyle=i\omega E_{1}A$
$\displaystyle+(i\omega-1)E_{2}\mid
A\mid^{2}A-\frac{E_{0}E_{3}(s+iw)}{s^{2}+w^{2}}A$ $\displaystyle+4i\omega
D_{1}\sin^{2}(\frac{q}{2})A+2iD\sin(q)\frac{\partial A}{\partial X_{2}}$
$\displaystyle+D\cos(q)\frac{\partial^{2}A}{\partial X_{1}^{2}}.$ (16)

Using the transformation $u_{k}=X_{k}-v_{g}T_{k}$ and $\tau_{k}=T_{k}$, we
obtain

$i\frac{\partial A}{\partial\tau_{2}}+\frac{P}{2}\frac{\partial^{2}A}{\partial
u_{1}^{2}}+Q\mid A\mid^{2}A+i\frac{R}{2}A=0.$ (17)

Equation (17) is the modified CGLE, which describes the evolution of modulated
pulses in the neuronal network.

The real dispersion coefficient $P$ is given by

$P=\frac{D\omega^{2}\cos q-D^{2}\sin^{2}q}{\omega^{3}}.$

The two complex dissipation coefficients $Q$ and $R$ are given by

$\displaystyle Q$ $\displaystyle=-\frac{E_{2}}{2\omega}+i\frac{E_{2}}{2},$
$\displaystyle R$
$\displaystyle=E_{1}+4D_{1}\sin^{2}\frac{q}{2}-\frac{E_{3}E_{0}}{s^{2}+\omega^{2}}+i\left(\frac{sE_{3}E_{0}}{w(s^{2}+\omega^{2})}\right).$

Here $Q_{r}$, $R_{r}$ and $Q_{i}$, $R_{i}$ are the real and imaginary parts of
$Q$, $R$ respectively.

The sign of $PQ_{r}$ plays a major role in determining the modulational
instability as the dispersion coefficient is real. According to the Benjamin-
Feir instability [30], positive and negative values of $PQ_{r}$ indicate that
plane pulses are unstable and stable, respectively. Clearly, this stability
condition is independent of wave propagation. Thus, in the positive domain of
$PQ_{r}$, we can find the nerve impulses for any wave carrier. We have shown
the variations of the coefficient, $PQ_{r}$ in Fig. 2 with respect to the wave
vector, $q$.

Figure 2: Variations of the coefficient $PQ_{r}$ with respect to the wave
vector $q$ for $\mu=0.003$ and $D=0.01$.

##  4 Traveling wave profiles and numerical results

In the following, we obtain the solution of Eq. (17) to understand the
dynamics of the envelope solitons by considering the purely real dissipation
term $R$ ($R_{i}=0$). We consider the solution of Eq. (17) in the following
form

$A(u_{k},\tau_{2})=\frac{A_{0}e^{\alpha}}{1+e^{(\alpha+\bar{\alpha})^{(}1+i\beta)}},$
(18)

where $\alpha=qu_{k}-\omega\tau_{2}$, $\beta=\gamma\pm\sqrt{2+\gamma^{2}}$ and
$\gamma=\frac{3Q_{r}}{2Q_{i}}$. Substituting Eq. (18) into Eq. (17), we obtain

$\displaystyle A$ $\displaystyle=A_{0}\frac{e^{-\alpha}+\cos 2\alpha\beta
e^{\alpha}}{2(\cosh 2\alpha+\cos 2\alpha\beta)}$
$\displaystyle+i\left(-A_{0}\frac{\sin 2\alpha\beta e^{\alpha}}{2(\cosh
2\alpha+\cos 2\alpha\beta)}\right).$ (19)

From Eq. (9), we obtain

$m=2(A_{r}\cos\theta-A_{i}\sin\theta),$ (20)

where $A_{r}$ and $A_{i}$ indicate the real and imaginary parts of $A$.
Finally, using Eqs. (19), (20) and the expression $v_{k}=\epsilon m_{k}$, we
obtain

$v_{k}=\epsilon
A_{0}\frac{\cos(\theta_{k}-2\alpha_{k}\beta)e^{\alpha_{k}}+\cos(\theta_{k})e^{-\alpha_{k}}}{2(\cosh(2\alpha_{k})+\cos(2\alpha_{k}\beta))}.$
(21)

Figure 3: Spatiotemporal evolution of the traveling wave profiles for the
slow-fast FHR system (1), using Eq. (21). A series of traveling pulses is
shown. The multi-bump solution and disappearance of bumps with the generation
of single solitary wave profiles and effects of perturbation on the impulses
are presented for the parameters $\mu=0.003$, $q=1.5$, $D=0.01$, $t=10$,
$A_{0}=0.1$ in (a) for $\epsilon=0.05$, (b) $\epsilon=0.1$, (c)
$\epsilon=0.25$ and (d) $\epsilon=0.7$. Figure 4: Traveling wave profiles of
analytical envelope solutions to system (1), using Eq. (21). Effects of small
parameter values $\mu$ on nerve impulses. Here $q=1.5$, $D=0.01$, $t=10$,
$A_{0}=0.1$, $\epsilon=0.7$. Panel (a) is for $\mu=0.000001$, panel (b) for
$\mu=0.1$, panel (c) for $\mu=0.4$ and panel (d) for $\mu=0.8$. Figure 5:
Traveling wave profiles of spatial envelope solutions to system (1), using Eq.
(21). Effects of small diffusion coefficient $D$ on the wave profiles for the
parameters $q=1.5$, $\mu=0.003$, $t=10$, $A_{0}=0.1$, $\epsilon=0.25$ for (a)
$D=0.001$, (b) $D=5$ and (c) $D=10$. Figure 6: Numerical results for traveling
wave pulses at different times $t$ (shown in the insets). Multi-bump solutions
in the form of (a) symmetric wave profiles with almost-same amplitudes at
three different times and (b) asymmetric wave profiles with different
amplitudes at the same three different times. The different colors denote the
different times shown in the insets at which the traveling pulses are plotted.
In both panels, we have used the parameters $\mu=0.003$, $q=1.5$, $D=0.01$,
$A_{0}=0.1$, with $\epsilon=0.1$ in panel (a) and $\epsilon=0.7$ in panel (b).

In Fig. 3, we show the effects of the perturbation parameter, $\epsilon$, on
the nerve impulses for fixed values of the remaining parameters. The constants
$A_{0}$ and $q$ assume small values [29]. The changes of the parameters do not
affect the dynamics. The value of $q$ is considered in the region mentioned in
Fig. 1(a). At lower perturbation $\epsilon$, we observe the envelope of a
multi-bump solution for fixed diffusive coupling $D=0.01$ changing in time,
where each of the bumps is unstable and transient [26]. As the perturbation
$\epsilon$ increases, there is not enough activity to generate new bumps. As a
result, there is no wave propagation. Finally, at higher perturbation
$\epsilon$, the nerve impulse is localized (see Fig. 3(d)) in the coupled
network, which is an asymmetric envelope soliton. The parameter $\epsilon$
plays a major role not only in wave profiles but also in their amplitudes. It
also affects the structural pattern of the wave profiles as depending on its
value, it can transform the pattern into single soliton pulses.

The disappearance of bumps changes the propagation activity. The individual
bump solution in the multi-bump wave profiles are transient and non-
propagating. The amplitude of the impulse increases as the perturbation
$\epsilon$ increases, which indicates the movement of more ions across the
cell’s membrane, leading to the generation of high-amplitude action
potentials. This phenomenon occurs due to the fluctuations in the sequence of
neuronal firing times [23, 41]. However, the nature of the traveling impulses
changes in time when the value of the slow-modulation parameter $\mu$ in the
slow subsystem changes for fixed higher values of the perturbation parameter
$\epsilon$. This represents single-soliton pulses and the amplitudes of the
single pulses decrease with the change in the values of $\mu$, as shown in
Fig. 4. The value of the slow-modulation parameter must be small in the
coupled dynamics to obtain higher amplitude traveling pulses.

Next, we study the dynamics of traveling-wave profiles for $\epsilon=0.25$ and
fixed time by changing the coupling strength, $D$, shown in Fig. 5. They
exhibit multi-pulses and the disappearance of both-end ripples changes the
propagation activity of the wave profiles. We can observe how multi-pulse wave
profiles emerge with surrounding standing pulses which collapses to single-
traveling pulses. The solitary pulse corresponds to the envelope of a multiple
bump solution. We also study the changes of the nerve impulse with respect to
time $t$ in Fig. 6. The amplitude of the envelope of a multi-bump is
independent of time (see Fig. 6 (a)). The form and amplitude of the asymmetric
envelope-soliton change with time (see Fig. 6 (b)), therefore it is
structurally unstable. The corresponding pulses are highly nonlinear envelope-
solitons.

Next, we are interested in finding the solution of Eq. (17), if the
dissipation term has both real and imaginary parts (i.e., $R_{i}\neq 0$). Let
the solution of Eq. (17) be $A_{new}=Ae^{i\delta\tau_{2}}$. From Eq. (17), we
get $\delta=-\frac{R_{i}}{2}$. Hence the solution is given by

$A_{new}=Ae^{-\frac{iR_{i}}{2}\tau_{2}}.$ (22)

##  5 Instability of plane waves

In the previous section, we noticed that Eq. (17) has both envelope of a
multi-bump and soliton wave profiles. In the following, we will find under
which conditions, plane waves are stable or unstable for small perturbations.
Amplitude modulated pulses emerge due to the instability of plane waves. We
obtain a plane wave profile of the form

$A(u_{k},\tau_{2})=A_{0}e^{i(au_{k}-f\tau_{2})},$ (23)

where $A_{0}$ is the plane wave amplitude and $a$ and $f$ are the wave number
and angular frequency, respectively. Substituting Eq. (23) into Eq. (17), from
the real part, we obtain

$f=\frac{P}{2}a^{2}+\frac{R_{i}}{2}-Q_{r}A_{0}^{2}$ (24)

and from the complex part,

$Q_{i}A_{0}^{2}+\frac{R_{r}}{2}=0.$ (25)

Equation (24) is known as the dispersion relation of plane waves. It is clear
from Eq. (24) that the wave angular frequency depends not only on the wave
number but also on the wave amplitude.

To study the instability of the plane wave, we consider a solution of the form

$A(u_{k},\tau_{2})=(A_{0}+a_{1}(u_{k},\tau_{2}))e^{i(au_{k}-f\tau_{2}+a_{2}(u_{k},\tau_{2}))},$
(26)

where the perturbation amplitude $a_{1}(u_{k},\tau_{2})$ is very small with
respect to the plane wave amplitude $A_{0}$. Using Eqs. (17), (26) and
neglecting the nonlinear terms of the perturbations $a_{1}$ and $a_{2}$, we
derive the equations

$\displaystyle-
A_{0}{a_{2}}_{\tau_{2}}+\frac{P}{2}{a_{1}}_{uu}-PA_{0}a{a_{2}}_{u}+2Q_{r}A_{0}^{2}a_{1}$
$\displaystyle=$ $\displaystyle 0,$ (27)
$\displaystyle{a_{1}}_{\tau_{2}}+\frac{P}{2}A_{0}{a_{2}}_{uu}+Pa{a_{1}}_{u}-R_{r}a_{1}$
$\displaystyle=$ $\displaystyle 0.$ (28)

Equations (27) and (28) describe the evolution of the perturbation. We
consider the solutions to Eqs. (27) and (28) in the form

$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle
a_{10}e^{i(\zeta_{1}u-\zeta_{2}\tau_{2})}+c.c.,$ (29) $\displaystyle a_{2}$
$\displaystyle=$ $\displaystyle
a_{20}e^{i(\zeta_{1}u-\zeta_{2}\tau_{2})}+c.c.,$ (30)

where $\zeta_{1}$ and $\zeta_{2}$ are the wave number of the perturbation and
its corresponding propagation frequency respectively and $c.c.$ stands for
complex conjugate. Generally, the wave number is real and the propagation
frequency complex.

Substituting the solutions (29) and (30) in Eqs. (27) and (28), we derive the
following linear homogeneous equations for $a_{10}$ and $a_{20}$

$\displaystyle\left(2Q_{r}A_{0}^{2}-\frac{P}{2}\zeta_{1}^{2}\right)a_{10}+iA_{0}(\zeta_{2}-Pa\zeta_{1})a_{20}$
$\displaystyle=$ $\displaystyle 0,$ (31)
$\displaystyle(-R_{r}-i(\zeta_{2}-Pa\zeta_{1}))a_{10}-\frac{P}{2}A_{0}\zeta_{1}^{2}a_{20}$
$\displaystyle=$ $\displaystyle 0,$ (32)

which can be written in the matrix form

$LL_{1}=0,$

where
$L=\left({\begin{array}[]{*{20}{c}}2Q_{r}A_{0}^{2}-\frac{P}{2}\zeta_{1}^{2}&iA_{0}(\zeta_{2}-Pa\zeta_{1})\\\
-R_{r}-i(\zeta_{2}-Pa\zeta_{1})&-\frac{P}{2}A_{0}\zeta_{1}^{2}\end{array}}\right)$
and $L_{1}=\left({\begin{array}[]{*{20}{c}}{a_{10}}\\\
a_{20}\end{array}}\right)$.

Equations (31) and (32) have the nontrivial solutions, if $L$ is a singular
matrix, i.e. if $\det(L)=0$, from which it results that

$(\zeta_{2}-Pa\zeta_{1})^{2}=\frac{P^{2}\zeta_{1}^{2}}{4}\left(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2}\right)+iR_{r}(\zeta_{2}-Pa\zeta_{1}).$
(33)

If $y=\zeta_{2}-Pa\zeta_{1}$, Eq. (33) becomes

$y^{2}-iR_{r}y-\frac{P^{2}\zeta_{1}^{2}}{4}\left(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2}\right)=0.$
(34)

Equation (33) is known as the dispersion relation of the perturbation.
Furthermore, for a fixed value of the wave number $\zeta_{1}$, the term
$\frac{Q_{r}}{P}$ plays a major role in controlling the dynamics of the
angular frequency $\zeta_{2}$. Solving Eq. (34), we obtain

$y=i\frac{R_{r}}{2}\pm\frac{\sqrt{P^{2}\zeta_{1}^{2}\left(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2}\right)-R_{r}^{2}}}{2}.$
(35)

Depending on the sign of the discriminant, we examine the following three
cases:

  * Case 1:

$P^{2}\zeta_{1}^{2}(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2})-R_{r}^{2}=0$,
i.e., $y=i\frac{R_{r}}{2}$.
The imaginary part of $\zeta_{2}$ is the same as the imaginary part of $y$,
which is given by $\frac{R_{r}}{2}$. From Fig. 7, it is clear that $R_{r}<0$,
which indicates the wave oscillates about its original value and after a
certain period of time the perturbations die out. Thus, the plane wave is
stable.

Figure 7: Effects of the wave vector $q$ on the real dissipation coefficient
$R_{r}$ for $\mu=0.003$ and $D=0.01$.

  * Case 2:

$P^{2}\zeta_{1}^{2}(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2})-R_{r}^{2}>0$,
i.e., $\Im(y)=i\frac{R_{r}}{2}$.
Similarly, the plane wave is also stable.

  * Case 3:

$P^{2}\zeta_{1}^{2}(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2})-R_{r}^{2}<0$,
i.e., $y=i\frac{R_{r}}{2}\pm
i\frac{\sqrt{R_{r}^{2}-P^{2}\zeta_{1}^{2}(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2})}}{2}.$
As we already know the plane wave is stable if $\Im(y)<0$, which gives

$P^{2}\zeta_{1}^{2}\left(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2}\right)>0.$

It is clear from the above expression that the plane wave solution is stable
if $P$ and $Q_{r}$ have opposite signs, i.e., if $PQ_{r}<0$. Similarly, we can
show that the plane wave solution is unstable, i.e., the perturbations grow
exponentially in time, if $PQ_{r}>0$.
In the case of modulational instability, the local growth rate or gain is
given by

$\displaystyle h$ $\displaystyle=\mid\Im y\mid=\mid\Im\zeta_{2}\mid$
$\displaystyle=\frac{1}{2}\Big{(}R_{r}+\sqrt{R_{r}^{2}-P^{2}\zeta_{1}^{2}(\zeta_{1}^{2}-\frac{4Q_{r}}{P}A_{0}^{2})}\Big{)}.$

The gain is maximum, i.e., the plane wave modulates itself, if the wave number
$\zeta_{1}=A_{0}\sqrt{\frac{Q_{r}}{P}}$ and the corresponding gain is given by

$\displaystyle
h=\frac{1}{2}\Big{(}R_{r}+\sqrt{R_{r}^{2}+3Q_{r}^{2}A_{0}^{4}}\Big{)}.$ (36)

It is clear from Eq. (36) that the amplitude of the plane wave and the
dissipation coefficient of CGLE (Eq. (17)) play an important role in
controlling the growth rate of modulational instability.

##  6 Conclusions and Discussion

In this paper, we studied nonlinear excitations in a locally diffusive network
of excitable slow-fast FitzHugh-Rinzel neurons. We showed that the dimension
of the model can be reduced by using perturbation theory and that it can
transform to an extended complex Ginzburg-Landau equation. The solitary
traveling wave profiles of the voltage dynamics are also discussed. Our
results show that the activities of the coupled neurons exhibit collective
dynamics that cannot be observed in the activity of single neurons.
Furthermore, we analyzed the propagation of wave profiles in an excitatory,
diffusively coupled network of slow-fast neurons, treating it as a spatially
extended medium. The model exhibits different neuro-computational features. We
showed how small amplitudes and short wavelengths affect the homogeneous
network and how they lead to traveling pulses. The wave speed was estimated
using perturbation theory. The stable traveling pulses exist for sufficiently
small $\epsilon$, however, the pulse is the envelope of the multiple traveling
wave profiles in which the multiple damping terms are transient and unstable.
The traveling wave profiles emerge via the appearance of multiple small-
amplitude oscillations at the beginning of the wave profile.

Our analytical treatment describes envelope solitary and multiple pulses [26]
observed in numerical simulations. Successive transitions from single pulses
to multi-bump solutions are described. The modulated traveling pulses and
their properties are observed and are influenced by the diffusive coupling and
the effects of small perturbations. Our analysis shows the appearance of
multiple small waves or series of pulses in the excitable network, that are
interesting to explore further as it may be possible to measure different
properties of traveling pulses in excitable media. The signature of these
types of time-dependent traveling-wave profiles can be effectively measured by
the variations in the parameters of certain slow-fast excitable systems [29,
30]. Neurons can connect electrically with other nearby neurons. This is a
type of weak, diffusive, coupling and is reminiscent of bursting-like activity
in $\beta$-cell islets in the pancreas that secrete insulin [28, 39]. From a
biological viewpoint, this shows that neurons “communicate” with each other
and participate in the collective processing of information exchange, a
portion of which is shared among them. As a result, the brain can effectively
process information not only in the temporal but also in the spatial domain
[27, 37]. Neuronal ensembles can process stimuli in different or similar ways.
The waves (action potentials) exhibit nonlinear-envelope solitons having an
up-and-down asymmetry in their amplitudes, as we also show in this paper. The
amplitude of the impulse increases as the perturbation $\epsilon$ increases,
which indicates the movement of more ions across the cell’s membrane, leading
to the generation of high-amplitude action potentials. This is caused by
fluctuations in the sequence of neuronal firing times [22, 25, 40]. In
extended spatial systems of neuronal populations, a transmembrane potential
difference can travel across neurons by means of traveling wave propagation.
Hence if one would want to measure propagation of nerve impulses in a neural
network, one would need to measure the transmembrane potential difference to
observe the phenomena discussed in this paper. Our work is theoretical and
supported by numerical results that confirm the analytical findings. Hence,
our findings are qualitative and important as they pave the way for
experimentalists to look out for them in experiments. Theoretical and
numerical results are important in understanding functional mechanisms in
neural computations [4, 5, 6]. A key challenge is to analyze the intrinsic
dynamics of neuronal activities and its characteristics [19, 20, 21].

## Appendix A

Equations (3) and (8) can be written in terms of the amplitudes $A$, $B$, $C$,
$F$, $G$ and $H$ in the forms

$\begin{array}[]{l}\left(\epsilon^{2}\frac{\partial^{2}A}{\partial
T_{1}^{2}}-\epsilon 2i\omega\frac{\partial A}{\partial
T_{1}}-\epsilon^{2}2i\omega\frac{\partial A}{\partial
T_{2}}-\omega^{2}A\right)e^{i\theta_{k}}+\Big{(}\epsilon^{2}\frac{\partial^{2}\bar{A}}{\partial
T_{1}^{2}}+\\\ \vspace{0.2cm}\epsilon 2i\omega\frac{\partial\bar{A}}{\partial
T_{1}}+\epsilon^{2}2i\omega\frac{\partial\bar{A}}{\partial
T_{2}}-\omega^{2}\bar{A}\Big{)}e^{-i\theta_{k}}+\left(-\epsilon^{2}4i\omega\frac{\partial
C}{\partial T_{1}}-\epsilon 4\omega^{2}C\right)\\\
e^{2i\theta_{k}}+\left(\epsilon^{2}4i\omega\frac{\partial\bar{C}}{\partial
T_{1}}-\epsilon
4\omega^{2}\bar{C}\right)e^{-2i\theta_{k}}+E_{0}\Big{(}Ae^{i\theta_{k}}+\overline{A}e^{-i\theta_{k}}+\\\
\vspace{0.2cm}\epsilon(B+Ce^{2i\theta_{k}}+\overline{C}e^{-2i\theta_{k}})\Big{)}+\epsilon^{2}\Big{(}E_{1}+E_{2}(A^{2}e^{2i\theta_{k}}+\bar{A}^{2}e^{-2i\theta_{k}}\\\
\vspace{0.2cm}+2A\bar{A})\Big{)}\left(-i\omega
Ae^{i\theta_{k}}+i\omega\bar{A}e^{-i\theta_{k}}\right)+\epsilon^{2}\frac{E_{2}}{3}\Big{(}A^{3}e^{3i\theta_{k}}+\bar{A}^{3}e^{-3i\theta_{k}}\\\
\vspace{0.2cm}+3A^{2}\bar{A}e^{i\theta_{k}}+3A\bar{A}^{2}e^{-i\theta_{k}}\Big{)}+\epsilon^{2}E_{3}(Fe^{i\theta_{k}}+\bar{F}e^{-i\theta_{k}})+o(\epsilon^{3})=\\\
\vspace{0.2cm}D\Big{(}(A+\epsilon\frac{\partial A}{\partial
X_{1}}+\epsilon^{2}\frac{\partial A}{\partial
X_{2}}+\frac{\epsilon^{2}}{2}\frac{\partial^{2}A}{\partial
X_{1}^{2}})e^{iq}e^{i\theta_{k}}+(\bar{A}+\epsilon\frac{\partial\bar{A}}{\partial
X_{1}}+\epsilon^{2}\frac{\partial\bar{A}}{\partial X_{2}}+\\\
\vspace{0.2cm}\frac{\epsilon^{2}}{2}\frac{\partial^{2}\bar{A}}{\partial
X_{1}^{2}})e^{-iq}e^{-i\theta_{k}}+\epsilon(B+\epsilon\frac{\partial
B}{\partial X_{1}}+(C+\epsilon\frac{\partial C}{\partial
X_{1}})e^{2iq}e^{2i\theta_{k}}+\\\
\vspace{0.2cm}(\bar{C}+\epsilon\frac{\partial\bar{C}}{\partial
X_{1}})e^{-2iq}e^{-2i\theta_{k}})-2(Ae^{i\theta_{k}}+\overline{A}e^{-i\theta_{k}}+\epsilon(B+Ce^{2i\theta_{k}}+\\\
\vspace{0.2cm}\overline{C}e^{-2i\theta_{k}}))+(A-\epsilon\frac{\partial
A}{\partial X_{1}}-\epsilon^{2}\frac{\partial A}{\partial
X_{2}}+\frac{\epsilon^{2}}{2}\frac{\partial^{2}A}{\partial
X_{1}^{2}})e^{-iq}e^{i\theta_{k}}+\\\
\vspace{0.2cm}(\bar{A}-\epsilon\frac{\partial\bar{A}}{\partial
X_{1}}-\epsilon^{2}\frac{\partial\bar{A}}{\partial
X_{2}}+\frac{\epsilon^{2}}{2}\frac{\partial^{2}\bar{A}}{\partial
X_{1}^{2}})e^{iq}e^{-i\theta_{k}}+\epsilon(B-\epsilon\frac{\partial
B}{\partial X_{1}}+\\\ \vspace{0.2cm}(C-\epsilon\frac{\partial C}{\partial
X_{1}})e^{-2iq}e^{2i\theta_{k}}+(\bar{C}-\epsilon\frac{\partial\bar{C}}{\partial
X_{1}})e^{2iq}e^{-2i\theta_{k}})\Big{)}+\epsilon^{2}D_{1}\\\
\vspace{0.2cm}\Big{(}-i\omega
Ae^{iq}e^{i\theta_{k}}+i\omega\bar{A}e^{-iq}e^{-i\theta_{k}}-2(-i\omega
Ae^{i\theta_{k}}+i\omega\bar{A}e^{-i\theta_{k}})-\\\ i\omega
Ae^{-iq}e^{i\theta_{k}}+i\omega\bar{A}e^{iq}e^{-i\theta_{k}}\Big{)}+o(\epsilon^{3})\hskip
71.13188pt{\rm(A1)}\end{array}$

and

$\begin{array}[]{l}\left(\epsilon\frac{\partial F}{\partial
T_{1}}+\epsilon^{2}\frac{\partial F}{\partial T_{2}}-i\omega
F\right)e^{i\theta_{k}}+\left(\epsilon\frac{\partial\bar{F}}{\partial
T_{1}}+\epsilon^{2}\frac{\partial\bar{F}}{\partial
T_{2}}+i\omega\bar{F}\right)e^{-i\theta_{k}}\\\
\vspace{0.2cm}+\epsilon^{2}\frac{\partial G}{\partial
T_{1}}+\left(\epsilon^{2}\frac{\partial H}{\partial T_{1}}-\epsilon 2i\omega
H\right)e^{2i\theta_{k}}+\left(\epsilon^{2}\frac{\partial\bar{H}}{\partial
T_{1}}+\epsilon 2i\omega\bar{H}\right)e^{-2i\theta_{k}}\\\
\vspace{0.2cm}+s\Big{(}Fe^{i\theta_{k}}+\overline{F}e^{-i\theta_{k}}+\epsilon(G+He^{2i\theta_{k}}+\overline{H}e^{-2i\theta_{k}})\Big{)}-\\\
E_{0}\Big{(}Ae^{i\theta_{k}}+\overline{A}e^{-i\theta_{k}}+\epsilon(B+Ce^{2i\theta_{k}}+\overline{C}e^{-2i\theta_{k}})\Big{)}+\\\
o(\epsilon^{3})=0.\hskip 179.25244pt{\rm(A2)}\end{array}$

#### Conflict of Interest

The authors have no conflicts to disclose.

#### Author Contributions

Arnab Mondal and Argha Mondal contributed equally to this work.

#### Data Availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## References

  * [1] A. Mondal, R. K. Upadhyay, A. Mondal, S. K. Sharma, Dynamics of a modified excitable neuron model: Diffusive instabilities and traveling wave solutions, Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (11) (2018) 113104.
  * [2] E. M. Izhikevich, Synchronization of elliptic bursters, SIAM Review 43 (2) (2001) 315–344.
  * [3] E. M. Izhikevich, Dynamical systems in neuroscience, MIT press, 2007.
  * [4] J. Ma, J. Tang, A review for dynamics in neuron and neuronal network, Nonlinear Dynamics 89 (3) (2017) 1569–1578.
  * [5] N. Iqbal, R. Wu, B. Liu, Pattern formation by super-diffusion in FitzHugh-Nagumo model, Applied Mathematics and Computation 313 (2017) 245–258.
  * [6] P. Muruganandam, M. Lakshmanan, Bifurcation analysis of the travelling waveform of FitzHugh-Nagumo nerve conduction model equation, Chaos: An Interdisciplinary Journal of Nonlinear Science 7 (3) (1997) 476–487.
  * [7] N. Corson, M. A. Aziz-Alaoui, Asymptotic dynamics of Hindmarsh-Rose neuronal system, Dynamics of Continuous, Discrete and Impulsive Systems, Series B: Applications and Algorithms (16) (2009) 535.
  * [8] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophysical Journal 1 (6) (1961) 445–466.
  * [9] J. Rinzel, Y. S. Lee, Nonlinear Oscillations in Biology and Chemistry (1986).
  * [10] J. Rinzel, Mathematical topics in population biology, morphogenesis, and neurosciences, Lecture Notes in Biomathematics 71 (1987) 267–281.
  * [11] R. G. Townsend, P. Gong, Detection and analysis of spatiotemporal patterns in brain activity, PLoS Computational Biology 14 (12) (2018) e1006643.
  * [12] J. G. Milton, P. H. Chu, J. D. Cowan, Spiral waves in integrate-and-fire neural networks, Advances in neural information processing systems (1993) 1001–1001.
  * [13] S. Kondo, T. Miura, Reaction-diffusion model as a framework for understanding biological pattern formation, Science 329 (5999) (2010) 1616–1620.
  * [14] S. R. Meier, J. L. Lancaster, J. M. Starobin, Bursting regimes in a reaction-diffusion system with action potential-dependent equilibrium, PloS One 10 (3) (2015) e0122401.
  * [15] P. Feng, J. Zhang, W. Wang, Spike-like solitary waves in incompressible boundary layers driven by a travelling wave, Chaos: An Interdisciplinary Journal of Nonlinear Science 26 (6) (2016) 063104.
  * [16] I. Belykh, E. De Lange, M. Hasler, Synchronization of bursting neurons: What matters in the network topology, Physical Review Letters 94 (18) (2005) 188101\.
  * [17] B. W. Connors, Y. Amitai, Generation of epileptiform discharges by local circuits in neocortex, Epilepsy: Models, Mechanisms and Concepts (1993) 388–424.
  * [18] J. W. Lance, Current concepts of migraine pathogenesis, Neurology 43 (6 Suppl 3) (1993) S11–5.
  * [19] B. W. Connors, M. J. Gutnick, Intrinsic firing patterns of diverse neocortical neurons, Trends in Neurosciences 13 (3) (1990) 99–104.
  * [20] E. M. Izhikevich, Which model to use for cortical spiking neurons?, IEEE Transactions on Neural Networks 15 (5) (2004) 1063–1070.
  * [21] R. Yafia, M. A. Aziz-Alaoui, Existence of periodic travelling waves solutions in predator prey model with diffusion, Applied Mathematical Modelling 37 (6) (2013) 3635–3644.
  * [22] G. S. Wig, B. L. Schlaggar, S. E. Petersen, Concepts and principles in the analysis of brain networks, Annals of the New York Academy of Sciences 1224 (1) (2011) 126–146.
  * [23] S. E. Folias, P. C. Bressloff, Stimulus-locked traveling waves and breathers in an excitatory neural network, SIAM Journal on Applied Mathematics 65 (6) (2005) 2067–2092.
  * [24] Z. P. Kilpatrick, P. C. Bressloff, Effects of synaptic depression and adaptation on spatiotemporal dynamics of an excitatory neuronal network, Physica D: Nonlinear Phenomena 239 (9) (2010) 547–560.
  * [25] D. J. Pinto, G. B. Ermentrout, Spatially structured activity in synaptically coupled neuronal networks: I. Traveling fronts and pulses, SIAM Journal on Applied Mathematics 62 (1) (2001) 206–225.
  * [26] Z. P. Kilpatrick, S. E. Folias, P. C. Bressloff, Traveling pulses and wave propagation failure in inhomogeneous neural media, SIAM Journal on Applied Dynamical Systems 7 (1) (2008) 161–185.
  * [27] E. A. Ermakova, E. E. Shnol, M. A. Panteleev, A. A. Butylin, V. Volpert, F. I. Ataullakhanov, On propagation of excitation waves in moving media: The FitzHugh-Nagumo model, PloS One 4 (2) (2009) e4454.
  * [28] J. A. Villacorta-Atienza, V. A. Makarov, Wave-processing of long-scale information by neuronal chains, PloS One 8 (2) (2013) e57440.
  * [29] F. M. Kakmeni Moukam, E. M. Inack, E. M. Yamakou, Localized nonlinear excitations in diffusive Hindmarsh-Rose neural networks, Physical Review E 89 (5) (2014) 052919.
  * [30] A. S. Tankou Tagne, C. N. Takembo, H. G. Ben-Bolie, P. Owona Ateba, Localized nonlinear excitations in diffusive memristor-based neuronal networks, PloS One 14 (6) (2019) e0214989.
  * [31] W. Malfliet, W. Hereman, The tanh method: I. Exact solutions of nonlinear evolution and wave equations, Physica Scripta 54 (6) (1996) 563.
  * [32] A.-M. Wazwaz, The tanh method: Exact solutions of the sine-Gordon and the sinh-Gordon equations, Applied Mathematics and Computation 167 (2) (2005) 1196–1210.
  * [33] R. A. Stefanescu, V. K. Jirsa, A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons, PLoS Computational Biology 4 (11) (2008) e1000219.
  * [34] J. Rinzel, W. C. Troy, Bursting phenomena in a simplified Oregonator flow system model, The Journal of Chemical Physics 76 (4) (1982) 1775–1789.
  * [35] J. Wojcik, A. Shilnikov, Voltage interval mappings for an elliptic bursting model, in: Nonlinear Dynamics New Directions, Springer International Publishing Switzerland, 2015, pp. 195–213.
  * [36] A. Mondal, S. K. Sharma, R. K. Upadhyay, A. Mondal, Firing activities of a fractional-order FitzHugh-Rinzel bursting neuron model and its coupled dynamics, Scientific Reports 9 (1) (2019) 1–11.
  * [37] C. A. Del Negro, C. F. Hsiao, S. H. Chandler, A. Garfinkel, Evidence for a novel bursting mechanism in rodent trigeminal neurons, Biophysical Journal 75 (1) (1998) 174–182.
  * [38] A. Mondal, A. Mondal, S. Kumar Sharma, R. Kumar Upadhyay, C. G. Antonopoulos, Spatiotemporal characteristics in systems of diffusively coupled excitable slow-fast FitzHugh-Rinzel dynamical neurons, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (10) (2021) 103122.
  * [39] M. Perez-Armendariz, C. Roy, D. C. Spray, M. V. Bennett, Biophysical properties of gap junctions between freshly dispersed pairs of mouse pancreatic beta cells, Biophysical Journal 59 (1) (1991) 76–92.
  * [40] S. Raghavachari, J. A. Glazier, Waves in diffusively coupled bursting cells, Physical Review Letters 82 (14) (1999) 2991.
  * [41] S. E. Folias, P. C. Bressloff, Breathing pulses in an excitatory neural network, SIAM Journal on Applied Dynamical Systems 3 (3) (2004) 378–407.