# Automatic Identification of MHD Modes in Magnetic Fluctuations Spectrograms
using Deep Learning Techniques

A. Bustos<EMAIL_ADDRESS>E. Ascasíbar  A. Cappa  R. Mayo-García
Department of Technology - CIEMAT, Av. Complutense 40, 28040 Madrid, Spain
National Fusion Laboratory - CIEMAT, Av. Complutense 40, 28040 Madrid, Spain

###### Abstract

The control and mitigation of MHD oscillations modes is an issue in fusion
science because they can contribute to the outward particle/energy flux and
can drive the device away from ignition conditions. It is then of general
interest to extract the mode information from large experimental databases in
a fast and reliable way. We present a software tool based on Deep Learning
that can identify these oscillations modes taking Mirnov coil spectrograms as
input data. It uses Convolutional Neural Networks that we trained with
manually annotated spectrograms from the TJ-II stellarator database. We have
tested several detector architectures, resulting in a detector AUC score of
0.99 on the test set. Finally, it is applied to find MHD modes in our
spectrograms to show how this new software tool can be used to mine large
databases.

###### keywords:

MHD Instabilities , Mirnov coil , Deep Learning , Data Mining

##  1 Introduction

The confinement of energetic/fast ions is one of the major issues in fusion
science, since those ions, either alphas produced by the fusion reactions or
proton/deuterons produced or accelerated by the different heating systems, are
the main source of plasma heating and current drive in reactor scenarios.
Their dynamic is strongly influenced by the magnetic field because the radial
particle drift depends on the product of the ion energy and other functions of
the magnetic field [1, 2]. Unfortunately, the presence of fast ions in a
magnetized plasma can destabilize magneto-hydrodynamic [MHD] oscillations that
interact with them and can expel them from the plasma bulk. This leads to an
important outward heat flux, either towards the vacuum vessel or any other
hardware in the device. Apart from the energy loss that drives the system away
from ignition conditions, the device itself can suffer physical damage and the
plasma may be polluted with impurities caused by sputtering processes. Among
the several types of MHD instabilities, one which is easily destabilized in
typical toroidal devices is the well-known shear Alfvén wave [3]. It is indeed
of great interest to understand both theoretically and experimentally the MHD
activity in order to control or mitigate the fast ion losses. Mirnov coils [4]
are diagnostics that can measure magnetic field fluctuations associated to the
Alfvén modes in specific spatial locations with a high sampling rate ($>1$
MHz).

A quantitative experimental study of these MHD oscillations requires the
identification of series of triplets $(t_{i},f_{i},s_{i})$ that define one of
those modes, being $i$ an integer index that runs though the mode lifetime.
Here, $t_{i}$ stands for time, $f_{i}$ is the dominant frequency at $t_{i}$
and $s_{i}$ is the mode intensity. This information is enough to characterize
the oscillation mode and to relate it to other plasma parameters. We have
attempted to extract this data from a few Mirnov coil spectrograms by means of
two classical techniques: standard image thresholding and Otsu’s thresholding
method. Both methods take a grayscale image with pixel intensity in the range
from 0 to 255 and produce a binary black and white image. Those white pixels
are identified as important oscillations in the spectrogram and the black
pixels to the background. In standard thresholding one simply defines a global
threshold to the pixel intensity and classifies the pixels in black or white
according ot this value. Otsu’s thresholding method is more sophisticated [5].
It applies a clustering method to the pixel intensity distribution to divide
it in two subsets in such way that the sum of the variances of each subset is
minimum.

In Fig. 1 we depict the application of these techniques to five spectrograms
taken from the TJ-II database up to 500 kHz. Each row of images corresponds to
one spectrogram, which is shown on the first column. The second column shows
the result of appliying standard thresholding with a threshold value of 127.
In the third and forth columns we apply Otsu’s method to the original
spectrogram and to the spectrogram smoothed by a 5x5 Gaussian filter. Finally,
in the last column we show the results of the mode identification obtained
with the method presented in this paper (Image Segmentator). Clearly the
classical methods produce very noisy results and they identify as modes a
large fraction of the background noise. One may think a possible improvement
would be to cut the spectrograms in the low frequency range where the noise is
more intense, but unfortunately some modes also live in this frequency range
(spectrograms d) and e)). The difficulties to establish objective conditions
to detect the modes and apply them to a large number of spectrograms motivates
us to try to solve this problem using Artificial Intelligence.

Figure 1: Application of standard thresholding methods to five different
Mirnov coil spectrograms.

Thus, we have developed a software tool based on Deep Learning that
automatically identifies modes in a Mirnov coil spectrogram. It is designed to
take as input a spectrogram (RGB image, for reasons explained in the next
sections) and produce a B/W image with the same resolution, where the white
pixels indicate the existence of MHD activity. In the context of supervised
learning, we have annotated a dataset of 500 spectrograms and trained a neural
network to identify the interesting modes in those images. To illustrate the
potential of the detector, we apply it to the whole dataset and build a
database of modes. Then we can obtain information of the distribution of
physical features of the population of modes, such as their frequency values,
duration, intensity, etc. The physics underlying a dataset of spectrograms can
now be faster derived, being possible to quickly and automatically identify
the modes characteristics and behavior. This is of utmost importance when the
analysis has to be made with a big dataset, like the TJ-II stellarator
database, that amounts up to more than 50,000 discharges. Moreover, the
methodology proposed in this work can be easily adapted to any other dataset
of spectrograms from any other fusion device worldwide.

The reminder of this paper is organized as follows: in Section 2 we briefly
review some basic concepts of computer vision and Deep Learning, explaining
the detection methodology and how we created the spectrograms dataset. In
Section 3 we describe the way the neural network is trained to perform the
mode identifications and in Section 4 we apply it to the whole dataset to
illustrate its applications. Finally, in Section 5 we present our conclusions
and next steps in this line of research.

##  2 Methodology

Deep Learning [6] comprises modern numerical tools and algorithms based on
deep neural networks. It is used in disparate fields to automatize human work.
Relevant examples of the use of Deep Learning are autonomous driving, face
recognition, recommender systems, anomaly detection, and treatment of medical
data among others. In the last years it has proven to be very useful in those
tasks and more, and can produce excellent results in problems that are
difficult to code or express mathematically. This is clearly the present case
because, as we mentioned in Sec. 1, it is very difficult to define objective
criteria of what regions of the spectrogram correspond to oscillation modes.

Deep Learning uses artificial neural networks with diverse architectures and
depths depending on the dataset and task assigned. There are usually several
millions of parameters to adjust in the neural network. In the paradigm of
supervised learning, the neural network uses a dataset of input and labeled
output events to minimize a cost function and then predict the output of
unlabeled data. The fit of the neural network parameters is known as the
_learning phase_ , an iterative procedure with high computational cost. The
proliferation of Graphical Processing Units [GPU’s] for scientific computation
has shortened drastically the training time for most Deep Learning problems,
spreading its use and pushing the R&D in new architectures and algorithms.

Starting from a Mirnov coil spectrogram, the mode detection problem on the
image can be treated as a pixel classification problem, where the probability
of being part of a mode is assigned to each pixel of the spectrogram. This is
compatible with an image segmentator architecture.

###  2.1 Image Segmentator

An image segmentator is a kind of neural network with an encoder-decoder
architecture that classifies the pixels of an input image in a predefined set
of classes. The encoder part takes the image and, though a set of
convolutional and pooling layers, reduces the image to a latent representation
in a lower dimension space. Recall that an RGB image has $3\times N_{x}\times
N_{y}$ dimensions, being $N_{x}$ and $N_{y}$ the image resolution in pixels.
The decoder part of the segmentator takes the latent representation and
creates an output image with the same resolution using multiple
deconvolutional and dense layers. Finally, the segmentator produces a series
of $N_{x}\times N_{y}$ arrays equal to the number of classes. Each matrix
element represents the probability that a pixel corresponds to a class in the
input image. This last part of the segmentator depends on the number of
classes we want to classify. For this work, a binary classificator is used,
since we just have two pixel classes: mode and no-mode. A comprehensive review
of these techniques and neural networks can be found in [7].

As we mentioned before, image segmentators have many numerical parameters
(weights), in the order of $10^{7}-10^{8}$. If we start from a random
initialization, we would require a huge annotated dataset of images and great
computing capabilities. Fortunately, it is possible to obtain the pre-trained
weights of several encoder architectures, saving time and resources. We have
tried three different encoders that are well established and known in the
field of computer vision, namely VGG [8], MobileNet [9] and Vanilla (alias for
simple, direct CNN encoder without fancy features).VGG and MobileNet are pre-
trained with large image databases of tens of thousands of images, whose
weights can be downloaded from public repositories. On the other hand, the
simpler Vanilla encoder has to be trained from random weights initialization.
The decoder part, responsible for the pixel classification, must be trained
with our spectrogram dataset. We have used two different decoders: FCN8 and
FCN32. The weights of the decoder part are more task-dependent so they are
randomly initialized. The number of free parameters to fit the neural network
depends mostly on the encoder, being around 70, 124, and 226 million
parameters for Vanilla, VGG and MobileNet respectively.

An implementation (in Python) of those networks can be found in [10] and it
was adapted to our purposes. All computations have been done using NVIDIA
Tesla P100 GPUs with 3,584 CUDA cores and 12 GB of memory located in the ACME
cluster at CIEMAT. The library versions are Tensorflow-gpu 1.15, CUDA 10.0 and
cuDNN 7.

###  2.2 Dataset Gathering

We have applied this mode detection method to the magnetic fluctuations
present in the TJ-II stellarator [11]. TJ-II is a medium size flexible heliac
with major radius equal to 1.5 m and minor radius of the order of 0.25 m,
depending on the magnetic configuration. TJ-II has several diagnostics that
measure plasma density, temperature, voltage, energy contained, magnetic
fluctuations, impurity concentration …in its database of $\sim 50,000$
discharges.

One of the TJ-II Mirnov coils, labeled MIR5C, has been extensively used in the
present study due to its particularly good signal-to-noise ratio. It is
located inside the vacuum vessel, around 4 cm away from the last close flux
surface of the plasma and protected by a thin-walled stainless steel tube. It
is positioned in a $\varphi$ = 298 degrees toroidal plane and it is oriented
so as to measure primarily the poloidal component of the magnetic field
created by the plasma. The corresponding measured voltage is digitized with 1
MHz sampling rate [12].

We firstly choose a dataset of $\sim 500$ TJ-II discharges to calculate the
Mirnov coil spectrograms, annotate and train our mode detector. We have made
sure that it contains images of a variety of TJ-II scenarios: discharges with
known Alfvén activity, discharges with no activity, discharges with different
heating methods and failed discharges. Failed discharges are discharges that
did not lead to a viable plasma for experimenting because of any reason. Thus
we made our dataset representative of the TJ-II database in terms of magnetic
fluctuations. The time interval used to calculate the spectrogram from the
Mirnov coil signal is 1 ms, more than two orders of magnitude shorter than the
typical TJ-II discharge and also shorter than the typical mode frequency
evolution times. With this we counted on 256 Fourier modes reaching up to 500
kHz, high enough to capture the typical MHD oscillation frequencies.

We then had to annotate, i.e. create the ground truth of all the 500
spectrograms to identify the modes activity and train the neural network. The
annotation process was done manually using any simple picture editor, freely
available in all operative systems. The intention was that the image
segmentator would mark and identify regions of the spectrograms in a similar
way as the humans do. We have annotated the images in a rough way, marking
approximately the surface of the image that contains a mode that we considered
interesting. This is a subjective task and probably two people would not
annotate an image in the same way. We intentionally tended to ignore MHD
activity in the low frequency region of the spectrogram because we focused on
the identification of shear Alfvén waves which usually appear above 50 kHz for
the typical plasma density, current and magnetic geometry of the TJ-II
stellarator.

All modes/structures that appeared in each spectrogram were annotated, not
only the Alfvén modes because it is not possible to identify them only with
the spectrogram information. The Alfvén modes will be identified later in a
post processing stage. We must mention that the annotation has been done by
all the authors, in order to introduce diversity and reduce human biases as
much as possible. This indeed would make the mode detector more flexible and
robust.

Finally we ended up with a dataset of 500 pairs of images. Each pair contains
an RGB spectrogram of the Mirnov coil and a B/W image that represents the
ground truth. In Fig. 2 we can see four examples of the dataset.

Figure 2: Four examples of manually annotated spectrograms. Notice the lack of
annotations in the last picture because no mode is visible on the spectrogram.

###  2.3 Performance Measurements

