# Spatiotemporal instabilities and pattern formation in systems of diffusively
coupled Izhikevich neurons

Argha Mondal${}^{1,2,^{*}}$, Chittaranjan Hens3, Arnab Mondal4, Chris G.
Antonopoulos2 1 School of Engineering, Amrita Vishwa Vidyapeetham, Amritapuri,
Kollam 690525, India
2 Department of Mathematical Sciences, University of Essex, Wivenhoe Park, UK
3 Physics and Applied Mathematics Unit, Indian Statistical Institute Kolkata,
WB, India
4 Department of Mathematics and Computing, Indian Institute of Technology
(Indian School of Mines), Dhanbad 826004, India

###### Abstract

Neurons are often connected, spatially and temporally, in phenomenal ways that
promote wave propagation. Therefore, it is essential to analyze the emergent
spatiotemporal patterns to understand the working mechanism of brain activity,
especially in cortical areas. Here, we present an explicit mathematical
analysis, corroborated by numerical results, to identify and investigate the
spatiotemporal, non-uniform, patterns that emerge due to instability in an
extended homogeneous 2D spatial domain, using the excitable Izhikevich neuron
model. We examine diffusive instability and perform bifurcation and fixed-
point analyses to characterize the patterns and their stability. Then, we
derive analytically the amplitude equations that establish the activities of
reaction-diffusion structures. We report on the emergence of diverse spatial
structures including hexagonal and mixed-type patterns by providing a
systematic mathematical approach, including variations in correlated
oscillations, pattern variations and amplitude fluctuations. Our work shows
that the emergence of spatiotemporal behavior, commonly found in excitable
systems, has the potential to contribute significantly to the study of
diffusively-coupled biophysical systems at large.

###### keywords:

Izhikevich model, bifurcation analysis, diffusive instabilities, amplitude
equations, spatiotemporal patterns

††journal: Journal of LaTeX Templates††footnotetext: Corresponding author
Email address<EMAIL_ADDRESS>
##  1 Introduction

The spatially structured patterns associated to wave propagation in coupled
neural populations can emerge due to high-low amplitude oscillations and
variations in the firing activities of single neurons in different network
topologies. One of the key challenges is to analyze the network dynamics and
characteristics [1, 2, 3, 4, 5]. The emergence of diverse spatial patterns is
often found to be associated with self-organization and competition within
neurons in a network. This may reflect the dynamics of signal processing and
functional connectivity in different brain areas [6, 7, 8]. Thus, it is
important to derive the spatial characteristics of the electrical activity,
investigate quantitatively and classify the various patterns observed [9, 10,
11, 12]. During wave propagation, single neurons are slightly depolarized,
however the probability of firing action potentials or generating high-
amplitude oscillations with spikes often increases depending on the synaptic
coupling strengths. These membrane currents depolarize the cell which produces
a voltage signal, composed of spikes. When a large-size network becomes
coherently depolarized during wave propagation, the local transmission between
interneurons in the network may be increased. As a result, various
spatiotemporal patterns such as spirals, traveling waves and hexagonal, mixed-
like patterns can be observed at different levels of activity [10, 13, 14, 15,
16]. Spatial patterns can also occur in Fitzhugh-Nagumo systems with spatial
interactions [17, 18, 19, 20, 21]. Wave propagation of neural activity is
often observed when a population of cortical cells is working during sensory
or motor event processing [7].

Here, we use a systematic, mathematical analysis based on amplitude equations
[22, 23] to investigate the formation of coherent structures of patterns. In
addition, rigorous numerical simulations are performed to identify the
diversity in wavy patterns. These patterns emerge due to instabilities in a
homogeneous medium (i.e., in a 2D spatial domain) in a reaction-diffusion
system of coupled Izhikevich neurons.

In our work, we present a comprehensive analysis to study the onset of pattern
formation and the dominant effects. Our work reveals how the analytical method
can be used to identify and analyze the spatiotemporal characteristics that
can help us in understanding the working mechanism of the collective dynamics.
The existence and stability of the spatial patterns are explored analytically
with the amplitude equations near the Turing bifurcation points. In [15], the
authors studied the spatiotemporal dynamics of pattern formation in a coupled
2D Morris-Lecar neural system and investigated the nonlinear responses of an
excitable, conductance-based, neuronal cable. They validated their numerical
results showing they are in good agreement with their analytical results using
amplitude equations and multiple-scale analysis, the same approach we adopt
here. Furthermore, they studied the formation of spatiotemporal patterns in
non-Turing regimes using the 2D Morris-Lecar model, whereas here we use
diffusively coupled Izhikevich neurons to study the formation of
spatiotemporal patterns in Turing regimes, following the same methodology.
This method presents a mathematical framework in the extended coupled
spatiotemporal systems to establish variety of structurally different patterns
with the modulation and stability analysis that are often found in excitable
coupled biophysical systems. Thus, our present work and the work in [15], show
clearly that the analytical approach used in our study here can be generalized
to study analytically the spatiotemporal instabilities and pattern formation
in other biophysically excitable models, including in Turing and non-Turing
regimes.

Particularly, in our work we explore the fluctuations and local nonlinear
excitations of diffusively coupled Izhikevich neurons [24, 25, 26] given by a
system of partial differential equations (PDEs). Performing a detailed
bifurcation analysis, we investigate the quiescent and oscillatory regimes and
explore the influence of coupling strengths on the system, which undergoes a
diffusion-driven instability, known as Turing instability [27, 28, 29, 30,
31]. We study how bursting regimes of a single excitable system can emerge in
a diffusively coupled system. The patterns associated with bursting range from
regular stable hexagons to distorted hexagons [21, 32, 33, 34]. To study the
modulation and stability of wave patterns (i.e., hexagons), we consider the
phase and amplitude of the oscillatory dynamics. We show that around the
Turing bifurcation points of the dynamics, spatiotemporal patterns emerge,
that are commonly found in extended biophysical systems [15, 35]. Thus, our
work has the potential to contribute significantly to the study of
diffusively-coupled biophysical systems at large.

##  2 The excitable single-neuron Izhikevich model

This work focuses on the complex dynamics of the Izhikevich model that
describes the electrical activity of a single neuron for a range of parameter
values. The model exhibits spike generation and a discontinuous resetting
process following the spikes. Particularly, its time evolution is described by
the set of ordinary differential equations [24, 25]

$\begin{array}[]{l}\frac{dv}{dt}=0.04v^{2}+5v+140-u+I=f_{1}(v,u),\\\
\frac{du}{dt}=a(bv-u)=f_{2}(v,u),\end{array}$ (1)

where $v$ is the membrane voltage and $u$ the recovery variable that measures
the activation of K+ and inactivation of Na+ ionic currents. When the voltage,
$v$ reaches its peak value, the following relation is applied after the spike
resetting constraint: if $v\geq{v_{peak}}=30$, then $v\leftarrow c$ and
$u\leftarrow u+d$. The resting potential depends on the parameter $b$ that
indicates the sensitivity of $u$ to the subthreshold fluctuations of the
voltage and, $a$ measures the timescale of the recovery variable, $u$. The
parameters $c$ and $d$ control the after-spike reset values of $v$ and $u$,
respectively. The function $(0.04v^{2}+5v+140)$ was derived using the spike
initiation dynamics of a cortical cell. The initial conditions are set to
$v=-63$ and $u=bv$, respectively [24, 25] and the stimulus currents are
injected to the neuron with the variable, $I$.

