# Estimation of Optical Aberrations in 3D Microscopic Bioimages
††thanks:  1https://github.com/kiraving/aberration

Kira Vinogradova  Center for Systems Biology Dresden
Max Planck Institute of Molecular Cell Biology and Genetics
Dresden, Germany
<EMAIL_ADDRESS>Eugene W. Myers  Center for Systems Biology Dresden
Max Planck Institute of Molecular Cell Biology and Genetics
Dresden, Germany
<EMAIL_ADDRESS>
###### Abstract

The quality of microscopy images often suffer from optical aberrations. These
aberrations and their associated point spread functions have to be
quantitatively estimated to restore aberrated images. The recent state-of-the-
art method PhaseNet, based on a convolutional neural network, can quantify
aberrations accurately but is limited to images of point light sources, e.g.
fluorescent beads. In this research, we describe an extension of PhaseNet
enabling its use on 3D images of biological samples. To this end, our method
incorporates object-specific information into the simulated images used for
training the network. Further, we add a Python-based restoration of images via
Richardson–Lucy deconvolution. We demonstrate that the deconvolution with the
predicted PSF can not only remove the simulated aberrations but also improve
the quality of the real raw microscopic images with unknown residual PSF. We
provide code1 for fast and convenient prediction and correction of
aberrations.

###### Keywords:

optics, bioimage analysis, 3D image regression, deep learning, image
restoration

##  I Introduction

An aberration is a deviation of a light beam from the trajectory described by
geometrical optics. Optical aberrations can be chromatic (color aberrations)
or monochromatic. Monochromatic aberrations make images look blurred or
distorted. In adaptive optics (AO), the branch of microscopy that studies
estimation and correction of aberrations, an estimation of monochromatic
aberrations is possible via finding the presence and amplitudes of the Zernike
polynomials [1] describing the PSF at a point in space. In terms of Deep
Learning, estimation of aberrations means predicting a 1D vector of amplitudes
for each Zernike mode.

AO methods aim at enhancing the image quality by correcting optical
aberrations, which need to be estimated first. There are two ways to estimate
aberrations: using a wavefront sensor [2] or computationally, where the latter
mostly relies on optimization algorithms [3, 4] or deep learning [5, 6, 7, 8,
9]. The advantage of sensorless approaches based on deep learning is their
simplified design of the optical system; however, they require a large number
of images for training.

State-of-the-art computational methods are limited to either trivial point
light sources [5, 3, 7, 8, 9], which reflect the point spread function (PSF)
by definition, or to objects with a simple well-defined shape like
microtubules [6]. The work [9] provides a comparison of six models, with five
of them being Convolutional Neural Networks (CNN) of different depths.
Remarkably, the authors demonstrated that the best performance on simulated
images of fluorescent beads was achieved by a rather small network with 7
convolutional layers combined with 5 max pooling and 3 dense layers. Some of
the state-of-the-art methods are designed only for two-dimentional images [8,
9], although the PSF affecting volumetric objects has three dimentions.

To the best of our knowledge, the PhaseNet publication [5] provides the only
publicly available datasets with experimentally introduced aberrations via a
point-scanning AO microscope with a spatial light modulator and a custom-built
widefield AO microscope with a deformable mirror. The PhaseNet datasets are
aberrated images of fluorescent beads that have a trivial spherical shape.

In this research, we explore the applicability of PhaseNet to objects of a
complex shape — the protoplasmic astrocytes in the hippocampus of Rattus
norvegicus [10] — for which we synthetically generate aberrated images.

##  II Methods

###  II-A PhaseNet

The core method used in this study, PhaseNet [5], is a CNN that takes a 3D
image as an input and outputs the amplitudes for each Zernike mode that it was
trained on. The CNN consists of five convolutional blocks (two 3D
convolutional layers and one 2D max-pooling layer) and three dense layers.

PhaseNet has been trained on synthetic images of 3D aberrated points. The
synthetic generator produced random aberrations with the amplitudes in a
predefined range. The training and validation datasets were formed using one
generator. PhaseNet has been tested on randomized synthetic data from the same
generator and on experimentally acquired 3D images of fluorescent beads on two
AO microscopes.

The model predicted the amplitudes for 11 nontrivial aberration modes, and the
visualization of the PSF and the wavefront reconstruction were shown. The
Python code for the model and data generation along with the data have been
made publicly available on Github by the authors. For details regarding the
data acquisition and CNN training, please refer to the original publication
[5].

Fig. 1 summarizes the training procedure and usage of PhaseNet.

Figure 1: Schematic of PhaseNet [5] algorithm. Training: the generator takes
the input parameters of the microscopy system (refractive index, numerical
aperture, etc.) and parameters for data generation (size of the point source,
noise to add, etc.) and produces the 3D PSF together with the ground truth
(aberration modes and their amplitudes). The CNN is trained in a supervised
manner on the generator’s output. Prediction: the trained CNN predicts the
vector of the aberration modes with their amplitudes (one float number per
mode) for an input 3D image of a bead. The predictions are used to visualize
the PSF and the wavefront.

###  II-B Synthetic Data Generator

The synthetic generator in PhaseNet produces random PSFs. The user defines the
nomenclature (ANSI and Noll are available), the list of aberration modes, the
range of amplitudes, the size of the generated PSF, and the output size (the
size of a crop of the PSF). The PSF can be cropped, and the cropping position
can be configured with the parameters “jitter” and “max jitter”. If “jitter”
is “false” (default), the crop is made in the center of the image; if it is
“true”, the crop is taken from a random part of the image. “Max jitter”
defines how far from the center a crop can be taken (by default, anywhere in
the image). The user can select the $z$ planes of the 3D image to train the
network on (for example, only a couple of central planes). It has been
demonstrated [5] that the best results are achieved when all $z$ planes are
used, which means setting “max jitter” to “None” (default) and training on
full-sized 3D images.

The parameters of the optical system (refractive index, light wavelength
$\lambda$, numerical aperture of the detection objective (NA), voxel size)
must be provided for correct calculation of the synthetic PSFs. The generator
convolves a point (sphere) with the synthetic PSF. Optionally, Gaussian noise
can be added after the convolution if its parameters (standard deviation,
mean, and signal-to-noise ratio (SNR)) are provided. For each of these
parameters, either a single float value or a range of values are accepted.
Also, Gaussian blur can be added.

###  II-C Noise Parameters

PhaseNet takes dataset-specific parameters, such as the noise parameters,
which are defined by the acquired images.

To retrieve the noise parameters, we used the “measure” function in Fiji
software [11] on selected areas of the image. The bright parts of the image
belonging to the astrocyte were considered the foreground, the dark areas —
background (noise). The generator accepts a range for each noise parameter,
therefore the mean intensity was measured on different parts of the astrocyte,
as well as the mean and standard deviation of the noise were measured twice.
By definition, the SNR is the ratio of the average power of a signal to the
average power of a background noise. The SNR was calculated as the foreground
mean divided by the noise mean.

###  II-D Data Generation with Phantoms

To include the information about the object structure, we passed a crop of the
raw unaberrated image into the generator, where it was convolved using
synthetic PSFs. In our experiments with synthetic data, the unaberrated image
is simply the downloaded raw image. In an experimental setup with acquisition
of real data using an AO microscope, an unaberrated image would be an image
without aberrations introduced by a deformable mirror (or by another hardware,
e.g. a spatial light modulator). It is important to provide such an
unaberrated image to ensure that the only PSF present is the residual PSF of
the microscope because it is present in all the data coming from this
microscope, whereas the introduced aberrations appear in images only when the
AO microscope is instructed to add them.

Our approach for generating images using an actual image aims to add prior
information about the object shape, the residual PSF of the microscope, and
other noise coming from the microscope besides Gaussian noise, which the
generator synthesizes.

We generated synthetically aberrated 3D images of astrocytes by cropping the
input image inside the generator, convolving the crop with random synthetic
PSFs, and adding random Gaussian noise. We refer to the input image that was
passed into the generator as a “phantom”. Training with two or more images
passed as phantoms is named “training with multiple phantoms” in this paper.

Generation of the images simulated the acquisition of a dataset with 11
Zernike modes. We limited our dataset to the Zernike polynomials up to the
fourth order because higher-order polynomials have a low impact on the image
quality. The trivial modes (piston, x-tilt, y-tilt, and defocus) were
excluded, as has been done in the original PhaseNet project.

In this research, we investigate whether the training with a phantom or
multiple phantoms is necessary for application of PhaseNet to complex objects.

###  II-E Restoration via Deconvolution

The goal of predicting aberrations is to precisely restore distorted images.
Restoration can be achieved using deconvolution, an operation inverse to
convolution, whose purpose is to restore a convolved image to its state before
the convolutional kernel was applied. In our case, the convolutional kernels
are PSFs (the synthetically introduced PSFs and the residual PSF of the
microscope), therefore deconvolution with predicted PSFs can be used for
restoring the images.

In our workflow, we incorporate the Python-based library Flowdec [12] which
implements the 3D deconvolution using the Richardson–Lucy algorithm [13, 14]
for restoration of aberrated images. Compared to commonly used
DeconvolutionLab2 [15], Flowdec has achieved the results faster. The
researchers provided a Jupyter notebook on their GitHub page comparing Flowdec
with DeconvolutionLab2 on a computer with 64G RAM, Intel Xeon CPU (x2), and
Nvidia GTX 1080 GPU (x2), which demonstrated that Flowdec was 55 times faster
on average.

We chose Flowdec because it is fast, able to take the PSFs generated by our
synthetic generator, and can be called directly from Python scripts or Jupyter
notebooks.

Our approach of training on phantoms and restoring images using Flowdec is
schematically shown in Fig. 2.

Figure 2: Schematic of training on phantoms followed by image restoration.
Training: the generator takes the input parameters of the microscopy system,
parameters for data generation, and phantoms (unaberrated 3D images of the
study samples). At each iteration, it convolves a PSF with a random 3D crop of
a randomly chosen phantom and outputs an aberrated image with the ground
truth, which are passed to the CNN for supervised training. Restoration: the
trained CNN predicts the vector of the aberration modes with their amplitudes
for an aberrated image of the sample. The CNN predictions can be used for
generating the PSF, which can then be used as a deconvolution kernel in
restoration. The restoration was done using the Flowdec library in Python. The
green parts of the scheme denote our modifications of PhaseNet.