In order to compare between several segmentator models and to give an
objective performance indicator we measure the performance of our segmentator
with two different quantities: the Intersection over Union [IoU] and the Area
Under the receiver operating characteristic Curve [AUC]. Here we must
emphasize that our dataset is unbalanced in terms of the classification
classes, and that the performance metrics must consider this. Some other
indicators, like the standard categorical accuracy, may lead to confusing
conclusions.

The IoU, also known as Jaccard Index, is a measure of similarity between two
images and takes values in $[0,1]$. The two images compared are in our case
the ground truth and the prediction of an image in the test set. IoU is
defined as the area of overlap divided by the area of union of the regions of
the image classified as modes. Then we can quantify the goodness of the mode
detector with the average of all IoU’s of all events in the test set.

The AUC is a more sophisticated measure [13], also in the $[0,1]$ range. The
receiver operating characteristic curve is created by plotting the true
positive rate [TPR] versus the false positive rate [FPR] at various detection
threshold settings. Here we understand as positive a pixel that belongs to a
mode and negative otherwise). This threshold is the minimum probability that
the segmentator must give to classify this pixel as positive. It is usually
fixed at 0.5, but if we vary this number between 0 and 1 we obtain a curve of
(TPR, FPR) points. The Area Under the Curve [AUC] is shown to be a good
measure of the segmentator quality. Values of AUC $\sim 0.95$ are generally
considered good enough, while AUC $\sim 0.99$ would be excellent.

##  3 Training the Model

The standard procedure in machine learning projects is to divide the dataset
in three subsets (train, validation and test sets) in order to separate the
learning and evaluation processes with objectivity and statistical
independence. The whole dataset is divided randomly so the subsets have the
same distribution, and their relative size is taken to be 70/15/15 for train,
validation and test.

The training set is used to iteratively adjust the weights of the image
segmentator so the loss function is minimized. We chose the categorical cross-
entropy of a binary classificator as the loss function. All along the learning
phase, the final model is the one model that minimizes the loss function on
the validation set. Then, we used the test set to evaluate the model using the
aforementioned performance metrics.

We trained the six combinations of encoders and decoders and compare them in
terms of the metrics IoU and AUC. We used an ADAM optimizer [14] to avoid
getting stuck in local minima of the loss function. Around 200 different
trainings are performed, varying not only the network architecture but also
other numerical parameters as the learning rate [$lr$], batch size [$bs$] and
the number of epochs. The learning rate is a parameter that controls the
amount that the weights are updated during the training phase.The batch size
is the number of images that are backpropagated simultaneously through the
network to smooth the loss function gradiant. One epoch denotes a complete
pass of the training set.

As an example of the convergence studies executed, in Table 1 we compare the
six architectures fixing $lr=0.0005,\,bs=8$ and training for 500 epochs. The
best model corresponds to a VGG encoder with FCN8 decoder. This is indeed the
best model found in all convergence studies, so it is used in the reminder of
this paper. We also show some additional information: the training time and
the categorical accuracy of the detector.

Encoder | Decoder | Tr. time [min] | IoU | AUC | Cat. Acc.
---|---|---|---|---|---
Vanilla | FCN32 | 100 | 0.233 | 0.976 | 0.989
VGG | FCN32 | 193 | 0.263 | 0.979 | 0.989
MobileNet | FCN32 | 391 | 0.260 | 0.973 | 0.988
Vanilla | FCN8 | 99 | 0.401 | 0.987 | 0.991
VGG | FCN8 | 187 | 0.427 | 0.991 | 0.991
MobileNet | FCN8 | 388 | 0.392 | 0.988 | 0.990
Table 1:  Training time and performance metrics for several image segmentator
architectures ($lr=0.0005,\,bs=8,$ num. epochs=500). We can see here that the
Categorical Accuracy is not the best performance indicator in this unbalanced
dataset because it does not vary much and it is very close to 1 in comparison
with the AUC and specially IoU.

Fig 3 shows the learning curves of the loss functions in the training and
validation sets as a function of the epoch during the learning stage. We can
see that the 500 epochs are almost enough to stabilize and that the loss in
the validation set does not increase as the train loss decreases. This means
that the training time is long enough and that no overfitting is present.
Training for longer time (2,000 epochs) does not give significant improvement.

The normalized confusion matrix at the pixel level is plotted in Fig. 4. It
provides more detailed information on the classification of the pixels.
Although an IoU = 0.427 and a True Positive Ratio = 0.579 may look small, we
will see in the next section that, with some data post processing, we have
built an effective mode detector. The reason for this is that the mode
detector produces noisy data, which must be smoothed and filtered before
identifying a region of the spectrogram as a mode. In Fig. 5 we can see some
examples of these noisy images.

Figure 3: Learning curves (loss function VS epoch) of the training and
validation sets. Figure 4: Normalized confusion matrix of the mode detector at
the pixel level. Figure 5: Input, ground truth and raw output of the mode
detector for shots #32912 (left) and #38940 (right) in the test set. Noisy
data appears in some regions and it must be cleaned.

##  4 Massive Analysis of Spectrograms and Database Creation

In this section we describe the use of all the 500 spectrograms of the dataset
as input images on the trained segmentator, the post processing of its output
and the application of descriptive statistics to characterize the oscillations
modes and locate the Alfvén activity. This post processing stage, including
the feed forward of the neural network on the complete dataset, takes around
60 minutes on a common desktop computer with no GPU.

###  4.1 Post processing

We first apply a blob detector to the segmentator output image (implemented in
the opencv library). This identifies connected regions of the spectrogram with
activity.In order to remove noise and find relevant modes, we disregarded all
blobs shorter than certain minimum duration, say 10 ms. Then, for each
mode/blob we calculated a list of time-frequency $(t_{i},f_{i})$ pairs that
represents it. Fixing a time interval in $t_{i}$, we calculated $f_{i}$ as the
average frequency weighted with the spectrogram values. This makes the
frequency estimator very robust with the possible fluctuations of the
upper/lower borders of the blob.

At this point we had a numerical description of every individual mode as a
frequency time series. This is the first time that this information was
extracted from the TJ-II database and can be used for multiple purposes. For
example, we can study the Alfvenic character of all modes in the database
using the correlation coefficient of the mode frequency with $n^{-1/2}$.

###  4.2 Dataset analysis

Once we had a functional mode detector, we applied it to all discharges we
have to show its features and possibilities. First we show some examples of
the detection of Alfvén modes using the correlation coefficient. In Fig. 6 we
can see three discharges with modes with correlation coefficient $\gtrapprox
0.9$. There are four plots of every TJ-II shot. Starting on the upper left
corner and following clockwise, we plot the spectrogram with the identified
modes; density and current time traces, and the evolution of the mode
frequencies; scatter plot of the mode frequency versus $n^{-1/2}$ and the
corresponding linear correlation coefficient; and finally a scatter plot of
the mode frequency and the plasma current.

Figure 6: Examples of Alfvén modes detected with the image segmentator after
some post processing. Since the plasma current is small and evolves slowly,
the $f_{\mathrm{Alfven}}\sim n^{-1/2}$ condition is nicely fulfilled.

Although in most detected modes the mentioned correlation coefficient is high,
we can also find some cases with MHD activity with frequency over 50 kHz in
which the plasma current changes abruptly and modifies the Shear Alfvén
Spectrum due to its impact on the rotational transform. In those cases the
direct comparison of frequency with $n^{-1/2}$ is misleading because the
plasma current variation must be taken into account [15]. Some examples of
this type of modes are shown in Fig. 7.

Figure 7: Examples in which, in spite the modes living in the Alfvenic
frequency range, the $f_{\mathrm{Alfven}}\sim n^{-1/2}$ condition is not
fulfilled without taking into account the strong evolution of the plasma
current along the discharge.

All mode information is stored in a non-relational database (MongoDB) for
further queries and analysis. In this way, any user can filter out the modes
according to his criteria and interests. For example, we can obtain the
distribution of several numerical features of the modes that were not
previously calculated in TJ-II.

In Fig. 8 we plot some data mining results obtained from the whole dataset of
detected modes. We can see a histogram of the duration of the 989 detected
modes (a). Most modes last for a few tens of milliseconds, with some of them
lasting more than 100 ms (recall that we neglect any mode shorter than 10 ms).
A similar histogram (b) showing all detected frequencies. Most modes
concentrate in frequencies between 100 and 350 kHz, with peaks around 130 and
250 kHz, are shown. In (c), the histograms of the correlation coefficient
between the frequency and $n^{-1/2}$ for each mode. This is very peaked around
1, meaning that most of the modes detected have Alfvén features and relatively
small and slow-evolving plasma current. In (d) we plot the distribution of the
mode duration restricted to modes with correlation coefficient larger than
0.9. Those examples show the potential of our method to extract information
from a big dataset of spectrograms and identify experimental features and/or
systematic behavior of the MHD oscillations.

(a)

(b)

(c)

(d)

Figure 8: Examples of data mining results: a) Distribution of the mode
duration in the dataset b) Distribution of the mode frequencies in the dataset
c)Distribution of the correlation coefficient between the mode frequency and
$n^{-1/2}$ and d) Distribution of the mode duration restricted to modes with
correlation coefficient between the frequency and $n^{-1/2}$ larger than 0.9.

###  4.3 Detector sensitivity

The sensitivity of the mode detector can be adjusted for future applications.
Recall that, given an input spectrogram, our binary segmentator estimates the
probability of each pixel to belong to a mode. The default value of the
acceptance probability is 0.5, but we can reduce it and increase the
sensitivity (and thus increase the number of false positives). If we increase
this threshold we change the behavior and operate the segmentator in a more
conservative way, rising the number of false negatives. A visual inspection of
some examples of the dataset can help to adjust this value depending on the
interests of the user.

In Fig. 9 we show four predictions of the same spectrogram using four
different sensitivity thresholds. We can see that as we decrease the
threshold, the area identified as mode increases and vice-versa. It is
important to notice that there is no need to train again the neural network to
use this functionality.

Figure 9: With thresholds larger than 0.5 we become more strict and
conservative, reducing the number of false positives but missing existing
modes. When this threshold is smaller, we increase the sensitivity of the
detector at the prize of, possibly, increasing the number of false positives.

##  5 Conclusions

In this work we have developed an analysis tool based on Deep Learning to
identify oscillation modes in Mirnov coil spectrograms measured in confined
plasma discharges. The target is to find series of frequency and time pairs
that locate oscillations modes in a spectrogram. We have trained such tool
with a manually annotated dataset of 500 spectrograms from the TJ-II
stellarator database. The detector turns out to work well, giving AUC scores
larger than 0.99 (excellent) and IoU larger than 0.4 (enough for our
purposes). To the authors’ knowledge, this is the first time that such
techniques are applied to this kind of dataset, and it has turned out to be a
valid tool that can save time to the experimental plasma physicists.

We can use it to effectively mine any set of discharges in TJ-II, creating a
database of modes with their physical characteristics. Once the segmentator is
trained, the database creation takes $\approx 7$ seconds per spectrogram on a
desktop computer. Then this information can be used to identify patterns,
physical behavior or check hypothesis.

It is worth to mention that the software developed and methods employed can be
applied to other facilities and diagnostics that produce images and/or
spectrograms like bolometry, Langmuir probes, heavy ion beam probe, etc.
Transfer learning techniques can be used here to speed up new trainings.

For the sake of completeness we also give some of the few examples where the
mode detector does not work properly. We observe in Fig. 10 modes that are
misdetected and modes that are merged by the detector and identified as a
single mode. Some of the mentioned errors might be avoided simply by
increasing the dataset size and train the neural network again. Usually, the
neural networks increase their performance with the dataset size, as opposed
to other simpler machine learning models that saturate after certain dataset
size. The problem with modes that crosses each other can be solved changing
the neural network architecture so the segmentator is able to resolve the
instance map, i.e. identify distinct modes. But this requires new annotation
and training processes and it is left for future work.

Finally we can say that, overall, the mode detector gives reasonable good
results and it is a potential tool for automatic analysis of spectrograms.

Figure 10: Some errors that sometimes appear in the mode detection. Starting
from the upper-left picture and following clockwise, we can see that the
segmentator missed parts of a long mode, missed a whole mode, merged many
modes in a very complex shape and detected ghost modes in a failed discharge.

## Acknowledgments

This work was partially funded by the Spanish Ministry of Science and
Innovation under contract number FIS2017-88892-P and
RTI2018-096006-B-I00(CODEC-OSE), and by the Comunidad de Madrid CABAHLA
project (S2018/TCS-4423). Simulations carried out made also use of the
computing facilities provided by the CYTED Network RICAP (517RT0529).