To study the behavior of the Izhikevich model (1) in the absense of diffusion
for different parameter sets and $I$ values, we perturb it around the fixed
point $E=(v_{0},u_{0})$ and compute the corresponding Jacobian matrix, $J$.
The equilibrium point $E$ can be deduced by solving the set of equations
$f_{1}(v_{0},u_{0})=0$ and $f_{2}(v_{0},u_{0})=0$ simultaneously for $v_{0}$
and $u_{0}$, which gives
$v_{0}=\frac{-(5-b)\pm\sqrt{(5-b)^{2}-0.16(140+I)}}{0.08}$ and $u_{0}=bv_{0}$.
The Jacobian of system (1) computed at $E=(v_{0},u_{0})$ is given by
$J=\left({\begin{array}[]{*{20}{c}}{{a_{11}}}&{{a_{12}}}\\\
{{a_{21}}}&{{a_{22}}}\end{array}}\right)$, where $a_{11}=0.08v_{0}+5$,
$a_{12}=-1$, $a_{21}=ab$ and $a_{22}=-a$. The fixed point, $E$ is locally
asymptotically stable if, trace$(J)=a_{11}+a_{22}<0$ and
$\det(J)=a_{11}a_{22}-a_{12}a_{21}>0$.

In the following, we consider specific types of biophysically plausible
parameter sets that give rise to Hopf bifurcations by varying $I$. This
amounts to finding spatial patterns using the method of amplitude equations
[22, 23] and identifying the onset of Hopf bifurcations. The parameter sets we
consider are the following: set I with $a=0.2$, $b=2$, $c=-56$, $d=-16$,
$I\in[-105.1,-103]$, set II with $a=-0.02$, $b=-1$, $c=-60$, $d=8$,
$I\in[78,80]$ and set III with $a=-0.026$, $b=-1$, $c=-45$, $d=0$,
$I\in[78,80]$. System (1) produces firings for these parameter sets with
varying $I$, ranging from quiescent to oscillatory behavior through chaotic
spiking (Fig. 1 (b), (c)), phasic to tonic spiking (Fig. 1 (d)-(f)) and phasic
bursting to regular bursting shown in Fig. 1 (g)-(i) [24, 25]. The numerical
simulations were performed using the fourth-order Runge-Kutta method with time
step $\delta t=0.01$. We have also checked that the results with a smaller
time step do not show any significant differences.

###  2.1 Bifurcation analysis

The bifurcation analysis of system (1) was performed using MatCont [36] in
Matlab, varying the bifurcation parameter $I$. At lower current stimulus
$I<-83.75$, the system has two fixed points which collide at $I=-83.75$ (see
Fig. 2(a), (b), where SN stands for saddle-node bifurcation) and vanish for
$I>-83.75$. The system changes its stability and becomes unstable at $I=-104$,
where an unstable limit cycle emerges, shown in Fig. 2(a), denoted by the
point HB for Hopf bifurcation. Particularly, a subcritical Hopf bifurcation
(HB) appears at $I=-104$. In Fig. 2(a), thick and dashed blue lines represent
the quiescent and oscillatory regions, respectively, and the red line, an
unstable limit cycle.

Figure 1: Temporal plots of the membrane voltage $v$ in system (1), where
panels (a)-(c) are for $I=-105.1$, $-103.9,-103$ (set I), (d)-(f) for $I=80$,
78.8, 78 (set II) and (g)-(i) for $I=80$, 78.9, 78 (set III).

In the following, we verify our numerical results performing an analytical
study. Particularly, the complex eigenvalues of $J$ computed at
$E=(v_{0},u_{0})$ are given by
$\gamma(I),\,\,\overline{\gamma(I)}=\alpha(I)\pm i\beta(I)$. Suppose that, for
a certain value $I=I_{0}$, the following conditions are met: (i)
$\alpha(I_{0})=0$, (ii) $\beta(I_{0})=\beta^{*}\neq 0$ and (iii)
${\left.{\frac{{d\alpha\left(I\right)}}{{dI}}}\right|_{I=I_{0}}}=\alpha^{*}\neq
0$. Then, system (1) undergoes a Hopf bifurcation at $I=I_{0}$. Conditions (i)
and (ii) are known as non-hyperbolicity conditions and condition (iii) as the
transversality condition. To verify analytically the Hopf bifurcation varying
$I$ for parameters in set I, we obtain $\alpha(I)=0.04(60+v_{0})$ and
$\beta(I)=0.04\sqrt{-v_{0}^{2}-130v_{0}-3975}$, where
$v_{0}=2.5\big{(}-15\pm\sqrt{-335-4I}\big{)}$. To show the existence of a Hopf
bifurcation, we find the value of $I=I_{0}$ at which conditions (i) - (iii)
are satisfied. Thus, solving the equation $\alpha(I)=0$ for $I$, we obtain
$I=I_{0}=-104$. Then, from condition (i) we have $\alpha(-104)=0$, from (ii)
$\beta\left({-104}\right)=0.6\neq 0$ and from (iii)
${\left.{\frac{{d\alpha\left(I\right)}}{{dI}}}\right|_{I=-104}}=0.0222\neq 0$.
Thus, we conclude that system (1) undergoes a Hopf bifurcation at $I=-104$,
which is in good agreement with the numerical simulations. The system
undergoes a super critical Hopf bifurcation at $I=79.070775$ and a saddle-node
bifurcation for parameters in set III (see Fig. 2(c), (d)). Interestingly, it
undergoes the same type of bifurcation for parameters in set II, and hence a
bifurcation plot for parameters in set II is not shown in Fig. 2.

Figure 2: Bifurcation diagrams of system (1), where $I$ is the bifurcation
parameter. Panels (a) and (b) are for parameters in set I and (c), (d) for
parameters in set III. Thick and dashed-blue curves represent stable and
unstable branches, i.e., stability analysis of the fixed points, green and red
curves represent stable and unstable limit cycles. Panels (b), (d) show the
nullclines for the emergence of saddle node bifurcations for three different
values of $I$. We note that system (1) undergoes the same type of bifurcation
for parameters in set II, and hence a bifurcation plot for this parameter set
is not shown. Note also that HB and HB1 stand for Hopf-bifurcation points and
SN and SN1 for saddle-node bifurcation points.

##  3 Diffusively coupled Izhikevich models on a 2D spatial domain

In the following, we consider the system of diffusively coupled Izhikevich
neurons