##  III Data

###  III-A Source Images

For this research, 3D images with known characteristics of the microscopy
setup (refractive index, $\lambda$, NA, voxel size) were taken from the open
database Cell Image Library (cellimagelibrary.org). The images were of rat
hippocampus at different ages taken at various resolutions [10]. We chose two
images [16, 17] on which the voxel size, fluorescent dye, and developmental
stage were identical. Further, we refer to them as the first [16] and the
second [17] raw image or dataset.

The images of astrocytes in a 4-week-old rat hippocampal area CA1
intracellularly injected with Lucifer Yellow [18] (${\lambda}=488$ nm) were
acquired with a single photon confocal microscope with $\mathrm{NA}=1.4$, oil
immersion medium, and voxel size of (0.068519, 0.068519, 0.2) ${\mu}$m in
$(x,y,z)$. For more details, please refer to [10, 16, 17].

###  III-B Training and Test Sets

The astrocyte images had a large size, and objects were not present in all
parts of the images. To ensure that the input images for PhaseNet contained an
object, a smaller 3D area was cropped out in both raw images. A 2D slice of
the second image [17] with a manual cropping is shown in Fig. 3. To speed up
the training, even smaller crops of (32, 32, 32) pixels (as in [5]) were made,
in the center or off-center, before they were passed to the CNN.

Figure 3: Cropping out the object in the 3D image [17] (the central 2D slice
shown in Fiji).

To the best of our knowledge, only images of fluorescent beads with
aberrations experimentally introduced by AO microscopes were publicly
available as the time of preparation of the manuscript. Acquisition of
bioimages of such kind is time-consuming and expensive because it requires an
AO microscope in the first place, as well as a suitable fixed fluorescent
biosample. Due to unavailability of bioimages with experimentally introduced
aberrations for our research, we synthesized different randomized training,
validation, and test data using our generation approach described in Section
II-D by passing the downloaded images of astrocytes as phantoms to the
generator. The generator was able to synthesize realistic images of
fluorescent beads according to our comparison of the synthetic images to the
real ones available in the publication [5]. Based on this, we assumed that the
generator could simulate realistic aberrations.

Images in the training and validation sets had random amplitudes for each of
11 aberration modes. The test sets consisted of 11 series of images for each
aberration mode in order to make the visualization of the predictions more
convenient and to compare them with the plots from the original PhaseNet
paper. The test images in the series for each aberration mode were convolved
with synthetic PSFs, which had only a random amplitude of the target mode
whereas other amplitudes were not present.

In addition to our synthetic data generation approach, we modified the
training procedure so that for validation the network can be provided with a
data generator that is different from the train set generator, or supplied
with an array of images. The latter functionality can be used in other
research, where real experimental images are available, for detection of
models’ overfitting to synthetic train data.

##  IV Results

###  IV-A Applying PhaseNet

We retrained PhaseNet with the new parameters of the microscope that determine
the synthetic PSFs. The network expects a certain pixel size, therefore the
pretrained weights from the original PhaseNet study cannot be reused for data
with objects of a different scale. In addition, in the training configuration,
we changed the noise parameters in the data generator according to the
microscope setup (Section II-C). The model was trained on synthetic PSFs
(point-like objects) without any other modifications except for changing the
microscope’s and the noise parameters, and the pixel size. The training
procedure is shown in Fig. 1A.

Then, synthetically aberrated images of astrocytes were generated by cropping
the image and convolving the crop with random synthetic PSFs (Section II-D).
Generation of the images simulated the acquisition of the dataset consisting
of 11 test series of images for each aberration mode. We tested retrained
PhaseNet on this new data and evaluated the predictions using the mean squared
error (MSE) averaged across all predictions. The model was unable to predict
aberrations, resulting in the MSE = 0.001718 on random crops. A possible
reason could be overfitting on the training object shape: the network learned
how aberrated point sources look like but could not generalize to nonspherical
objects in the test set. Alternatively, the cause could be the inability of
PhaseNet to understand the aberrations on complex shapes, which is
investigated below.

For visualization purposes, the results for one of the aberration modes (mode
8 in ANSI nomenclature, horizontal coma) are plotted in Fig. 4. In the plot,
blue stars on the diagonal line are the desired ground truth predictions for
the present aberration mode, yellow dots are the predicted values. The inset
in the bottom right corner shows statistics of predictions for all modes
(minimum, maximum, median, the first and third quartiles, and outliers).
Ideally, the variation of predictions for the modes, except the present one,
should be minimal and centered around zero $\mu$m.

Figure 4: Results of PhaseNet trained on points and tested on astrocytes. Blue
stars on the diagonal line are the desired ground truth predictions for the
present aberration mode (mode 8, horizontal coma). Yellow dots indicate the
predicted values. The inset shows the statistics of predictions for all the
modes.

###  IV-B Training with Phantom Image

To investigate whether PhaseNet can learn the aberrations on an object of a
nontrivial shape, we trained the model with the image [16] as the phantom. In
this experiment, the generator was instructed to crop the center of the
phantom, meaning that the object on the images remained the same and only the
introduced aberrations varied. The results (Fig. 5) proved that the model can
learn the representations of aberrations on a nontrivial object, such as a
fragment of an astrocyte, if the training data includes the information about
the object.

To understand whether the model overfits on the object that was present in the
supplied phantom, the model was tested on the random crops (“jitter” = “true”,
see Section II-B) from the first image [16] around the center and from the
second image [17] (unseen by the model during training). A significant drop in
performance from MSE = 0.000027 on the central crop to 0.000299 on the random
crops around the center and 0.002186 on the unseen data indicated that the
model strongly overfits on the object if its shape does not vary.

Figure 5: Results of training on the phantom with central cropping and testing
on the central crop convolved with random PSFs. Blue stars on the diagonal
line are the desired ground truth predictions for the present aberration mode
(mode 8, horizontal coma). Yellow dots indicate the predicted values. The
inset shows the statistics of predictions for all modes.

###  IV-C Training with Phantom and Random Cropping

We explored whether the model can learn the aberrations if the shape of the
object varies. The variety of the training data was enlarged by the random
off-center cropping (“jitter” = “true”, see Section II-B). This experiment
proves (Fig. 6) that the model can perform well on objects of varying shapes
if it has seen these shapes during the training (Table I, row “Random crop”).

Figure 6: Results of the training on random crops of the phantom and testing
on random crops. Blue stars on the diagonal line are the desired ground truth
predictions for the present aberration mode (mode 8, horizontal coma). Yellow
dots indicate the predicted values. The inset shows the statistics of
predictions for all modes.

###  IV-D Training with Multiple Phantoms

We supplied the generator with two phantoms: one used in the previous
experiments and another cropped from the second image. Our generator randomly
selected a phantom from the list of images at each iteration of generating a
synthetically aberrated image (at each training step). The phantoms were
randomly cropped. Testing on random crops from the first and second images
with random aberrations gave almost the same MSE on both sets: 0.000129 on the
first one and 0.000131 on the second. Clearly, this model performed better on
the second dataset than the previous CNNs because it had seen the objects from
this dataset during the training.

To further prove the idea that the model overfits on the shapes of the objects
from the training set, we tested it on another region from the second image
that did not overlap with the crops used for training and validation. As a
result, we observed a significant drop in the MSE to 0.001554. Nevertheless,
training on more various data led to a lower MSE compared to previous
experiments.

Additionally, we tried training on the central region of the first image,
validation on the central region of the second image, and testing on another
region of the second image. Overfitting was observed after 400 epochs, which
supports our assumption that the network overfits on the seen shapes.

The results of all experiments are summarized in Table I. The models were
trained for 50,000 epochs with 5 steps per epoch and the batch size of 2
images. Training of one model on a laptop with Intel Core i7-8850H CPU, single
NVIDIA Quadro P3200 GPU with 8 Gb memory, and 64 Gb RAM lasted 30 hours.

TABLE I: Comparison of results | Mean MSE of testing on:
---|---
Training on | Central crop | Random crop | Unseen set
Points (III-A) | naa | naa | 0.001718
Central crop (III-B) | 0.000027 | 0.000299 | 0.002186
Random crop (III-C) | 0.000040 | 0.000045 | 0.001793
Two phantomsb (III-D) | 0.000090 | 0.000129 | 0.001554
aSets with central and random crops were unseen.
bTraining with multiple phantoms, two were provided.

###  IV-E Restoration of the Aberrated Images

As shown in Fig. 2, our workflow includes the restoration of aberrated 3D
images. An aberrated image is cropped in three dimensions and passed to the
model which predicts the aberration modes with their amplitudes. The
predictions are used for generating the PSF. The full-size aberrated 3D input
image is then deconvolved with the PSF image. An example of a restoration of
the synthetically aberrated image is shown in Fig. 7.

Figure 7: Restoration of an aberrated 3D image using deconvolution with the
predicted 3D PSF. 2D slices are shown.

###  IV-F Restoration of the Source Images

It is possible to improve the quality of the initial source images by the
deconvolution of the initial raw image with the PSF predicted by our model
trained on the image as the phantom. To do so, the raw image has to be cropped
to the size accepted by the network. The network takes the crop without
additional modifications and outputs the vector of amplitudes for the residual
PSF. The PSF is generated based on the prediction and used as the
deconvolutional kernel.

The raw source image [16] and its restoration are shown on Fig. 8. The
restoration of the manually cropped area (cropping is shown on Fig. 3) of the
second source image [17] is visualized on Fig. 9.

Figure 8: Restoration of the raw 3D image using deconvolution with the
predicted 3D PSF. The central 2D slices are shown.

Figure 9: Restoration of the cropped source 3D image (Fig. 3) using
deconvolution with the predicted 3D PSF. The central 2D slice is shown.

##  V Discussion

The images in this study are of objects of a complex structure (astrocytes),
and crops from the images contain smaller objects varying in shape and size.
Due to unavailability of such bioimages with aberrations experimentally
introduced by an AO microscope, the aberrated images of astrocytes were
synthetically generated.