## References

  * [1] P. Helander and D. J. Sigmar, Collisional Transport in Magnetized Plasmas (Cambridge University Press, Cambridge, 2001).
  * [2] R. Balescu, Transport Processes in Plasmas: Neoclassical Transport Theory (Elsevier Science Ltd, Netherlands, 1988).
  * [3] G.Vlad, F.Zonca, and S.Briguglio, Rivista del Nuovo Cimento 22, (1999).
  * [4] S. Mirnov and I. Semenov, Soviet Physics JETP 33, 6, (1971).
  * [5] N. Otsu, IEEE Transactions on Systems, Man and Cybernetics 9, 62 (1979).
  * [6] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning (MIT Press, USA, 2016), http://www.deeplearningbook.org.
  * [7] S. Minaee et al., ArXiv abs/2001.05566, (2020).
  * [8] K. Simonyan and A. Zisserman, CoRR abs/1409.1556, (2015).
  * [9] A. G. Howard et al., (2017), cite arxiv:1704.04861.
  * [10] D. Gupta and R. J. Wala, Image Segmentation Keras : Implementation of Segnet, FCN, UNet, PSPNet and other models in Keras., https://github.com/divamgupta/image-segmentation-keras, 2020.
  * [11] R. Jimenez-Gomez et al., Fusion Sci. Technol. 51, 20 (2007).
  * [12] R. Jiménez-Gómez et al., Nuclear Fusion 51, 033001 (2011).
  * [13] A. P. Bradley, Pattern Recognition 30, 1145 (1997).
  * [14] D. P. Kingma and J. Ba, 3rd International Conference for Learning Representations (2015), cite arxiv:1412.6980.
  * [15] B. Sun, M. Ochando, and D. López-Bruna, Nuclear Fusion 55, 093023 (2015).


# Spatiotemporal characteristics in systems of diffusively coupled excitable
slow-fast FitzHugh-Rinzel dynamical neurons

Arnab Mondal  Department of Mathematics and Computing, Indian Institute of
Technology (Indian School of Mines), Dhanbad 826004, India  Argha Mondal
Corresponding author<EMAIL_ADDRESS>School of Engineering, Amrita
Vishwa Vidyapeetham, Amritapuri, Kollam 690525, India  Department of
Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester CO4 3SQ,
UK  Sanjeev Kumar Sharma  Department of Mathematics and Computing, Indian
Institute of Technology (Indian School of Mines), Dhanbad 826004, India
Ranjit Kumar Upadhyay  Department of Mathematics and Computing, Indian
Institute of Technology (Indian School of Mines), Dhanbad 826004, India  Chris
G. Antonopoulos  Department of Mathematical Sciences, University of Essex,
Wivenhoe Park, Colchester CO4 3SQ, UK

###### Abstract

In this paper, we study an excitable, biophysical system that supports wave
propagation of nerve impulses. We consider a slow-fast, FitzHugh-Rinzel neuron
model where only the membrane voltage interacts diffusively, giving rise to
the formation of spatiotemporal patterns. We focus on local, nonlinear
excitations and diverse neural responses in an excitable 1- and 2-dimensional
configuration of diffusively coupled FitzHugh-Rinzel neurons. The study of the
emerging spatiotemporal patterns is essential in understanding the working
mechanism in different brain areas. We derive analytically the coefficients of
the amplitude equations in the vicinity of Hopf bifurcations and characterize
various patterns, including spirals exhibiting complex geometric
substructures. Further, we derive analytically the condition for the
development of antispirals in the neighborhood of the bifurcation point. The
emergence of broken target waves can be observed to form spiral-like profiles.
The spatial dynamics of the excitable system exhibits 2- and multi-arm spirals
for small diffusive couplings. Our results reveal a multitude of neural
excitabilities and possible conditions for the emergence of spiral-wave
formation. Finally, we show that the coupled excitable systems with different
firing characteristics, participate in a collective behavior that may
contribute significantly to irregular neural dynamics.

Keywords: FitzHugh-Rinzel neurons, Diffusion, Spatiotemporal characteristics,
Amplitude equations, Hopf bifurcation, Ginzburg-Landau equation,
Synchronization, Antispirals, Target waves

In this paper, we study the slow-fast FitzHugh-Rinzel biophysical excitable
neuron model. We demonstrate the appearance of diverse dynamical behavior
depending on spike-bursting, firing responses. We also discuss analytically
and numerically the bifurcation analysis. To understand the effects of
diffusion on spatiotemporal pattern formation in coupled FitzHugh-Rinzel
neurons, we consider a FitzHugh-Rinzel system where only the membrane voltage
interacts diffusively. We compute the synchronization index to study the
collective behavior of the coupled systems and describe various patterns,
including target waves, multi-arm spirals and antispirals. Finally, we derive
analytically the amplitude equations to validate the existence of antispiral
patterns in the vicinity of Hopf bifurcations.

##  1 Introduction

The complex dynamics of excitable neurons is related to ionic concentrations
inside and outside neural membranes and play a major role in generating action
potentials. This provides a better way to understand brain functioning in
normal and pathophysiological states [1, 2, 3, 4, 5, 6]. Neurons receive
incoming sensory inputs, encode them into different biophysical variables and
produce relevant outputs [3]. Biophysical mechanisms are dynamical in nature
and can be studied using concepts from dynamical systems theory. In extended
biophysical systems of coupled excitable cells, neural populations produce a
transmembrane potential difference that travels across neurons by means of
wave propagation [7, 8, 9, 6]. This may cause the appearance of synchronous
activity in groups of neighboring neurons. Such types of synchronized behavior
are playing a major role in signal processing in neural populations [10].

Here, we consider a slow-fast, FitzHugh-Rinzel (FHR) neuron model where only
the membrane voltage interacts diffusively, giving rise to the formation of
spatiotemporal patterns. The FHR model exhibits diverse bursting activities.
In this context, bursting refers to rapid changes in membrane voltage
oscillations which involve simultaneous changes between active phases (high
amplitude oscillations or bursts of spikes) and silent phases (low
amplitude/subthreshold oscillations or amplitude death, i.e., quiescence) [2].
We study a dynamical system of three coupled partial differential equations
(PDEs) for 2-dimensional spatial pattern formation based on the FHR model,
where the recovery and/or slow variables (slow-fast dynamics) interact with
the membrane voltage variable. We seek to explore spatial instabilities
considering 1- and 2-dimensional diffusion processes and the nature of
spatiotemporal patterns. It is the symmetry breaking nature of spatial systems
that gives rise to Turing-like structures [11]. Only one spatial variable,
i.e., the membrane potential, determines the scale of the patterns observed.
In our system, there exists only one spatially distributed scale, and thus the
spatial instability occurs due to the dynamics rather than due to differences
in spatial scales [12].

The mathematical analysis for the emergence of spatial structures is important
to understand a wide range of biophysical and pathological phenomena [13, 14,
15, 16, 17, 18]. The work here is motivated by our earlier work in [12] and
other previously studied diffusively coupled biophysical excitable systems,
such as those in [19, 20, 21, 22, 23, 24, 25, 26]. To the best of our
knowledge, a clear analytical study describing the dynamics of a diffusively
coupled, slow-fast, neuron model in which only the membrane voltage is
spatially distributed, has not yet been deeply explored with respect to
pattern formation and emergence of spirals. Often biophysical phenomena
observed in experiments can be reproduced by solving numerically a nonlinear
system of coupled PDEs that describes them. However, the underlying dynamical
behavior of real biophysical systems generally involves a large number of
variables.

The excitable FHR model studied here is an extended version of the FitzHugh-
Nagumo (FHN) model [27]. Excitable, biophysical media represented by the slow-
fast dynamics of electrically coupled FHR neurons [28, 29, 30, 31] give rise
to wave propagation when the system becomes excited above a threshold, called
a traveling wave, that travels through the nerve cells. The responses
propagate along the axons, allowing the consideration of the system as a
spatially distributed neural medium. The ideal biophysical nerve membrane
model that produces bursting was formulated and studied numerically in [29,
30]. The FHR model offers a wider view of the original FHN model by taking
into consideration a special type of bursting, known as elliptic bursting. As
a consequence, the FHR model with different neurocomputational properties
exhibits various responses at different levels of biophysical plausibility,
similar to real neurons.

Here, we consider diffusively coupled identical FHR neurons described by a
system of three coupled PDEs. The model of the single deterministic nerve cell
is studied theoretically and numerically. Different parameter regimes that
correspond to qualitatively different biophysical behaviors are investigated.
The self-sustained and self-organized behavior are investigated for a
localized stimulus current with different synaptic couplings. The possible
oscillatory behavior is considered and analytical conditions are derived. We
find that weak stimulus injected to the system when it resides at the stable
steady-state, produces nontrivial firing responses. Consequently, the stimulus
perturbations can be applied to the system deterministically, leading to the
emergence of different spatiotemporal phenomena.

Particularly, we explore the dynamics in a continuous, 1-dimensional piece of
spatially extended neural cable [32, 8, 33] and in a square, considered as a
2-dimensional piece of spatially extended neural tissue. In the 1-dimensional
case, neurons are coupled along a chain of excitable cells by means of a
1-dimensional diffusion term. We find suitable parameter regimes for which the
system shows different fluctuations and pattern formations in the vicinity of
two Hopf bifurcations for a certain range of coupling strengths. In the
2-dimensional case, we employ the amplitude equations and elucidate phenomena
of spatiotemporal pattern formation. Our analysis is based on the complex
Ginzburg-Landau equation (CGLE) [34, 35, 36, 37, 38] for the derivation of the
amplitude equations using multiple space and time scales. The method provides
a Taylor series expansion of the original nonlinear equations with many power
operators and depends on the expansion of the linear and nonlinear terms of a
small perturbation parameter close to the onset of instability. The
coefficients of CGLE obtained here, may have their own values which might
relate to certain dynamical behavior in the diffusion process under different
experimental conditions. We also derive analytically the stability of the
coupled system.

We report on the presence of various emerging spatiotemporal patterns for
different electrical impulses in the 2-dimensional system. The slow-fast
dynamics supports the formation of spirals due to excitabilities of the slow
time-scale in the single-neuron model. Various spiral patterns emerge in the
diffusive system when it is not in a stable, stationary state. Spiral and
target waves are known as ordered waves, which are often observed in the
extended excitable system. Generally, suitable periodic forcing plays a major
role in developing target waves [22, 39, 40, 41]. Spiral or antispiral waves
behave asymptotically as plane waves far from their cores [42, 43] and a
variety of complex patterns can manifest, such as target waves, spirals-
antispirals, hexagons, spot-stripes, mixed patterns, etc., [44, 7, 45, 46, 47,
48]. Spirals are particular types of patterns that rotate around a central
point, known as the rotor. Rotors show powerful rhythmic activity by sending
rotating-type, robust, wavy patterns outwards. The head of a stable rotating
free spiral wave moves around a circular core. The motion of the spiral’s tip
controls the dynamics of spiral waves. Local excitatory connections play
essential role in the dynamics and emergence of spirals during certain brain
functions in the cortex [49]. Antispiral waves observed near the Hopf
bifurcation and the direction of wave-propagation are determined by the
competition between waves and their surrounding bulk oscillations. The
direction of the wave-propagation (either outwards or inwards) plays a role in
characterizing spiral or antispiral waves. Spiral patterns may pinpoint to a
mechanism that changes the irregular dynamics of cortical neurons to rhythmic
behavior [50]. Particularly, the mechanism changes the frequency of
oscillations and amplitudes, as well as the spatial coherence activity. Our
results suggest that spiral waves can emerge in a 2-dimensional spatially
extended system with dynamics around a Hopf bifurcation. We report on the
formation of single- and multi-arm spiral patterns. Spiral waves emerge when
the target waves break-up and are associated with spatial heterogeneity in the
system [46].

Spatiotemporal patterns can often be observed in coupled excitable
oscillators, cardiovascular systems and in neural systems [41, 22, 51, 39, 52,
40]. It is important to understand the characteristics and formation of these
spirals. It has been observed that stable spiral patterns can emerge in an
oscillatory system around a Hopf bifurcation associated with the spatially
uniform steady-state solution. The spiral waveforms have been investigated in
bursting media in [53]. These results may be related to the formation and
development of spiral patterns in neural systems, especially in cortical areas
in the brain. These types of spatiotemporal patterns may also be related to
the characteristics of the generation of activity recorded in EEG signals, to
epileptic seizures, to migraine-related issues and to the transmission of
visual images in the cortex [41, 49, 53, 4, 54, 9]. Our results show that
diffusively coupled FHR neurons with different firing characteristics,
depending on control parameters, participate in the collective behavior of the
system.