$\begin{array}[]{l}\frac{{\partial v}}{{\partial
t}}=f_{1}(v,u)+D_{11}\Big{(}\frac{\partial^{2}v}{\partial
x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\Big{)},\\\ \frac{{\partial
u}}{{\partial t}}=f_{2}(v,u)+D_{22}\Big{(}\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\Big{)},\end{array}$ (2)

on a 2D spatial domain $\Omega$, extending the single deterministic Izhikevich
model (1). In this framework, $v\equiv v(x,y,t)$, $u\equiv u(x,y,t)$ and $t$,
$x$ and $y$ represent the time, $x$ and $y$ spatial coordinates, respectively.
Therefore, the solution $(v,u)$ represents the activity of the Izhikevich
neurons in $x$ and $y$ at time $t$. Parameters $D_{11}$ and $D_{22}>0$
quantify the diffusion coupling strengths. The initial conditions are given by
the known functions $v(x,y,t=0)$ and $u(x,y,t=0)$ for $x,y\in\Omega$, where
$\Omega$ is the bounded-square or spatial domain. Moreover, we consider zero-
flux boundary conditions $\frac{{\partial v}}{{\partial n}}=\frac{{\partial
u}}{{\partial n}}=0$, $x,y\in\partial\Omega$ and $t>0$, where $n$ is the
outward normal to the boundary of $\Omega$, denoted by $\partial\Omega$. The
zero-flux boundary conditions imply that the cell membranes are impermeable to
ionic movements at the boundaries [8].

For diffusive systems such as system (2), the Turing condition [27, 29]
implies that its steady-state solution is unstable around a fixed point and
stable for the system with no diffusion, such as system (1). Next, we perturb
system (2) around the uniform steady-state condition and linearize it around
the nontrivial fixed point, obtaining the characteristic equation, assuming
the particular solution $\left({\begin{array}[]{*{20}{c}}v\\\
u\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{{v_{0}}}\\\
{{u_{0}}}\end{array}}\right)+\varepsilon\left({\begin{array}[]{*{20}{c}}{{v_{k}}}\\\
{{u_{k}}}\end{array}}\right){e^{\lambda_{k}t+ikr}}+c.c.+o\left({{\varepsilon^{2}}}\right)$,
where $c.c$. stands for complex conjugate, $\lambda_{k}$ is the wave length
and $k$ the wave number in the $r$ direction, where $r=(x,y)$ is the
directional vector in two dimensions. The Jacobian matrix, $J_{D}$, of Eq. (2)
computed at the fixed point $E=(v_{0},u_{0})$ is then given by
$J_{D}=J-k^{2}D=\left({\begin{array}[]{*{20}{c}}{{a_{11}-D_{11}k^{2}}}&{{a_{12}}}\\\
{{a_{21}}}&{{a_{22}-D_{22}k^{2}}}\end{array}}\right),$ where
$D=\left({\begin{array}[]{*{20}{c}}{{D_{11}}}&{{0}}\\\
{{0}}&{{D_{22}}}\end{array}}\right)$. The stability conditions of system (2)
at $E=(v_{0},u_{0})$ are given by
trace$(J_{D})=a_{11}+a_{22}-D_{11}k^{2}-D_{22}k^{2}<0$ and
$\det(J_{D})=(a_{11}-D_{11}k^{2})(a_{22}-D_{22}k^{2})-a_{12}a_{21}>0$. The
eigenvalues of $J_{D}$ are given by
$\lambda=\frac{1}{2}\Big{(}\mbox{trace}(J_{D})\pm\sqrt{\mbox{trace}^{2}(J_{D})-4\det(J_{D})}\Big{)}$,
which is also known as the dispersion relation. If we consider the stable
fixed point of system (1), then trace$(J_{D})<0$ as $a_{11}+a_{22}<0$, $D>0$
and $k^{2}>0$. Thus, the conditions for the Turing instability become
$\mbox{trace}(J)=a_{11}+a_{22}<0$, $\det(J)=a_{11}a_{22}-a_{12}a_{21}>0$ and
$\det(J_{D})=\big{(}a_{11}-D_{11}k^{2}\big{)}\big{(}a_{22}-D_{22}k^{2}\big{)}-a_{12}a_{21}<0$.

Figure 3: Plots of the dispersion relation for system (2) for various
couplings: (a) $I=-105.1$, $D_{11}=0.1$, where the green, red and blue curves
correspond to $D_{22}=4$, 6.68117 and 9 (set I), (b) $I=80$, $D_{22}=0.01$,
where the green, red and blue lines correspond to $D_{11}=1.1$, 1.45612 and
1.9 (set III), respectively. Panels (c) and (d) show plots of Eq. (25) for (c)
$I=-105.1$, $D_{11}=0.1$, $D_{22}=9$ (set I) and (d) $I=80$, $D_{11}=1.46$,
$D_{22}=0.01$ (set III). Note that $\lambda$ are the eigenvalues of the
Jacobian matrix $J_{D}$, thus Re$(\lambda)$ are their real parts and that
$\zeta_{1}$, $\zeta_{2}$, $\zeta_{3}$ are the modes discussed in Sec. 3.2. (e)
Bifurcation diagram of model (2) for parameter set I with $I=-105.1$,
$D_{11}=0.1$ and $D_{22}=9$. $H_{0}$: hexagonal patterns with $\varphi=0$,
$H_{\pi}$: hexagonal patterns with $\varphi=\pi$. Solid and dashed lines
indicate stable and unstable states.

To appreciate the Turing instability, we have computed numerically and present
in Fig. 3(a), (b) the dispersion relation for various couplings. The function
$\det(J_{D})$ is quadratic in $k^{2}$, say $H(k^{2})$, which can be written as
$H(k^{2})=D_{11}D_{22}k^{4}+(-a_{22}D_{11}-a_{11}D_{22})k^{2}+\det(J)$. The
minimum value of $H(k^{2})$ is attained at the critical wave number $k=k_{T}$,
where

$\displaystyle k_{T}^{2}=\frac{a_{22}D_{11}+a_{11}D_{22}}{2D_{11}D_{22}}.$ (3)

Note that, setting $H(k_{T}^{2})=0$ and using Eq. (3), we obtain the critical
wave number $k_{T}=\left(\frac{\det(J)}{\det(D)}\right)^{\frac{1}{4}}$. In the
next section, we derive the amplitude equations of system (2) in the vicinity
of the Turing bifurcation point, where the system exhibits critical slowing
down. Here, we only consider constant diagonal diffusion couplings as this
leads to interesting phenomena in excitable spatial systems. In the
neighborhood of a Turing bifurcation point, the spatial symmetry of the system
is destroyed and a stationary form arises in both time and oscillatory
dynamics that can be described by the amplitude equations, discussed in the
next section.

###  3.1 Amplitude Equations

The characteristics of different pattern formations around the bifurcation
point can be analyzed further by studying the amplitude equations of system
(2). Particularly, these equations can be derived by using a multiple-scale
analysis [22, 23, 37]. Perturbing the system around the fixed point
$E=\left({{v_{0}},{u_{0}}}\right)$ by considering $v=\tilde{v}+{v_{0}}$ and
$u=\tilde{u}+{u_{0}}$, we obtain the set of equations

$\begin{array}[]{l}\frac{{\partial\tilde{v}}}{{\partial
t}}={a_{11}}\tilde{v}+{a_{12}}\tilde{u}+0.04{{\tilde{v}}^{2}}+D_{11}{\nabla^{2}}\tilde{v},\\\
\frac{{\partial\tilde{u}}}{{\partial
t}}={a_{21}}\tilde{v}+{a_{22}}\tilde{u}+D_{22}{\nabla^{2}}\tilde{u},\end{array}$
(4)

where $\tilde{v}$ and $\tilde{u}$ are small perturbations to $v_{0}$ and
$u_{0}$, respectively. Equation (4) can be written in the compact form

$\frac{{\partial P}}{{\partial t}}=LP+N,$ (5)

where $L$ and $N$ are the linear operator and nonlinear term, respectively,
that we discuss below. We Taylor-expand the bifurcation parameter $I$, which
is close to the bifurcation threshold, $I_{T}$, as
$I-I_{T}=\varepsilon{I_{1}}+{\varepsilon^{2}}{I_{2}}+{\varepsilon^{3}}{I_{3}}+o\left({{\varepsilon^{4}}}\right)$,
where $\left|\varepsilon\right|<1$. We can write $P$ and $N$ as

$\displaystyle P$
$\displaystyle=\left({\begin{array}[]{*{20}{c}}{{\tilde{v}}}\\\
{{\tilde{u}}}\end{array}}\right)=\varepsilon\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)+{\varepsilon^{2}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{2}}}\\\
{{\tilde{u}_{2}}}\end{array}}\right)+{\varepsilon^{3}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{3}}}\\\
{{\tilde{u}_{3}}}\end{array}}\right)+o\left({{\varepsilon^{3}}}\right),$ (14)
$\displaystyle N$
$\displaystyle={\varepsilon^{2}}{n_{2}}+{\varepsilon^{3}}{n_{3}}+o\left({{\varepsilon^{3}}}\right).$
(15)

In the last equation, $n_{2}$ and $n_{3}$ represent the second and third order
of $\varepsilon$ and, are given by
${n_{2}}=\left({\begin{array}[]{*{20}{c}}{0.04{\tilde{v}_{1}}^{2}}\\\
0\end{array}}\right)$ and
${n_{3}}=\left({\begin{array}[]{*{20}{c}}{0.08\tilde{v}_{1}\tilde{v}_{2}}\\\
0\end{array}}\right)$, respectively. The linear operator $L$ can be written in
the form

$\begin{array}[]{l}L={L_{T}}+\left({I-I_{T}}\right){L_{1}}+{\left({I-I_{T}}\right)^{2}}{L_{2}}+o\left({{{\left({I-I_{T}}\right)}^{3}}}\right)={L_{T}}+\left({I-I_{T}}\right)M+o\left({{{\left({I-I_{T}}\right)}^{2}}}\right),\end{array}$
(16)