For correction of substantially aberrated bioimages, we recommend to test our
method on experimentally aberrated data (e.g. acquired by an AO microscope
with a deformable mirror or using a spatial light modulator). Nevertheless,
our approach demonstrates convincing improvements of the quality of the source
raw images where a small residual PSF of the microscopic system should be
present.

The results and conclusions are valid for biological samples of a structure
and complexity similar to the ones of the protoplasmic astrocytes in the
hippocampus of Rattus norvegicus. However, training PhaseNet without a phantom
may work on approximately spherical objects or on samples that contain beads
by design or spherical structures by nature.

We noticed that in the configuration file of the “single mode model” on the
official Github page with the code for the PhaseNet paper [5], the range for
generating synthetic PSFs was set to (–0.075 $\mu$m, 0.075 $\mu$m), and
decided to train our models in the same range.

We used the learning rate of 0.0003 in training our models and observed that
it leads to better results than with a learning rate of 0.0001 suggested in
the PhaseNet publication. To determine the number of training epochs, we tried
training for 10,000, 50,000, and 70,000 epochs. The decrease in the loss
function between 50,000 and 60,000 epochs was insignificant; between 60,000
and 70,000, it reached a plateau.

##  VI Conclusions

This research describes a method for estimation of monochromatic aberrations
on 3D biological images and image restoration by correction of aberrations via
deconvolution with the estimated PSFs.

PhaseNet trained on synthetically generated PSFs that represent only the
effects of aberrations on point light sources may not be directly applicable
to images of biological samples with complex shapes. Providing the generator
with an unaberrated image of the sample (phantom) that represents the
information about the shape of the object allows applying PhaseNet to
nontrivial biological samples.

We suggest that PhaseNet trained on synthetically aberrated images of
astrocytes overfits on the shapes of the seen objects and cannot perform well
on the completely unseen data. To ensure a good performance, we propose to
include the information about the object’s shape by supplying the generator
with multiple phantoms representing the variety of shapes.

We modified the data generator from the PhaseNet study so that it can take a
list of images as phantoms. The training procedure was modified in the way
that the model can be trained with a user-defined validation dataset.

Additionally, we include the guidelines for image restoration using the
model’s predictions. Moreover, we demonstrate that our approach is suitable
for restoration of raw microscopic images without additionally introduced
aberrations if the parameters of the microscope (refractive index, $\lambda$,
NA, voxel size) are known.

## Acknowledgments

We thank Coleman W. Broaddus, Alexandr Dibrov, and Debayan Saha (MPI-CBG) for
scientific discussions at different stages of the research. We thank Coleman
W. Broaddus, Tom Burke (MPI-CBG), and Fabio Cunial (Broad Institute) for
providing feedback on the manuscript.

## References

  * [1] F. Zernike, “Beugungstheorie des schneidenver-fahrens und seiner verbesserten form, der phasenkontrastmethode,” _physica_ , vol. 1, no. 7-12, pp. 689–704, 1934.
  * [2] P. Pozzi, D. Wilding, O. Soloviev, H. Verstraete, L. Bliek, G. Vdovin, and M. Verhaegen, “High speed wavefront sensorless aberration correction in digital micromirror based confocal microscopy,” _Opt. Express_ , vol. 25, pp. 949–959, 01 2017.
  * [3] I. Vishniakou and J. D. Seelig, “Differentiable model-based adaptive optics for two-photon microscopy,” _Optics Express_ , vol. 29, no. 14, pp. 21 418–21 427, 2021.
  * [4] L. Murray, C. Dainty, J. Coignus, and F. Felberer, “Wavefront correction of extended objects through image sharpness maximisation - art. no. 60181a,” _Proceedings of SPIE - The International Society for Optical Engineering_ , vol. 6018, 12 2005.
  * [5] D. Saha, U. Schmidt, Q. Zhang, A. Barbotin, Q. Hu, N. Ji, M. J. Booth, M. Weigert, and E. W. Myers, “Practical sensorless aberration estimation for 3d microscopy with deep learning,” _Optics express_ , vol. 28, no. 20, pp. 29 044–29 053, 2020.
  * [6] Q. Hu, M. Hailstone, and M. J. Booth, “Efficient and versatile aberration correction through sensorless adaptive optics,” in _Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XXIX_ , vol. 11966\. SPIE, 2022, pp. 21–26.
  * [7] A. J.-W. Whang, Y.-Y. Chen, T.-H. Yang, C.-T. Lin, Z.-J. Jian, and C.-H. Chou, “Zernike coefficient prediction technique for interference based on generation adversarial network,” _Applied Sciences_ , vol. 11, no. 15, p. 6933, 2021.
  * [8] M. Liu, D. N. Lopez, and G. C. Spalding, “Experimental implementation of wavefront sensorless real-time adaptive optics aberration correction control loop with a neural network,” in _Emerging Topics in Artificial Intelligence 2020_ , G. Volpe, J. B. Pereira, D. Brunner, and A. Ozcan, Eds., vol. 11469, International Society for Optics and Photonics. SPIE, 2020, p. 114691S. [Online]. Available: https://doi.org/10.1117/12.2569647
  * [9] Q. Tian, C. Lu, B. Liu, L. Zhu, X. Pan, Q. Zhang, L. Yang, F. Tian, and X. Xin, “Dnn-based aberration correction in a wavefront sensorless adaptive optics system,” _Opt. Express_ , vol. 27, no. 8, pp. 10 765–10 776, Apr 2019\. [Online]. Available: http://opg.optica.org/oe/abstract.cfm?URI=oe-27-8-10765
  * [10] E. A. Bushong, M. E. Martone, and M. H. Ellisman, “Maturation of astrocyte morphology and the establishment of astrocyte domains during postnatal hippocampal development,” _International Journal of Developmental Neuroscience_ , vol. 22, no. 2, pp. 73–86, 2004.
  * [11] J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid _et al._ , “Fiji: an open-source platform for biological-image analysis,” _Nature methods_ , vol. 9, no. 7, pp. 676–682, 2012.
  * [12] E. Czech, B. A. Aksoy, P. Aksoy, and J. Hammerbacher, “Cytokit: A single-cell analysis toolkit for high dimensional fluorescent microscopy imaging,” _bioRxiv_ , 2018. [Online]. Available: https://www.biorxiv.org/content/early/2018/12/14/460980
  * [13] W. H. Richardson, “Bayesian-based iterative method of image restoration,” _JoSA_ , vol. 62, no. 1, pp. 55–59, 1972.
  * [14] L. B. Lucy, “An iterative technique for the rectification of observed distributions,” _The astronomical journal_ , vol. 79, p. 745, 1974.
  * [15] D. Sage, L. Donati, F. Soulez, D. Fortun, G. Schmit, A. Seitz, R. Guiet, C. Vonesch, and M. Unser, “Deconvolutionlab2: An open-source software for deconvolution microscopy,” _Methods_ , vol. 115, pp. 28–41, 2017.
  * [16] Bushong, Eric and Martone, Maryann and Ellisman, Mark, “Ccdb:1065, rattus norvegicus, protoplasmic astrocyte. cil. dataset.” https://doi.org/doi:10.7295/W9CCDB1065, 2002, accessed: 2022-20-07.
  * [17] ——, “Ccdb:1073, rattus norvegicus, protoplasmic astrocyte. cil. dataset.” https://doi.org/doi:10.7295/W9CCDB1073, 2002, accessed: 2022-20-07.
  * [18] P. V. Belichenko and A. Dahlström, “Studies on the 3-dimensional architecture of dendritic spines and varicosities in human cortex by confocal laser scanning microscopy and lucifer yellow microinjections,” _Journal of neuroscience methods_ , vol. 57, no. 1, pp. 55–61, 1995.


# Long-range interaction induced collective dynamical behaviors

K. Sathiyadevi1, V. K. Chandrasekar1, D. V. Senthilkumar2 and M. Lakshmanan3
1Centre for Nonlinear Science & Engineering, School of Electrical &
Electronics Engineering, SASTRA Deemed University, Thanjavur -613 401, Tamil
Nadu, India.
2School of Physics, Indian Institute of Science Education and Research,
Thiruvananthapuram-695551, India.
3Centre for Nonlinear Dynamics, School of Physics, Bharathidasan University,
Tiruchirappalli - 620 024, Tamil Nadu, India<EMAIL_ADDRESS>
###### Abstract

Long-range interacting systems are omnipresent in nature. We investigate here
the collective dynamical behavior in a long-range interacting system
consisting of coupled Stuart-Landau limit cycle oscillators. In particular, we
analyze the impact of a repulsive coupling along with a symmetry breaking
coupling. We report that the addition of repulsive coupling of sufficient
strength can induce a swing of the synchronized state which will start
disappearing with increasing disorder as a function of the repulsive coupling.
We also deduce analytical stability conditions for the oscillatory states
including synchronized state, solitary state, two-cluster state as well as
oscillation death state. Finally, we have also analyzed the effect of power-
law exponent on the observed dynamical states.

##  1 Introduction

Complex dynamical systems exhibit a variety of dynamical patterns through
distinct interesting transitions due to the nature of coupling architectures.
Among the various coupling geometries, investigations on coupled networks with
global (all to all) coupling have become an active area of research in recent
times, because of their prevalence in many natural systems [1, 2, 3, 4].
Moreover, such a type of network architecture is often used to study the
collective properties of chemical reactions [5, 6], physical systems [7, 8, 9,
10], electronic devices [11] and biological systems [12, 13]. Globally coupled
dynamical systems exhibit diverse dynamical patterns such as synchronization,
clustering and many more states including one of the emergent hybrid dynamical
behavior called chimera state [14, 15, 16]. In addition, the network of
coupled dynamical systems exhibit different transitions with coexisting
distinct dynamical states. But still the underlying mechanism behind such
transitions, and multistabilities among the distinct dynamical behaviors,
remain unclear. Further it is also known that globally (mean-field) coupled
systems can act as long-range interacting systems under certain limitations
[17]. Studies on long-range interaction have also been an active area of
research over the last several decades [18, 19, 20, 21]. Generally, any
complex interaction decaying as a power law at large distances, without any
specific length, can be considered as long range interaction. One can find
such systems in astrophysics, plasma physics, atomic physics, nuclear physics
and hydrodynamics [22, 23, 24, 25]. Due to the widespread nature of long-range
interacting systems, it is essential to understand the dynamical behavior of
an ensemble of oscillators with this kind of interactions. The effect of such
interactions can be well understood with mean-field (global) coupling [26,
27].