The paper is organized as follows: In Sec. 2, we discuss the single FHR model
and perform a bifurcation analysis to reveal the multitude of dynamical
behaviors as a function of the external current. In Sec. 3, we study,
analytically and numerically, the effects of diffusion on coupled FHR neurons
arranged in a 1-dimensional configuration. We also discuss the synchronization
index that we use to quantify the collective behavior in the spatially
extended systems. In Sec. 4, we study the diffusive properties of the dynamics
of coupled FHR neurons, arranged in a 2-dimensional spatially extended mesh.
We describe analytically the diffusive properties using the amplitude equation
and present the results that reveal complex structures for positive diffusive
couplings, such as target waves, spiral patterns, and two- and multi-arm
spiral waves. Finally, in Sec. 5, we conclude our work and discuss it in the
framework of other studies in the field.

##  2 The excitable FHR model

The FHR model is the extended version of the original FHN model [27].
Incorporating slow dynamics, it can exhibit a wide range of different firing
patterns in certain parameter regimes. The FHR model was introduced by
FitzHugh and Rinzel [28, 29, 30, 31] and is also known as the elliptic
bursting model [31]. It is described by the system of ordinary differential
equations

$\displaystyle\dot{u}$ $\displaystyle=f(u)-v+w+I,$ $\displaystyle\dot{v}$
$\displaystyle=\delta(a+u-bv),$ (1) $\displaystyle\dot{w}$
$\displaystyle=\mu(c-u-w),$

where $u$ represents the membrane voltage of the nerve cell and, $v$ and $w$,
the recovery and slow variables, respectively. Particularly, $v$ provides a
slower negative feedback in the system. Here,

$f(u)=u-\frac{u^{3}}{3}$ (2)

is the cubic nonlinearity in the membrane voltage and allows for regenerative,
self-excitation via positive feedback. $I$ indicates a constant external
current stimulus and $\delta$, $a$, $b$, $\mu$ and $c$ are the parameters of
the system. Parameter $\mu$ is a small time-scale number that controls the
slow variable $w$ and $c$ plays a role similar to $a$ in the FHN model [27].
The decrease of $a$ and $c$ gives rise to longer intervals between consecutive
bursting activities, with the system showing relatively fixed times in burst
duration. Additionally, with the increase of $a$, the interburst intervals
become shorter and the oscillations change to tonic spiking [28, 31]. In our
work, we have used $\delta=0.08$, $a=0.7$, $b=0.8$, $\mu=0.002$, $c=-0.775$
and we will be varying $I$.

To understand the dynamics of system (2), we study its local and global
stability properties performing a bifurcation analysis. To this end, its
nontrivial fixed point is denoted by $E=(u_{0},v_{0},w_{0})$, where
$v_{0}=(a+u_{0})/b$ and $w_{0}=c-u_{0}$. In this framework, $u_{0}$ can be
obtained from the real solution to the cubic equation

$u_{0}^{3}+Pu_{0}+Q=0,$ (3)

where $P=3/b$ and $Q=-3(I-a/b+c)$. Its discriminant is given by
$\Delta=-4P^{3}-27Q^{2}<0$, provided that $b>0$, which implies that Eq. (3)
has only one real root. Particularly, the real solution to Eq. (3) is given by

$u_{0}=-\frac{1}{3}\bigg{(}\Delta_{2}+\frac{\Delta_{0}}{\Delta_{2}}\bigg{)},$
(4)

where $\Delta_{0}=-3P$,
$\Delta_{2}=\Bigg{(}\bigg{(}\Delta_{1}+\sqrt{\Delta_{1}^{2}-4\Delta_{0}^{3}}\bigg{)}/2\Bigg{)}^{\frac{1}{3}}$
and $\Delta_{1}=27Q$. This results to system (2) admitting the unique
nontrivial fixed point $E=(u_{0},v_{0},w_{0})$, where $u_{0}$,
$v_{0}=(a+u_{0})/b$ and $w_{0}=c-u_{0}$ are given by Eq. (4).

The Jacobian matrix, $J$, of system (2), computed at the fixed point $E$, is
given by

$J=\begin{pmatrix}a_{11}&a_{12}&a_{13}\\\ a_{21}&a_{22}&a_{23}\\\
a_{31}&a_{32}&a_{33}\end{pmatrix},$

where $a_{11}=1-u_{0}^{2}$, $a_{12}=-1$, $a_{13}=1$, $a_{21}=\delta$,
$a_{22}=-b\delta$, $a_{23}=0$, $a_{31}=-\mu$, $a_{32}=0$ and $a_{33}=-\mu$.
The characteristic equation is given by

$\lambda^{3}+a_{1}\lambda^{2}+a_{2}\lambda+a_{3}=0,$ (5)

where $a_{1}=-1+b\delta+\mu+u_{0}^{2}$,
$a_{2}=\delta-b\delta+b\mu\delta+b\delta u_{0}^{2}+\mu u_{0}^{2}$ and
$a_{3}=\mu\delta+b\mu\delta u_{0}^{2}$. Analyzing the stability of the system
in a neighborhood of the fixed point $E$ using the Routh-Hurwitz criterion
[55], implies that $a_{1}>0$, $a_{2}>0$, $a_{3}>0$ and that the second order
minor $a_{1}a_{2}-a_{3}>0$. Next, we prove a theorem that shows for which
values of $a$, $b$ and $c$, the fixed point $E$ is globally asymptotically
stable when $I=0$.

######  Theorem 2.1.

The fixed point $E=(u_{0},v_{0},w_{0})$ of system (2) is globally
asymptotically stable in the exterior of the ellipsoid

$\frac{1}{3}\bigg{(}u^{2}-\frac{3}{2}\bigg{)}^{2}+b\bigg{(}v-\frac{a}{2b}\bigg{)}^{2}+\bigg{(}w-\frac{c}{2}\bigg{)}^{2}=\frac{3b+a^{2}+bc^{2}}{4b}$

for applied current stimulus $I=0$.

###### Proof.

We consider the Lyapunov function $L:\mathbb{R}^{3}\rightarrow\mathbb{R}$

$L(u,v,w)=u^{2}+\frac{1}{\delta}v^{2}+\frac{1}{\mu}w^{2},$

where $\delta$ and $\mu$ are positive parameters appearing in system (2).
Then,

  * •

$L(u,v,w)\geq 0$ for all $(u,v,w)\in\mathbb{R}^{3}$,

  * •

$L(u,v,w)=0$ if and only if $(u,v,w)=(0,0,0)$, and

  * •

all sublevel sets of $L$ are bounded, i.e., $L(u,v,w)\rightarrow\infty\mbox{
as }(u,v,w)\rightarrow\infty$.

This shows that $L(u,v,w)$ is positive definite. Then, the time-derivative of
$L$ is given by

$\displaystyle\frac{dL(u,v,w)}{dt}$
$\displaystyle=2\bigg{(}u\dot{u}+\frac{1}{\delta}v\dot{v}+\frac{1}{\mu}w\dot{w}\bigg{)}$
$\displaystyle=-2\Bigg{(}\frac{1}{3}\bigg{(}u^{2}-\frac{3}{2}\bigg{)}^{2}+b\bigg{(}v-\frac{a}{2b}\bigg{)}^{2}+\bigg{(}w-\frac{c}{2}\bigg{)}^{2}-\frac{3b+a^{2}+bc^{2}}{4b}\Bigg{)}.$

It follows that $\frac{dL}{dt}$ is negative outside the ellipsoid

$\frac{1}{3}\bigg{(}u^{2}-\frac{3}{2}\bigg{)}^{2}+b\bigg{(}v-\frac{a}{2b}\bigg{)}^{2}+\bigg{(}w-\frac{c}{2}\bigg{)}^{2}=\frac{3b+a^{2}+bc^{2}}{4b}$

and that

$\frac{dL(0,0,0)}{dt}=0.$

Thus, the fixed point $E$ is globally asymptotically stable for $I=0$ in the
exterior of the specified ellipsoid. ∎

Following [28, 31], we consider the three sets of external currents $I$: (a)
$I=0.2$ (set 1), (b) $I=0.43$ (set 2) and (c) $I=0.5$ (set 3) around the two
Hopf bifurcation points of system (2) that we discuss next.

###  2.1 Bifurcation analysis of the single FHR model

Here, we study the bifurcation properties of the single FHR model (2)
considering the current, $I$, as the bifurcation parameter. We use MatCont, a
Matlab package for the numerical continuation and bifurcation analysis of
continuous, parameterized dynamical systems [56]. We plot the results in Fig.
1.

Particularly, for $I<0.137$, system (2) has a stable focus node (solid cyan
line), exhibiting a quiescent state (Fig. 1). The system changes its stability
and a stable limit cycle appears around $I=0.137$, which corresponds to the
supercritical Hopf bifurcation (HB1). As $I$ increases further, the steady
state appears again through a second supercritical Hopf bifurcation (HB2) at
$I=3.16298$. For $I>3.16298$, the system has a stable focus node (solid cyan
line), which corresponds to the quiescent state shown in yellow in the fourth
inset. For $0.137<I<3.16298$, the system exhibits different types of firing
patterns, shown in blue, red and green in the insets in Fig. 1, such as
elliptic-type bursting, mixed-mode oscillations and tonic spiking. With
further increase in $I$, the system shows a quiescent state, shown as the
yellow curve in Fig. 1 for $I=3.8$. In the range $0.137<I<3.16298$, the limit
cycle changes its stability through period-doubling bifurcations at
$I=0.1925$, $I=0.4683$ and $I=3.1075$. The system has an unstable focus for
$0.137<I<0.6$ and a saddle node for $0.7<I<2.6$ (dotted cyan line). Increasing
$I$ even more, the system has an unstable focus node up to $I\approx 3.16298$.
The existence of period-doubling bifurcations in the system is one of the
routes to chaos. The system is chaotic where both the fixed point and limit
cycle are unstable. As an example, we calculated the Lyapunov exponents (LEs)
for $I=0.4$, where both the fixed point and limit cycle are unstable and found
that they converge approximately to the values
$(0.000121,-0.004023,-0.522634)$, ordered in descending order. As the largest
LE is positive, it confirms the system is chaotic. Further on, in Fig. 1,
solid and cyan dotted lines denote the quiescent and oscillatory regions,
respectively. The stable and unstable limit cycles are shown in solid magenta
and dotted black lines, respectively.

Figure 1: Bifurcation diagram of the FHR system (2), where solid and dotted
cyan lines denote the quiescent and oscillatory regions, respectively. The
stable and unstable limit cycles are shown in solid magenta and dotted black
lines. Different oscillatory regimes and bifurcations as a function of the
bifurcation parameter $I$ are also shown in the upper insets, where we present
plots of the time-series $u$ as a function of time $t$ for the parameter sets
1, 2 and 3, shown in blue (elliptic-type bursting), red (mixed-mode
oscillations), green (tonic spiking) and with further increase of $I$, it
shows quiescent state. Note that HB1 and HB2 correspond to the two critical
points $I=0.137$ and $I=3.16298$, respectively, where Hopf bifurcations occur
in system (2).

In the following, we verify our numerical results with respect to the critical
bifurcation points HB1 and HB2 using an analytical approach [57]. The complex
eigenvalues of the linearized FHR system of system (2) around the fixed point
$E$ are given by $\gamma(I)=\alpha(I)\pm i\beta(I)$, where $i$ is the
imaginary unit of the complex numbers. Suppose for a certain value of $I$, say
$I\equiv I_{0}$, the following conditions are met: (i) $(I_{0})=0$, (ii)
$\beta(I_{0})=\beta^{*}\neq 0$ and (iii)
$\frac{d\alpha(I)}{dI}\Big{|}_{I=I_{0}}=\alpha^{*}\neq 0$. Then, system (2)
undergoes a Hopf bifurcation at $I\equiv I_{0}$. Conditions (i) and (ii) are
known as the non-hyperbolicity conditions and condition (iii), as the
transversality condition. To verify analytically the Hopf bifurcation at
$I_{0}$, we consider $I$ as the bifurcation parameter, with $\alpha(I)$ and
$\beta(I)$ given by

$\displaystyle\alpha(I)$
$\displaystyle=0.31133-0.33333u_{0}^{2}+(-68.6643+172.1666u_{0}-83.3333u_{0})/A^{\frac{1}{3}}-0.000333A^{\frac{1}{3}}$