where
${L_{i}}=\frac{1}{{i!}}\frac{{{\partial^{i}}L}}{{\partial{I^{i}}}},\,i=1,2,\;{L_{T}}=\left({\begin{array}[]{*{20}{c}}{{a_{11}^{T}}+D_{11}{\nabla^{2}}}&{{a_{12}^{T}}}\\\
{{a_{21}^{T}}}&{{a_{22}^{T}}+D_{22}{\nabla^{2}}}\end{array}}\right)$ and
$M=\left({\begin{array}[]{*{20}{c}}{{m_{11}}}&0\\\ 0&0\end{array}}\right)$
with $m_{11}=\frac{-0.08}{0.08v_{0}+5-b}$.

Next, we introduce the multiple-scale [23, 37] $\frac{\partial}{{\partial
t}}=\varepsilon\frac{\partial}{{\partial{T_{1}}}}+{\varepsilon^{2}}\frac{\partial}{{\partial{T_{2}}}}+o\left({{\varepsilon^{2}}}\right)$,
where ${T_{1}}=\varepsilon t$ and ${T_{2}}={\varepsilon^{2}}t$. Substituting
Eqs. (14) - (16) into Eq. (5), we obtain

$\begin{array}[]{l}\left({\varepsilon\frac{\partial}{{\partial{T_{1}}}}+{\varepsilon^{2}}\frac{\partial}{{\partial{T_{2}}}}}\right)\left[{\varepsilon\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)+{\varepsilon^{2}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{2}}}\\\
{{\tilde{u}_{2}}}\end{array}}\right)+{\varepsilon^{3}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{3}}}\\\
{{\tilde{u}_{3}}}\end{array}}\right)}\right]=\\\
\left({{L_{T}}+\left({I-I_{T}}\right)M}\right)\left[{\varepsilon\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)+{\varepsilon^{2}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{2}}}\\\
{{\tilde{u}_{2}}}\end{array}}\right)+{\varepsilon^{3}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{3}}}\\\
{{\tilde{u}_{3}}}\end{array}}\right)}\right]+{\varepsilon^{2}}{n_{2}}+{\varepsilon^{3}}{n_{3.}}\end{array}$
(17)

Equating the coefficients of $\varepsilon$, $\varepsilon^{2}$ and
$\varepsilon^{3}$ in both sides of Eq. (17), we obtain the equations

${L_{T}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)=0,$ (18)
${L_{T}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{2}}}\\\
{{\tilde{u}_{2}}}\end{array}}\right)=\frac{\partial}{{\partial{T_{1}}}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)-{I_{1}}M\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)-{n_{2}}\buildrel\Delta\over{=}\left({\begin{array}[]{*{20}{c}}{{F_{\tilde{v}}}}\\\
{{F_{\tilde{u}}}}\end{array}}\right),$ (19)

and

${L_{T}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{3}}}\\\
{{\tilde{u}_{3}}}\end{array}}\right)=\frac{\partial}{{\partial{T_{1}}}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{2}}}\\\
{{\tilde{u}_{2}}}\end{array}}\right)+\frac{\partial}{{\partial{T_{2}}}}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)-{I_{1}}M\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{2}}}\\\
{{\tilde{u}_{2}}}\end{array}}\right)-{I_{2}}M\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)-{n_{3}}.$ (20)

Solving Eq. (18), we obtain
${\tilde{v}_{1}}=\frac{{a_{11}^{T}}D_{22}-{a_{22}^{T}}D_{11}}{2D_{11}{a_{21}^{T}}}={f}$
and ${\tilde{u}_{1}}=1$. We can write

$\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{1}}}\\\
{{\tilde{u}_{1}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{{f}}\\\
1\end{array}}\right){\sum\limits_{j=1}^{3}{{\Theta_{j}}}e^{i{k_{j}}.r}+c.c},$

where $(\tilde{v}_{1},\tilde{u}_{1})$ is the linear combination of the
eigenvectors that corresponds to the zero eigenvalue of the linear operator
$L_{T}$. In this context, $\Theta_{j}$ denotes the amplitude of the mode
$e^{i{k_{j}}.r}$. We use the Fredholm solvability condition [33, 37] to find
the nontrivial solution of Eq. (19). The right hand-side of Eq. (19) and the
eigenvectors corresponding to the zero eigenvalue of $L_{T}^{+}$ must be
orthogonal. $L_{T}^{+}$ indicates the adjoint operator of $L_{T}$ and the
eigenvectors of $L_{T}^{+}$ are given by
$\left({\begin{array}[]{*{20}{c}}{{1}}\\\ g\end{array}}\right)e^{ik_{j}r}$,
where $j=1,2,3$ and
$g=\frac{{a_{22}^{T}}D_{11}-{a_{11}^{T}}D_{22}}{2D_{22}{a_{21}^{T}}}$.
Applying the orthogonality condition,
$\left({\begin{array}[]{*{20}{c}}{{1}}\,\,\,\,{{g}}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}{{F_{\tilde{p}}^{i}}}\\\
F_{\tilde{q}}^{i}\end{array}}\right)=0$, we obtain the equations

$\begin{array}[]{l}\left({f+g}\right)\frac{{\partial{\Theta_{1}}}}{{\partial{T_{1}}}}={I_{1}}{{m_{11}}{f}}{\Theta_{1}}-2{h_{1}}{{\bar{\Theta}}_{2}}{{\bar{\Theta}}_{3}},\\\
\left({f+g}\right)\frac{{\partial{\Theta_{2}}}}{{\partial{T_{1}}}}={I_{1}}{{m_{11}}{f}}{\Theta_{2}}-2{h_{1}}{{\bar{\Theta}}_{1}}{{\bar{\Theta}}_{3}},\\\
\left({f+g}\right)\frac{{\partial{\Theta_{3}}}}{{\partial{T_{1}}}}={I_{1}}{{m_{11}}{f}}{\Theta_{3}}-2{h_{1}}{{\bar{\Theta}}_{1}}{{\bar{\Theta}}_{2}},\end{array}$

where $F_{\tilde{v}}^{i}$, $F_{\tilde{u}}^{i}$ are the coefficients of
$e^{ik_{i}.r}$ in $F_{\tilde{v}}$, $F_{\tilde{u}}$ and
${h_{1}}=-0.04{f_{1}}^{2}$. Solving Eq. (19), we obtain