On the other hand, studies on synchronization and its control have been
receiving much attention due to their potential applications in many physical,
biological, social and even in man-made systems [28]. Particularly, a large
scale neuronal synchronization causes several brain disorders such as
Parkinson’s disease, tremor and epilepsy [29, 30, 31, 32]. Though many
investigations have been made to understand their intricacies in different
systems, still it demands a deep knowledge to control the mechanism of
synchronized states. Therefore a continuous effort is still on to explore many
aspects associated with the phenomenon of synchronization [33, 34]. Out of
numerous investigations on synchronized states, Daido and Nakanishi had found
a novel interesting behavior called a “swing of the synchronized state” in a
large population of globally coupled Stuart-Landau oscillators where the
inhomogeneity among the synchronized regions emerges by introducing diffusion
in a globally coupled system in the presence of nonisochronicity parameter.
The investigation by Daido and Nakanishi revealed that the synchronized state
is destabilized while increasing the coupling strength which again gets
stabilized upon increasing the coupling strength further. The destabilization
of stable synchronized state and again re-stabilization of destabilized
synchronized state is called the swing of the synchronized state [35, 36].
Later, such an interesting dynamical behavior was also found in a globally
coupled system of Stuart-Landau oscillators with symmetry breaking coupling in
the presence of nonisochronicity parameter by Premalatha et al [37]. These
authors reported that the symmetry breaking coupling facilitates the widening
of the dynamical regime by twice as compared with the symmetry preserving
coupling, and an increase of disorderliness when increasing the
nonisochronicity parameter. Besides the above, recent investigations revealed
that the trade-off between repulsive and attractive couplings can exhibit
diverse collective dynamical states, including spontaneous symmetry breaking,
along with distinct multistablities in a single network [38, 39, 40]. In the
present study, we explore whether the addition of repulsive coupling along
with symmetry breaking coupling can induce the swing by mechanism and, if so,
what is its influence on the degree of disorder among the coupled systems.
Indeed, we find the emergence of swing of the synchronized state while
increasing the strength of the repulsive coupling up to a sufficient value
beyond which it will disappear with increasing disorderliness at strong
repulsive coupling.

In this report, we investigate the dynamical transitions in a network of long-
range coupled Stuart-Landau oscillators with the addition of repulsive link in
the symmetry breaking coupling. Initially, we show the emergence of a swing of
synchronized state for sufficient strength of the repulsive coupling in the
absence of power-law exponents and increasing disorder with the increase of
the strength of the repulsive coupling. Further, the swing of the synchronized
state disappears at strong repulsive coupling strength and the systems attain
a completely disordered state. In the swing of the synchronized state the
dynamical transitions are observed through the synchronized state to solitary
state and cluster state and finally re-emergence of the synchronized state. We
have also deduced the analytical stability conditions for the oscillatory
states of synchronized states, solitary state and cluster states. Finally, we
obtain the stability conditions for oscillation death state. We also identify
the dynamical transitions for different power-law exponents and report that
the system becomes completely desynchronized when increasing the power-law
exponent to appropriate higher values.

The rest of the paper is organized as follows. In Sec. 2, we introduce our
model of long-range coupled Stuart-Landau oscillators. The numerical results
for swing by mechanism and global behavior of the coupled system are discussed
in the absence of power-law exponents (i.e. the system is all-to-all coupled)
in Sec. 3. The corresponding analytical results are discussed in Sec. 4. We
analyze the dynamical transitions by varying the power-law exponents in Sec.
5. Finally, we summarize our findings in Sec. 6.

##  2 The Model

Stuart-Landau limit cycle oscillator is a general, paradigmatic model, which
can be used to model a variety of nonlinear dynamical systems near the Hopf
bifurcation point [41, 42, 43, 44]. We consider a one-dimensional ring network
of coupled Stuart-Landau oscillators and investigate its dynamical behavior by
including a repulsive link along with a symmetry breaking coupling. The
proposed coupling is a symmetry breaking coupling which breaks the rotational
symmetry of the coupled system explicitly. The governing equations of the
network can be represented as

$\displaystyle\quad\dot{z}_{j}=f(z_{j})+\frac{\kappa}{\eta(\alpha)}\,\bigg{[}\sum_{j\neq
k}\frac{1}{|j-k|^{\alpha}}\big{(}Re({z_{k}}-{z}_{j})-i\,q\,Im({z_{k}}-{z}_{j})\big{)}\bigg{]},\,\,j,k=1,2,...,N,$
(1)

where $f(z_{j})=(\lambda+i\omega)z_{j}-(1-ic)|z_{j}|^{2}z_{j}$, the state
variables being $z_{j}=x_{j}+iy_{j}\in C$. In the absence of coupling the
considered system preserves symmetry (rotational invariance), that is the
system equation remains invariant under the transformation $z_{j}\rightarrow
z_{j}e^{i\theta}$. In contrast, the rotational invariance is not preserved
upon introducing the couplings and hence these types of couplings are known as
symmetry breaking couplings [45, 46]. Here $N$ is the total number of the
oscillators in the network and $\eta(\alpha)=\sum_{k\neq
j}^{N}{|j-k|^{-\alpha}}$ is the normalization constant. When the strength of
the interaction lies between $0\leq\alpha\leq 1$, the interaction is said to
be of long range. On the otherhand $\alpha>1$ is considered as the short
range. The edge of long-range interaction is achieved through $\alpha=0$, when
the system becomes mean-field (globally) coupled whereas
$\alpha\rightarrow\infty$ indicates that systems are coupled only to the
nearest neighbors [47, 48, 49, 50]. Here, $\kappa=((N-1)/N)\epsilon$, where
$\epsilon$ is the coupling strength and $q$ is the strength of the repulsive
coupling. $\lambda$ and $\omega$ correspond to the bifurcation parameter and
the natural frequency of the system, respectively. c is the nonisochronicity
parameter. Throughout the manuscript, the parameter values are fixed as
$\lambda=1.0$, $\omega=2.0$, $c=2.5$ and $N=100$. The initial state for the
state variables $x_{i}$ and $y_{i}$ are chosen randomly between -1 to 1 and
Runge-Kutta fourth order integration scheme is used for all our simulations
with a time step 0.01.

##  3  Dynamical transitions for power-law exponent $\alpha=0$

###  3.1  Swing of synchronized states

Figure 1:  Spatio-temporal plots for $q=0.1$ as a function of $\epsilon$: (a)
synchronized state for $\epsilon=1.05$, (b) solitary state for
$\epsilon=1.35$, (c) two cluster state for $\epsilon=1.65$ and (d)
synchronized state for $\epsilon=2.05$. The corresponding frequency plots are
shown in the lower panels (e)-(h). The insets in (e)-(h) denote the snapshots
of the variables $x_{i}$. Figure 2:  Phase-portrait trajectories of the
representative oscillators for (a) synchronized state, (b) solitary state and
(c) two-cluster oscillatory state. Dashed line denotes the phase-portrait
trajectory of the synchronized group and solid line in Fig. 2(b) represents
the solitary oscillators, while the solid line in Fig. 2(c) denotes second
cluster. Filled circles denote the Poincaré points corresponding to the
synchronized group whereas Poincaré points for the quasiperiodic oscillations
of the solitary oscillators are denoted by filled diamonds. The parameter
values are the same as in Figs. 1(a)-(c).

The emergence of swing by mechanism due to diffusive coupling has been
revealed both in the cases of symmetry preserving and symmetry breaking
couplings, which induce inhomogeneity among the synchronized regions resulting
in the swing of synchronized states [37]. In order to understand the impact of
the additional repulsive coupling on the emergence of the swing of the
synchronized states in the ensemble of Stuart-Landau oscillators with symmetry
breaking coupling, at first, we have plotted the space-time evolution as a
function of the coupling strength $\epsilon$ for the repulsive coupling
$q=0.1$ (see Figs. 1(a)-1(d)). Further to distinguish the distinct dynamical
states, we have estimated the frequencies of the oscillators corresponding to
each of the observed dynamical states, which are depicted in Figs. 1(e)-1(h),
while the snapshots of the dynamical states are shown in their respective
insets. From Figs. 1(a) and 1(e), it is clear that a completely synchronized
state (SYN) emerges for the coupling strength $\epsilon=1.05$. In this state,
all the oscillators oscillate with the same frequency as is evident from Fig.
1(e), while the inset in Fig. 1(e) clearly illustrates that each oscillator in
the network oscillates with the same amplitude. Increasing the coupling
strength further to $\epsilon=1.35$, the system exhibits a solitary state (SS)
as shown in Figs. 1(b) and 1(f). A very few of the oscillators makes an
excursion away from the synchronized group and oscillate with their own
frequency and amplitude [51, 52], is evident from Fig. 1(f). In contrast to
[52], the solitary oscillators in Fig. 1(f) do not have identical mean
frequencies. In fact, their mean frequencies are distributed between four
different levels. Emergence of two-cluster (2C) state is observed on further
increasing the coupling strength to $\epsilon=1.65$ (see Figs. 1(c) and 1(g)).
It is to be noted that in the 2C state, the oscillators in the network
oscillate in two different groups with different amplitudes but with the same
frequency. By increasing the coupling strength further, again the system re-
enters into a completely synchronized state (see Figs. 1(d) and 1(h)),
characterized by identical amplitude and frequency. Thus, it becomes evident
that the phenomenon of swing of synchronized states, that is the re-emergence
of the synchronized state after being collapsed in a certain range of
parameters, persists even with the addition of repulsive coupling. For a
clearer insight on the nature of the observed dynamical states, we have
depicted the phase-space dynamics of the representative oscillators for each
of the observed states along with their Poincaré points in Fig. 2. The
representative oscillator from the synchronized group is denoted by a dashed
line, which always exhibits periodic oscillations. The filled circles
correspond to the Poincaré points of the synchronized group of oscillators.
The Poincaré points of the solitary oscillators are represented by filled
diamonds. The representative oscillators from the synchronized state
$(x_{2},x_{15})$ exhibiting periodic oscillations are evident from Fig. 2(a)
for $\epsilon=1.05$. In the case of the solitary state, the oscillator from
the synchronized group $(x_{2})$ exhibits periodic oscillations as in Fig.
2(a), while the solitary oscillator $(x_{15})$ exhibits quasi-periodic
oscillations. The closed loop of the Poincaré points confirms the quasi-
periodic nature of the solitary oscillator. The representative oscillators
$(x_{2},x_{18})$ from the two-cluster state (see Fig. 2(c)) indicate that the
oscillators from the distinct clusters follow different trajectories with
different amplitudes.