and

$\beta(I)=\frac{118.9304-298.2023u_{0}^{2}+144.338u_{0}^{4}}{A^{\frac{1}{3}}}-0.000577A^{\frac{1}{3}},$

where

$\displaystyle A$
$\displaystyle=10^{5}\bigg{(}931.04315-3529.5225u_{0}^{2}+3873.75u_{0}^{4}-1250u_{0}^{6}$
$\displaystyle+10.392\sqrt{-67.14+25.33u_{0}^{2}+21.64u_{0}^{4}+606.6u_{0}^{6}-375.4u_{0}^{8}}\bigg{)},$
$\displaystyle u_{0}$ $\displaystyle=\frac{4.27494}{B}-0.292402B,$
$\displaystyle B$
$\displaystyle=\bigg{(}99-60I+1.414\sqrt{6463-5940I+1800{I^{2}}}\bigg{)}^{\frac{1}{3}}.$

Thus, to find the critical points $I_{0}$ of system (2) where Hopf
bifurcations occur, we find the values $I\equiv I_{0}$, where conditions
(i)-(iii) hold. To this end, solving the equation $\alpha(I)=0$, we obtain
$I=0.137$ and $I=3.16298$. This implies that

$\displaystyle\alpha(0.137)$ $\displaystyle=\alpha(3.16298)=0,$
$\displaystyle\beta(0.137)$ $\displaystyle=\beta(3.16298)=-0.2792,$
$\displaystyle\frac{d\alpha(I)}{dI}\bigg{|}_{I=0.137}$ $\displaystyle=0.443,$
$\displaystyle\frac{d\alpha(I)}{dI}\bigg{|}_{I=3.16298}$
$\displaystyle=-0.443.$

Consequently, system (2) undergoes two Hopf bifurcations at the critical
points $I=0.137$ (HB1) and $I=3.16298$ (HB2), which are in good agreement with
the numerical results shown in Fig. 1.

In the next section, we study the diffusive properties of coupled FHR neurons
arranged in a 1-dimensional configuration to investigate the collective
dynamics and emergence of biologically plausible patterns around the two Hopf
bifurcations.

##  3 Diffusively coupled FHR neurons arranged in a 1-dimensional
configuration

Here, we consider a 1-dimensional configuration of excitable, diffusively
coupled FHR neurons through the membrane voltage $u$. The corresponding system
is governed by the system of PDEs

$\displaystyle\frac{\partial u}{\partial t}$
$\displaystyle=f(u)-v+w+I+D\frac{\partial^{2}u}{\partial x^{2}},$
$\displaystyle\frac{\partial v}{\partial t}$ $\displaystyle=\delta(a+u-bv),$
(6) $\displaystyle\frac{\partial w}{\partial t}$ $\displaystyle=\mu(c-u-w),$

where $t$ is the time and $u(x,t=0)$, $v(x,t=0)$ and $w(x,t=0)$ for
$x\in\Psi$, the initial conditions. $\Psi$ is the space of the 1-dimensional
configuration and the function $f(u)$ is given by Eq. (2). We consider that
only the membrane voltage $u$ diffuses to neighboring neurons, modeled by the
1-dimensional diffusion term $D\frac{\partial^{2}u}{\partial x^{2}}$, where
$D\geq 0$ is the diffusion coefficient or diffusion coupling or constant
synaptic coupling strength. By 1-dimensional configuration we mean that the
dimension of $\Psi$ is 1, thus $u$, $v$ and $w$ are functions of only one
spatial coordinate, i.e., of $x$. Furthermore, we assume the zero-flux
boundary condition,

$\frac{\partial u}{\partial n}=\frac{\partial v}{\partial n}=\frac{\partial
w}{\partial n}=0$ (7)

for $x\in\partial\Psi$, where $\partial\Psi$ is the boundary of $\Psi$. In
this framework, $\frac{\partial}{\partial n}$ denotes the directional
derivative along the outward normal $n$ to the boundary $\partial\Psi$. The
reason for considering a zero-flux boundary condition is that it makes the
1-dimensional membrane impermeable at the two edges, i.e., no ions flow inside
or outside the boundaries, and thus the membrane acts as an isolated cable [8,
12]. Next, we prove an important theorem.

######  Theorem 3.1.

If the fixed point $E$ of the FHR system (2) is locally asymptotically stable,
then the fixed point of the corresponding 1-dimensional diffusive model (3) is
also locally asymptotically stable.

###### Proof.

We consider the particular solution of the linearized model

$\begin{pmatrix}u\\\ v\\\ w\end{pmatrix}=\begin{pmatrix}u_{0}\\\ v_{0}\\\
w_{0}\end{pmatrix}+\epsilon\begin{pmatrix}u_{k}\\\ v_{k}\\\
w_{k}\end{pmatrix}e^{\lambda_{k}t+ikx}+c.c.+o(\varepsilon^{2})$

of system (3), where $c.c$. stands for “complex conjugate” terms,
$\lambda_{k}$ is the wave length, $k>0$ the wave number along the $x$
direction and $x$ the directional vector in $\Psi$. The Jacobian matrix of
system (3) at the fixed point $E=(u_{0},v_{0},w_{0})$ is given by

$J=\begin{pmatrix}a_{11}-Dk^{2}&a_{12}&a_{13}\\\ a_{21}&a_{22}&a_{23}\\\
a_{31}&a_{32}&a_{33}\end{pmatrix}$

and its characteristic equation by

$\lambda_{k}^{3}+b_{1}\lambda_{k}^{2}+b_{2}\lambda_{k}+b_{3}=0.$ (8)

The coefficients of Eq. (8) are given by $b_{1}=a_{1}+Dk^{2}$,
$b_{2}=a_{2}+b\delta Dk^{2}+\mu Dk^{2}$ and $b_{3}=a_{3}+b\mu\delta Dk^{2}$,
where $a_{i}$, $i=1,2,3$ are the coefficients of Eq. (5). We analyze the
stability of the fixed point $E=(u_{0},v_{0},w_{0})$ using the Routh-Hurwitz
criterion [55], resulting in the fixed point $E$ being stable if $b_{1}>0$,
$b_{2}>0$, $b_{3}>0$ and if the second order minor $b_{1}b_{2}-b_{3}$ is
positive. If we consider the stable fixed point of system (2), then $b_{1}>0$,
$b_{2}>0$, $b_{3}>0$ and $b_{1}b_{2}-b_{3}>0$ as $a_{1}>0$, $a_{2}>0$,
$a_{3}>0$, $a_{1}a_{2}-a_{3}>0$, $D>0$ and $k^{2}>0$. Thus, the fixed point
$E$ of the 1-dimensional diffusive model (3) is locally asymptotically stable.
∎

It is worth it to note that if the fixed point $E$ of system (2) is unstable,
then the fixed point of the corresponding 1-dimensional diffusive model (3)
can be made stable by increasing appropriately the value of the diffusion
coefficient [58]. For example, if we consider set 1, the fixed point
$E=(-0.939127,-0.298909,0.164127)$ of system (2) is unstable as
$a_{1}=-0.0520405<0$, $a_{2}=0.0743373>0$, $a_{3}=0.000272891>0$ and
$a_{1}a_{2}-a_{3}=-0.00414144<0$. Interestingly, for the diffusive model (3),
the same fixed point is stable if $b_{1}=-0.0520405+Dk^{2}>0$,
$b_{2}=0.0743373+0.066Dk^{2}>0$, $b_{3}=0.000272891+0.000128Dk^{2}>0$ and
$b_{1}b_{2}-b_{3}=-0.00414144+0.0707747Dk^{2}+0.066D^{2}k^{4}>0$, i.e., when
$Dk^{2}>0.05563$, meaning that $D>0.05563/k^{2}$, where $k^{2}>0$. Thus, the
fixed point of the corresponding diffusive model (3) can be made stable by
increasing appropriately the diffusive coupling, $D$. One can prove
analytically similar results for sets 2 and 3, where the corresponding
conditions for the stability of the diffusive model (3) are given by
$Dk^{2}>0.255616$ and $Dk^{2}>0.315046$, respectively.

In Fig. 4(a), the dashed blue, brown and green curves delineate the boundaries
of the stable and unstable regions of sets 1, 2 and 3, and are given by

$\displaystyle D$ $\displaystyle=\frac{0.05563}{k^{2}},$ (9) $\displaystyle D$
$\displaystyle=\frac{0.255616}{k^{2}},$ (10) $\displaystyle D$
$\displaystyle=\frac{0.315046}{k^{2}},$ (11)

respectively. The solid blue, brown and green lines for sets 1, 2 and 3,
respectively, were derived numerically and separate the stable from the
unstable regions in system (3). Particularly, the critical $D$ values where
these lines are located, were computed using the PDEPE toolbox in Matlab and
the bisection method. The method solves system (3) numerically for different
$D$ values until it converges to the critical $D$ at which the transition from
unstable to quiescent (i.e., stable) dynamics occurs. For each $D$, we checked
visually the dynamics based on plots such as those in Fig. 2. The method
resulted at $D=6$ for set 1 (solid blue line), at $D=0.98$ for set 2 (solid
brown line) and at $D=0.9$ for set 3 (solid green line).

Figure 2: Spatiotemporal plots of the system of coupled FHR neurons (3)
arranged in a 1-dimensional configuration for different values of the
diffusion coefficient, $D$. Panels (a), (b) and (c) are for set 1 with
$D=0.0001$, $D=1$ and $D=8$, respectively. Panels (d), (e) and (f) are for set
2 with $D=0.0001$, $D=0.1$ and $D=1$, respectively, and panels (g), (h) and
(i) are for set 3 with $D=0.0001$, $D=0.01$ and $D=1$, respectively. The color
bars encode the values of the membrane voltages $u(x,t)$ of diffusively
coupled, FHR neurons in system (3). Figure 3: Temporal evolution of
arbitrarily chosen neurons of the system of coupled FHR neurons (3) arranged
in a 1-dimensional configuration for different values of the diffusion
coefficient, $D$. The values of $D$ and sets of parameters are the same as
those in the corresponding panels in Fig. 2. The yellow, black, green, blue
and red curves depict the temporal evolution of the arbitrarily chosen FHR
oscillators in the coupled systems.

Next, we investigate numerically the effects of the diffusive coupling, $D$,
on the formation of spatiotemporal patterns for the three sets of parameters.
In the numerical simulations, we consider system (3), which is a system of
coupled FHR neurons arranged in a 1-dimensional configuration, spaced at equal
distances $\Delta x=0.1$, i.e., we consider $N_{x}=500$ spatial points and
also a time-step $\Delta t=0.01$. Moreover, we consider the zero-flux boundary
condition in Eq. (7) and initial conditions in the vicinity of the fixed point
$E$ of system (3). Hence depending on the value of $D$, $E$ can be either
unstable or stable. We use the PDEPE toolbox in Matlab to solve numerically
system (3) and obtain its spatial behavior. The individual neurons for sets 1,
2 and 3 exhibit different voltage responses $u$ over time, i.e., firing
patterns, which result in various collective dynamics that we study in Subsec.
3.1 via the synchronization index $R$ of Eq. (12).

First, for parameters in set 1, system (3) for $D=0$ (which is essentially
system (2)) exhibits mixed-mode type oscillations (elliptic bursting) as can
be observed in blue in the left inset in Fig. 1. For smaller diffusion
coupling, e.g., $D=0.0001$, the system exhibits inhomogeneous irregular
firing, shown in Fig. 2(a). However, as $D$ increases, it exhibits different
firing patterns with bursting characteristics, known as elliptic bursting, a
kind of mixed-mode type oscillations. It produces from typical bursting to
rapid changes in membrane voltage oscillations, which involve simultaneous
spike changes between low- and high-amplitude oscillations. This can be
appreciated in Figs. 2(b) and 3(b). The yellow horizontal lines in Fig. 2(b)
indicate the peaks of the generated action potentials, shown in Fig. 3(b). For
higher diffusion coefficient $D=0.6$, system (3) for set 1 becomes completely
synchronized, confirmed by the synchronization index, $R=1$ that we discuss in
Subsec. 3.1. At sufficiently higher diffusion coefficient $D=8$, all neurons
show oscillation death as shown in Fig. 2(c). In panels (a)-(c) in Fig. 3, we
show the corresponding time-series of arbitrary neurons for $D=0.0001$,
particularly of neurons 1 and 8.