$\begin{array}[]{l}\left({\begin{array}[]{*{20}{c}}{{\tilde{v}_{2}}}\\\
{{\tilde{u}_{2}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{{V_{0}}}\\\
{{U_{0}}}\end{array}}\right)+\sum\limits_{j=1}^{3}{\left({\begin{array}[]{*{20}{c}}{{V_{j}}}\\\
{{U_{j}}}\end{array}}\right)e^{i{k_{j}}.r}}+\sum\limits_{j=1}^{3}{\left({\begin{array}[]{*{20}{c}}{{V_{jj}}}\\\
{{U_{jj}}}\end{array}}\right)e^{i2{k_{j}}.r}}+\left({\begin{array}[]{*{20}{c}}{{V_{12}}}\\\
{{U_{12}}}\end{array}}\right)e^{i({k_{1}}-{k_{2}}).r}+\\\
\left({\begin{array}[]{*{20}{c}}{{V_{23}}}\\\
{{U_{23}}}\end{array}}\right)e^{i({k_{2}}-{k_{3}}).r}+\left({\begin{array}[]{*{20}{c}}{{V_{31}}}\\\
{{U_{31}}}\end{array}}\right)e^{i({k_{3}}-{k_{1}}).r}+c.c,\end{array}$ (21)

where, the coefficients are given in the Appendix. Again, using the
orthogonality condition, we get

$\begin{array}[]{l}\vspace{0.2cm}\left({f+g}\right)\left({\frac{{\partial{U_{1}}}}{{\partial{T_{1}}}}+\frac{{\partial{\Theta_{1}}}}{{\partial{T_{2}}}}}\right)={m_{11}}f\left({{I_{1}}{U_{1}}+{I_{2}}{\Theta_{1}}}\right)+H\left({{{\bar{Q}}_{2}}{{\bar{\Theta}}_{3}}+{{\bar{Q}}_{3}}{{\bar{\Theta}}_{2}}}\right)-\,\left[{{\sigma_{1}}{{\left|{{\Theta_{1}}}\right|}^{2}}+{\sigma_{2}}\left({{{\left|{{\Theta_{2}}}\right|}^{2}}+{{\left|{{\Theta_{3}}}\right|}^{2}}}\right)}\right]{\Theta_{1}},\\\
\vspace{0.2cm}\left({f+g}\right)\left({\frac{{\partial{U_{2}}}}{{\partial{T_{1}}}}+\frac{{\partial{\Theta_{2}}}}{{\partial{T_{2}}}}}\right)={m_{11}}f\left({{I_{1}}{U_{2}}+{I_{2}}{\Theta_{2}}}\right)+H\left({{{\bar{Q}}_{1}}{{\bar{\Theta}}_{3}}+{{\bar{Q}}_{3}}{{\bar{\Theta}}_{1}}}\right)-\,\left[{{\sigma_{1}}{{\left|{{\Theta_{1}}}\right|}^{2}}+{\sigma_{2}}\left({{{\left|{{\Theta_{2}}}\right|}^{2}}+{{\left|{{\Theta_{3}}}\right|}^{2}}}\right)}\right]{\Theta_{2}},\\\
\left({f+g}\right)\left({\frac{{\partial{U_{3}}}}{{\partial{T_{1}}}}+\frac{{\partial{\Theta_{3}}}}{{\partial{T_{2}}}}}\right)={m_{11}}f\left({{I_{1}}{U_{3}}+{I_{2}}{\Theta_{3}}}\right)+H\left({{{\bar{Q}}_{1}}{{\bar{\Theta}}_{2}}+{{\bar{Q}}_{2}}{{\bar{\Theta}}_{1}}}\right)-\,\left[{{\sigma_{1}}{{\left|{{\Theta_{1}}}\right|}^{2}}+{\sigma_{2}}\left({{{\left|{{\Theta_{2}}}\right|}^{2}}+{{\left|{{\Theta_{3}}}\right|}^{2}}}\right)}\right]{\Theta_{3}}.\end{array}$
(22)

Furthermore, the amplitude ${\Lambda_{j}}$, $j=1,2,3$ can be expanded as
${\Lambda_{j}}=\varepsilon{\Theta_{j}}+{\varepsilon^{2}}{U_{j}}+o({\varepsilon^{3}})$.
Using the expression for $\Lambda_{j}$ and Eq. (22), we obtain the amplitude
equations that correspond to $\Lambda_{j}$,

$\begin{array}[]{l}{\tau_{0}}\frac{{\partial{\Lambda_{1}}}}{{\partial
t}}=\mu{\Lambda_{1}}+h{{\bar{\Lambda}}_{2}}{{\bar{\Lambda}}_{3}}-\left[{\sigma_{1}^{\prime}{{\left|{{\Lambda_{1}}}\right|}^{2}}+\sigma_{2}^{\prime}\left({{{\left|{{\Lambda_{2}}}\right|}^{2}}+{{\left|{{\Lambda_{3}}}\right|}^{2}}}\right)}\right]{\Lambda_{1}},\\\
{\tau_{0}}\frac{{\partial{\Lambda_{2}}}}{{\partial
t}}=\mu{\Lambda_{2}}+h{{\bar{\Lambda}}_{1}}{{\bar{\Lambda}}_{3}}-\left[{\sigma_{1}^{\prime}{{\left|{{\Lambda_{2}}}\right|}^{2}}+\sigma_{2}^{\prime}\left({{{\left|{{\Lambda_{1}}}\right|}^{2}}+{{\left|{{\Lambda_{3}}}\right|}^{2}}}\right)}\right]{\Lambda_{2}},\\\
{\tau_{0}}\frac{{\partial{\Lambda_{3}}}}{{\partial
t}}=\mu{\Lambda_{3}}+h{{\bar{\Lambda}}_{1}}{{\bar{\Lambda}}_{2}}-\left[{\sigma_{1}^{\prime}{{\left|{{\Lambda_{3}}}\right|}^{2}}+\sigma_{2}^{\prime}\left({{{\left|{{\Lambda_{1}}}\right|}^{2}}+{{\left|{{\Lambda_{2}}}\right|}^{2}}}\right)}\right]{\Lambda_{3}}.\end{array}$
(23)

All the parameters involved in Eqs. (22) and (23) are discussed in the
Appendix.

###  3.2 Amplitude Stability

Next, we study the stability of the patterns obtained near the bifurcation
point. The amplitude $\Lambda_{j}$ can be decomposed into the mode
$\zeta_{j}=|\Lambda_{j}|$ and a corresponding phase-angle $\varphi_{j}$. Then,
substituting ${\Lambda_{j}}={\zeta_{j}}e^{i{\varphi_{j}}}$ into Eq. (23), and
separating the real and imaginary parts, we obtain the following system of
ordinary differential equations for the real variables

$\begin{array}[]{l}{\tau_{0}}\frac{{\partial\varphi}}{{\partial
t}}=-h\frac{{\zeta_{1}^{2}\zeta_{2}^{2}+\zeta_{1}^{2}\zeta_{3}^{2}+\zeta_{2}^{2}\zeta_{3}^{2}}}{{{\zeta_{1}}{\zeta_{2}}{\zeta_{3}}}}\sin\varphi,\\\
{\tau_{0}}\frac{{\partial{\zeta_{1}}}}{{\partial
t}}=\mu{\zeta_{1}}+h{\zeta_{2}}{\zeta_{3}}\cos\varphi-\sigma_{1}^{\prime}\zeta_{1}^{3}-\sigma_{2}^{\prime}\big{(}\zeta_{2}^{2}+\zeta_{3}^{2}\big{)}{\zeta_{1}},\\\
{\tau_{0}}\frac{{\partial{\zeta_{2}}}}{{\partial
t}}=\mu{\zeta_{2}}+h{\zeta_{1}}{\zeta_{3}}\cos\varphi-\sigma_{1}^{\prime}\zeta_{2}^{3}-\sigma_{2}^{\prime}\big{(}\zeta_{1}^{2}+\zeta_{3}^{2}\big{)}{\zeta_{2}},\\\
{\tau_{0}}\frac{{\partial{\zeta_{3}}}}{{\partial
t}}=\mu{\zeta_{3}}+h{\zeta_{1}}{\zeta_{2}}\cos\varphi-\sigma_{1}^{\prime}\zeta_{3}^{3}-\sigma_{2}^{\prime}\big{(}\zeta_{1}^{2}+\zeta_{2}^{2}\big{)}{\zeta_{3}},\end{array}$
(24)

where $\varphi={\varphi_{1}}+{\varphi_{2}}+{\varphi_{3}}$. Here, we are
interested in the stable steady-state. From the first equation in Eq. (24), we
have that $\varphi=0$ or $\pi$, as $h\zeta_{i}\neq 0$. The state corresponding
to $\varphi=0$ and $\varphi=\pi$ is stable when $h>0$ and $h<0$, respectively.
Then, the amplitudes of Eq. (24) are given by

$\begin{array}[]{l}{\tau_{0}}\frac{{\partial{\zeta_{1}}}}{{\partial
t}}=\mu{\zeta_{1}}+h{\zeta_{2}}{\zeta_{3}}-\sigma_{1}^{\prime}\zeta_{1}^{3}-\sigma_{2}^{\prime}\big{(}\zeta_{2}^{2}+\zeta_{3}^{2}\big{)}{\zeta_{1}},\\\
{\tau_{0}}\frac{{\partial{\zeta_{2}}}}{{\partial
t}}=\mu{\zeta_{2}}+h{\zeta_{1}}{\zeta_{3}}-\sigma_{1}^{\prime}\zeta_{2}^{3}-\sigma_{2}^{\prime}\big{(}\zeta_{1}^{2}+\zeta_{3}^{2}\big{)}{\zeta_{2}},\\\
{\tau_{0}}\frac{{\partial{\zeta_{3}}}}{{\partial
t}}=\mu{\zeta_{3}}+h{\zeta_{1}}{\zeta_{2}}-\sigma_{1}^{\prime}\zeta_{3}^{3}-\sigma_{2}^{\prime}\big{(}\zeta_{1}^{2}+\zeta_{2}^{2}\big{)}{\zeta_{3}}.\end{array}$
(25)

System (25) admits different kinds of solutions such as stationary states,
stripe patterns, hexagonal patterns and mixed states. To investigate pattern
formation, we perform a linear stability analysis. Depending on the parameters
$\mu$, $\sigma_{1}^{\prime}$, $\sigma_{2}^{\prime}$ and $h$, different
structural patterns can emerge in system (2). As the Jacobian matrix of system
(25) is given by

$\left({\begin{array}[]{*{20}{c}}{\mu-3\sigma_{1}^{\prime}\zeta_{1}^{2}-\sigma_{2}^{\prime}\big{(}{\zeta_{2}^{2}+\zeta_{3}^{2}}\big{)}}&{\left|h\right|{\zeta_{3}}-2\sigma_{2}^{\prime}{\zeta_{1}}{\zeta_{2}}}&{\left|h\right|{\zeta_{2}}-2\sigma_{2}^{\prime}{\zeta_{1}}{\zeta_{3}}}\\\
{\left|h\right|{\zeta_{3}}-2\sigma_{2}^{\prime}{\zeta_{1}}{\zeta_{2}}}&{\mu-3\sigma_{1}^{\prime}\zeta_{2}^{2}-\sigma_{2}^{\prime}\big{(}{\zeta_{1}^{2}+\zeta_{3}^{2}}\big{)}}&{\left|h\right|{\zeta_{1}}-2\sigma_{2}^{\prime}{\zeta_{2}}{\zeta_{3}}}\\\
{\left|h\right|{\zeta_{2}}-2\sigma_{2}^{\prime}{\zeta_{1}}{\zeta_{3}}}&{\left|h\right|{\zeta_{1}}-2\sigma_{2}^{\prime}{\zeta_{2}}{\zeta_{3}}}&{\mu-3\sigma_{1}^{\prime}\zeta_{3}^{2}-\sigma_{2}^{\prime}\big{(}{\zeta_{1}^{2}+\zeta_{2}^{2}}\big{)}}\end{array}}\right),$

we can identify the following cases:

  * 1.

The stationary state is given by ${\zeta_{1}}={\zeta_{2}}={\zeta_{3}}=0$ and
is stable for $\mu<0$ and unstable for $\mu>0$.

  * 2.

The hexagonal pattern exists when
$\zeta={\zeta_{1}}={\zeta_{2}}={\zeta_{3}}=\frac{{\left|h\right|\pm\sqrt{{h^{2}}+4\left({{\sigma_{1}^{\prime}}+2{\sigma_{2}^{\prime}}}\right)\mu}}}{{2\left({{\sigma_{1}^{\prime}}+2{\sigma_{2}^{\prime}}}\right)}}$,
with $\varphi=0$ or $\varphi=\pi$ and
$\mu>{\mu_{1}}=\frac{{-{h^{2}}}}{{4\left({{\sigma_{1}^{\prime}}+2{\sigma_{2}^{\prime}}}\right)}}$.
The hexagonal pattern $H_{\pi}$ (when $\varphi=\pi$) is stable for
$\mu<{\mu_{4}}=\frac{{2{\sigma_{1}^{\prime}}+{\sigma_{2}^{\prime}}}}{{{{\left({{\sigma_{2}^{\prime}}-{\sigma_{1}^{\prime}}}\right)}^{2}}}}{h^{2}}$
only, whereas, $H_{0}$ is always unstable, when $\varphi=0$.

  * 3.

Finally, the mixed state exists when $\sigma_{2}^{\prime}>\sigma_{1}^{\prime}$
and $\mu>\sigma_{1}^{\prime}\zeta_{1}^{2}=\mu_{3}$, where
${\zeta_{1}}=\frac{|h|}{\sigma_{2}^{\prime}-\sigma_{1}^{\prime}}$,
${\zeta_{2}}=\sqrt{\frac{\mu-\sigma_{1}^{\prime}\zeta_{1}^{2}}{\sigma_{1}^{\prime}+\sigma_{2}^{\prime}}}$,
${\zeta_{3}}=\sqrt{\frac{\mu-\sigma_{1}^{\prime}\zeta_{1}^{2}}{\sigma_{1}^{\prime}+\sigma_{2}^{\prime}}}$,
and is always unstable.

Figure 4: Pattern formation in system (2) of diffusively coupled, 2D
Izhikevich neurons. First row of panels: set I for $I=-105.1$, $D_{11}=0.1$
and $D_{22}=9$: at (a) $t=100$, (b) $t=150$ and (d) $t=200$. Second row of
panels: set III for $I=80$, $D_{11}=1.46$ and $D_{22}=0.01$: at (e) $t=60$,
(f) $t=63$ and (h) $t=65$. Panels (c) and (g) show zoom-in versions of panels
(b) and (f), respectively. Note that in all panels, the colorbar encodes the
values of the solution $u(x,y)$ of system (2).

In the next section, we present numerical results that corroborate the
analytical results derived in this section using the amplitude equations.

###  3.3 Numerical results for diffusively coupled Izhikevich neurons on a 2D
spatial domain

Here, we analyze the spatiotemporal evolution of patterns exhibited by system
(2) for different diffusion couplings. We consider an extended medium with 2D
spatial diffusion and $N\times N$ coupled, equally spaced, identical
Izhikevich neurons given by system (2), where $N=200$. We solve numerically
system (2) using the forward Euler method with a central difference
representation, which is a finite-difference scheme, with step sizes $\Delta
x=\Delta y=\Delta=0.5$ and $\delta t=0.001$ for the spatial mesh and time-
integration step size, respectively. The initial conditions are considered
with appropriate periodic perturbations as in the case of system (1). We
denote the variations in the membrane voltage of a single neuron at the
spatial coordinate $(i,j)$ at time $t$ as $v_{i,j}(t)$. The partial derivative
is then approximated by
$\frac{{{\partial^{2}}v}}{{\partial{x^{2}}}}+\frac{{{\partial^{2}}v}}{{\partial{y^{2}}}}=\frac{1}{{{\Delta^{2}}}}\big{(}{{v_{i-1,\,j}}+{v_{i+1,\,j}}+{v_{i,\,j-1}}+{v_{i,\,j+1}}-4{v_{i,\,j}}}\big{)}$.
We note that here we want to investigate the characteristics of pattern
formation in the context of the diffusive system (2), where the diffusion
couplings account for the synaptic couplings among the Izhikevich neurons on
the 2D spatial domain.

In the following, we consider positive diffusion couplings $D_{11}$ and
$D_{22}$ such that a Turing instability can manifest in the dynamics of system
(2). Then, we identify the conditions for the emergence and stability of
stationary, hexagonal and mixed-state patterns, using the amplitude equations.
The conditions for the emergence and stability properties of these
spatiotemporal patterns depend on the parameter values of the system, among
which are the diffusion couplings $D_{11}$ and $D_{22}$. Hence any pair of
$D_{11}$ and $D_{22}$ values that satisfy the above conditions, could be a
choice.

The spatial patterns observed are organized into simple hexagons and mixed-
state arrangements of spot-like formations, which are temporarily stationary.
The patterns that arise around the spatially uniform steady-state solution
occur at a finite wave number. For $I=-105.1$, $D_{11}=0.1$ and $D_{22}=9$ of
set I, we obtain $h=-0.3079$, $\sigma_{1}^{\prime}=1.0675$,
$\sigma_{2}^{\prime}=2.9011$, $\mu_{1}=-0.0034$, $\mu_{3}=0.0301$,
$\mu_{4}=0.1420$ and $\mu=0.0106$. Clearly, $\mu>0$, which indicates the
stationary state is unstable. Here, $\mu>\mu_{1}$ and $h<0$, which pinpoint to
the existence of hexagonal patterns, $H_{\pi}$ [21, 23, 33, 34]. Also, since
$\mu<\mu_{4}$, it confirms the existence of stable hexagons. We perform the
linear stability analysis for hexagonal patterns and pay particular attention
to the destabilizing modes that are moved with respect to the patterns.
Summarize the above analysis, we can show our results as a bifurcation diagram
in Fig. 3(e). The $H_{\pi}$ patterns are stable in the region $\mu<\mu_{4}$.
The $H_{0}$ patterns are always unstable. We have validated the existence of
spatial modulation of hexagonal patterns by solving Eq. (25) and present them
in Fig. 4 (a)-(d) for long integration times. A zoom-in, shown in Fig. 4(c),
confirms their appearance, which is shown in small patches. The hexagonal
patterns, $H_{\pi}$, are observed on the whole domain. Interestingly, we have
shown both the low-amplitude oscillations (deep-blue regions) and high
amplitude oscillations (small yellow patches) to validate the time-independent
hexagonal patterns. The numerical values for the existence of stable hexagons
are given by $\zeta_{1}=\zeta_{2}=\zeta_{3}=0.000000071\approx 0$, which are
also verified by the analytical values discussed in Subsec. 3.2. The spatial
instabilities lead to spatially modulated hexagons and other mixed states
characterized by mixed-type patterns. Similarly, for $I=80$, $D_{11}=1.46$ and
$D_{22}=0.01$ of set III, we find $h=0.0058$,
$\sigma_{1}^{\prime}=-0.0002534$, $\sigma_{2}^{\prime}=-0.00023657$,
$\mu_{1}=0.0116$, $\mu_{3}=-30.2297$, $\mu_{4}=-88.6813$ and $\mu=0.0118$.
Here, the stationary state is unstable as $\mu>0$. The above numerical values
confirm the existence of unstable mixed-type patterns (as
$\sigma_{2}^{\prime}>\sigma_{1}^{\prime}$ and $\mu>\mu_{3}$) and unstable
hexagons, $H_{0}$, (as $\mu>\mu_{1}$ and $h>0$), shown in Fig. 4(e)-(h).
However, as time increases, the patterns become distorted and mixed-state
patterns emerge. The coexistence of ordered and disordered states are shown in
Fig. 4(f). With further increase in time, most of the neurons’ activity leads
to a steady state, shown in Fig. 4(h). Finally, plots of the corresponding
time-series are shown in Fig. 3(c)-(d), using Eq.(25).

##  4 Conclusions

In this work, we introduced a mathematical methodology to analyze and
characterize the propagation of wavy patterns in diffusively coupled
Izhikevich neurons on a 2D spatial domain. Our approach focused on the
development of a theory including diffusive, Turing instabilities. The non-
uniform, spatiotemporal patterns exhibit collective dynamics that can be found
in neural recordings for particular working mechanisms in brain functioning
[3, 4, 8, 9, 12, 35, 38].

First, we described the oscillatory behavior of a single Izhikevich neuron and
presented the corresponding bifurcation analysis. Next, we explored
diffusively coupled Izhikevich neurons on a 2D configuration. We characterized
the spatial instability and derived a three-mode normal form. Later, this
normal form was used to study pattern formation in the diffusive system. We
explained theoretically and showed numerically the emergence of hexagonal
patterns. Even for suitable parameter regimes at which stable hexagons exist,
persistent periodic and irregular dynamical behavior can be found [32]. We
focused on instabilities that lead to the emergence of various patterns in the
Turing domain. Sometimes, the oscillatory short-wave instabilities can form
stable, periodically modulated, hexagons. Furthermore, the spatial
instabilities lead to spatially modulated hexagons and other mixed states
characterized by mixed-type patterns. Consequently, we have been able to
verify the observed spatiotemporal activities with our analytical method.
Thus, our present study and the work in [15], show clearly that the
mathematical approach used here can be generalized to study analytically
spatiotemporal instabilities and pattern formation in excitable extended
neural models, including in Turing and non-Turing regimes. The diffusively
coupled system showed dynamical coherent structures and the single-neuron
model showed highly variable complex firings. This may be useful in future
work to investigate the spatial scale of neural features with various
oscillatory signals as the low-amplitude oscillations affect the high-
amplitude oscillations [4]. Our theoretical approach explored some of the
mechanisms of pattern formation pertinent to systems of spiking neurons [3,
5]. The observed hexagonal patterns might reflect the rhythmic activities of
neural dynamics such as information processing in locally connected diffusive
networks. The reported results have potential to contribute significantly to
the study of diffusively-coupled neural systems in general.

Appendix
The coefficients of Eq. (21) are given by
$\left({\begin{array}[]{*{20}{c}}{{V_{0}}}\\\
{{U_{0}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{{z_{v0}}}\\\
{{z_{u0}}}\end{array}}\right)\Big{(}{\left|{{\Theta_{1}}}\right|^{2}}+{\left|{{\Theta_{2}}}\right|^{2}}+{\left|{{\Theta_{3}}}\right|^{2}}\Big{)},\,\,{V_{j}}={f}{U_{j}}$,
$\left({\begin{array}[]{*{20}{c}}{{V_{jj}}}\\\
{{U_{jj}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{{z_{v1}}}\\\
{{z_{u1}}}\end{array}}\right)\Theta_{j}^{2},\,\,\left({\begin{array}[]{*{20}{c}}{{V_{jk}}}\\\
{{U_{jk}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{{z_{v2}}}\\\
{{z_{u2}}}\end{array}}\right){\Theta_{j}}{\bar{\Theta}_{k}}$, where
${z_{v0}}=\frac{2{a_{22}^{T}}h_{1}}{A},\,\,{z_{u0}}=\frac{-2{a_{21}^{T}}h_{1}}{A}$,
${z_{v1}}=\frac{\big{(}{a_{22}^{T}}-4D_{22}k_{T}^{2}\big{)}h_{1}}{B},\,\,{z_{u1}}=\frac{-{a_{21}^{T}}h_{1}}{B}$,
${z_{v2}}=\frac{2\big{(}{a_{22}^{T}}-3D_{22}k_{T}^{2}\big{)}h_{1}}{C},\,\,{z_{u2}}=\frac{-2{a_{21}^{T}}h_{1}}{C}$,
$A={a_{11}^{T}}{a_{22}^{T}}-{a_{12}^{T}}{a_{21}^{T}}\mbox{ and
}B=\Big{(}{a_{11}^{T}}-4D_{11}k_{T}^{2}\Big{)}\Big{(}{a_{22}^{T}}-4D_{22}k_{T}^{2}\Big{)}-{a_{12}^{T}}{a_{21}^{T}},C=\Big{(}{a_{11}^{T}}-3D_{11}k_{T}^{2}\Big{)}\Big{(}{a_{22}^{T}}-3D_{22}k_{T}^{2}\Big{)}-{a_{12}^{T}}{a_{21}^{T}}$.
The coefficients of Eq. (22) are given by
$H=0.08{f}^{2},\,\,{\sigma_{1}}=-0.08(z_{v0}+z_{v1})f$ and
${\sigma_{2}}=-0.08(z_{v0}+z_{v2})f$. Finally, the coefficients of Eq. (23)
are given by
${\tau_{0}}=\frac{f+g}{{{I_{T}}{m_{11}}{f}}},\,\,\mu=\frac{{I-I_{T}}}{{{I_{T}}}},\,\,h=\frac{{H}}{{{I_{T}}{m_{11}}{f}}},\,\,\sigma_{1}^{\prime}=\frac{{{\sigma_{1}}}}{{{I_{T}}{m_{11}}{f}}}$
and $\sigma_{2}^{\prime}=\frac{{{\sigma_{2}}}}{{{I_{T}}{m_{11}}{f}}}$.

## References

  * [1] B. Ermentrout, Neural networks as spatio-temporal pattern-forming systems, Reports on Progress in Physics 61 (4) (1998) 353.
  * [2] E. M. Izhikevich, Dynamical Systems in Neuroscience, MIT Press, 2007.
  * [3] A. Keane, P. Gong, Propagating waves can explain irregular neural dynamics, Journal of Neuroscience 35 (4) (2015) 1591–1605.
  * [4] R. G. Townsend, P. Gong, Detection and analysis of spatiotemporal patterns in brain activity, PLOS Computational Biology 14 (12) (2018) e1006643.
  * [5] M. Song, M. Kang, H. Lee, Y. Jeong, S.-B. Paik, Classification of spatiotemporal neural activity patterns in brain imaging data, Scientific Reports 8 (1) (2018) 1–12.
  * [6] J.-Y. Wu, X. Huang, C. Zhang, Propagating waves of activity in the neocortex: what they are, what they do, The Neuroscientist 14 (5) (2008) 487–502.
  * [7] X. Huang, W. Xu, J. Liang, K. Takagaki, X. Gao, J.-Y. Wu, Spiral wave dynamics in neocortex, Neuron 68 (5) (2010) 978–990.
  * [8] 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).
  * [9] S. J. Schiff, X. Huang, J.-Y. Wu, Dynamical evolution of spatiotemporal patterns in mammalian middle cortex, Physical Review Letters 98 (17) (2007) 178102\.
  * [10] L. Yang, M. Dolnik, A. M. Zhabotinsky, I. R. Epstein, Turing patterns beyond hexagons and stripes, Chaos 16 (3) (2006) 037114.
  * [11] H. Jiang, Z. Hou, Formation of spiral waves with substructure in a bursting media, Chaos 25 (12) (2015) 123105.
  * [12] R. G. Townsend, S. S. Solomon, S. C. Chen, A. N. J. Pietersen, P. R. Martin, S. G. Solomon, P. Gong, Emergence of complex wave patterns in primate cerebral cortex, Journal of Neuroscience 35 (11) (2015) 4657–4662.
  * [13] A. E. Tikhomirova, V. A. Volpert, Nonlinear dynamics of endothelial cells, Applied Mathematics Letters 20 (2) (2007) 163–169.
  * [14] Y. Chen, G. Rangarajan, M. Ding, Stability of synchronized dynamics and pattern formation in coupled systems: Review of some recent results, Communications in Nonlinear Science and Numerical Simulation 11 (8) (2006) 934–960.
  * [15] A. Mondal, S. K. Sharma, R. K. Upadhyay, M. A. Aziz-Alaoui, P. Kundu, C. Hens, Diffusion dynamics of a conductance-based neuronal population, Physical Review E 99 (4) (2019) 042307.
  * [16] J. Ma, C.-N. Wang, W.-Y. Jin, Y. Wu, Transition from spiral wave to target wave and other coherent structures in the networks of Hodgkin-Huxley neurons, Applied Mathematics and Computation 217 (8) (2010) 3844–3852.
  * [17] B. Ambrosio, M. A. Aziz-Alaoui, Synchronization and control of coupled reaction-diffusion systems of the FitzHugh-Nagumo type, Computers & Mathematics with Applications 64 (5) (2012) 934–943.
  * [18] B. Ambrosio, M. A. Aziz-Alaoui, V. L. E. Phan, Large time behaviour and synchronization of complex networks of reaction-diffusion systems of FitzHugh-Nagumo type, IMA Journal of Applied Mathematics 84 (2) (2019) 416–443.
  * [19] T. Carletti, H. Nakao, Turing patterns in a network-reduced FitzHugh-Nagumo model, Physical Review E 101 (2) (2020) 022203.
  * [20] M. Kuznetsov, A. Kolobov, A. Polezhaev, Pattern formation in a reaction-diffusion system of Fitzhugh-Nagumo type before the onset of subcritical Turing bifurcation, Physical Review E 95 (5) (2017) 052208.
  * [21] N. Iqbal, R. Wu, B. Liu, Pattern formation by super-diffusion in fitzhugh–nagumo model, Applied Mathematics and Computation 313 (2017) 245–258.
  * [22] M. Ipsen, F. Hynne, P. G. Sørensen, Amplitude equations for reaction-diffusion systems with a Hopf bifurcation and slow real modes, Physica D 136 (1-2) (2000) 66–92.
  * [23] H. Zhao, X. Huang, X. Zhang, Turing instability and pattern formation of neural networks with reaction-diffusion terms, Nonlinear Dynamics 76 (1) (2014) 115–124.
  * [24] E. M. Izhikevich, Simple model of spiking neurons, IEEE Transactions on Neural Networks 14 (6) (2003) 1569–1572.
  * [25] E. M. Izhikevich, Which model to use for cortical spiking neurons?, IEEE Transactions on Neural Networks 15 (5) (2004) 1063–1070.
  * [26] W. W. Teka, R. K. Upadhyay, A. Mondal, Spiking and bursting patterns of fractional-order izhikevich model, Communications in Nonlinear Science and Numerical Simulation 56 (2018) 161–176.
  * [27] X. Tang, Y. Song, Bifurcation analysis and turing instability in a diffusive predator-prey model with herd behavior and hyperbolic mortality, Chaos, Solitons & Fractals 81 (2015) 303–314.
  * [28] X. Tang, J. Li, Chemotaxis induced Turing bifurcation in a partly diffusive bacterial and viral diseases propagation model, Applied Mathematics Letters 100 (2020) 106037.
  * [29] A. M. Turing, The chemical basis of morphogenesis, Bulletin of Mathematical Biology 52 (1-2) (1990) 153–197.
  * [30] Q. Ouyang, H. L. Swinney, Transition from a uniform state to hexagonal and striped Turing patterns, Nature 352 (6336) (1991) 610–612.
  * [31] S. Kondo, T. Miura, Reaction-diffusion model as a framework for understanding biological pattern formation, Science 329 (5999) (2010) 1616–1620.
  * [32] F. Sain, H. Riecke, Instabilities and spatio-temporal chaos in hexagon patterns with rotation, Physica D: Nonlinear Phenomena 144 (1-2) (2000) 124–141.
  * [33] S. Yuan, C. Xu, T. Zhang, Spatial dynamics in a predator-prey model with herd behavior, Chaos: An Interdisciplinary Journal of Nonlinear Science 23 (3) (2013) 033102.
  * [34] T. Zhang, Y. Xing, H. Zang, M. Han, Spatio-temporal dynamics of a reaction-diffusion system for a predator–prey model with hyperbolic mortality, Nonlinear Dynamics 78 (1) (2014) 265–277.
  * [35] H. Mhatre, A. Gorchetchnikov, S. Grossberg, Grid cell hexagonal patterns formed by fast self-organized learning within entorhinal cortex, Hippocampus 22 (2) (2012) 320–334.
  * [36] A. Dhooge, W. Govaerts, Y. A. Kuznetsov, H. G. E. Meijer, B. Sautois, New features of the software MatCont for bifurcation analysis of dynamical systems, MCMDS 14 (2) (2008) 147–175.
  * [37] S. Yan, D. Jia, T. Zhang, S. Yuan, Pattern dynamics in a diffusive predator-prey model with hunting cooperations, Chaos, Solitons & Fractals 130 (2020) 109428.
  * [38] C. Liu, J. Wang, H. Yu, B. Deng, X. Wei, H. Li, K. A. Loparo, C. Fietkiewicz, Dynamical analysis of parkinsonian state emulated by hybrid izhikevich neuron models, Communications in Nonlinear Science and Numerical Simulation 28 (1-3) (2015) 10–26.