Further to quantify the dynamical transitions and to characterize the swing of
synchronized states, we have calculated the standard deviation as a function
of $\epsilon$ for three different values of the strength of the repulsive
coupling. The standard deviation is calculated using the relation,

$\displaystyle\sigma=\langle(\overline{|x_{j}-\overline{x_{j}}|^{2}})^{1/2}\rangle_{t},$
(2)

where the bar denotes the average over $1\leq j\leq N$ and $<\cdot>_{t}$
represents the time average.

Figure 3:  Standard deviation as a function of $\epsilon$ (after leaving the
transients time, $5\times 10^{4}$ units) for (a) $q=0.1$ (b) $q=0.3$ and (b)
$q=0.7$. SYN, SS and DS denote the synchronized state, solitary and
desynchronized states. 2C and 3C are the two and three cluster oscillatory
states. OD is the oscillation death state.

Initially, the dynamical transitions will be analyzed for a smaller value of
the strength of the repulsive coupling. The standard deviation is shown in
Fig. 3 as a function of the coupling strength $\epsilon$ for three different
values of the repulsive coupling. Fig. 3(a) corresponds to the strength of the
repulsive coupling $q=0.1$. It is evident from the null value of the standard
deviation in the range of $\epsilon\in(0.08,1.09)$ that all the oscillators in
the network evolve in synchrony with each other. Further, the finite value of
the standard deviation in a narrow range of $\epsilon\in(1.09,1.78)$ indicates
that all the oscillators in the network evolve independently (asynchronously).
The oscillators in this range of $\epsilon$ actually exhibit either solitary
state or two cluster state. The null value of the standard deviation beyond
$\epsilon=1.78$ elucidates that all the oscillators are in complete synchrony.
Thus the value of the standard deviation as a function of the coupling
strength $\epsilon$ corroborates the re-emergence of the synchronized state
after a short range of excursion away from the synchronized state, thereby
illustrating the phenomenon of the swing of the synchronized state.

Similarly, Figs. 3(b) and 3(c) are plotted for two distinct values of the
strength of the repulsive coupling (i.e. $q=0.3$ and $q=0.7$). It is clear
that increasing $q$ increases the range of disorder among the coupled
oscillators as a function of the coupling strength $\epsilon$ (see Figs.
3(a)-3(c)). The synchronized state re-emerges for $q=0.3$ as a function of the
coupling strength $\epsilon$ as indicated by the value of the standard
deviation in Fig. 3(b). It is also to be noted that the synchronized state
characterized by the null value of the standard deviation does not re-emerge
for a larger value of the repulsive coupling strength (see Fig. 3(b) for
$q=0.7$) thereby elucidating the disappearance of the swing by mechanism for
larger ranges of the strength of the repulsive coupling. In Figs. 3(a)-3(c),
the shaded regions denote the desynchronized (DS) or three-cluster (3C)
(disordered) states, whose spread increases with the strength of the repulsive
coupling. Also, we observed the oscillation death state at strong coupling
strengths (see Fig.3(c)) which will be detailed in the following.

Figure 4:  Spatio-temporal evolutions for $q=0.7$, as a function $\epsilon$,
(a) desynchronized (DS) state for $\epsilon=1.3$ (b) three cluster (3C) state
for $\epsilon=2.8$, (c) oscillation death state for $\epsilon=4.0$. The
corresponding snapshots of the variables $x_{i}$ are shown in the lower panels
(d)-(f).

To understand the dynamical transitions at even stronger repulsive coupling
strengths, we have plotted space-time plots and snapshots for $q=0.7$ for
three different values of $\epsilon$ in Fig. 4. From the figure, it is evident
that there is a transition from desynchronized state to oscillation death
state through the three-cluster state as a function of the coupling strength
$\epsilon$ for the chosen value of the strength of the repulsive coupling. The
coupled Stuart-Landau oscillators evolve independently as shown in Figs. 4(a)
and 4(d) for $\varepsilon=1.3$. The asynchronous oscillators get segregated
into three clusters for further increase in the attractive coupling as
depicted in Figs. 4(b) and 4(e) for $\epsilon=2.8$. Finally, the oscillators
randomly populate the upper and lower branches of the inhomogeneous steady
state as shown in Figs. 4(c) and 4(f) for $\epsilon=4.0$. The global dynamical
transitions will be delineated in the following.

###  3.2 Global dynamical behavior of coupled oscillators

Figure 5:  Two parameter phase diagram for globally coupled Stuart-Landau
oscillators in $(q,\epsilon)$ space. SYN, SS and DS are the synchronized,
solitary and desynchronized states, respectively. 2C and 3C denote the two-
cluster and three-cluster oscillatory states. The oscillation death is denoted
by OD. $R_{1}$, $R_{2}$, $R_{3}$ and $R_{4}$ are the multistability regions
between SYN-SS, SYN-2C, SYN-OD and DS-3C, respectively.

To unravel the global dynamical behavior of the coupled Stuart-Landau
oscillators, we have depicted a two parameter phase diagram in the
$(q,\epsilon)$ space in Fig. 5. The coupled Stuart-Landau oscillators evolve
in synchrony for lower values of attractive $(\epsilon)$ and repulsive $(q)$
couplings. For lower values of the strength of repulsive coupling ($q<0.5$),
the dynamical system admits the re-emergence of synchronized state via the
solitary and two-cluster states as a function of the coupling strength,
thereby corroborating the existence of the phenomenon of the swing of the
synchronization in a rather large range of the parameters. Further, there is a
coexistence of the synchronized state with the solitary state and two-cluster
state in the regimes $R_{1}$ and $R_{2}$, respectively. For a larger value of
$\epsilon$ (approximately around $\epsilon=4.5$), and $q$ around $0.45$, there
is a coexistence of the oscillation death state (indicated by the checked
regimes) with the synchronized state in the region $R_{3}$. The complete
synchronization region is obtained by choosing nearly the same initial states
for all the oscillators. For larger values of the repulsive coupling strength,
the coupled system exhibits transition from homogeneous oscillatory state
(synchronized state) to inhomogeneous steady state (oscillation death state)
via asynchronous and three-cluster states as a function of $\epsilon$. In the
oscillation death state the oscillators populate both the upper and lower
branches of the inhomogeneous steady state randomly. In addition, there is a
coexistence of desynchronized state with the three-cluster state in the region
denoted by $R_{4}$. Thus it is evident that smaller values of the repulsive
coupling strength favors the onset of the phenomenon of the swing of the
synchronized state. On the other hand, stronger repulsive coupling favors the
independent evolution of the oscillators thereby preventing them from
entrainment with each other in a large range of the attractive coupling
$\epsilon$ as is evident from Fig. 5. In addition, the strong repulsive
coupling favors the onset of the inhomogeneous death state in a large range of
$\epsilon$.

The numerical boundaries of each regions were distinguished by using the
notions of $D$-factor [53] and the strength of incoherence $(SI)$ [37, 54].
Initially, the oscillatory and oscillation death state are identified through
the $D$-factor. Then the distinct oscillatory dynamical behaviors are
classified from the $SI$, which takes the values $0$ and $1$ for synchronized
and desynchronized states, respectively. For the two-cluster (2C), three
cluster (3C) and solitary state (SS), the values of $SI$ lies between
$0<SI<1$. Further, 2C, 3C and SS states are distinguished by using the removal
of discontinuities [54]. The analytical stability conditions for the observed
dynamical states are detailed in the next section.

##  4 Analytical results

In the following, we will deduce the analytical stability conditions for the
observed dynamical states.

###  4.1 Stability conditions for distinct oscillatory states

To start with, we consider the solitary state which is characterized by a very
few oscillators making excursion away from the majority of the synchronized
oscillators. Hence, the equation corresponding to the coupled Stuart-Landau
oscillators system (1), can be decomposed as,

$\displaystyle\frac{dz_{s}}{dt}=f(z_{s})+\epsilon[pRe(Z_{s})+\bar{p}Re(Z_{d})-Re(z_{s})]-iq\epsilon[pIm(Z_{s})+\bar{p}Im(Z_{d})-Im(z_{s})],$
(3)
$\displaystyle\frac{dz_{d}}{dt}=f(z_{d})+\epsilon[pRe(Z_{s})+\bar{p}Re(Z_{d})-Re(z_{d})]-iq\epsilon[pIm(Z_{s})+\bar{p}Im(Z_{d})-Im(z_{d})],$
(4)

where $z_{s}$ represents the state of synchronized oscillators and $z_{d}$
represents the desynchronized (solitary) oscillators. Further, we have
$f(z_{s})=(1+i\omega)z_{s}-(1-ic)|z_{s}|^{2}z_{s}$ and
$f(z_{d})=(1+i\omega)z_{d}-(1-ic)|z_{d}|^{2}z_{d}$, $\bar{p}=1-p$,
$s=1,2,...l_{1}$, $d=1,2,...l_{2}$, $l_{1}+l_{2}=N$. Here $N$ is the total
number of oscillators in the network, $p=l_{1}/N$, $\bar{p}=l_{2}/N$,
$p+\bar{p}=1$. Also $Z_{s}=\frac{1}{l_{1}}\sum_{s=1}^{l_{1}}z_{s}$ and
$Z_{d}=\frac{1}{l_{2}}\sum_{d=1}^{l_{2}}z_{d}$.