Next, for parameters in set 2, system (3) for $D=0$ (i.e., system (2)) shows
bursting activity as evidenced in the second inset from the left in Fig. 1.
For weak coupling $D=0.0001$, it gives rise to inhomogeneity, which results in
a train of irregular spikes, shown in Figs. 2(d) and 3(d). As the diffusion
coefficient is increased to $D=0.1$, the system exhibits mixed-mode
oscillations with different spike-numbers in single bursts, shown in Figs.
2(e) and 3(e). For even higher diffusion coefficient, e.g., for $D=1$, the
system converges to a quiescent state and the dynamics shows a fixed point
(Figs. 2(f) and 3(f)).

Last, we consider the case of tonic spiking for parameters in set 3. For
$D=0.0001$, 0.01 and 1, spatiotemporal patterns and corresponding time-series
of arbitrary neurons are shown in Figs. 2 (g)-(i) and 3 (g)-(i), respectively.
In this case, system (3) exhibits elliptic type bursting which comprises
complex firing activities, more complex than bursting and spiking activities
as it involves spikes and bunches of small-amplitude oscillations. This might
be the reason set 1 needs higher diffusive coupling for the neurons to
synchronize, compared to sets 2 and 3.

With the systematic change in $D$, the dynamics of system (3) transits from
the regime of inhomogeneous instability to a uniform steady-state through the
formation of regular structures of trains of tonic spiking at intermediate $D$
values. For sets 2 and 3 and intermediate diffusion coefficients, $D=0.1$ and
$D=0.01$, all neurons become synchronized. This synchronization property is
verified by the synchronization index $R$ that we discuss in Subsec. 3.1.

The fixed point $E$ of system (2) is unstable for the values of the diffusion
coefficient $D$ in panels (a), (b), (d), (e), (g), (h) in Figs. 2 and 3. The
stabilizing mechanism responsible for the stable dynamics in panels (c), (f),
(i) in Figs. 2 and 3 results from the stability properties of the fixed point
$E$ of systems (2) and (3). As we discussed before, the fixed point $E$ of
system (3) can be made stable by increasing appropriately the value of the
diffusion coefficient. The diffusion coefficient $D$ in panels (c), (f), (i)
in Figs. 2 and 3 is equal to 8, 1, 1, respectively, and for these values, the
fixed point $E$ of system (3) is stable. Hence starting with initial
conditions in the vicinity of the stable fixed point results in the solutions
remaining in its vicinity. This is shown in Fig. 3(c), (f), (i) for the three
sets of parameters. However, that is not the case for smaller $D$ values where
the fixed point $E$ of system (3) is unstable as can be seen in panels (a),
(b), (d), (e), (g), (h) in Figs. 2 and 3. In these cases, solutions that start
in the vicinity of the unstable fixed point $E$ drift away in time (see panels
(a), (b), (d), (e), (g), (h) in Fig. 3). Hence to stabilize system (3) in the
vicinity of the fixed point $E$, one must find appropriately big values of the
diffusion coefficient, $D$, for which solutions that start in the vicinity of
$E$, remain in its vicinity for $t>0$. Consequently, the diffusion
coefficient, $D$, plays a major role in the stabilization mechanism in system
(3) as it determines the stability properties of the fixed point $E$.

###  3.1 Synchronization index

To quantify the collective behavior in systems of coupled oscillators (in our
case of diffusively coupled FHR neurons), a synchronization measure based on
mean field theory, can be defined. The synchronization index, $R$, ranges in
$[0,1]$ [45, 48], where 0 indicates complete desynchronization (i.e., all
oscillators are completely desynchronized) and 1 complete synchronization
(i.e., all oscillators are completely synchronized). The synchronization
index, $R$, is defined by

$R=\frac{\langle F^{2}\rangle-\langle
F\rangle^{2}}{\frac{1}{N_{x}}\sum_{i=1}^{N_{x}}(\langle
u_{i}^{2}\rangle-\langle u_{i}\rangle^{2})},$ (12)

where

$F=\frac{1}{N_{x}}\sum_{i=1}^{N_{x}}u_{i}$

is the average of the membrane voltages $u_{i}$ over all spatial points and
$N_{x}$ the number of spatial points. The notation $\langle\cdot\rangle$
denotes the mean value of the argument over time.

The results for the synchronization index, $R$, as a function of the diffusion
coefficient, $D$, are shown in Fig. 4(b) for sets 1, 2 and 3, in blue, brown
and green, respectively. For sets 1 and 2, the coupled neurons are
desynchronized for weak diffusive coupling, $D$, as $R$ is close to 0. This
indicates that all neurons exhibit almost desynchronized oscillations as shown
in the spatial plots in panels (a) and (d) in Fig. 3. However, as $D$
increases, all neurons of the three sets exhibit complete synchronization,
evidenced in panels (b), (c), (e) and (f) in Fig. 3. Interestingly, $R$ is
almost equal to 1 for weak coupling $D$ for set 3, which can be appreciated in
Fig. 4(b).

Figure 4: Boundaries of stable and unstable regions and synchronization index,
$R$, for the sets of parameters 1, 2 and 3. (a) The dashed blue, brown and
green curves (see Eqs. (9), (10), (11)) indicate the boundaries of the stable
and unstable regions of model (3) for the parameter sets 1, 2 and 3,
respectively. The solid blue, brown and green lines are located at $D=6$ (set
1), $D=0.98$ (set 2) and $D=0.9$ (set 3), are derived numerically with the
help of the bisection method and, separate the stable from the unstable
regions. Here, “us” stands for unstable region. (b) The synchronization index,
$R$, as a function of the diffusion coefficient, $D$, for the sets of
parameters 1, 2 and 3, shown in blue, brown and green, respectively. Note that
in panel (b), the line segments connect the numerically computed values of $R$
plotted by small filled circles.

##  4 Diffusively coupled FHR neurons arranged in a 2-dimensional
configuration

In this section, we investigate different types of pattern formations for
diffusively coupled FHR neurons in a square domain, $\Psi$. In this case, the
mathematical model is given by the system of coupled PDEs

$\displaystyle\frac{{\partial u}}{{\partial t}}$
$\displaystyle=f(u)-v+w+I+D\Bigg{(}\frac{{\partial^{2}u}}{{\partial
x^{2}}}+\frac{{\partial^{2}u}}{{\partial y^{2}}}\Bigg{)},$
$\displaystyle\frac{{\partial v}}{{\partial t}}$
$\displaystyle=\delta\left({a+u-bv}\right),$ (13)
$\displaystyle\frac{{\partial w}}{{\partial t}}$
$\displaystyle=\mu\left({c-u-w}\right),$

where we use similar initial and boundary conditions as in Sec. 3, adapted to
the case. Thus, here $u\equiv u(x,y,t)$, $v\equiv v(x,y,t)$ and $w\equiv
w(x,y,t)$, where $x$ and $y$ are the spatial coordinates and $t$ is the time.
System (13) can be written in the compact form

$\frac{{\partial X}}{{\partial t}}=H(X,I)+D_{1}\nabla^{2}X,$ (14)

where

$X=\begin{pmatrix}X_{1}\\\ X_{2}\\\ X_{3}\end{pmatrix}=\begin{pmatrix}u\\\
v\\\ w\end{pmatrix},\;H=\begin{pmatrix}u-\frac{{u^{3}}}{3}-v+w+I\\\
\delta({a+u-bv})\\\ \mu({c-u-w})\end{pmatrix}$

and

$D_{1}=\begin{pmatrix}D&0&0\\\ 0&0&0\\\ 0&0&0\end{pmatrix}.$

Equation (14) will prove useful in the next, where we discuss the amplitude
equations that we use to study analytically the formation of spatiotemporal
patterns in the vicinity of the two Hopf bifurcation points, HB1 and HB2.

###  4.1 Amplitude equations

The dynamics of system (13) in the vicinity of a Hopf bifurcation can be
described by the amplitude equations, also known as CGLE. It’s general form is
given by [35, 36, 38]

$\frac{\partial W}{\partial
t}=W-(1+i\alpha){\big{|}W\big{|}}^{2}W+(1+i\beta)\nabla^{2}W,$ (15)

where $i$ denotes the imaginary unit of the complex numbers, i.e., $i^{2}=-1$
and $\nabla^{2}$ is the 2-dimensional Laplacian operator. As the real
parameters $\alpha$ and $\beta$ vary, the complex amplitude, $W$, exhibits
rich dynamics. In the following, we seek to find the values of $\alpha$ and
$\beta$ in the vicinity of the two Hopf bifurcations, HB1 and HB2 (as
discussed in Subsec. 2.1), i.e., around the bifurcation points $I=0.137$ and
$I=3.16298$.

The Jacobian matrix $J$ in this case has two complex conjugate, imaginary
eigenvalues $\lambda_{1,2}=\pm\omega i$. The right eigenvectors that
correspond to the eigenvalues $\pm i\omega$ are denoted by $U_{1}$,
$U_{2}=\bar{U}_{1}$ and the left by $U_{1}^{+}$, $U_{2}^{+}=\bar{U}_{1}^{+}$.
For simplicity, we consider $U=U_{1}=\bar{U}_{2}$ and
$U^{+}=U_{1}^{+}=\bar{U}_{2}^{+}$. Furthermore, the right and left
eigenvectors can be normalized according to $U_{i}^{+}U_{j}=\delta_{ij}$,
where $\delta_{ij}$ is the Kronecker delta, i.e.,
$U^{+}U=\bar{U}^{+}\bar{U}=1$ and $U^{+}\bar{U}=\bar{U}^{+}U=0$. In Subsec.
2.1, we have shown that the two Hopf bifurcations, HB1 and HB2, occur at
$I=0.137$ and $I=3.16298$, respectively. Interestingly, the eigenvalues of $J$
for both are the same and are given by $\lambda_{1,2}=\pm 0.279302i$ and
$\lambda_{3}=-0.0036$ as $J$ contains the term $u_{0}^{2}$ and its value is
the same for both bifurcation points (i.e., $u_{0}=\pm 0.96829$). The right
and left eigenvectors that correspond to the eigenvalue
$\lambda_{1}=0.279302i$ are
$U=(0.963143,0.0600545-0.262111i,-0.00005+0.0069i)^{T}$ and
$U^{+}=(0.5198-0.1159i,-0.0108+1.8586i,-0.4018-1.8639i)$, respectively. These
again refer to both bifurcation points as the Jacobian matrix $J$ and its
eigenvalues are the same for both.

To compute $\alpha$ and $\beta$ in Eq. (15), we start with the unscaled
Ginzburg-Landau equation [35, 36]

$\dot{W}=\sigma_{1}IW-g\big{|}W\big{|}^{2}W+d\nabla^{2}W$ (16)

and seek to find the values of $\sigma_{1}$, $g$ and $d$ using the values in
Table 1 in [35, 36]. Particularly, parameters $\sigma_{1}$, $g$ and $d$ are
given by

$\displaystyle\sigma_{1}$ $\displaystyle=U^{+}H_{XI}U+U^{+}H_{XX}(U,l_{001}),$
$\displaystyle g$
$\displaystyle=-\Big{(}U^{+}H_{XX}(U,l_{110})+U^{+}H_{XX}(\bar{U},l_{200})+\frac{1}{2}U^{+}H_{XXX}(U,U,\bar{U})\Big{)},$
$\displaystyle d$ $\displaystyle=U^{+}D_{1}U,$

respectively, where

$\displaystyle H_{{XX}}\left({\xi,\,\eta}\right)$
$\displaystyle=\sum\limits_{i,j=1}^{3}{\left.{\frac{{\partial^{2}H}}{{\partial
X_{i}\partial X_{j}}}}\right|_{E=(u_{0},v_{0},w_{0})}\xi_{i}\eta_{j}},$
$\displaystyle H_{{XXX}}\left({\xi,\,\eta,\,\zeta}\right)$
$\displaystyle=\sum\limits_{i,j,k=1}^{3}{\left.{\frac{{\partial^{3}H}}{{\partial{X}_{i}\partial{X}_{j}\partial{X}_{k}}}}\right|}_{E=(u_{0},v_{0},w_{0})}\xi_{i}\eta_{j}\zeta_{k},$
$\displaystyle l_{001}$ $\displaystyle=(0.4571,0.5714,-0.4571)^{T},$
$\displaystyle l_{110}$ $\displaystyle=(0.8212,1.0265,-0.8212)^{T},$
$\displaystyle l_{200}$
$\displaystyle=(-0.2395-2.1447i,-0.3070-0.0009i,0.0077-0.0008i)^{T}.$

In our case,

$H=\begin{pmatrix}u-\frac{u^{3}}{3}-v+w+I\\\ 0.08(0.7+u-0.8v)\\\
0.002(-0.775-u-w)\end{pmatrix},$