If we choose $p\approx 1$, then $\bar{p}\approx 0$ and hence the synchronized
oscillators are decoupled from the solitary oscillators (see Eq. (3) and Eq.
(4)) [35, 36, 37, 55, 56]. As a consequence, the solution of the
synchronization state can be deduced as $z_{s}=e^{i(\omega+c)t}$. Now, upon
substituting $z_{s}$, Eq. (4) can be rewritten as

$\frac{dz_{d}}{dt}=f(z_{d})+\epsilon[Re(e^{i(\omega+c)t})-Re(z_{d})]-iq\epsilon[Im(e^{i(\omega+c)t})-Im(z_{d})].$
(5)

The above equation corresponding to the solitary oscillator is also found to
admit the solution of the synchronized state $z_{d}=z_{s}=e^{i(\omega+c)t}$,
which may be stable or unstable depending on the strength of the repulsive
coupling. Thus it is evident that the above synchronization manifold is one of
the possible dynamical states of the coupled Stuart-Landau oscillators.

In order to unravel the other possible collective dynamical behaviors, we
apply the ansatz $z_{d}(t)=Z(t)e^{i(\omega+c)t}$, where $Z(t)$ is the
amplitude of the desynchronized state. Inserting this into Eq. (5), we obtain

$\displaystyle\frac{dZ}{d\tau}=(1-ic)Z-(1-ic)|Z|^{2}Z+\frac{1}{2}(1-Z)\epsilon\bar{q}+\frac{1}{2}(1-\bar{Z})\epsilon\bar{q}e^{-2i(\omega+c)t},$
(6)

where $\bar{q}=1-q$. Assuming that $\omega+c>>1$, we can averge Eq. (6) over
the fast oscillation period $[0,2\pi/(\omega+c)]$. This yields

$\displaystyle\frac{dZ}{d\tau}=(1-ic)Z-(1-ic)|Z|^{2}Z+\frac{1}{2}(1-Z)\epsilon\bar{q}.$
(7)

Now, Eq. (7) has the following three fixed points

$\displaystyle Z_{1}=1,$ (8)
$\displaystyle|Z_{2,3}|^{2}=\frac{(1+c^{2})-\epsilon\bar{q}\mp\sqrt{(1-c^{2}-\epsilon\bar{q})^{2}-(1+c^{2})\epsilon^{2}\bar{q}^{2}}}{2(1+c^{2})}.$
(9)

The stability of the fixed point $Z_{1}=1$ corroborates the stability of the
complete synchronization manifold and the other two fixed points determine the
stability of desynchronized (solitary) oscillator which exist when,

$\displaystyle\epsilon<\epsilon_{SN}=\frac{1-c^{2}(q-1)\sqrt{(1+c^{2})^{3}(q-1)^{2}}+q}{(1+c^{2})(q-1)^{2}}.$
(10)

Carrying out the linear stability analysis of Eq. (7), we find that the
determinant of the Jacobian matrix for the fixed points $Z_{2}$ and $Z_{3}$
can be expressed as

$det(J_{Z}{{}_{{}_{2,3}}})=\frac{c^{6}-1-d_{1}c^{4}+d_{2}c^{2}-2\epsilon(q-1)\mp\sqrt{\Delta}d_{3}}{2(c^{4}-1)},$
(11)

where $d_{1}=\epsilon^{2}(q-1)^{2}-1-2\epsilon(q-1)$,
$d_{2}=\epsilon^{2}(q-1)^{2}-1$, $d_{3}=1+c^{2}-\epsilon$ and
$\Delta=(c^{2}-1)^{2}(1+c^{4}-c^{2}(\epsilon^{2}(q-1)^{2}-2(1+\epsilon)(q-1))+2\epsilon(q-1))$.
The value of the determinant is negative for $Z_{2}$, and positive for
$Z_{3}$, which imply that the fixed point $Z_{2}$ turns out to be a saddle
whereas the other fixed point $Z_{3}$ is either an unstable node or a focus
for $|c|\leq 1$. Further the unstable node or focus turns out to be a stable
attractor when $|c|>1$ and $\epsilon>\epsilon_{HB}$, where $\epsilon_{HB}$ is
given by

Figure 6:  (a) Analytical stability curves for distinct dynamical states in
$(q,\epsilon)$ space. Dashed and dotted curves represent the saddle node and
Hopf bifurcation curves. The solid curve is the numerically obtained boundary.
The lines connected by filled triangles and open circles denote the boundary
for synchronized and oscillation death states, respectively. $R_{1}$, $R_{2}$
and $R_{3}$ are the multistability regions between SYN-SS, SYN-2C and SYN-OD,
respectively. (b) The largest Lyapunov exponent $\Lambda_{max}$ in the
$(q,\epsilon)$ space. The dotted line in Fig. (b) denotes the stable boundary
of the synchronized state.

$\displaystyle\epsilon_{HB}=\frac{2(2-\sqrt{(5+2c^{2}-c^{4})(q-1)^{2}}-2q)}{(1+c^{2})(q-1)^{2}},$
(12)

which is obtained by equating the trace of the Jacobian matrix of Eq. (7) to
zero $(tr(J_{Z}{{}_{{}_{3}}})=0)$, so that

$\displaystyle
tr(J_{Z}{{}_{{}_{3}}})=\frac{1}{(c^{4}-1)}[2\Delta+(c^{2}-1)^{2}\epsilon(q-1)].$
(13)

The Hopf bifurcation and saddle node curves are denoted by dotted and dashed
lines, respectively, in Fig. 6(a). The lower boundary, indicated by the solid
(green) line, separating the completely synchronized state and the
multistability state $R_{1}$ is obtained by solving Eq. (7) numerically. In
the region $R_{1}$, the system exhibits solitary state, where a small number
of individual oscillators excurs away from the synchronized group and
oscillates with quasi-periodic oscillations. The solitary state further gets
destabilized via $\epsilon_{HB}$ and gives birth to two-cluster (2C)
oscillatory state, which is denoted by the region $R_{2}$. In the 2C state,
the oscillators in the network split into two groups and oscillate with
different amplitudes. On the other hand the cluster state (see region $R_{2}$)
is stabilized via $\epsilon_{HB}$ which is further destabilized via
$\epsilon_{SN}$ (dashed line).

In addition, from the linear stability criterion for the solution $Z=1$ given
in Eq. (8) one can show that the synchronized state is stable in the entire
region of $(q,\epsilon)$ space. But it is known that the repulsive coupling
can destabilize the synchronized region at strong repulsive coupling strengths
in the parametric space. In order to find the analytical condition for the
synchronized manifold and to find the limit of stable solitary and two-cluster
states as a function of the repulsive coupling strength, we use the master
stability approach [57]-[59]. In the synchronized manifold, $x_{k}=x$, and
$y_{k}=y$, $\forall\,k=1,2...,N$, where $x(t)$ and $y(t)$ are the solutions
for the uncoupled system (1). For an infinitesimal perturbation to the
synchronized solution, Eq. (1) can be reduced as

$\displaystyle{\dot{\bm{\eta}_{1}}}=(1-3x^{2}-y^{2}-2cxy){\bm{\eta}_{1}}-(\omega+cx^{2}+3cy^{2}+2xy){\bm{\eta}_{2}}+\epsilon{\bf
G}{\bm{\eta}_{1}},$
$\displaystyle{\dot{\bm{\eta_{2}}}}=(\omega+3cx^{2}+cy^{2}-2xy){\bm{\eta}_{1}}+(1-x^{2}-3y^{2}+2cxy){\bm{\eta}_{2}}-q\epsilon{\bf
G}{\bm{\eta}_{2}},$ (14)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad{k=1,2,...N}$

Here $\bm{{\eta_{1}}}=(\eta_{1}^{x},\eta_{2}^{x},...,\eta_{N}^{x})^{T}$ and
$\bm{{\eta_{2}}}=(\eta_{1}^{y},\eta_{2}^{y},...,\eta_{N}^{y})^{T}$. We note
that $\eta_{j}^{x}$ and $\eta_{j}^{y}$ are the deviations of $x_{j}$ and
$y_{j}$ from the synchronized solution $(x,y)$. ${\bf G}$ is the coupling
matrix which is defined as

${\mathbf{G}}=\frac{1}{N}\left({\begin{array}[]{cccc}1-N&1&\cdots&1\\\
1&1-N&\cdots&1\\\ \vdots&\vdots&\ddots&\vdots\\\
{1}&{1}&\cdots&{1-N}\end{array}}\right)$

Then the perturbation from the synchronized manifold is defined as
${\beta_{jk}}$ $(j=1,2)$, where $\beta_{jk}={\bm{\eta_{j}^{T}V_{k}}}$.
$\bm{V_{k}}$ is the eigenvector of the coupling matrix G. The diagonalized
variational equation from Eq. (14) can be written as,

$\displaystyle{\dot{\beta}_{1k}}=(1-3x^{2}-y^{2}-2cxy)\beta_{1k}-(\omega+cx^{2}+3cy^{2}+2xy)\beta_{2k}+\epsilon\lambda_{k}\beta_{1k},$
$\displaystyle{\dot{\beta}_{2k}}=(\omega+3cx^{2}+cy^{2}-2xy)\beta_{1k}+(1-x^{2}-3y^{2}+2cxy)\beta_{2k}-q\epsilon\lambda_{k}\beta_{2k},$
(15)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad{k=1,2,...N}$

where $\lambda_{k}^{\prime}$s are the eigenvalues of the coupling matrix G.
The stable region of the synchronized state is obtained from the eigenvalues
of the coupling matrix G, which turn out to be $\lambda_{0}=0$,
$\lambda_{k}=-1$, $k=1,2,...N-1$. The eigenvalue $\lambda_{0}=0$ corresponds
to the perturbation parallel to the synchronization manifold, while the other
eigenvalues correspond to the perturbations transverse to the synchronization
manifold. The transverse eigenmodes should be damped out to have a stable
synchronization manifold. By substituting the eigenvalues in Eq. (15) and by
finding the largest Lyapunov exponents ($\Lambda_{max}$), one can obtain the
stable boundaries of the synchronized region which separates the boundary
between synchronized and desynchronized states. Whenever the largest Lyapunov
exponents acquires negative value the synchronized manifold is stable. The
line connecting the filled triangles in Fig. 6 (a) denote the boundary between
the synchronized and desynchronized states and the corresponding largest
Lyapunov exponents in $(q,\epsilon)$ space is shown in Fig. 6 (b). Finally, we
will deduce the analytical stability condition for oscillation death state in
the following.

###  4.2 Stability condition for oscillation death state

The system of coupled Stuart-Landau oscillators also exhibits oscillation
death state at strong coupling strengths, where the system suppresses the
oscillations and ultimately attains stable inhomogeneous steady states. Here,
we observed the total population of the network splits into either upper or
lower branches of inhomogeneous states. The corresponding dynamical equations
can be written as,

$\displaystyle\frac{dz_{d1}}{dt}=f(z_{d1})+\epsilon[pRe(z_{d1})+\bar{p}Re(z_{d2})-Re(z_{d1})]$
$\displaystyle\qquad\qquad\qquad\qquad\quad\quad\quad\quad-
iq\epsilon[pIm(z_{d1})+\bar{p}Im(z_{d2})-Im(z_{d1})],$ (16)
$\displaystyle\frac{dz_{d2}}{dt}=f(z_{d2})+\epsilon[pRe(z_{d1})+\bar{p}Re(z_{d2})-Re(z_{d2})]$
$\displaystyle\qquad\qquad\qquad\qquad\quad\quad\quad\quad-
iq\epsilon[pIm(z_{d1})+\bar{p}Im(z_{d2})-Im(z_{d1})],$ (17)

where $f(z_{d1})=(1+i\omega)z_{d1}-(1-ic)|z_{d1}|^{2}z_{d1}$ and
$f(z_{d2})=(1+i\omega)z_{d2}-(1-ic)|z_{d2}|^{2}z_{d2}$. The fixed point
solutions for $p=0.5$, of Eq. (16) and Eq. (17) can be expressed as

$\displaystyle x_{1,2}^{*}$ $\displaystyle=$
$\displaystyle\mp\frac{((1+q)\epsilon-
f_{0})}{\sqrt{2}c(\epsilon-1)-\omega}\sqrt{\frac{1}{(1+c^{2})^{2}(1+q)\epsilon}(f_{1}\epsilon^{2}+f_{2}\epsilon+f_{3}),}$
$\displaystyle y_{1,2}^{*}$ $\displaystyle=$
$\displaystyle\frac{c(\epsilon-1)-\omega}{((1+q)\epsilon-d)}\,\,x_{1,2}^{*}\,,$
(18)

where $f_{0}=\sqrt{h_{1}\epsilon^{2}+h_{2}(1-q)\epsilon-4h_{3}}$,
$f_{1}=q(1+q)+c^{2}(3q-1)$,
$f_{2}=(c^{2}+q)f_{0}+1+5c^{2}+q(1-3c^{2})+2c\omega(1-c^{2})-4cq\omega$,
$f_{3}=h_{4}f_{0}+2c^{3}\omega-4c^{2}-6c\omega-2\omega^{2}(1-c^{2})$,
$h_{1}=q^{2}+4c^{2}q+2q+1$, $h_{2}=4c(c+\omega)$,
$h_{3}=c^{2}+2c\omega+\omega^{2}$ and $h_{4}=1-c^{2}-2c\omega$. In the
oscillation death state, each oscillator in the network settles down in one of
the inhomogeneous steady states. Further using the linear stability analysis,
we can find the stability conditions for oscillation death state. The
inhomogeneous steady state becomes stable at

$\epsilon=\frac{1}{4q}(q-1)+\sqrt{1+q(2+4c^{2}+16\omega(c+\omega)+q)}.$ (19)

The stable boundary (line connected by unfilled circles) of OD region is shown
in Fig. 6. The stable boundary of OD state agrees well with the numerical
boundary. In addition, there is a coexistence of OD with SYN state in the
region $R_{3}$ (see Fig. 6). From the stable boundaries of the SYN, SS and 2C
states, the swing of synchronized state is strongly evident. The obtained
stable boundaries of the synchronized state, solitary and cluster states in
Fig. 6 are in good agreement with the numerical boundaries between the
dynamical states in Fig. 5.

##  5  Dynamical transitions for varying power-law exponent $(\alpha\geq 0)$

Previously, we have analysed the dynamical transitions for power-law exponent
$\alpha=0$. In this section, we study the dynamical behavior while varying the
power-law exponent $(\alpha)$ [60, 61], by considering the ring network of
coupled Stuart-Landau oscillators (1) for $N=101$.

Figure 7:  Standard deviation as a function of coupling strength $\epsilon$
(after leaving the transients time, $5\times 10^{4}$ units) for (a) $q=0.1$
(b) $q=0.3$ and (b) $q=0.7$. The open diamonds, circles, squars and triangular
points connecting by lines denote the power-law exponent
$\alpha=0.0,~{}\alpha=0.5,~{}\alpha=1.0$ and $\alpha=2.0$, respectively.

Figs. 7(a)-7(c) are plotted for three different values of the repulsive
coupling strengths $q=0.1$, $q=0.3$ and $q=0.7$, as a function of coupling
strength $\epsilon$. The different lines connected by points, namely diamonds,
circles, squares and triangles, respectively represent the power-law exponents
$\alpha=0.0,~{}\alpha=0.5,~{}\alpha=1.0$ and $\alpha=2.0$. Fig. 7(a) clearly
illustrates the emergence of the swing of synchronized state, for $\alpha=0$.
Increasing $\alpha$ to 0.5 and 1.0 give rise to a widening of the disorder
region among the synchronized region. Upon increasing $\alpha$ to still higher
values $(\alpha=2.0)$, the system exhibits completely desynchronized
(disordered) state for the entire range of coupling strength $\epsilon$.
Similar dynamical transitions are observed for $q=0.3$ and $q=0.7$ in Figs.
7(b) and 7(c). From Figs. 7, we have identified the fact that increasing the
power-law exponent can increase the disorderliness among the dynamical
regions.

Figure 8:  Standard deviation as a function $\alpha$ by fixing $\epsilon=1.0$
for different values of $q$ (after leaving the transients time, $5\times
10^{4}$ units). The open diamonds, circles, squares and triangular points
connected by lines denote the power-law exponents for
$q=0.1,~{}q=0.3,~{}q=0.5$ and $q=0.7$, respectively.

In order to show this clearly we have also plotted the dynamical transitions
as a function of power-law exponent by fixing the coupling strength at
$\epsilon=1.0$ for four different strengths of repulsive coupling $q$ (see
Fig. 8). The different lines connected by points such as diamonds, circles,
squares and triangles represent different values of the strength of repulsive
coupling, namely $q=0.1,~{}q=0.3,~{}q=0.5$ and $q=0.7$, respectively. For all
the values of $q$, the system is synchronized upto a critical value of
$\alpha$ (i.e. $\alpha\approx 1.2$), which disappears while increasing the
power-law exponent $\alpha$ to further higher values and the system becomes
completely desynchronized (disordered).

##  6 Conclusion

Long-range interacting systems are very rich and intriguing due to their
omnipresent nature in real world applications. Recently, many investigations
have been carried out for an in-depth understanding of long-range interacting
systems, including galaxies, quantum spin models and cold atom models. In this
article, we have investigated the collective dynamical behavior in a long
range interacting system of coupled Stuart-Landau limit cycle oscillators with
a combination of attractive and repulsive couplings as a prototypical example.
The dynamical behaviors were found as a function of the strength of the
repulsive coupling. Interestingly, we find the emergence of the swing of the
synchronized state, where the stable synchronized state is destabilized by the
emergence of solitary or cluster states which is again stabilized while
increasing the coupling strength. This re-emergence of synchronized state is
called the swing by mechanism. We have also found increasing disorder among
the coupled Stuart-Landau oscillators when increasing the strength of the
repulsive coupling. Finally, we have deduced the analytical boundary condition
for the synchronized state, solitary and cluster states and also for
oscillation death state. We have also found that the analytically obtained
boundaries match perfectly with the numerical boundaries. The dynamical
transition is also observed as a function of power-law exponent and we find
that the system exhibits completely desynchronized (disordered) state while
increasing the power-law exponent. The observed results may shed light for the
control of complete synchronization in long-range interacting systems.

## Acknowledgments

KS sincerely thanks the CSIR for a fellowship under SRF Scheme
(09/1095(0037)/18-EMR-I). The work of VKC forms part of a research project
SERB-DST Fast Track scheme for young scientists under Grant No.
YSS/2014/000175 and the CSIR Project under Grant No.03(14444)/18/EMR II. DVS
is supported by the CSIR EMR Grant No. 03(1400)/17/EMR-II. The work of ML is
supported by DST-SERB Distinguished Fellowship.

## References