thus

$\frac{\partial^{2}H}{\partial
X_{1}^{2}}\bigg{|}_{E}=\begin{pmatrix}-2u_{0}\\\ 0\\\ 0\end{pmatrix}\mbox{ and
}\;\frac{\partial^{3}H}{\partial X_{1}^{3}}\bigg{|}_{E}=\begin{pmatrix}-2\\\
0\\\ 0\end{pmatrix},$

with all other terms in $H_{XX}$ and $H_{XXX}$ being equal to zero. With this
in mind, we obtain $\sigma_{1}=0.4432-0.0988i$, $g=0.3642+2.1015i$ and
$d=D(0.5006-0.1117i)$. Next, using the transformation

$\displaystyle W$ $\displaystyle=\sqrt{I}W^{\prime},$ $\displaystyle t$
$\displaystyle=t^{\prime}/I,$ $\displaystyle x$
$\displaystyle=x^{\prime}/\sqrt{I},$ (17) $\displaystyle y$
$\displaystyle=y^{\prime}/\sqrt{I},$

we obtain

$\frac{\partial W}{\partial t}=\sigma_{1}W-g|W|^{2}W+d\nabla^{2}W,$ (18)

which is known as the scaled complex Ginzburg-Landau equation. Transformation
(17) renders CGLE independent of the distance from the bifurcation point. The
dynamics of system (13) can be described by the solution of Eq. (18), except
for a short initial decay of transient eigenmodes. Finally, introducing the
new transformation

$\displaystyle W$
$\displaystyle=\sqrt{\frac{\Re(\sigma_{1})}{\Re(g)}}W^{\prime}e^{i\frac{\Im(\sigma_{1})}{\Re(\sigma_{1})}t^{\prime}},$
$\displaystyle t$ $\displaystyle=\frac{t^{\prime}}{\Re(\sigma_{1})},$
$\displaystyle x$
$\displaystyle=\sqrt{\frac{\Re(d)}{\Re(\sigma_{1})}}x^{\prime},$
$\displaystyle y$
$\displaystyle=\sqrt{\frac{\Re(d)}{\Re(\sigma_{1})}}y^{\prime},$

we obtain the dimensionless CGLE (15), that is the equation

$\frac{\partial W}{\partial
t}=W-(1+i\alpha){\left|{W}\right|}^{2}W+(1+i\beta)\nabla^{2}W,$ (19)

where $\alpha=\frac{\Im(g)}{\Re(g)}=5.7706$ and
$\beta=\frac{\Im(d)}{\Re(d)}=-0.2230$. In this framework, $\Re(*)$ and
$\Im(*)$ denote the real and imaginary parts of the argument, respectively.
Clearly, $\alpha+\beta=5.5476>0$, which confirms the existence of antispiral
patterns in the vicinity of the two Hopf bifurcations [42], HB1 and HB2, at
the critical points $I=0.137$ and $I=3.16298$, respectively.

###  4.2 Results

Here, we focus on the diffusively coupled FHR model (13) and study numerically
the emergence of complex structures in a square spatial domain $\Psi$ for
$D>0$, considering $I$ as the varying parameter. The method of amplitude
equations was employed in Subsec. 4.1 to show the existence of antispiral
patterns in the vicinity of the two Hopf bifurcations, HB1 and HB2. In the
following, we discuss specific patterns that emerge around HB1.

To do so, we construct a spatially extended medium with 2-dimensional spatial
diffusion using an $N_{x}\times N_{y}=100\times 100$ mesh grid. We solve
numerically the resulting system using the forward Euler method with central
differences, a finite-difference scheme, with spatial step size $\Delta=\Delta
x=\Delta y=1.25$ and integration time-step $\Delta t=0.1$. The initial
conditions are considered with appropriate periodic perturbations from the
initial conditions of the system in Sec. 2, using similar boundary conditions
as in the 1-dimensional case, adapted to the case. We denote by $u_{i,j}(t)$
the membrane voltage of a neuron at the node $(i,j)$ on the grid at time $t$.
Consequently, the sum of the spatial, second-order partial derivatives
$\frac{{{\partial^{2}}u}}{{\partial{x^{2}}}}+\frac{{{\partial^{2}}u}}{{\partial{y^{2}}}}$
can be approximated by

$\frac{{{\partial^{2}}u}}{{\partial{x^{2}}}}+\frac{{{\partial^{2}}u}}{{\partial{y^{2}}}}=\frac{1}{{{\Delta^{2}}}}\left({{u_{i-1,\,j}}+{u_{i+1,\,j}}+{u_{i,\,j-1}}+{u_{i,\,j+1}}-4{u_{i,\,j}}}\right).$

We note that in the following, we will study the spatiotemporal
characteristics in the context of the nonlinear, coupled, diffusive system
(13), where the coupling indicates the synaptic, diffusive coupling among
neurons on the nodes of the $N_{x}\times N_{y}$ grid. Furthermore, we
construct a square $N_{x}\times N_{y}$ domain $\Psi$ (grid) using nearest
neighbor connectivities and explore the evolution of target waves and spirals
for fixed values of the diffusion coefficient, $D$, using long-time numerical
integration, i.e., up to $t=20000$.

Figure 5: Formation of target waves and spiral patterns resulting from the
system of diffusively coupled FHR neurons (13) arranged in a 2-dimensional
space for set 1 with $D=0.25$ at $t=3450$, 3500, 20000, 40000 in panels (a) to
(d), set 2 with $D=0.3$ at $t=2000$, 3000, 5000, 40000 in panels (e) to (h)
and set 3 with $D=0.009$ at $t=3500$, 5000, 20000, 40000 in panels (i) to (l).
The color bars encode the values of the membrane voltages $u(x,y,t)$ of the
diffusively coupled FHR neurons of system (13).

We start with set 1, for which the single-neuron model (2) with no diffusion,
exhibits mixed-mode oscillations, often known as elliptic type bursting [28].
For $D=0.25$, system (13) gives rise to the emergence of target waves at
$t=3450$ shown in Fig. 5(a). We also find that these target waves are sparse
as shown in the same panel. The reason might be that most of the neurons
cannot activate neighboring neurons in short time intervals. However, as time
increases to $20000$, small spiral-like activity can be seen, bounded by
target waves [41]. The developed target waves become denser, occupying the
whole extend of the mesh grid, as shown in panels (b) and (c) in Fig. 5. The
coexistence of spirals and target waves can be observed at $t=40000$ in Fig.
5(d), which suggests the patterns in the spatial domain are stable.

Next, we consider set 2 with $D=0.3$, where individual neurons exhibit mixed-
mode oscillations. At the initial stage, multi-arm spiral waves emerge from
broken waves, which can be observed in Fig. 5(e). As the time increases to
$t=3000$, target waves and various wavefronts emerge as can be observed in
Fig. 5(f). The spiral-like patterns coexist with target waves until $t=5000$,
shown in Fig. 5(g), and the spirals are not stable. The reason may be the
excitability of neurons in the arm segments of small spirals being high,
attracting other small segments of spiral arms and forming more stable
wavefronts. Interestingly, at even higher times, for example at $t=40000$,
only spirals can be observed, shown in Fig. 5(h), with the centers of the
wavefronts frequently changing over time. A group of small-size spiral
profiles emerge that change into multi-arm spiral waves. Our numerical results
show the emergence of different wavefronts that develop and break over time,
with spirals emerging thereafter.

Multi-arm spiral waves emerge for set 3 for small diffusion coefficient
values. For example, at $D=0.009$, target waves appear to cover partially the
grid. Initially, we observe the formation of distinct clusters consisting of
small curls on the blue backdrop of the mesh grid, that can be seen in panels
(i) and (j) in Fig. 5. However, as time increases, for example at $t=20000$
and $t=40000$, multi-arm spirals emerge, shown in panels (k) and (l) in Fig.
5. This shows that the small curls collide with neighboring curls that extend
to the whole mesh. With further increase in time, multi-arm spirals emerge.

Next, we discuss the effect of diffusive coupling on the formation of target
waves and spirals for long integration times. At coupling strength $D=0.05$,
the formation of two-arm spiral waves [41] can be observed at $t=20000$, shown
in Fig. 6(a) for set 1. Only a few spirals are visible. However, as $D$
increases slightly, for $D=0.09$ and $D=0.1$, antispiral waves emerge, shown
in panels (b) and (c) in Fig. 6. The results are verified using Eq. (19). It
follows from the study of the amplitude equations that the antispirals exist
in the vicinity of the two supercritical Hopf points and that they depend on
$\alpha$ and $\beta$. Solving the system numerically, we can see the
antispirals for suitable values of the diffusive coupling, $D$, depending on
$\alpha$ and $\beta$. These multi-arm antispiral waves with an increasing
number of arms extend to the whole spatial grid. For even higher diffusion
values, for example for $D=0.245$, double spirals can be observed, shown in
Fig. 6(d). Periodic forcing and coupling strength play a major role in
promoting target waves by suppressing spiral waves [40, 41]. Increasing
slightly the diffusive coupling to $D=0.25$, we observe the same route of
transition, i.e., two spirals that collide and generate target wave patterns,
shown in the first panel in Fig. 5.

When considering set 2 and for smaller couplings, for example for $D=0.1$,
multi-arm spirals emerge. For even higher $D$, for example for $D=0.25$ and
$D=0.3$, the multi-arm spiral waves dissolve and, target waves and spirals
emerge, shown in panels (e) to (g) in Fig. 6. With the further increase of $D$
to 0.9, only spiral patterns can be observed, shown in Fig. 6(h).

For set 3, with increasing $D$, multi-arm spirals emerge, occupying a greater
extend on the spatial grid, observed in panels (i) to (l) in Fig. 6. We cross-
checked the appearance of similar spatiotemporal patterns on a bigger,
$N_{x}\times N_{y}=500\times 500$ mesh grid and we obtained similar results,
with no significant changes. Finally, the same wavy patterns emerged when
running simulations for even longer times, i.e., up to $t=50000$.

Figure 6: Snapshots of pattern formation of the system of diffusively coupled
FHR neurons (13) arranged in a 2-dimensional configuration. Panel (a)
$D=0.05$, (b) $D=0.09$, (c) $D=0.1$, (d) $D=0.245$ for set 1, panel (e)
$D=0.1$, (f) $D=0.25$, (g) $D=0.3$, (h) $D=0.9$ for set 2, panel (i)
$D=0.005$, (j) $D=0.007$, (k) $D=0.009$ and panel (l) $D=0.01$ for set 3. The
color bars encode the values of the membrane voltages $u(x,y,t)$ of the
diffusively coupled FHR neurons of system (13).

##  5 Conclusions and Discussion

In this paper, we considered the biophysically motivated, slow-fast, excitable
FitzHugh-Rinzel neuron model, which provides a variety of neuronal responses.
It exhibits a diverse repertoire of firing patterns for certain fixed sets of
parameters and different external current stimuli. We studied the FitzHugh-
Rinzel model theoretically and numerically in different dynamical regimes.
Further, we discussed its local and global stability. We performed a
bifurcation analysis and demonstrated the qualitative changes between stable
steady and oscillatory states. We incorporated 1-dimensional diffusion to
construct a 1-dimensional chain of neurons and investigated several complex
behaviors that pertain to synchronization. The stability analysis was
performed for the diffusively coupled system and we showed how the stability
changes from an unstable regime to a stable with the increase of the diffusion
coefficient. The diffusively coupled system changes its dynamical
characteristics as the coupling strength changes. Next, we studied the
diffusive system on a square spatial domain. The amplitude equation,
describing the onset of spirals close to the Hopf bifurcations, was derived
and numerical simulations were provided. As the parameters change, we showed
the appearance of target and spiral waves for three sets of parameters that
correspond to three distinct oscillatory regimes.

We explored emerging target waves and spirals as dynamical features and
verified them analytically using amplitude equations. Our investigation shows
that the spatial dynamics of the slow-fast model exhibits two- and multi-arm
spiral waves for low diffusion couplings. We found that multi-arm spiral waves
with increasing number of arms occupy the 2-dimensional mesh grid and that
they are unstable. Bursting in the activity of single neurons is a bio-
physiological phenomenon, however, bursting in the activity of neural
populations might be pathological [53, 8]. Interestingly, these emergent
patterns may be relevant to the synchronized activities of neural populations,
particularly related to neurological diseases [8]. The propagation of neuronal
impulses can be relevant to brain functioning [15, 41, 4, 40]. In a weakly
coupled system of pancreatic $\beta$-cells, bursting varies, where pancreatic
$\beta$-cells secrete insulin [47]. The analysis of the mechanisms underlying
spatial profiles of activity in the neural tissue is important in
understanding a wide range of biophysical and pathological phenomena [49, 3,
53]. Moreover, the form of the nerve impulse propagation is relevant to
certain brain pathologies [49, 50, 18, 6].