## References

  * [1] De Monte S, d’Ovidio F and Mosekilde E 2003 Coherent regimes of globally coupled dynamical systems Phys. Rev. Lett. 90, 5
  * [2] Sailer X, Zaks M, and Schimansky-Geier L 2005 Collective dynamics in an ensemble of globally coupled FHN systems Fluctuation and Noise letters 5 L299-L304
  * [3] Moon J-Y, l Lee U, Blain-Moraes S, Mashour G A 2015 General relationship of global topology, local dynamics, and directionality in large-scale brain networks PLoS Comput. Biol. 11 e1004225
  * [4] Pikovsky A and Rosenblum M 2015 Dynamics of globally coupled oscillators: Progress and perspectives Chaos 25 097616
  * [5] Kuramoto Y 1981 Rhythms and turbulence in populations of chemical oscillators Physica A: Statistical Mechanics and its Applications 106 128-143.
  * [6] Kuramoto Y, Chemical Oscillations, Waves and Turbulence (Springer, Berlin, 1984).
  * [7] Chavanis P H and Rieutord M 2003 Statistical mechanics and phase diagrams of rotating self-gravitating fermions Astronomy and Astrophysics 412 1.
  * [8] Vega H J and Sánchez 2002 Statistical mechanics of the self-gravitating gas: I. Thermodynamic limit and phase diagrams Nuclear Physics B 625 407-518.
  * [9] Bouchet F and Sommeria J 2002 Emergence of intense jets and Jupiter’s great red spot as maximum-entropy structures Journal of Fluid Mechanics 464 165-207.
  * [10] Bouchet F and Barré J Statistical mechanics of systems with long range interactions J. phys. conf. series 31 18-26.
  * [11] Hadley P, Beasley M R, and Wiesenfeld K 1988 Phase locking of Josephson-junction series arrays Phys. Rev. B 38 8712
  * [12] Glass L 2001 Synchronization and rhythmic processes in physiology Nature (London) 410 277; Mosekilde E, Maistrenko Y, and Postnov D 2002 Chaotic Synchronization: Applications to Living Systems Chaotic Synchronization (World Scientific, Singapore, 2002);
  * [13] Lerch H P, Stange P, Mikhailov A S, and Hess B 2002 Mutual synchronization of molecular turnover cycles in allosteric enzymes III. Intramolecular Cooperativity J. Phys. Chem. B 106 3237
  * [14] Böhm F, Zakharova A, Schöll E, Lüdge K 2015 Amplitude-phase coupling drives chimera states in globally coupled laser networks Phys. Rev. E 91 040901
  * [15] Schmidt L and Krischer K 2015 Clustering as a prerequisite for chimera states in globally coupled systems Phys. Rev. Lett. 114 034101
  * [16] Schmidt L and Krischer K 2015 Chimeras in globally coupled oscillatory systems: From ensembles of oscillators to spatially continuous media. Chaos 25 064401
  * [17] Miritello G, Pluchino A, and Rapisarda A 2009 Phase transitions and chaos in long-range models of coupled oscillators Euro. Phys. Lett 85 10007
  * [18] Dauxois T, Ruffo S, Arimondo E and Wilkens M 2002 Dynamics and thermodynamics of systems with long range interactions, Lect. Notes Phys. (Springer), 602
  * [19] Gupta S, and Ruffo S 2017 The world of long-range interactions: A bird’s eye view, Int. J. Modern phys. A 32 1741018
  * [20] Bouchet F, Gupta S 2010 Mukamel D, Thermodynamics and dynamics of systems with long-range interactions Physica A: Statistical Mechanics and its Applications 389 4389-4405
  * [21] Nardini C, Gupta S, Ruffo S, Dauxois T, Bouchet F 2012 Kinetic theory for non-equilibrium stationary states in long-range interacting systems, Journal of Statistical Mechanics: Theory and Experiment 2012 L01002
  * [22] Campa A, Dauxois T, and Ruffo S 2009 Statistical mechanics and dynamics of solvable models with long-range interactions Physics Reports 480 57-159
  * [23] Chavanis P H 2006 Phase transitions in self-gravitating systems Int. J. of Modern Phys. B 20 3113-3198
  * [24] Elskens Y, Escande D F 2002 Microscopic Dynamics of Plasmas and Chaos, IOP Publishing, Bristol.
  * [25] Miller J 1990 Statistical mechanics of Euler equations in two dimensions Phys. Rev. Lett. 65 2137
  * [26] Maródi M, d’Ovidio F, and Vicsek T 2002 Synchronization of oscillators with long range interaction: Phase transition and anomalous finite size effects Phys. Rev. E 66 011109
  * [27] Viana R L, Batista A M, Batista C A S, de Pontes J C A, Silva F A S and Lopes S R 2012 Bursting synchronization in networks with long-range coupling mediated by a diffusing chemical substance Commun. Nonlinear Sci. Numer. Simulat., 17 2924-2942
  * [28] Pikovsky, A., Rosenblum, M., and Kurths, J. 2001 Synchronization: A Universal Concept in Nonlinear Sciences (The city: Cambridge University Presss)
  * [29] Hammond C, Bergman H, and Brown P 2007 Pathological synchronization in Parkinson’s disease: networks, models and treatments Trends Neurosci.  30 357-364
  * [30] Popovych O V and Tass P A 2014 Control of Abnormal Synchronization in Neurological Disorders Front Neurol. 5 268\.
  * [31] Zheng Y, Wang G,Li K, Bao G and Wang J 2014 Epileptic seizure prediction using phase synchronization based on bivariate empirical mode decomposition Clinical Neurophys. 125 1104-1111
  * [32] Sabesan S, Narayanan K, Prasad A, Spanias A, Sackellares J C and Iasemidis L D 2003  Biomedical sciences instrumentation 39 129-135
  * [33] Bhowmick S K, Ghosh D and Dana S K 2011 Synchronization in counter-rotating oscillators Chaos 21 033118
  * [34] Koluda P, Perlikowski P, Czolczynski K and Kapitaniak T 2014 Synchronization of two self-excited double pendula The European Physical Journal Special Topics 223 613-629
  * [35] Daido H and Nakanishi K 2006 Diffusion induced inhomogeneity in globally coupled oscillators: Swing-by mechanism Phys. Rev. Lett. 96 054101
  * [36] Daido H and Nakanishi K 2007 Aging and clustering in globally coupled oscillators textitPhys. Rev. E 75 056206
  * [37] Premalatha K, Chandrasekar V K, Senthilvelan M, and Lakshmanan M 2015 Impact of symmetry breaking in networks of globally coupled oscillators Phys. Rev. E 91 052915
  * [38] Sathiyadevi K, Karthiga S, Chandrasekar V K, Senthilkumar D V, Lakshmanan M 2017 Spontaneous symmetry breaking due to the trade-off between attractive and repulsive couplings Phys. Rev. E 95 042301
  * [39] Sathiyadevi K, Chandrasekar V K, Senthilkumar D V, Lakshmanan M 2018 Distinct collective states due to trade-off between attractive and repulsive couplings Phys. Rev. E 97 032207
  * [40] Mishra A, Hens C R, Bose M, Roy P K, and Dana S K 2015 Chimeralike states in a network of oscillators under attractive and repulsive global coupling Phys. Rev. E 92 062920
  * [41] Thompson M C and Gal P L 2004 The Stuart–Landau model applied to wake transition revisited European Journal of Mechanics B/Fluids  23 219-228
  * [42] Moon J-Y, Lee U, Blain-Moraes S and Mashour G A 2015 General relationship of global topology, local dynamics, and directionality in large-scale brain networks PLoS Comput. Biol. 11 e1004225
  * [43] Frasca M, Bergner A, Kurths J, and Fortuna L 2012 Bifurcation in a star-like network of Stuart-Landau oscillators. Int. J. of Bifur. and Chaos 22 1250173
  * [44] Punetha N, Varshney V, Sahoo S, Saxena G, Prasad A, Ramaswamy R 2018 Dynamical effects of breaking rotational symmetry in counter-rotating Stuart-Landau oscillators Phys. Rev. E 98 022212
  * [45] Zakharova A, Schneider I, Yuliya Kyrychko Y, and Schöll E 2013 Time delay control of symmetry-breaking primary and secondary oscillation death Euro. Phys. Lett 104(5) 50004
  * [46] Schneider I, Kapeller M, Loos S, Zakharova A, Fiedler B Schø”ll E 2015 Stable and transient multicluster oscillation death in nonlocally coupled networks Phys. Rev. E 92 052915
  * [47] Xiang Li 2006 Phase synchronization in complex networks with decayed long-range interactions Physica D: Nonlinear Phenomena 223 242-247
  * [48] Iubini S, Cintio P D, Lepri S, Livi R, and Casetti L 2017 Heat transport in oscillator chains with long-range interactions coupled to thermal reservoirs 2018 Phys. Rev. E 97 032102
  * [49] D Bagchi 2017 Thermal transport in the Fermi-Pasta-Ulam model with long-range interactions Phys. Rev. Lett. 95 032102
  * [50] Christodoulidi H, Tsallis C, and Bountis T 2014 Fermi-Pasta-Ulam model with long-range interactions: Dynamics and thermostatistics Euro. Phys. Lett 108 40006
  * [51] Maistrenko Y, Penkovsky B and Rosenblum M 2014 Solitary state at the edge of synchrony in ensembles with attractive and repulsive interactions Phys. Rev. E 69 060901(R)
  * [52] Jaros P, Brezetsky S, Levchenko R, Dudkowski D, Kapitaniak T, Maistrenko Y 2018 Solitary states for coupled oscillators with inertia Chaos 28 011103
  * [53] Majhi S, Muruganandam P, Ferreira F, Ghosh D, and Dana S K 2018 Asymmetry in initial cluster size favors symmetry in a network of oscillators Chaos 28 081101
  * [54] Gopal R, Chandrasekar V K, Venkatesan A, and Lakshmanan M 2014 Observation and characterization of chimera states in coupled dynamical systems with nonlocal coupling Phys. Rev. E 89 052914
  * [55] Sheeba J H, Chandrasekar V K, and Lakshmanan M 2011 General coupled-nonlinear-oscillator model for event-related (de) synchronization Phys. Rev. E 84 036210
  * [56] Chandrasekar V K, Gopal R, Venkatesan A, and Lakshmanan M 2014 Mechanism for intensity-induced chimera states in globally coupled oscillators Phys. Rev. E 90 062913
  * [57] Pecora L M and Carroll T L 1998 Master stability functions for synchronized coupled systems Phys. Rev. Lett. 80 2109
  * [58] Sun J, Bollt E M, and Nishikawa T 2009 Master stability functions for coupled nearly identical dynamical systems  Euro. Phys. Lett. 85 60011
  * [59] Acharyya S and Amritkar R E 2015 Synchronization of nearly identical dynamical systems: Size instability Phys. Rev. E 92 052902
  * [60] Banerjee T, Dutta P S, Zakharova A and Schöll E 2016 Chimera patterns induced by distance-dependent power-law coupling in ecological networks Phys. Rev. E 94 032206
  * [61] Gupta A, Banerjee T and Dutta P S 2017 Increased persistence via asynchrony in oscillating ecological populations with long-range interaction Phys. Rev. E 96 042202