Finally, we discussed in detail how nonlinearities in the biophysical,
excitable model change the distribution of single-cell characteristics into
different types of patterns. Spiral waves emerge frequently in cortical areas
with limited lifespan. They can modify cortical activities affecting the
oscillation frequency and spatial coherence-like activity. The emergence of
spiral waves during sleep-like states varies greatly and, can also organize
and modulate the pathological patterns during epilepsy. Spiral waves have been
observed to play a major role in organizing irregular dynamics in cortical
neurons and rhythmic behavior [49]. Our results reveal a multitude of neural
excitabilities and possible conditions for the emergence of spiral-wave
formation in diffusively coupled FitzHugh-Rinzel systems with different firing
characteristics.

## Data availability

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

## References

  * [1] B. Ermentrout.  Neural networks as spatio-temporal pattern-forming systems.  Rep. Prog. Phys., 61(4):353–430, 1998.
  * [2] E. M. Izhikevich.  Neural excitability, spiking and bursting.  Int. J. Bifurc. Chaos, 10(06):1171–1266, 2000.
  * [3] E. M. Izhikevich.  Dynamical Systems in Neuroscience.  MIT press, 2007.
  * [4] A. Keane and P. Gong.  Propagating waves can explain irregular neural dynamics.  J. Neurosci., 35(4):1591–1605, 2015.
  * [5] J. Ma and J. Tang.  A review for dynamics in neuron and neuronal network.  Nonlin. Dyn., 89(3):1569–1578, 2017.
  * [6] R. G. Townsend and P. Gong.  Detection and analysis of spatiotemporal patterns in brain activity.  PLOS Comput. Biol., 14(12):e1006643, 2018.
  * [7] S. Kondo and T. Miura.  Reaction-diffusion model as a framework for understanding biological pattern formation.  Science, 329(5999):1616–1620, 2010.
  * [8] S. R. Meier, J. L. Lancaster, and J. M. Starobin.  Bursting regimes in a reaction-diffusion system with action potential-dependent equilibrium.  PLOS One, 10(3), 2015.
  * [9] J. G Milton, P. H. Chu, and J. D. Cowan.  Spiral waves in integrate-and-fire neural networks.  In Advances in neural information processing systems, pages 1001–1006, 1993.
  * [10] I. Belykh, E. De Lange, and M. Hasler.  Synchronization of bursting neurons: what matters in the network topology.  Phys. Rev. Lett., 94(18):188101, 2005.
  * [11] A. M. Turing.  The chemical basis of morphogenesis.  Phil. Trans. R. Soc. Lond. B, 237:37–72, 1952.
  * [12] A. Mondal, S. K. Sharma, R. K. Upadhyay, M. A. Aziz-Alaoui, P. Kundu, and C. Hens.  Diffusion dynamics of a conductance-based neuronal population.  Phys. Rev. E, 99(4):042307, 2019.
  * [13] S. Kondo.  The reaction-diffusion system: a mechanism for autonomous pattern formation in the animal skin.  Genes Cells, 7(6):535–541, 2002.
  * [14] Jr. R. Erichsen and L. G. Brunnet.  Multistability in networks of Hindmarsh-Rose neurons.  Phys. Rev. E, 78(6):061917, 2008.
  * [15] M. Grace and M. T. Hütt.  Regulation of spatiotemporal patterns by biological variability: general principles and applications to dictyostelium discoideum.  PLOS Comput. Biol., 11(11):e1004367, 2015.
  * [16] Z. P. Kilpatrick and P. C. Bressloff.  Effects of synaptic depression and adaptation on spatiotemporal dynamics of an excitatory neuronal network.  Physica D, 239(9):547–560, 2010.
  * [17] A. Roxin, H. Riecke, and S. A. Solla.  Self-sustained activity in a small-world network of excitable neurons.  Phys. Rev. Lett., 92(19):198101, 2004.
  * [18] M. Song, M. Kang, H. Lee, Y. Jeong, and S. Paik.  Classification of spatiotemporal neural activity patterns in brain imaging data.  Sci. Rep., 8(1):1–12, 2018.
  * [19] B. Ambrosio, M. A. Aziz-Alaoui, and V. L. E. Phan.  Global attractor of complex networks of reaction-diffusion systems of Fitzhugh-Nagumo type.  Discrete Contin. Dyn. Syst. B, 23(9):3787, 2018.
  * [20] B. Ambrosio, M. A. Aziz-Alaoui, and V. L. E. Phan.  Large time behaviour and synchronization of complex networks of reaction-diffusion systems of FitzHugh-Nagumo type.  IMA J Appl Math, 84(2):416–443, 2019.
  * [21] G. Gambino, M. C. Lombardo, G. Rubino, and M. Sammartino.  Pattern selection in the 2d FitzHugh-Nagumo model.  Ric. Mat., 68(2):535–549, 2019.
  * [22] M. Jun, Y. Li-Jian, W. Ying, and Z. Cai-Rong.  Spiral wave in small-world networks of Hodgkin-Huxley neurons.  Commun. Theor. Phys., 54(3):583, 2010.
  * [23] M. Kuznetsov, A. Kolobov, and A. Polezhaev.  Pattern formation in a reaction-diffusion system of Fitzhugh-Nagumo type before the onset of subcritical turing bifurcation.  Phys. Rev. E, 95(5):052208, 2017.
  * [24] X. Liao, Q. Xia, Y. Qian, L. Zhang, G. Hu, and Y. Mi.  Pattern formation in oscillatory complex networks consisting of excitable nodes.  Phys. Rev. E, 83(5):056204, 2011.
  * [25] Anotida Madzvamuse, Andy HW Chung, and Chandrasekhar Venkataraman.  Stability analysis and simulations of coupled bulk-surface reaction–diffusion systems.  Proc. R. Soc. A, 471(2175):20140546, 2015.
  * [26] Q. Wang, M. Perc, Z. Duan, and G. Chen.  Delay-enhanced coherence of spiral waves in noisy Hodgkin-Huxley neuronal networks.  Phys. Lett. A, 372(35):5681–5687, 2008.
  * [27] R. FitzHugh.  Impulses and physiological states in models of nerve membrane.  Biophys. J., pages 445–466, 1961.
  * [28] E. M Izhikevich.  Synchronization of elliptic bursters.  SIAM Rev., 43(2):315–344, 2001.
  * [29] J. Rinzel.  A formal classification of bursting mechanisms in excitable systems.  In Mathematical Topics in Population Biology, Morphogenesis and Neurosciences, pages 267–281. Springer, 1987.
  * [30] J. Rinzel and W. C. Troy.  Bursting phenomena in a simplified oregonator flow system model.  J. Chem. Phys, 76(4):1775–1789, 1982.
  * [31] J. Wojcik and A. Shilnikov.  Voltage interval mappings for an elliptic bursting model.  In Nonlin. Dyn., pages 195–213. Springer, 2015.
  * [32] M. Bär and M. Eiswirth.  Turbulence due to spiral breakup in a continuous excitable medium.  Phys. Rev. E, 48(3):R1635, 1993.
  * [33] M. Perc.  Fluctuating excitability: a mechanism for self-sustained information flow in excitable arrays.  Chaos Solitons Fractals, 32(3):1118–1124, 2007.
  * [34] I. S. Aranson and L. Kramer.  The world of the complex Ginzburg-Landau equation.  Rev. Mod. Phys., 74(1):99, 2002.
  * [35] M. Ipsen, F. Hynne, and P. G. Sørensen.  Amplitude equations for reaction-diffusion systems with a Hopf bifurcation and slow real modes.  Physica D, 136(1-2):66–92, 2000.
  * [36] M. Ipsen, L. Kramer, and P. G. Sørensen.  Amplitude equations for description of chemical reaction-diffusion systems.  Phys. Rep., 337(1-2):193–235, 2000.
  * [37] Y. Kuramoto.  Springer series in synergetics.  In Chemical Oscillations, Waves and Turbulence, volume 19, page 130\. Springer Berlin, 1984.
  * [38] E. P. Zemskov, V. K. Vanag, and I. R. Epstein.  Amplitude equations for reaction-diffusion systems with cross diffusion.  Phys. Rev. E, 84(3):036216, 2011.
  * [39] J. Ma, C. Wang, W. Jin, and Y. Wu.  Transition from spiral wave to target wave and other coherent structures in the networks of Hodgkin-Huxley neurons.  Appl. Math. Comput., 217(8):3844–3852, 2010.
  * [40] X. Wu and J. Ma.  The formation mechanism of defects, spiral wave in the network of neurons.  PLOS one, 8(1), 2013.
  * [41] B. Hu, J. Ma, and J. Tang.  Selection of spiral waves in a regular network of neurons.  PLOS One, 8(7), 2013.
  * [42] E. M. Nicola, L. Brusch, and M. Bär.  Antispiral waves as sources in oscillatory reaction-diffusion media.  J. Phys. Chem. B, 108(38):14733–14740, 2004.
  * [43] Yunfan Gong and David J Christini.  Antispiral waves in reaction-diffusion systems.  Phys. Rev. Lett., 90(8):088302, 2003.
  * [44] Yanmei Kang, Yaqian Chen, Yuxuan Fu, Zuolei Wang, and Guanrong Chen.  Formation of spiral wave in Hodgkin-Huxley neuron networks with gamma-distributed synaptic input.  Commun Nonlinear Sci Numer Simul, 83:105112, 2020.
  * [45] Yong Liu, Ying Xu, and Jun Ma.  Synchronization and spatial patterns in a light-dependent neural network.  Commun Nonlinear Sci Numer Simul, 89:105297, 2020.
  * [46] J. Ma, B. Hu, C. Wang, and W. Jin.  Simulating the formation of spiral wave in the neuronal system.  Nonlin. Dyn., 73(1-2):73–83, 2013.
  * [47] S. Raghavachari and J. A. Glazier.  Waves in diffusively coupled bursting cells.  Phys. Rev. Lett., 82(14):2991, 1999.
  * [48] C. Wang, M. Lv, A. Alsaedi, and J. Ma.  Synchronization stability and pattern selection in a memristive neuronal network.  Chaos, 27(11):113108, 2017.
  * [49] X. Huang, W. Xu, J. Liang, K. Takagaki, X. Gao, and J. Wu.  Spiral wave dynamics in neocortex.  Neuron, 68(5):978–990, 2010.
  * [50] S. J. Schiff, X. Huang, and J. Wu.  Dynamical evolution of spatiotemporal patterns in mammalian middle cortex.  Phys. Rev. Lett., 98(17):178102, 2007.
  * [51] G. P. Kremmydas and A. V. Holden.  Spiral-wave meandering in reaction-diffusion models of ventricular muscle.  Chaos Solitons Fractals, 13(8):1659–1669, 2002.
  * [52] F. X. Witkowski, L. J. Leon, P. A. Penkoske, W. R. Giles, M. L. Spano, W. L. Ditto, and A. T. Winfree.  Spatiotemporal evolution of ventricular fibrillation.  Nature, 392(6671):78–82, 1998.
  * [53] H. Jiang and Z. Hou.  Formation of spiral waves with substructure in a bursting media.  Chaos, 25(12):123105, 2015.
  * [54] Suyu Liu and Qingyun Wang.  Transition dynamics of generalized multiple epileptic seizures associated with thalamic reticular nucleus excitability: A computational study.  Commun Nonlinear Sci Numer Simul, 52:203–213, 2017.
  * [55] Edmund X DeJesus and Charles Kaufman.  Routh-Hurwitz criterion in the examination of eigenvalues of a system of nonlinear ordinary differential equations.  Phys. Rev. A, 35(12):5288, 1987.
  * [56] A. Dhooge, W. Govaerts, Yu. A. Kuznetsov, H. G. E. Meijer, and B. Sautois.  New features of the software MatCont for bifurcation analysis of dynamical systems.  MCMDS, 14(2):147–175, 2008.
  * [57] W. M. Liu.  Criterion of Hopf-bifurcations without using eigenvalues.  J. Math. Anal. Appl., 182:250–256, 1994.
  * [58] B. Dubey and J. Hussain.  Modelling the interaction of two biological species in a polluted environment.  J. Math. Anal. Appl., 246:58–79, 2000.