††thanks: Corresponding author.

# Dynamical robustness of network of oscillators

Soumen Majhia Physics Department, University of Rome Tor Vergata, Via della
Ricerca Scientifica 1, 00133 Rome, Italy
Physics and Applied Mathematics Unit, Indian Statistical Institute, 203 B. T.
Road, Kolkata 700108, India  Biswambhar Rakshita Department of Mathematics,
Amrita School of Physical Sciences, Coimbatore, Amrita Vishwa Vidyapeetham,
India  Amit Sharmaa Department of Physics, University Institute of Sciences,
Chandigarh University, Mohali 140413, India  Jürgen Kurths  Potsdam Institute
for Climate Impact Research - Telegraphenberg A 31, Potsdam, 14473, Germany
Humboldt University Berlin, Department of Physics, Berlin, 12489, Germany
Dibakar Ghosh<EMAIL_ADDRESS>Physics and Applied Mathematics Unit,
Indian Statistical Institute, 203 B. T. Road, Kolkata 700108, India

###### Abstract

Most complex systems are nonlinear, relying on emergent behavior resulting
from many interacting subsystems, which are often characterized by oscillatory
dynamics. Having collective oscillatory behavior is an essential requirement
for an appropriate functioning of various real-world systems. Complex networks
have proven to be exceptionally efficient in elucidating the topological
structures of both natural and artificial systems, as well as describing
diverse processes taking place over them. Remarkable advancements have been
achieved in recent years in comprehending the emergent dynamics atop complex
networks. Specifically, among other processes, a large body of works intend to
explore the dynamical robustness of complex networks, which is the networks’
ability to withstand dynamical degradation in the network constituents while
maintaining collective oscillatory dynamics. Indeed, various physical and
biological systems are recognized to undergo a decline in their dynamic
activities, whether occurring naturally or influenced by environmental
factors. The impact of such damages on network performance can be significant,
and the system’s robustness is indicative of its capability to maintain
fundamental functionality in the face of dynamic deteriorations, often called
aging. This review offers a comprehensive excerpt of notable research
endeavors that scrutinize how networks sustain global oscillation under a
growing number of inactive dynamical units. We present the contemporary
research dedicated to the theoretical understanding and the enhancement
mechanisms of the dynamical robustness in complex networks. Our emphasis lies
on various network topologies and coupling functions, elucidating the
persistence of networked systems. We cover variants of system characteristics
from heterogeneity in network connectivity to heterogeneity in the dynamical
units. Finally we discuss challenges ahead in this potential field and open
areas for future studies.

a These authors equally contributed to the manuscript

###### pacs:

89.75.-k, 89.75.Hc, 05.45.Xt

Keywords: Complex networks, coupled oscillators, dynamical robustness, aging
transition

###### Contents

  1. I Introduction
  2. II Dynamical robustness for different network structures
    1. II.1 Dynamical robustness of globally coupled networks
      1. II.1.1 Diffusive coupling
      2. II.1.2 Time-delay diffusive coupling
      3. II.1.3 Mean-field diffusive coupling
      4. II.1.4 Weighted conjugate coupling
      5. II.1.5 Attractive and repulsive interaction
      6. II.1.6 Dynamical robustness in the absence of inactive oscillators
    2. II.2 Dynamical robustness of complex networks
      1. II.2.1 Homogeneously coupled complex networks
      2. II.2.2 Dynamical robustness of heterogeneous networks
      3. II.2.3 Weighted complex networks
      4. II.2.4 Correlated networks
    3. II.3 Dynamical robustness of multiplex networks
    4. II.4 Dynamical robustness of long-range connectivity networks
  3. III Dynamical robustness of quantum oscillators
  4. IV Dynamical robustness of biological networks
    1. IV.1 Spatial metapopulation networks
    2. IV.2 Neuronal networks
  5. V Enhancement and maintenance of dynamical robustness
    1. V.1 Controlled diffusion method
    2. V.2 Mean-field feedback method
    3. V.3 Addition of active oscillators
    4. V.4 Self-feedback delay
    5. V.5 Asymmetric interactions
    6. V.6 Low-pass filtering mechanism
  6. VI Machine learning techniques to predict aging transition
  7. VII Conclusions and future perspectives

##  I Introduction

Many natural systems are nonlinear and rarely isolated, and hence
understanding of such complex systems requires system-level interpretation.
The complexity of many social, biological, and physical systems emanates from
the intricacy in the patterns of interaction among their constituents. Within
this paradigm, the field of network science has been developed as the ideal
platform for providing tools for modeling and analysis of complex systems
barabasi2016network ; strogatz2001exploring ; newman2003structure ;
boccaletti2006complex ; albert2002statistical ; watts1998collective ;
barabasi1999emergence . These networked systems often function through the
emergent behavior of many interacting units, where each unit exhibits
oscillatory dynamics. Consequently, interacting oscillatory systems constitute
an efficient framework to model many complex systems. In such a framework, the
inherent dynamics of an individual node can be modeled as a system of
nonlinear differential equations, and different coupling functions can
describe the interactions among the nodes stankovski2017coupling . Even simple
nonlinear systems, when connected with each other through rather simple
coupling functions, can generate complex collective dynamics. With various
types of coupling schemes and network structures, a wide variety of collective
dynamics have been explored. For example, synchronization processes are
ubiquitous in nature and has been studied extensively in populations of
locally interacting elements in the context of physical, biological, chemical,
and technological systems boccaletti2002synchronization ;
pecora1990synchronization ; mirollo1990synchronization ;
arenas2008synchronization ; ghosh2022synchronized . Among various types of
partial synchronization, the coexistence of coherent and incoherent patterns,
commonly referred to as the chimera state panaggio2015chimera ;
majhi2019chimera ; abrams2004chimera ; bera2017chimera ;
montbrio2004synchronization , has been in the focus over the last decade
across a diverse array of systems. Coupled oscillators, therefore, play a key
role in many areas of science and technology, and their dynamics help us to
explore various emergent behaviors of living and non-living systems.

Among all the research perspectives concerning dynamics of networks, studies
focused on network robustness, which refers to the ability to withstand even
strong perturbations, holds considerable significance from various aspects.
This phenomenon of network robustness can be conceptualized in two distinct
ways: structural and dynamical robustness. i) Structural robustness addresses
the endurance of network activities when faced with structural perturbations,
which could involve the removal of links (bond percolation) or nodes (site
percolation) within a network cohen2010complex ; albert2000error ;
callaway2000network ; cohen2000resilience ; dorogovtsev2008critical ;
buldyrev2010catastrophic ; vespignani2010fragility ; gao2016universal . ii) On
the other hand, dynamical robustness of complex networks, in general, is
defined as the network’s capacity to sustain its dynamic activity despite
local disturbances. Specifically, throughout this review, we uncover dynamical
robustness of networks of coupled oscillators which refers to the network’s
ability to sustain its global dynamical activity even when a portion of its
dynamical components are functionally degraded. In the literature, this is,
popularly presented in terms of aging transition. Aging transition is an
emergent collective phenomenon in networked systems comprising of self-
oscillatory and non self-oscillatory nodes. This transition occurs when the
network shifts from a globally oscillatory state to a oscillation quenched
state as the proportion of inactive (aged) oscillators exceeds a critical
threshold.

Possessing collective oscillatory behavior is an essential prerequisite for
the regular functioning of many complex systems. Examples include circadian
rhythms bernard2007synchronization , biological pacemaker cells
miake2002biological , cardiac and respiratory systems Jalife_1998_Chaos  etc.
Oscillation plays important roles in other several dynamic processes within
both single cell and multicellular processes kruse2005oscillations ;
chen2010modeling . Additionally, coupled oscillator models are applicable to
electric power-grid networks, where components such as power sources need to
be synchronized to the same frequency Menck_2014_Nat.Com ;
filatrella2008analysis . Robust oscillatory dynamics is thus a fundamental
characteristic of such systems winfree1980geometry ; strogatz2001exploring .
Despite being regularly subjected to internal and external disturbances, these
systems can maintain their rhythmic activities to a certain degree. If a
limited number of units below a specific threshold fail to generate
oscillatory behavior, the remaining units can compensate, enabling the entire
system to resiliently preserve its proper functioning. However, if a
substantial number of units transit to an inactive quenched state, it can
significantly impede their functions, potentially resulting in a partial
collapse or even complete failure of the system in question.

Keeping this fundamental and inevitable context in mind, in the last two
decades, researchers from the nonlinear dynamics community have been working
on the dynamical robustness of a network of coupled oscillators Daido_2004_PRL
; Daido_2007_PRE ; tanaka2012dynamical . It is defined as the ability of a
network to sustain its collective macroscopic oscillation when a few of its
nodes fail to produce rhythmic dynamics due to local degradation
tanaka2012dynamical . Daido and Nakanishi Daido_2004_PRL  in their pioneering
work laid the mathematical framework for studying dynamical robustness. They
examined a situation in which oscillating nodes gradually transform into
stable points. If the quantity of nodes transitioning to stable points exceeds
a critical threshold, it could disrupt the typical oscillatory patterns in
these systems, leading to an abrupt phase shift toward a universally non-
oscillatory state. They explored a scenario, where oscillatory nodes
progressively transition to fixed points. If the number of nodes that shift to
fixed points surpasses a critical threshold, the usual oscillatory behavior of
these systems could be disrupted, resulting in a sudden phase transition
towards a globally non-oscillatory state. This phenomenon, characterized by
such abrupt and catastrophic emergence, is termed an aging transition. The
authors demonstrated that in a global network of diffusively coupled
oscillators, an aging transition can be characterized by a universal scaling
law of an order parameter involving inactivation fraction and the strength of
coupling. After their initial work a series of research articles have been
published. Owing to its widespread relevance, aging transition has been
studied in diverse models with different coupling functions and network
structures. Among the many significant attempts made along this topic, Pazó et
al. pazo2006universal  studied aging transition in an ensemble of globally
coupled Morris-Lecar model which exhibits a saddle-node bifurcation on an
invariant circle. A similar scaling law is established for an ensemble of
excitable and oscillatory units. Tanaka et al. tanaka2012dynamical  have
explored the aging transition in a complex network and have shown that scale-
free networks are highly resilient to random inactivation but extremely
vulnerable to targeted inactivation of low-degree oscillators with respect to
the dynamical robustness. Their finding is not in agreement with the
structural robustness of a scale-free network where high degree nodes play a
crucial role. Dynamical robustness has been studied in the context of
metapopulation dynamics by Kundu et al. Rakshit_2017_PRE . Their results
reveal how the network topology plays a crucial role in metapopulation
survivavbility. Thakur et al. thakur2014time  studied the influence of time-
delayed coupling on the nature of the aging transition in globally coupled
Stuart-Landau oscillators. Their findings divulge that time delay in the
coupling does not favor dynamical robustness. In Ref. morino2011robustness
the authors demonstrated aging transition in a multi-layer network of active
and inactive units, while contemplating with various interlayer coupling
functions.

As mentioned above, such global dynamical degradation in the form of aging
transition can have outright impacts on the substantive system performance.
Thus for many of real-world and man-made systems, continued dynamic
oscillatory activity of the components are often extremely crucial for
maintaining proper functioning. For instance, in case of neuronal activity
lisman2008neural , from physiological processes like cell necrosis within
organs gurtner2007progress  to cardiac and respiratory systems jalife1998self
, robust oscillatory dynamics is quite necessary. Due to such high practical
importance many researchers proposed remedial measures to enhance the
dynamical robustness against aging or deterioration of individual units. Among
the notable attempts in this regard, Liu et al. liu2016enhancing  proposed an
efficient method to enhance dynamical persistence by introducing an additional
parameter that controls the diffusion rate. In Ref. kundu2018resumption
authors established that by adding a linear positive mean-field feedback term,
network’s dynamical robustness can be improved substantially. Effectiveness of
self-feedback delay to increase the dynamical resilience has been studied by
Sharma and Rakshit sharma2021enhancement . The increasing body of vast
literature on aging transition and dynamical robustness is itself a testimony
of its relevance in various fields of science and engineering. In this report,
we intend to provide an exhaustive overview on aging transition by integrating
prevailing knowledge achieved in the last two decades. Thus the relevant
results and methodology related to dynamical robustness will be more generally
accessible for researchers in diverse communities of science and technology.
There are also several open challenging problems emphasized here.

##  II Dynamical robustness for different network structures

Based on the network topologies, we discuss the dynamical robustness by taking
different types of networks. Most of the results are explored for globally
connected network and complex networks. We also discuss the results on
dynamical robustness for multiplex, time-varying and long-ranged networks.

###  II.1 Dynamical robustness of globally coupled networks

Initially, we address the outcomes regarding dynamical robustness concerning
instantaneous diffusive coupling, followed by an exploration of time-delay
diffusive coupling topology. Subsequently, we introduce findings related to
weighted conjugate coupling and interactions characterized as attractive-
repulsive. Finally, we delve into dynamical robustness in scenarios where
inactive oscillators are absent.

####  II.1.1 Diffusive coupling

Figure 1: The normalized order parameter $R$ (Eq. (3)) is plotted as a
function of $p$ to illustrate the aging transition in the globally coupled
Stuart–Landau limit-cycle system (Eq. (1)). Different (increasing) values of
the coupling strength $K$ are chosen ranging from $K=1$ to $K=8$, with $a=2$
and $b=1$. Reprinted figure with permission from Ref.Daido_2004_PRL

Diffusive coupling represents the predominant form of coupling observed in
numerous real-world systemsrosenblum2003synchronization , consequently
garnering significant attention when investigating aging transitions. We first
discuss the aging transition phenomena in a system of $N$ globally coupled
Stuart–Landau oscillators Daido_2004_PRL . The mathematical form of all-to-all
diffusively coupled network is,

$\dot{z}_{j}(t)=(\alpha_{j}+i\omega-|z_{j}(t)|^{2})z_{j}(t)+\frac{K}{N}\sum_{l=1}^{N}(z_{l}-z_{j}),$
(1)

where $j=1,2,...,N$ and $z=x+iy$ is a complex variable. Here $\omega$
signifies the internal frequency of each oscillator, and $\alpha_{j}$ is the
bifurcation parameter of the $j$-th oscillator, denoting its proximity to a
Hopf bifurcation point. When $\alpha_{j}>0$, each isolated Stuart-Landau
oscillator exhibits a stable sinusoidal oscillation, while it converges to the
stable trivial fixed point $z_{j}=0$ for $\alpha_{j}<0$. The second term on
the right-hand side of Eq.(1) denotes the presence of diffusive coupling with
a strength represented by $K$.

The framework introduced by Daido and Nakanishi Daido_2004_PRL  facilitates
the study of robustness in a manner where an active oscillator with
$\alpha_{j}=a>0$ transitions to an inactive state with $\alpha_{j}=-b<0$, with
a parameter $p$ representing the fraction of inactive oscillators. To
simplify, we can designate the group of oscillators as active ones, denoted by
$j=1,2,...,N(1-p)$, and the remaining oscillators as inactive, denoted by
$j=N(1-p)+1,N(1-p)+2,...,N$. This implies that $0<p<1$ signifies the fraction
of inactive oscillators in the entire networked system. The degree of
macroscopic oscillation of the whole network is then measured by the order
parameter $|\bar{Z}(p)|$,

$\bar{Z}(p)=\frac{1}{N}\sum_{j=1}^{N}z_{j},$ (2)

and subsequently by the normalized order parameter $R$, which is defined as
follows:

$R=\frac{|\bar{Z}(p)|}{|\bar{Z}(0)|}.$ (3)

As one increases the inactivation parameter $p$, the order parameter gradually
diminishes and at a critical value $p=p_{c}$, the loss of global oscillation
takes place and leads to a transition in the mean-field dynamics called aging
transition. Figure 1 portrays the normalized order parameter $R$ as a function
of $p$. It depicts the gradual loss of global oscillation for various
strengths of coupling constant $K$. For decreasing $K$ the critical value of
the inactivation parameter $p_{c}$ monotonically increases until it reaches
unity at a threshold value of $K=K_{c}$, below which $p_{c}$ remains at unity.

Synchronized activities among the oscillators permit us to assume that within
each group of active and inactive nodes oscillators behave identically. With
this presumption in mind, we assign $z_{j}$ as $A$ for the active ensemble and
$z_{j}$ as $I$ for the inactive ensemble of oscillators. Consequently, Eq. (1)
transforms into the subsequent interconnected systems:

$\displaystyle\dot{A}(t)$ $\displaystyle=$
$\displaystyle(a-Kp+i\omega-|A(t)|^{2})A(t)+KpI(t),$ $\displaystyle\dot{I}(t)$
$\displaystyle=$ $\displaystyle(-b-Kq+i\omega-|I(t)|^{2})I(t)+KqA(t),$

where $q=1-p.$ A linear stability analysis of the reduced system leads to an
analytical formula of the critical point $p_{c}$ as

$\displaystyle p_{c}=\frac{a(K+b)}{(a+b)K}.$ (5)

In the limiting case, $\lim_{K\to\infty}p_{c}=\frac{a}{a+b}$. Consequently,
one can conclude that the dynamical robustness is stronger when the active
oscillators have larger amplitude of oscillation. We can derive the scaling
property of the order parameter near the critical point $p_{c}$ Daido_2004_PRL
. The scaling law is represented as $|\bar{Z}|\propto(p_{c}-p)^{\beta}$, where
the critical exponent $\beta$ varies based on the coupling strength in the
following manner,

$\displaystyle\beta=\begin{cases}1/2&\text{for $K<K_{c}$}\\\ 1&\text{for
$K=K_{c}$}\\\ 3/2&\text{for $K>K_{c}$}.\end{cases}$

The critical exponent $\beta$ increases by increasing the coupling strength
$K$.

Until recently, the investigation of dynamical robustness was primarily
focused on coupled Stuart–Landau oscillators, which exhibit typical sinusoidal
oscillations. However, many natural systems can be modeled by networks of non-
sinusoidal oscillators, such as the Van der Pol oscillator. Rakshit et al.
rakshit2020abnormal  have studied the aging transition in a network of
globally coupled Van der Pol oscillators and highlighted how it differs from
the Stuart-Landau model. Their investigation uncovers distinct pathways to
aging transition in networks of Van der Pol oscillators compared to typical
sinusoidal oscillators like Stuart–Landau oscillators. Unlike sinusoidal
oscillators, where the order parameter smoothly undergoes a second-order phase
transition, they observed an unconventional phase transition characterized by
the abrupt emergence of unbounded trajectories at a critical point. Through
detailed bifurcation analysis, they elucidated this abnormal phase transition,
demonstrating that it is driven by the boundary crisis of a limit-cycle
oscillator, paving the way for an unconventional and discontinuous aging
transition path.

####  II.1.2 Time-delay diffusive coupling

Dynamical systems characterized by temporal delays, are prevalent in nature.
They manifest due to finite signal propagation time and memory effects in a
diverse range of natural phenomena sismondo1990synchronous ; yeung1999time ;
boutle2007nino ; dhamala2004enhancement ; lakshmanan2011dynamics ;
reddy1998time , spanning physical, chemical, engineering, economic, and
biological domains, including their respective networks. Several dynamical
systems can be portrayed by delay differential equations with single constant
delay mackey1977oscillation ; ikeda1980optical ; fort1999time , discrete
delays saad1998comparative ; asea1999time ; kuang1993delay , distributed delay
haldane1934contribution ; baylor1974electrical , state-dependent delay
insperger2008criticality ; hartung2006functional , and time-dependent delay
kye2004synchronization ; senthilkumar2007delay ; ghosh2010generalized ;
majhi2016synchronization ; ghosh2010projective .

Let us now introduce a time delay in the signal transmission and observe how
it affects the dynamical robustness of networked systems. For simplicity, we
contemplate with globally (all-to-all) coupled systems of oscillators. The
system of $N$ all-to-all coupled Stuart-Landau oscillators subject to linear
time-delayed coupling can then be represented by the following set of
equations,

$\begin{array}[]{lcl}\dot{z_{j}}=(\alpha_{j}+i\omega-|{z_{j}}|^{2}){z_{j}}+\dfrac{\kappa^{\prime}}{N}\sum\limits_{k=1,k\neq
j}^{N}[{z_{k}}(t-\tau)-{z_{j}}(t)];\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ j=1,2,\cdots,N,\end{array}$ (6)

with $\kappa^{\prime}=2\kappa$ as the interaction strength and $\tau$ as the
parameter describing the time delay in the signal transmission.

Choosing a network of $N=500$ dynamical units, we depict in Fig. 2(a) how the
normalized order parameter $R$ (Eq. (3)) changes with respect to increasing
inactivation represented by the parameter $p$, for various values of $\tau$.
We keep $\kappa^{\prime}=5$ fixed and start with the non-delay case (i.e.,
$\tau=0$). Then we choose a higher time-delay with $\tau=0.01$ and observe a
slightly faster aging transition in the system. More significant change takes
place when we increase $\tau=0.03$ where the critical inactivation ratio
$p_{c}$ evidently decreases. This scenario remains valid as we increase the
delay to $\tau=0.05$ and $\tau=0.07$. In fact, with increasing $\tau$ this
scenario becomes more pronounced and the aging transition occurs faster. So,
for a fixed interaction strength, the critical inactivation ratio $p_{c}$
decreases with increasing $\tau$. Thus, the aging takes place faster and hence
the robustness of the networked system decreases due to the introduction of
time-delay.

Going further, one can perform a linear stability analysis of this reduced
system of equations around the origin (as in Ref. thakur2014time ), ending up
with a characteristic equation for the eigenvalues, from which the aging
islands in the $(\kappa,\leavevmode\nobreak\ \tau)$ parameter plane can be
determined. Subsequently, for a better perception of the robustness subject to
time-delay, along with varying the time-delay, we consider simultaneous change
in the interaction strength, and plot the phase diagram in the
$(\kappa,\leavevmode\nobreak\ \tau)$ parameter plane for different values of
the inactivation ratio $p$ (in Fig. 2(b)). Specifically, $p=0.0,0.2,$ and
$0.4$ are chosen and it is clear that the aging island expands for raising the
ratio of the inactive oscillators. The aging islands expand in both directions
of $\kappa$ and $\tau$ reflecting the fact that along with the interaction
strength and the inactive elements, the time-delay also carries the capability
of suppressing oscillations to the trivial fixed point and hence to decrease
the robustness of the system.

Figure 2: (a) The normalized order parameter $R$ for the system (6) as a
function of the proportion $p$ of inactive oscillators for different values of
the time-delay $\tau$ for $N=500$ dynamical units, whenever $\omega=10$ and
$\kappa^{\prime}=5$ with $a=2,b=1$. (b) Aging islands in the $(\kappa,\tau)$
parameter plane for various values of $p$. Reprinted figure with permission
from Ref.thakur2014time

Further, the study by Rahman et al. rahman2017aging  examines a globally
connected network comprising active and inactive oscillators with distributed-
delay coupling. It establishes conditions for aging transition, derived for
both uniform and gamma delay distributions. The findings suggest that in the
case of a uniform distribution, increasing the width of the delay
distribution, while maintaining the same mean delay, enables an aging
transition to occur with a smaller coupling strength and a lower proportion of
inactive elements. Provided the coupling strength falls within a specific
range and the mean time delay is substantial, it may be feasible to achieve an
aging transition for any proportion of inactive oscillators for a gamma
distribution.

####  II.1.3 Mean-field diffusive coupling

The study of mean-field coupling is extensively studied due to its prevalence
in numerous natural occurrences within the realms of physics, biology and
engineering. Effects of mean-field diffusion have previously been explored in
synchronization shiino1989synchronization , multistable dynamics of synthetic
genetic networks garcia2004modeling ; ullner2007multistability ;
koseska2007inherent ; ullner2008multistability , oscillation suppression
processes mirollo1990amplitude ; de2003coherent ; sharma2012amplitude ;
ghosh2015revival , and also in the dynamics of chimera-death banerjee2015mean
. Recently, the dynamical robustness in the coupled Stuart-Landau oscillators
through mean-field diffusion has been reported in bera2020 . The mathematical
model of the coupled oscillators is given by

$\displaystyle\dot{z}_{j}$ $\displaystyle=$
$\displaystyle(\alpha_{j}+i\omega-|z_{j}|^{2})z_{j}+\epsilon(Q\bar{z}-z_{j}),$
(7)

where $\bar{z}$ is mean value which represent by
$\frac{1}{N}\sum_{l=1}^{N}z_{l}$. In this coupling scheme, the control
parameter $Q$ associated with the mean-field interaction describes the influx
of the mean- field in the dynamical units. The parameter $Q$ essentially
controls the rate of mean-field interaction in the diffusive coupling among
the interacting oscillatory systems. When $Q=0$, there is no interaction
between the oscillators and they behave like uncoupled ones subjected to self-
feedback, while $Q=1$ maximizes the interaction with the mean field
PRE_2012_Amit . Linear stability analysis around the origin enables us to
analytically derive the critical value $p_{c}$ for Eq. (7) as bera2020 ,

$\displaystyle
p_{c}=(b+\epsilon)\left(\frac{1}{a+b}+\frac{a-\epsilon}{Q\epsilon(a+b)}\right),\leavevmode\nobreak\
\leavevmode\nobreak\ with\leavevmode\nobreak\ \leavevmode\nobreak\
Q>1-\frac{a}{\epsilon},$ (8)

where $b-a-Q\epsilon+2\epsilon>0$. The above expression for $p_{c}$ clearly
indicates that lowering the mean-field parameter $Q$ has a negative effect on
the dynamical robustness of the networked system. This result is in agreement
with the previous study that the mean-field control parameter plays an
important role in the suppression of oscillations PRE_2012_Amit . As
$Q\rightarrow 0$, the effect of the mean-field interaction in coupling
decreases, which in turn hinders the collective macroscopic oscillation.
Figure 3 depicts the normalized order parameter $R$ (Eq. (3)) as a function of
the inactivation parameter $p$ for various values of $Q$. It is discernible
from the plot for the maximum mean-field scenario (with $Q=1$) that the order
parameter sharply decreases with increasing $p$, and eventually vanishes at
the critical value $p_{c}\sim 0.75$, indicating the occurrence of an aging
transition. The value of $Q$ is decreased subsequently and the order parameter
is plotted for $Q=0.9$, for which aging transition takes place at a lower
value of $p$, and hence the dynamical robustness of the coupled system
reduces. Similarly for a lower value of $Q=0.8$, the robustness decreases even
more. This trend remains intact for even lower values of the mean-field
parameter $Q=0.7$ and $Q=0.6$. These numerical results are in line with the
analytically obtained critical values of the inactivation ratio as well. Thus,
these plots demonstrate that a lower mean-field density forces the entire
system to dynamically collapse for lower values of the critical inactivation
ratio, and hence leads to the decrement of dynamical robustness. Similar
results are observed for networks of delay-coupled systems as well bera2020 .

Figure 3: The influences of the mean-field density parameter for the system
(7) on the aging transition scenarios has shown for a certain coupling value.
The order parameter $R$ versus the fraction of inactivation ratio p
corresponding to the various mean-field density parameter $Q=0.6$, $0.7$,
$0.8$, $0.9$ and $Q=1.0$ with fixed coupling strength $\epsilon=2.0.$
Reprinted figure with permission from Ref.bera2020 .

####  II.1.4 Weighted conjugate coupling

Apart from the traditional coupling via the similar variables among the
dynamical systems, coupling through dissimilar or conjugate variables is also
natural in a number of experimental scenarios in which the units are coupled
by feeding the output of one into the other. Examples include the experiments
by Kim et al. kim2005scaling  of coupled semiconductor laser systems. Later,
in studying the phenomenon of oscillation quenching saxena2012amplitude ;
koseska2013oscillation , this form of coupling is utilized in detail
karnatak2007amplitude ; punetha2018dynamical , as coupling through dissimilar
variables naturally breaks the rotational symmetry of diffusively coupled
systems. Conjugate coupling has further been shown to be effective in
enhancing coherence zhao2018enhancing  and inducing explosive death
zhao2018explosive .

Here we consider a globally coupled network of Stuart-Landau limit cycle
oscillators with weighted conjugate coupling Ponrasu2019 , for which the
governing equation can be represented as,

$\displaystyle\dot{z}_{j}$ $\displaystyle=$
$\displaystyle(\alpha_{j}+i\omega-|z_{j}|^{2})z_{j}+\frac{\epsilon
m_{(1,2)}}{N}\sum_{k=1}^{N}[(Img(z_{k})-\beta Re(z_{j}))+i(Re(z_{k})-\beta
Img(z_{j}))],$ (9)

where $\epsilon$ is the overall coupling strength, and $\beta$ is the feedback
control parameter ($0\leq\beta\leq 1$). Here $Re(z)$ and $Img(z)$ are the real
and imaginary parts of $z$, respectively. $\alpha_{j}$ is the bifurcation
parameter for Stuart-Landau oscillators, specifically we choose $\alpha_{j}=2$
for the active set of oscillators, and $\alpha_{j}=-1$ for the inactive ones.
The intra-group coupling strengths of active and inactive oscillators are
$m_{1}$ and $m_{2}$. When $m_{1}=m_{2}$, the interaction is symmetric, while
for $m_{1}\neq m_{2}$ the interaction becomes asymmetric.

For the stability analysis, the global network is divided into two groups of
active (A) and inactive (I) oscillators. The $A=A_{r}+iA_{im}$, and
$I=I_{r}+iI_{im}$, where $A_{r},I_{r}$ and $A_{im},I_{im}$ are the real and
imaginary variables which satisfy the following equations,

$\displaystyle\dot{A}_{r}$ $\displaystyle=$ $\displaystyle
aA_{r}-A_{im}\omega-
A_{r}A^{2}_{im}-A^{3}_{r}+m_{1}\epsilon[(1-p)A_{im}+pI_{im}-\beta A_{r}],$
$\displaystyle\dot{A}_{im}$ $\displaystyle=$ $\displaystyle
aA_{im}+A_{r}\omega-
A_{im}A^{2}_{r}-A^{3}_{im}+m_{1}\epsilon[(1-p)A_{r}+pI_{r}-\beta A_{im}],$
$\displaystyle\dot{I}_{r}$ $\displaystyle=$ $\displaystyle
bI_{r}-I_{im}\omega-
I_{r}I^{2}_{im}-I^{3}_{r}+m_{2}\epsilon[(1-p)A_{im}+pI_{im}-\beta I_{r}],$
$\displaystyle\dot{I}_{im}$ $\displaystyle=$ $\displaystyle
bI_{im}+I_{r}\omega-
I_{im}I^{2}_{r}-I^{3}_{im}+m_{2}\epsilon[(1-p)A_{r}+pI_{r}-\beta I_{im}].$
(10)

Figure 4: (a-d) The order parameter $R$ for the system (9) varies with the
inactivation ratio $p$ across various coupling strengths, and (e-h) the phase
diagram in ($\epsilon-p$) space for the natural frequency $w=1$ for
$(m_{1},m_{2})=(1,1),(1,2),(2,1)$ and $(1,3)$ (left to right). Dynamics
regions: HSS for homogenous steady state corresponding to the aging transition
(black), IHSS for inhomogeneous steady state (orange), and OS for oscillatory
state (yellow).

The determination of the critical value $p_{c}$ involves conducting a linear
stability analysis on Eq. (II.1.4) at the origin, denoted as $(\bf A,\bf
I)=(\bf 0,\bf 0)$. For this purpose, the 4th order characteristic equation is
deduced from the Jacobian, and the Hopf and pitchfork bifurcation points using
the Routh-Hurwitz stability criterion are computed (for more details see
Ponrasu2019 ).

For the numerical analysis, first we consider the symmetric case
($m_{1}=m_{2}=1$) in the conjugate coupled system of Stuart-Landau oscillators
with the natural frequency $\omega=1$ and the control parameter $\beta=1$. In
Figure 4(a), we illustrate the variation of the order parameter $R$ as a
function of the inactivation ratio $p$ across various coupling strengths
$\epsilon$. The critical value $p_{c}$ of aging transition to the homogeneous
steady state (HSS) decreases with increasing the coupling strength $\epsilon$.
To show the effect of $\epsilon$ on the aging transition, we portray the
corresponding phase diagram in the $p-\epsilon$ parameter plane in Fig. 4(e),
which confirms the fall of $p_{c}$ with coupling strength for HSS. To untangle
the impact of asymmetry parameter on the aging transition, we assume three
cases with the parameter sets $(m_{1},m_{2})=(1,2),(2,1)$ and $(1,3)$. The
aging transition in terms of $R$ and the corresponding phase diagrams for
these three cases are respectively depicted in Figs. 4(b-d) and 4(f-h). As the
impact of inactive (active) oscillators on the active (inactive) oscillators
intensifies with rising values of $m_{2}(m_{1})$, it either enhance or shrink
the HSS region in the parameter plane.

Figure 5: Phase diagrams for the system (9) in the $(\epsilon,\omega)$ space
for $p=0.8$ for conjugate coupled Stuart-Landau oscillators for coupling
strength for intra groups (a) $m_{1}=m_{2}=1$, (b) $m_{1}=1,m_{2}=2$,
(c)$m_{1}=2,m_{2}=1$ and (d) $m_{1}=1,m_{2}=3$.

The impact of the natural frequency on the aging transition is shown in the
($\epsilon-\omega$) parameter plane for different combinations of
($m_{1},m_{2}$) at fixed inactivation ratio $p=0.8$ in Fig. 5(a-d). In the
parameter plane, dynamical transitions occur through three bifurcations
(similar to Fig. 4) among the oscillatory regime (OS), HSS and the
inhomogeneous steady state (IHSS). The variation in the spread of the aging
transition island as a function of $\omega$ and $\epsilon$ for different
combinations of ($m_{1},m_{2}$) is depicted in Fig. 5. Increasing the coupling
strength of inactive oscillators through $m_{2}$ enlarges the aging transition
region (HSS), as a result of the influence of the inactive oscillators over
the active group of oscillators. On the other hand, with the increase of
$m_{1}$ the aging transition region, i.e., the HSS region shrinks in the
parameter plane due to the influence of active oscillators over the inactive
ones in the globally conjugate coupled Stuart-Landau oscillators in Figs.
5(a-d).

####  II.1.5 Attractive and repulsive interaction

Up to this point, we have explored the concept of aging transition,
predominantly focusing on scenarios where interactions between dynamical nodes
are attractive (positive). However, real-world systems often exhibit a more
intricate nature, involving mixed coupling with both attractive (positive) and
repulsive (negative) connections majhi2020perspective . It is important to
emphasize that the coexistence of positive and negative couplings can
simultaneously occur in diverse physical, ecological, and biological systems
giron2016synchronization ; daido1987population ; dixit2020static . In recent
years, numerous studies have delved into the emergent dynamics within an
ensemble of oscillators featuring attractive-repulsive coupling. These
investigations have showcased the emergence of manifold collective dynamics,
including chimera states, solitary states, extreme events, amplitude (or
oscillation) death, and anti-phase synchrony in a network of oscillatory nodes
experiencing both attractive and repulsive interactions majhi2020perspective ;
hens2013oscillation ; dixit2020static ; sharma2021explosive ;
chowdhury2021antiphase ; chowdhury2020effect . The interplay between
attractive and repulsive couplings can lead to suppression of oscillatory
activities among coupled oscillators. As a result, the phenomenon of aging
transition exhibits qualitative differences in the presence of competitive
attractive-repulsive interactions. In this direction the very first work was
carried out by Bera bera2020additional . The mathematical form of a $N$
coupled network of Stuart-Landau oscillators having attracting and repulsive
interactions is,

$\displaystyle\dot{z}_{j}$ $\displaystyle=$
$\displaystyle(\alpha_{j}+i\omega-|z_{j}|^{2})z_{j}+\frac{K}{N}\sum_{k=1}^{N}C_{jk}G(z_{k},z_{j})-\frac{\epsilon}{N}\sum_{k=1}^{N}B_{jk}H(z_{k},z_{j}),$
(11)

where $j=1,2,\dots,N$. Here $C_{jk}$ and $B_{jk}$ represent the adjacency
matrices of attractive and mean repulsive interactions, respectively. The
functional forms of attractive and repulsive interactions are defined by
$G(z_{k},z_{j})=(z_{k}-z_{j})$ and $H(z_{k},z_{j})=(z_{k}+z_{j})$. In this
context, the parameters $K$ and $\epsilon$ represent the strengths of
attractive and repulsive interactions, respectively.

To determine the critical threshold $p_{c}$ of the aging parameter, we set
$z_{j}=A$ for all active oscillators and $z_{j}=I$ for all inactive
oscillators. Consequently for a globally connected network Eq. (11) reduced
to,

$\displaystyle\dot{A}$ $\displaystyle=$
$\displaystyle(a+i\omega-|A|^{2}-2\epsilon-Kp-p\epsilon)A+(K-\epsilon)pI,$
$\displaystyle\dot{I}$ $\displaystyle=$
$\displaystyle(-b+i\omega-|I|^{2}-\epsilon-Kq-p\epsilon)I+(K-\epsilon)qA,$

where $q=1-p$. Employing a linear stability analysis of this reduced model
around the origin, we derive the critical value $p_{c}$ as

$\displaystyle p_{c}$ $\displaystyle=$
$\displaystyle\frac{(a-2\epsilon)(b+K+\epsilon)}{(a+b)(K-\epsilon)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for $\epsilon+K\geq a$}.$ (13)

When $\epsilon+K\leq a$, the critical parameter $p_{c}$ remains at unity. In
the limiting case as $\epsilon$ approaches 0, $\lim_{\epsilon\to
0}p_{c}=\frac{a(b+K)}{K(a+b)}$, aligning with the critical aging parameter for
a diffusively coupled global network. Additionally,
$\lim_{K\to\infty}p_{c}=\frac{a-2\epsilon}{a+b}$. Consequently, it can be
inferred that an increase in the repulsive interaction strength $\epsilon$
leads to a decrease in dynamical robustness.

We plot the order parameter $R$ as a function of aging parameter $p$ for
various scenarios. Figs. 6(a) and 6(b) depict the transition scenario for
$\epsilon=0$ and $\epsilon=0.2$ respectively, considering a spectrum of
attractive coupling strengths. The general trend of aging transition remains
same for both settings (absence and presence of repulsive interaction) for
increasing attractive coupling strength $K$. But, for the latter case of
$\epsilon=0.2$, the aging occurs for smaller values of the inactivation ratio
$p$ compared to the former case of $\epsilon=0$, implying a decline of the
robustness. This illustration demonstrates that the inclusion of repulsive
interaction significantly reduces the dynamical robustness. In Fig. 6(c), we
portray, how the transition takes place for increasing values of $\epsilon$,
with a fixed attractive coupling strength $K=0.2$. The figure clearly shows
that for increasing $\epsilon$, the aging transition keeps occurring for lower
values of $p$. Thus, Fig. 6(c) further confirms that repulsive interactions
negatively impacts dynamical robustness.

Figure 6: The order parameter $R$ as a function of aging parameter $p$, for
the system (11). (a) In absence of repulsive interactions ($\epsilon$=0 ) for
different values of attractive coupling strengths $K$. (b) In presence of
repulsive interactions ($\epsilon=0.2$) for different values of $K$. (c) For
different values of repulsive coupling strength $\epsilon$ at a fixed value of
$K=2$. Figure 7: The average amplitude (A) for the system (14) varying with
the inactivation ratio $p$ for various values of (a) attractive coupling
strength $k_{a}$ at fixed repulsive coupling $k_{r}=1$ and for different
repulsive coupling strength $k_{r}$ at fixed attractive coupling strength (b)
$k_{a}=5$ ,(c) $k_{a}=7$ and (d) $k_{a}=9$ for $N=500$ coupled Stuart-Landau
oscillators with $\omega=5$. Reprinted figure with permission from
Ref.sharma2022 .

——————————————————————————————

In contrast to the above study based upon the assumption of attractive and
repulsive interactions in both real and imaginary variables of Stuart-Landau
oscillators, our next analysis delves into the robustness of an alternative
form of attractive and repulsive interactions, as outlined in sharma2022 .
Consider $N$ coupled Stuart-Landau oscillators with attractive and repulsive
interactions given by,

$\displaystyle\dot{z}_{j}(t)$ $\displaystyle=$
$\displaystyle(\alpha_{j}+i\omega-|z_{j}(t)|^{2})z_{j}(t)+\frac{k_{a}}{\lambda_{j}}\sum_{l=1}^{N}B_{jl}(\mbox{Re}(z_{l})-\mbox{Re}(z_{j}))-\frac{k_{r}}{\lambda_{j}}\sum_{l=1}^{N}B_{jl}(\mbox{Img}(z_{l})-\mbox{Img}(z_{j})),$
(14)

Figure 8: The ($p,k_{r}$) parameter plane’s phase diagram is depicted for a
network of Stuart-Landau oscillators (14) that are coupled in an attractive-
repulsive manner for (a) when the attractive coupling strength is set to
$k_{a}=3$, and (b) when it is adjusted to $k_{a}=7$. The ($k_{a},k_{r}$)
parameter plane exhibits a phase diagram for two inactivation ratio (c)
$p=0.3$ and (d) $p=0.6$. The regions corresponding to homogeneous aging
transition (HAT), inhomogeneous aging transition (IHAT), and oscillatory state
(OS) are indicated by light blue, orange, and white colors, respectively.
Reprinted figure with permission from Ref.sharma2022 . Figure 9: Average
amplitude ($A$) and variance ($\rho$) for the system (14) as a function of
aging parameter $p$ for repulsive coupling (a) $k_{r}=2$ and (d) $k_{r}=3.6$
at a fixed attractive coupling $k_{a}=7$. The time series of the $x_{j}$
variable for the oscillator is depicted at various values of the inactivation
ratio (b) $p=0.4$, (c) $p=0.6$ at fixed $k_{r}=2$, and (e) $p=0.1$, (f)
$p=0.3$. Reprinted figure with permission from Ref.sharma2022 .

Here $k_{a}$, and $k_{r}$ indicate the magnitude of attractive and repulsive
interactions in coupled oscillators, respectively. $B_{jl}$ is the connection
matrix, and $\lambda_{j}$ is the average degree of the $j$th node of the
coupled oscillators. Two order parameters, namely, the average amplitude
$A_{amp}=\frac{1}{N}\sum_{i=1}^{N}(\langle x_{i,max}\rangle_{t}-\langle
x_{i,min}\rangle_{t})$ and the variance
$\rho=\sigma^{2}(\langle(x_{i})\rangle_{t})$ are used to identify the aging
transition. These order parameters are employed to differentiate between
homogeneous and inhomogeneous steady states. Specifically, when
$A_{amp}=\rho=0$, all oscillators converge to HSS, while $A_{amp}=0$ and
$\rho\neq 0$ indicate the presence of oscillators in IHSS.

In Fig. 7, the average amplitude ($A_{amp}$) is plotted with respect to $p$
for different combinations of the attractive coupling $k_{a}$ and the
repulsive coupling strength $k_{r}$. Figure 7(a) suggests that as we augment
$k_{a}$, the critical inactivation ratio $p_{c}$ decreases, indicating a
faster occurrence of the phase transition due to aging. This observation
aligns with earlier studies suggesting that increased attractive couplings
serve to dampen oscillations. Conversely, it is evident from Fig. 7(b) that
repulsive interactions enhance robustness to large extent. Yet, in Figs.
7(c-d), it becomes evident that as the repulsive interaction strength $k_{r}$
reaches higher values, a sudden phase transition unfolds in the order
parameter $A$. This transition is characterized by its catastrophic nature,
with the magnitude of the order parameter providing no forewarning of the
impending collapse of the system. Consequently, our observation suggests that
depending on parameter values, repulsive interactions have the potential to
bolster dynamical robustness in coupled oscillators, but concurrently, they
might instigate an abrupt suppression of macroscopic oscillations within the
network.

Subsequently, we delve into the mechanism behind this sudden catastrophic
phase transition by examining the 2-dimensional parameter space. The phase
diagram on the $k_{r}-p$ plane is presented in Figs. 8(a-d) for various values
of attractive coupling strength. This figure delineates three distinct
regions: OS (oscillatory state), IHAT (inhomogeneous aging transition state),
and HAT (homogeneous aging transition state). With an increase in $k_{a}$, we
note an expansion in both the IHAT and the HAT regions. Additionally, it is
observable that for higher values of $p$, the transition from homogeneous
steady-state to an inhomogeneous one occurs with variations in $k_{r}$. This
result is also apparent in Fig. 9, where we have graphed $A_{amp}$ and $\rho$
as functions of $p$, alongside the corresponding time series near the aging
transition point. In Fig. 9(a), a typical aging transition is evident as the
order parameters $A$ and $\rho$ exhibit smooth functions of $p$ and converge
to zero at a critical value $p=p_{c}$. The corresponding time-series of the
oscillators on both sides of the aging transition point $p_{c}$ is depicted in
Figs. 9(b-c). However, Figs. 9(d-f) portray an aging transition that is
qualitatively distinct from the previous one. In this case, a discontinuity in
the order parameters is noticeable at the aging transition point $p_{c}$.
Additionally, the corresponding time series indicates an aging transition
through an inhomogeneous steady-state, where oscillators settle into three
different steady states after the aging transition. Therefore, we can infer
that the abrupt discontinuous jump in the order parameter is a consequence of
a bifurcation from the oscillatory state to inhomogeneous steady states.

####  II.1.6 Dynamical robustness in the absence of inactive oscillators

Thus far, whatever we have discussed above is the scrutinization of the
phenomenon of dynamical robustness in ensembles of dynamical systems achieved
by increasing the proportion of inactive oscillatory systems. Now, we move
aside from this specification of the dynamical units, and demonstrate that
even an ensemble of co- and counter-rotating oscillatory systems
prasad2010universal ; bhowmick2011synchronization  can exhibit similar aging
transition upon raise in the fraction of counter-rotating oscillators. The
impact of the mean-field feedback in the symmetry preserving as well as the
symmetry breaking coupling is explored in Ref. sathiyadevi2019aging . The
network becomes dynamically vulnerable completely and hence the regime of
global oscillation transits to aging through a Hopf bifurcation, whereas the
transition from aging to oscillation death occurs via a pitchfork bifurcation.

The dynamical evolution of a network of $N$ globally interacting Stuart-Landau
oscillators with symmetry preserving coupling is described as,

$\begin{array}[]{lcl}\dot{z_{j}}=(\lambda+i\omega_{j}-|{z_{j}}|^{2})z_{j}+\dfrac{K}{N}\sum\limits_{k=1}^{N}(\alpha
z_{k}-z_{j});\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
j=1,2,\cdots,N,\end{array}$ (15)

where $\alpha$ refers to the strength of the mean-field feedback. Moreover,
$\omega_{j}$ is the frequency of the $j$-th dynamical unit. Specifically, for
a frequency $+\omega$, the system rotates in a counter-clockwise direction,
and in the clockwise direction whenever the frequency is $-\omega$. Analogous
to the studies with inactive oscillators, we divide the entire system into two
groups. To be precise, the frequency $\omega_{j}=\omega$ for $j\in
1,2,\dots,N-Np$ and frequency $\omega_{j}=-\omega$ for $j\in N-Np+1,\dots,N$,
and consequently the parameter $p$ characterizes the proportion of the
counter-rotating oscillators in the networked system.

Figure 10:  Normalized order parameter $R$, for the system (15) with respect
to the fraction $p$ of the counter-rotating oscillators for (a) $\alpha=1.0$
and (b) $\alpha=0.95$, with fixed $\omega=5$, $\lambda=1$ and $N=100$. Four
different values of the interaction strength are chosen as $K=1,4,7$, and
$10$. The corresponding phase diagrams in the $(K,p)$ parameter plane are
plotted in (c) and (d), respectively. The oscillatory and aging regimes are
respectively represented by OS and AG. Reprinted figure with permission from
Ref.sathiyadevi2019aging .

Further choosing the network size $N=100$, frequency $\omega=5.0$, and
defining the normalized order parameter $R$ similarly as before, we plot $R$
with respect to the fraction $p$ of counter-rotating oscillators in Fig. 10,
for two different values of the mean-field feedback strength $\alpha=1.0$ and
$\alpha=0.95$. The dynamical outcomes being symmetric for both the ranges
$p\in(0,0.5)$ and $p\in(0.5,1)$, we effectively confine ourselves to the range
$p\in(0,0.5)$ in the plots. In Fig. 10(a), we encounter that whenever
$\alpha=1.0$, with $K=1$, $R$ remains non-zero finite for the entire range of
$p$, depicting oscillatory regime of all the individual units and hence the
whole network. However, with higher coupling strengths $K=4,7,$ and $10$, the
order parameter $R$ transits from non-zero to zero value and hence aging
transition takes place through Hopf bifurcations at $p_{HB}=0.28,0.20,$ and
$0.19$, respectively. Thus, increasing coupling strength decreases the
critical value of $p$ at which aging occurs via Hopf bifurcation. In order to
perceive the impact of the mean-field feedback on the dynamical robustness, we
depict similar transitions for a lower feedback strength $\alpha=0.95$, with
the same set of interaction strengths (in Fig. 10(b)). Similarly as before,
the entire system shows oscillatory nature for $K=1$. For increasing $K=4,7,$
and $10$, the aging transition occurs at the Hopf bifurcation points
$p_{HB}=0.24,0.15,$ and $0.10$ respectively. Thus, even for a tiny decrement
in the feedback strength $\alpha$, the onset of aging transition takes place
for smaller values of the fraction $p$. For a better understanding, we further
portray this scenario through the cooresponding phase diagrams in the $(K,p)$
parameter planes, in the lower panel of Fig. 10. As can be witnessed, for both
values of the mean-field feedback strength ( Figs. 10(c) and 10(d)), the
networked system displays oscillatory solutions for a small interaction
strength and smaller proportion of the counter-rotating oscillators. However,
as $K$ and $p$ surpass certain values, aging transitions take place. The
regimes of oscillation and aging are classified as OS and AG, respectively.
More importantly, for smaller feedback strength, the aging region expands in
the phase diagram, as expected from earlier observations.
Our investigation thus demonstrates how networked dynamical systems comprising
of co- and counter-rotating oscillators instead of active-inactive
decomposition, can also give rise to the phenomenon of aging transition and,
hence, a potential of the study of dynamical robustness of networked systems
from a different perspective.

###  II.2 Dynamical robustness of complex networks

In the above, we have examined the dynamical robustness of globally coupled
networks under various coupling schemes. In this section, we turn our
attention to aging transitions within different complex network topologies. In
the past two decades, network science has experienced significant advancements
due to the discovery of various topological structures in many real-world
networks barabasi2016network ; newman2003structure ; boccaletti2006complex ;
strogatz2001exploring ; albert2002statistical ; estrada2012structure ;
costa2011analyzing . The complexity of a network structure can be
characterized by the connectivity properties of interaction pathways (links)
among network components (nodes). In terms of the degree distribution (the
probability distribution of node degrees across the network), complex networks
mainly fall into two categories: homogeneous and heterogeneous networks.
Homogeneous networks, exemplified by random graphs erdHos1960evolution  and
small-world models watts1998collective , exhibit a binomial or Poisson degree
distribution where node degrees cluster around the mean degree. In contrast,
heterogeneous networks like scale-free networks display a heavy-tailed degree
distribution that approximately follows a power-law distribution
barabasi1999emergence . We will first discuss the results obtained for
homogeneous complex networks, followed by the study of heterogeneous networks.
Then, we will present the findings for weighted complex networks, and finally,
we will discuss the dynamical robustness within a correlated network topology.

####  II.2.1 Homogeneously coupled complex networks

Tanaka et al. tanaka2012dynamical  extended the work of Daido and Nakanishi
Daido_2004_PRL  to complex network structures. They considered a homogeneously
coupled network having a Poisson degree distribution. In particular they
explored an Erdös-Rényi random graph erdHos1960evolution  by considering a
network of diffusively coupled oscillators as follows:

$\dot{z_{j}}(t)=(\alpha_{j}+i\omega-|z_{j}(t)|^{2})z_{j}(t)+\frac{K}{N}\sum_{l=1}^{N}A_{jl}(z_{l}-z_{j}).$
(16)

Using the system reduction techniques proposed by Daido and Nakanishi
Daido_2004_PRL , the critical value $p_{c}$ of the inactivation parameter can
be analytically calculated. For a homogeneous network, the degree of a node
can be approximated by the average degree of the network. With this
assumption, we can conclude that each active oscillator is connected to
neighboring $(1-p)\langle k\rangle$ oscillators and each inactive oscillator
is to $p\langle k\rangle$ neighboring oscillators. Here $\langle k\rangle$
represents the mean degree of the network. When we designate $z_{j}$ as ’A’
for the active group and ’I’ for the inactive group of oscillators, Eq. (16)
in its original form transforms then into the following coupled system:

$\displaystyle\dot{A}(t)$ $\displaystyle=$
$\displaystyle(a-Kpd+i\omega-|A(t)|^{2})A(t)+KpdI(t),$
$\displaystyle\dot{I}(t)$ $\displaystyle=$
$\displaystyle(-b-K(1-p)d+i\omega-|I(t)|^{2})I(t)+K(1-p)dA(t).$

A linear stability analysis of the equilibrium point $(A,I)=(0,0)$ leads to

$\displaystyle p_{c}^{hom}=\frac{a(Kd+b)}{(a+b)Kd}$ $\displaystyle\text{for
$K>K_{c}^{hom}$},$ (18)

where $d$ represents the link density, defined as $\langle k\rangle/N$, and
$K_{c}^{hom}=a/d$ provides the critical coupling strength, below which
$p_{c}^{hom}=1$. Here we consider $\alpha_{j}=a$ for an active oscillator and
$\alpha_{j}=-b$ for an inactive oscillator. From the Eq. (18), we can conclude
that for a fixed $K$, the robustness increases with the decrease of link
density $d$.

####  II.2.2 Dynamical robustness of heterogeneous networks

To understand how robust a networked dynamical system is under increasing
inactivation of its dynamical units whenever the degree distribution of the
network is heterogeneous, we present a detailed comparative analysis of aging
transition between homogeneous and heterogeneous networks tanaka2012dynamical
, while considering a network of $N$ dynamical systems as described by the
networked system Eq. (16). In order to derive the desired expression for the
critical inactivation ratio $p_{c}$ for heterogeneous networks under random
inactivation, we follow the degree-weighted mean field approximation
pastor2001epidemic . Then the original system (Eq. (16) can be approximated as

$\begin{array}[]{lcl}\dot{z_{j}}=(\alpha+i\omega-|z_{j}|^{2})z_{j}+\dfrac{K{k_{j}}}{N}\Big{[}(1-p)M_{A}(t)+pM_{I}(t)-z_{j}\Big{]},\\\
\end{array}$ (19)

where

$\begin{array}[]{lcl}M_{A}(t)=\dfrac{\sum_{j\in
S_{A}}{k_{j}}z_{j}(t)}{\sum_{j\in S_{A}}{k_{j}}},\leavevmode\nobreak\
\leavevmode\nobreak\ M_{I}(t)=\dfrac{\sum_{j\in
S_{I}}{k_{j}}z_{j}(t)}{\sum_{j\in S_{I}}{k_{j}}}\end{array}$ (20)

are the degree-weighted mean fields for the active and inactive sets of
dynamical units, with ${k_{j}}(j=1,2,...,N)$ being the degree of the $j$-th
node in the network.
Based upon the fact that the oscillators display phase synchronization with
frequency $\Omega$, let us now assume that the state variables can be written
as $z_{j}(t)=r_{j}(t)e^{i(\Omega t+\theta)}$, $r_{j}$ being the amplitude and
$\theta$ being the phase shift. On replacing this in Eq. (19), we get

$\small{\begin{array}[]{lcl}\dot{r}_{j}=\Big{(}\alpha_{j}-\dfrac{Kk_{j}}{N}-r_{j}^{2}\Big{)}r_{j}+\dfrac{Kk_{j}}{N}\Big{[}(1-p)R_{A}(t)+pR_{I}(t)\Big{]},\\\
\end{array}}$ (21)

where

$\begin{array}[]{lcl}R_{A}(t)=\dfrac{\sum_{j\in
S_{A}}{k_{j}}r_{j}(t)}{\sum_{j\in S_{A}}{k_{j}}},\leavevmode\nobreak\
\leavevmode\nobreak\ R_{I}(t)=\dfrac{\sum_{j\in
S_{I}}{k_{j}}r_{j}(t)}{\sum_{j\in S_{I}}{k_{j}}}.\end{array}$ (22)

Further supposing that in the stationary oscillatory regime, $R_{A}(t)$ and
$R_{I}(t)$ are time-independent, the phase transition from the oscillatory
($R_{A},R_{I}>0$) regime to the non-oscillatory ($R_{A}=R_{I}=0$) regime
eventuates due to the change in the stability of the equilibrium point at the
origin. The stability is governed by the following Jacobian matrix

$J_{0}=\begin{pmatrix}\dfrac{\partial G_{A}(R_{A},R_{I})}{\partial
R_{A}}&\dfrac{\partial G_{A}(R_{A},R_{I})}{\partial R_{I}}\\\ \dfrac{\partial
G_{I}(R_{A},R_{I})}{\partial R_{A}}&\dfrac{\partial
G_{I}(R_{A},R_{I})}{\partial
R_{I}}\end{pmatrix}\Biggr{\rvert}_{R_{A}=R_{I}=0},$

in which

$\begin{array}[]{lcl}G_{A}(R_{A},R_{I})=\dfrac{\sum_{j\in
S_{A}}{k_{j}}r_{j}^{*}(R_{A},R_{I})}{\sum_{j\in S_{A}}{k_{j}}},\\\ \\\
G_{I}(R_{A},R_{I})=\dfrac{\sum_{j\in
S_{I}}{k_{j}}r_{j}^{*}(R_{A},R_{I})}{\sum_{j\in S_{I}}{k_{j}}}.\end{array}$
(23)

The stationary amplitude $r_{j}^{*}$ is a positive real solution of the
following equation,

$\begin{array}[]{lcl}r_{j}^{3}-\Big{(}\alpha_{j}-\dfrac{Kk_{j}}{N}\Big{)}r_{j}-\dfrac{Kk_{j}}{N}\Big{(}(1-p)R_{A}+pR_{I}\Big{)}=0.\end{array}$
(24)

Eq. (24) has only one positive real root if

$\begin{array}[]{lcl}\alpha_{j}-\dfrac{Kk_{j}}{N}<0,\leavevmode\nobreak\
\forall j\in S_{A}.\end{array}$ (25)

Now we differentiate Eqns. (23) and (24) with respect to $R_{A}$ and we find
the $(1,1)$-th entry of $J_{0}$ as,

$\begin{array}[]{lcl}\dfrac{\partial G_{A}}{\partial
R_{A}}\Bigr{\rvert}_{R_{A}=R_{I}=0}=\dfrac{\dfrac{(1-p)}{N}K}{\sum_{j\in
S_{A}}{k_{j}}}\Big{[}\sum\limits_{j\in
S_{A}}\dfrac{{k_{j}}^{2}}{(Kk_{j}/N)-\alpha_{j}}\Big{]}\\\ \\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \simeq\dfrac{1}{d}\Big{(}\dfrac{1}{N}\sum\limits_{j\in
S_{A}}\dfrac{{d_{j}}^{2}}{d_{j}-\alpha_{j}/K}\Big{)},\par\\\ \end{array}$ (26)

in which $d_{j}=k_{j}/N\leavevmode\nobreak\ (j=1,2,...,N)$ is the normalized
degree of the $j$-th node, where the approximation $\sum_{j\in
S_{A}}{k_{j}}\simeq(1-p)dN^{2}$ is used. Also, in the large $N$ limit, the
link density is $d=\langle k\rangle/(N-1)$.

If we now define

$\begin{array}[]{lcl}H(K,\alpha)=\dfrac{1}{N}\sum\limits_{j=1}^{N}\dfrac{{d_{j}}^{2}}{d_{j}-\alpha/K},\end{array}$
(27)

then we can write

$\begin{array}[]{lcl}\dfrac{\partial G_{A}}{\partial
R_{A}}\Bigr{\rvert}_{R_{A}=R_{I}=0}\simeq(1-p)H(K,a)/d.\end{array}$ (28)

Similarly, we obtain

$\begin{array}[]{lcl}\dfrac{\partial G_{A}}{\partial
R_{I}}\Bigr{\rvert}_{R_{A}=R_{I}=0}\simeq pH(K,a)/d,\\\ \dfrac{\partial
G_{I}}{\partial R_{A}}\Bigr{\rvert}_{R_{A}=R_{I}=0}\simeq(1-p)H(K,-b)/d,\\\
\dfrac{\partial G_{I}}{\partial R_{I}}\Bigr{\rvert}_{R_{A}=R_{I}=0}\simeq
pH(K,-b)/d.\\\ \end{array}$ (29)

Thus we arrive at

$J_{0}=\begin{pmatrix}(1-p)H(K,a)/d&pH(K,a)/d\\\ \leavevmode\nobreak\
\leavevmode\nobreak\ (1-p)H(K,-b)/d&\leavevmode\nobreak\ \leavevmode\nobreak\
pH(K,-b)/d\end{pmatrix}.$

Figure 11:  (a) Critical inactivation ratio $p_{c}$ as a function of the
interaction strength $K$ in networked systems possessing $d\sim 0.08$
($N=3000$ and $\langle k\rangle=240$). The dashed and solid black curves
respectively stand for the analytically obtained results in Eqs. (18) and
(30). Red triangles and blue diamonds represent the numerical outcomes. The
error bars correspond to the variance for $10$ network realizations. (b) The
critical ratio $p_{c}$ with respect to the link density $d$ in networked
system with $N=3000$ and $K=30$. Reprinted figure with permission from
Ref.tanaka2012dynamical .

The equilibrium point at $R_{A}=R_{I}=0$ changes its stability near the phase
transition, providing us with the following critical inactivation ratio

$\begin{array}[]{lcl}{p_{c}}^{het}=\dfrac{H(K,a)-d}{H(K,a)-H(K,-b)}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{for}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ K>{K_{c}}^{het},\end{array}$ (30)

in which the definition of $H$ follows from Eq. (27). The result is true if
$K>a/d_{min}$ where $d_{min}=k_{min}/N$.

Figure 11 portrays a comparative analysis between the critical inactivation
ratios obtained for the homogeneous and the heterogeneous networks. In Fig.
11(a), we depict the critical ratio $p_{c}$ with respect to the interaction
strength $K$, with a fixed linked density of $d\sim 0.08$ for a network of
size $N=3000$ and the mean-degree $\langle k\rangle=240$. The dashed and solid
black curves respectively stand for the analytically obtained outcomes in Eqs.
(LABEL:sf2) and (30), whereas the red triangles and blue diamonds represent
the numerical results. As can be witnessed, the analytically obtained
expressions are in sufficiently good agreement with the numerical results.
$p_{c}$ values expectedly start decreasing with increasing $K$, after the
respective certain critical interaction strengths ${K_{c}}^{hom}$ and
${K_{c}}^{het}$. The critical ratio $p_{c}$ for the homogeneous network is
smaller than that for the heterogeneous one for the entire range of the
coupling strength $K\in[0,50]$. Similar results are shown in Fig. 11(b), but
this time for varying link density $d\in[0,0.1]$ and fixed coupling strength
$K=30$. Qualitatively the same scenario remains valid here as far as the
evolution of $p_{c}$ are concerned.

This analytical study for random failures has further been generalized to
targeted attacks yuan2017robustness . The study presents a universal formula
for the critical fraction of inactive units, applicable to both random
failures and targeted attacks on networked systems. It examines the impact of
targeting nodes based on their degrees, starting with identical oscillators
and homogeneous edge weights. The theory is subsequently extended to networks
with heterogeneous edge weights and non-identical oscillators. The analytical
findings are confirmed through extensive numerical simulations. On the other
hand, Tanaka et al. tanaka2014dynamical  presented a general formula for the
critical inactivation ratio for interacting heterogeneous oscillators. This is
done while assuming different values of the units’ intrinsic parameters,
instead of choosing the same fixed parameter value for all the active and all
the inactive elements. The study demonstrates that increasing heterogeneity in
the oscillator components of networks leads to improved dynamical robustness,
as evidenced by the comparison of critical values for networks with various
extents of heterogeneity. This observation is further validated for networks
of Morris-Lecar neuronal systems communicating through electrical synapses.

####  II.2.3 Weighted complex networks

The majority of the early research on dynamical robustness focused on networks
without taking weightings into account. But in reality, many complex networks
are indeed weighted, and the connection strength of the nodes highly
influences the network dynamics barrat2004architecture ; song2005highly ;
scannell1999connectional . Both the topology and the strength of a network’s
connections have an impact on its dynamics. In particular, studies have shown
that in complex networks, the presence of both degree heterogeneity and weight
heterogeneity tends to impede full synchronization. However, synchronization
can be notably enhanced and made unaffected by these types of heterogeneity,
when the weight distribution is appropriately integrated with the degree
distribution motter2005network ; motter2005weighted ; motter2005enhancing .

He et al. Zhan_2013_Physica-A  first studied dynamical robustness in a
weighted complex network. Their investigation highlights the roles of high and
low degree nodes on dynamical robustness in weighted complex networks and
generalizes the works of Tanaka et al. tanaka2012dynamical . Study of
dynamical robustness in a weighted complex network has been carried out for
$N$ diffusively coupled Stuart-Landau oscillators as expressed below,

$\dot{z_{j}}(t)=(\alpha_{j}+i\omega-|z_{j}(t)|^{2})z_{j}(t)+K\sum_{l=1}^{N}W_{jl}A_{jl}(z_{l}-z_{j}).$
(31)

Here $A=(A_{jl})$ and $W=(W_{jl})$ are the adjacency and weight matrix of the
network, respectively. The weight matrix $W$ is defined as
$W_{jl}=\frac{1}{k_{j}^{\beta}}$, where $k_{j}$ is the degree of the $j$th
node and $\beta$ is a tunable parameter, whose value determines whether the
network is weighted ($\beta\neq 0$) or unweighted ($\beta=0$). Dynamical
robustness has been studied following the mathematical framework proposed by
Daido and Nakanishi Daido_2004_PRL .

Figure 12: The critical values $p_{c}$ (for the system (31)) are plotted
against $K$ for random inactivation (circles), targeted inactivation of high-
degree nodes (stars), and that of low-degree nodes (squares) in a weighted
heterogeneous networks. $K_{c}^{R}$, $K_{c}^{H}$, and $K_{c}^{L}$ represent
their corresponding critical coupling strengths $K$ at $p_{c}=1$. (a)
$\beta=0.2$, (b) $\beta=0.6$, (c) $\beta=1.0$, and (d) $\beta=1.4$. For
$\beta<1$, an intersection point ($K_{cin},p_{cin}$) exists, whereas for
$\beta>1$, it vanishes and the curve of target on the low-degree nodes is
always above that of target on the high-degree nodes. In (c),
$K_{c}^{R}=K_{c}^{H}=K_{c}^{L}=K_{cin}$. Reprinted figure with permission from
Ref.Zhan_2013_Physica-A .

In Fig. 12, the critical values $p_{c}$ are depicted against $K$ for different
$\beta$ values in a weighted heterogeneous network. The graph includes data
for random inactivation (circles), targeted inactivation of high-degree nodes
(stars), and targeted inactivation of low-degree nodes (squares). Three
distinct critical couplings at $p_{c}=1$ are indicated as $K_{c}^{R}$,
$K_{c}^{H}$, and $K_{c}^{L}$. In Fig. 12(a-b) we observe that all the critical
curves intersect at a point, denoted by $(K_{cin},p_{cin})$. Clearly for
$0<\beta<1$, $p_{c}^{L}<p_{c}^{H}$, while $K<K_{cin}$ but for $K>K_{cin}$ it
becomes completely opposite. This implies the crucial role of low-degree nodes
in impacting the dynamical robustness which can only happen in a weakly
weighted network and for weak coupling. For $\beta=1.0$, $p_{cin}=1$ and
$K_{c}^{R}=K_{c}^{H}=K_{c}^{L}=K_{cin}$ as shown in Fig. 12(c). From Fig.
12(d) it is quite obvious that for $\beta>1$, $K_{c}^{H}<K_{c}^{R}<K_{c}^{L}$.
By applying the heterogeneous mean-field approximation, we can derive the
critical value of the inactivation parameter as

$\displaystyle p_{c}^{het}=\frac{F(K,a)-\langle k\rangle N}{F(K,a)-F(K,-b)}$
$\displaystyle\text{for $K>K_{c}^{hetR}$},$ (32)

where
$F(K,\alpha)=\sum_{j=1}^{N}\frac{Kk_{j}^{2-\beta}}{-\alpha+Kk_{j}^{1-\beta}}$,
and $K_{c}^{hetR}=max_{j\in S_{A}}\\{ak_{j}^{\beta-1}\\}$. Here we consider
$\alpha_{j}=a$ for an active oscillator and $\alpha_{j}=-b$ for an inactive
oscillator.

The phenomenon reported in the Ref. tanaka2012dynamical , where heterogeneous
networks are more susceptible to the failure of low-degree nodes rather than
high-degree nodes, is observed specifically in weakly weighted networks and
under weak coupling conditions. However, when considering the impact of
weighted connections, we discover that the susceptibility to high-degree node
failures is more widespread and occurs across a broader parameter range. This
finding indicates a strong alignment between dynamical and structural
robustness. In the case of the unweighted version of a heterogeneous network,
high-degree nodes (hubs) are impacted by numerous neighbors, whereas low-
degree nodes are influenced by only a small number of neighboring nodes,
rendering them relatively isolated. Consequently, active low-degree nodes can
sustain relatively high dynamical activity in comparison to active high-degree
nodes. Consequently, carrying out targeted inactivation on low-degree nodes
(rather than high-degree nodes) leads to a notable decrease in the network’s
dynamical activity. This aligns precisely with the observations made by Tanaka
et al. in Ref.tanaka2012dynamical . Nonetheless, as the overall coupling
strength K increases, the exchange of information and dynamics among these
networked nodes becomes more effective, eliminating the isolating effect
previously experienced by low-degree nodes. As the weighted coupling scale
$\beta$ increases, the coupling on each node (influenced by the weight matrix
W) becomes more evenly distributed, diminishing the system’s heterogeneity
associated with the adjacency matrix $A$. A notable instance is $\beta=1$,
where all nodes experience the same level of input signal intensity.
Consequently, the distinctive isolating effect of low-degree nodes diminishes
as well. Under these conditions, it is plausible to infer that the typical
influential behavior exhibited by high-degree nodes extends across a broader
parameter range, aligning with the findings of the structural robustness
analysis.

Recently, Ray et al. Ray_2020_EPL  delved into aging transitions within a
weighted heterogeneous network. In their investigation, weights are randomly
selected from a uniform distribution $[0,w]$, where no connection exists
between node degrees and the weights. Their findings reveal a direct link
between weight heterogeneity and the aging transition point $p_{c}$. It was
observed that dynamical robustness diminishes as the mean value $\bar{w}$ of
the weight distribution increases. Heightened heterogeneity among connection
weights correlates with reduced dynamical robustness. Moreover, the analytical
expression of the critical value $p_{c}$ depends on both the mean weights and
the network’s average degree.

####  II.2.4 Correlated networks

Up to this point, we have explored the dynamical robustness of complex
networks with structures characterized by degree distributions, which
represent the probability distribution of the number of connections per node
throughout the network. In particular, we have examined the dynamical
robustness of these networks in relation to both homogeneously and
heterogeneously connected networks featuring various degree distributions.
Nonetheless, it is crucial to recognize that the degree distribution alone
does not completely define the network’s topology. Networks with identical
degree distributions can exhibit diverse network structures. These
distinctions can be quantified by examining network assortativity concerning
node degrees (degree-degree correlations), assessing the clustering
coefficient, and considering various network
characteristicsnewman2002assortative ; mcgraw2005clustering . Network
assortativity, in particular, measures the correlation between a node’s degree
and the degrees of its neighboring nodes. In assortative networks, there is a
positive correlation, meaning that nodes tend to link with others of similar
degrees. Conversely, in disassortative networks, there is a negative
correlation, indicating that high-degree nodes are more inclined to connect
with low-degree nodes. Here, we examine the influence of network assortativity
on the dynamical robustness of coupled oscillator networks, as explored in the
study sasai2015robustness .

The assortativity coefficient $r$ is determined by computing the Pearson
correlation coefficient of the degrees between pairs of connected nodes, and
is calculated as follows,

$\displaystyle r=\frac{1}{\sigma^{2}_{q}}\sum_{j}\sum_{k}jk(E(j,k)-Q(j)Q(k)).$
(33)

Here $Q$ represents the probability distribution of the remaining degree,
which quantifies the probability that a node in the end of a randomly chosen
edge has $k$ edges except for the chosen one. The distribution of $Q$ is
derived from the degree distribution $P(k)$ as
$Q(k)=\frac{(k+1)P(k+1)}{\sum_{j}jP(j)}$. The term
$\sigma^{2}_{q}=\sum_{k}k^{2}Q(k)-(\sum_{k}Q(k))^{2}$ is the variance of the
distribution $Q(k)$, and $E(j,k)$ represents the joint probability
distribution of the remaining degrees of two vertices. This distribution
exhibits symmetry in undirected graphs and adheres to sum rules
$\sum_{j}\sum_{k}E(j,k)=1$ and $\sum_{j}E(j,k)=Q(k)$. The assortativity
coefficient $r$ in Eq (33) can be rewritten as follows,

$\displaystyle
r=\frac{4M\sum_{m}j_{m}k_{m}-[\sum_{m}(j_{m}+k_{m})]^{2}}{2M\sum_{m}(j^{2}_{m}+k^{2}_{m})-[\sum_{m}(j_{m}+k_{m})]^{2}},$
(34)

here $M$ represents number of edges in the network, $m\in{1,2,..,M}$ is the
index of edges, and $j_{m}$ and $k_{m}$ represent the degrees of the two nodes
$j$, and $k$ connected by the edge $m$. The assortativity coefficient $r$
varies within the range of $-1$ to $1$. A value of $r$ greater than $0$
signifies assortative networks, $r=0$ indicates uncorrelated networks, and an
$r<0$ implies disassortative networks.

We investigate the dynamical robustness by altering $r$ through two
approaches: Greedy edge rewiring (GER) xulvi2004reshuffling  and Stochastic
edge rewiring (SER) newman2002assortative . We start with an uncorrelated
network, characterized by $r=0$, and then perform edge reshuffling without
permitting self-loops or overlaps. We then randomly select two existing edges
of the network, given by the connected node pairs $(v_{1},w_{1})$ and
$(v_{2},w_{2})$. The remaining degrees of these node pairs are denoted as
$(j_{1},k_{1})$ and $(j_{2},k_{2})$ respectively. i) The degree of the
connected node pairs is used to guide the edge rewiring in the GER approach.
We arrange the remaining degrees, namely, $j_{1}$, $j_{2}$, $k_{1}$, and
$k_{2}$ in descending order and assign them new labels, $l_{1}$, $l_{2}$,
$l_{3}$, and $l_{4}$, respectively, ensuring that $l_{1}\geq l_{2}\geq
l_{3}\geq l_{4}$. There are three possible ways to partition the four nodes
into two pairs of connected nodes as depicted in Figure 13(a). To enhance the
assortativity coefficient, we opt for Case I to establish connections between
nodes with more analogous degrees if the current state is Case II or III.
Conversely, to reduce the assortativity coefficient, we employ Case III,
establishing an edge between nodes with the highest and lowest degrees if the
current state is Case I or II. With repeated employment of edge rewiring in a
greedy manner we can monotonically increase or decrease $r$ until it no longer
changes. ii) In the Stochastic Edge Rewiring (SER) technique, we iteratively
rewire edges in a stochastic manner, altering the network’s assortativity. We
aim to construct a network that adheres to a predefined joint probability
distribution E(j, k) for the remaining degrees. This process is implemented
using a numerical method based on newman2002assortative .

Figure 13: Network rewiring method for changing the network assortativity $r$.
(a) The greedy edge-rewiring (GER) method xulvi2004reshuffling . (b) The
stochastic edge-rewiring (SER) method. Reprinted figure with permission from
Ref.sasai2015robustness . Figure 14: The critical value $p_{c}$ is plotted as
a function of the assortativity coefficient $r$ for random and targeted
inactivation, for the system (35). We start from a Erdös-Rényi uncorrelated
network and employ (a) GER method and (b) SER method, while altering the
assortativity coefficient $r$. We start from the BA scale-free network with
$r=0$ and employ (c) GER method and (d) SER method, and alter the
assorttativity coefficient $r$. Reprinted figure with permission from
Ref.sasai2015robustness .

To study the effect of network assortivity on dynamical robustness, we
consider the network model consisting of $N$ diffusively coupled Stuart-Landau
oscillators which is described by,

$\displaystyle\dot{z}_{j}$ $\displaystyle=$
$\displaystyle(\alpha_{j}+i\omega-|z_{j}|^{2})z_{j}+\frac{K}{N}\sum_{k=1}^{N}A_{jk}(z_{k}-z_{j}).$
(35)

First, we investigate the dynamical robustness of an uncorrelated Erdős-Rényi
random graph erdHos1960evolution , with degrees centered around the average
degree. Subsequently, we modify the network to make it assortative or
disassortative by applying the edge-rewiring algorithms outlined above.
Figures 14(a) and 14(b) display the critical value $p_{c}$ as a function of
$r$ for the networks generated using the GER and SER methods, respectively. In
both panels, for random inactivation, the value of $p_{c}$ remains almost
constant, regardless of the $r$ value. This is due to the fact that the number
of inactive oscillators in the vicinity of each oscillator node is not
influenced by the value of $r$. The targeted inactivation of high-degree and
low-degree oscillator nodes both exhibit a monotonic increase in the value of
$p_{c}$ with $r$, as depicted in Figs. 14(a) and 14(b). This outcome can be
attributed to the fact that in more assortative networks, the amplitudes of
the active oscillators, which play a dominant role in determining the order
parameter, are larger. Consequently, it can be concluded that network
assortativity has a positive impact on the dynamical robustness of oscillator
networks when facing targeted inactivation.

Next, we explore the impacts of assortativity on dynamical robustness of
correlated networks with power-law degree distributions barabasi1999emergence
. The critical fraction $p_{c}$ plotted versus $r$ for the GER and SER
techniques are shown in Figs. 14(c) and 14(d), respectively. In assortative
networks, it is evident from the figure that for all types of inactivation,
the value of $p_{c}$ consistently increases as $r$ progresses from $0$. In
assortative networks, connections tend to be formed between high-degree nodes
and between low-degree nodes. Consequently, when targeting the inactivation of
high-degree nodes, low-degree active oscillators that are connected to a few
inactive oscillators can sustain large oscillation amplitudes. Likewise, when
targeting the inactivation of low-degree nodes, high-degree active oscillators
connected to a few inactive oscillators can also maintain large oscillation
amplitudes. These nodes, which preserve substantial oscillation amplitudes,
contribute significantly to the high value of $p_{c}$, signifying a highly
robust oscillatory behavior. Thus, assortativity plays a positive role in
enhancing dynamical robustness. In contrast, for disassortative network as $r$
is decreased from $0$, the dynamical robustness increases after a slight
downward trend for the GER method, as shown in Fig. 14(c), but it gradually
decreases until $r=-0.5$ for the SER method, as shown in Fig. 14(d). In
summary, we conclude that network assortativity enhances dynamical robustness,
while the impact of network disassortativity on dynamical robustness depends
on the specific edge-rewiring methods employed.

###  II.3 Dynamical robustness of multiplex networks

Recent studies have further confirmed that the functions emerging within a
single network can exert a notable impact on other networks. Specifically, a
node within one network is often found to be a component of another network.
From ecological pilosof2017multilayer  and climate systems
donges2011investigating , via physical and transportation systems to social
networks szell2010multirelational , these interconnections have been
identified across various contexts. This underscores the effectiveness of an
interdependent gao2012networks ; buldyrev2010catastrophic , particularly
multilayer (multiplex) boccaletti2014structure ; kivela2014multilayer ,
network of networks architecture in describing numerous systems and scenarios.
Thus, a multilayer network boccaletti2014structure ; kivela2014multilayer ;
de2013mathematical ; pilosof2017multilayer ; vaiana2020multilayer  is simply
defined to be a network with a set of nodes, together with the concept of
layers in addition to nodes and links. These layers are often present in order
to represent different types or aspects of interactions (links). In
particular, a special type of multilayer network having same number of nodes
in each layer, in which each node in each layer posseses exactly one
connection with a node (its replica) in another layer is known to be a
multiplex network nicosia2013growing ; gomez2013diffusion ;
verbrugge1979multiplexity ; granell2013dynamical ; majhi2019emergence ;
cardillo2013emergence .

Let us here consider a multilayer network consisting of $L$ layers, each layer
comprised of $N$ globally coupled nodes. Casting the dynamics of the nodes of
the network by Stuart-Landau oscillatory systems, the time evolution of the
entire multilayer network is then governed by the following Eq.
morino2011robustness ,

$\begin{array}[]{lcl}\dot{z_{j}}^{[l]}=\\{\alpha_{j}+i\omega-|{z_{j}}^{[l]}|^{2}\\}{z_{j}}^{[l]}+\dfrac{\epsilon}{N}\sum\limits_{k=1}^{N}({z_{k}}^{[l]}-{z_{j}}^{[l]})+\sigma
H({z_{j}}^{[1]},{z_{j}}^{[2]},\dots,{z_{j}}^{[L]}),\par\end{array}$ (36)

where ${z_{j}}^{[l]}$ represents the state of the $j$-th node in the $l$-th
layer, $j=1,2,\dots,N;\leavevmode\nobreak\ l=1,2,\dots,L$. Here $\epsilon$ is
the intralayer coupling strength and $\sigma$ corresponds to the interlayer
interaction strength with $H(\cdot)$ being the function characterizing the
interlayer connections.

Then, three different interlayer coupling functions are assumed, namely mean-
field, chain and diffusive interlayer interactions which are realized by the
following forms

$H({z_{j}}^{[1]},{z_{j}}^{[2]},\dots,{z_{j}}^{[L]})=\begin{cases}\frac{1}{L}\sum\limits_{m=1}^{L}{z_{j}}^{[m]}&\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{(case I)},\\\
{z_{j}}^{[l-1]}/2&\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{(case
II)},\\\
\frac{1}{L}\sum\limits_{m=1}^{L}({z_{j}}^{[m]}-{z_{j}}^{[l]})&\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{(case III)}.\end{cases}$

Next, we modify the order parameter for the multilayer network as
$R=\dfrac{1}{NL}|\sum\limits_{l=1}^{L}\sum\limits_{j=1}^{N}{z_{j}}^{[l]}|$
characterizes the intensity of the global oscillation in the entire multilayer
networked system. We then set the value of the system parameters as
$\alpha_{j}=2$ for an active oscillator, $\alpha_{j}=-1$ for an inactive
oscillator, and $\omega=3$, while keeping the network size $N=3000$ and the
coupling strength fixed at $\epsilon=8$.

Figure 15: (a) An exemplary multiplex network with $N=5$ nodes in each of
$L=2$ layers. The red and blue circles respectively stand for the active and
inactive oscillators. The ratio of the inactive units is $p=0.4$. The black
solid and green dashed lines indicate all-to-all intralayer connections and
local interlayer connections, respectively. (b) The order parameter $R$ with
respect to the fraction of inactive elements $p$, for three different
scenarios: single-layer network, multi-layer network ($L=2$, case II), and
multi-layer network ($L=2$, case III). Reprinted figure with permission from
Ref.morino2011robustness .

In Fig. 15(a), an exemplary multi-layer network is portrayed with $N=5$ and
$L=2$, in which the fraction of inactive units is chosen as $p=0.4$. Figure
15(b) depicts the variation in the order parameter $R$ as a function of the
inactivation ratio $p$ for the three different network setups. To be precise,
the first plot (in red) corresponds to the single-layer case, i.e., when
$L=1$. As witnessed, starting from the normalized unit value, $R$
monotonically decreases and eventually drops down to zero near $p\sim 0.74$
implying an aging transition around $p_{c}\sim 0.74$. Next, $R$ is shown (in
green) for the bilayer framework ($L=2$), precisely for case II, while keeping
the interlayer coupling strength $\sigma$ at $\sigma=1.5$. This time, the
critical inactivation ratio $p_{c}\sim 0.94$ for which the aging transition
occurs is much higher than that of the earlier case (single-layer network).
This implies that the robustness of the bilayer network for case II is higher
than that of the single-layer network. Finally, $R$ is plotted (in blue) for
the multilayer network ($L=2$) under case III, for a much longer $\sigma=8$.
In contrast to the previous case, now the aging transition takes place earlier
than the single-layer network. Specifically, the critical ratio $p_{c}\sim
0.72$ of the inactive elements is smaller than that for the single-layer case.
Thus, as compared to the single-layer network formulation, the robustness of
the multilayer network can be higher or lower depending on the functional form
of the interlayer coupling.

###  II.4 Dynamical robustness of long-range connectivity networks

As far as the interactions among the constituents of a complex system are
concerned, we have, thus far, confined ourselves only on the short-range
direct communications. However in networked systems, interactions arise not
only from the direct connections between nodes but also from indirect long-
range communications facilitated by numerous other existing pathways that link
the nodes. Long-range connectivity campa2009statistical ; gupta2010slow ;
joyce2014scaling ; teles2015temperature ; bountis2016mechanical ;
bouchet2010thermodynamics  is omnipresent in complex systems and hence has
recently emerged as one of the flourishing areas of research. In particular,
researchers have studied the presence of long-range interactions in various
networks characterized by a power-law decay. The examples include biological
networks raghavachari1995spatially , Rydberg atoms saffman2010quantum ,
hydrodynamic interaction uchida2011generic , plasmas levin2008collisionless ,
nuclear spins alvarez2015localization , and in climate, called teleconnections
donges2009complex . From synchronization marodi2002synchronization ;
chowdhury2010synchronization ; rakshit2020synchronization ;
rakshit2021neuronal  to chimera state banerjee2016chimera  and oscillation
quenching sathiyadevi2019long , various phenomena have been examined in
networks subject to long-range communications among the nodes.

Here, let us first describe long-range interaction in networks through the
adjacency matrices corresponding to different paths. Let $G=(V,E)$ be a
network consisting of $N$ nodes in which $V$ and $E$ represent the collections
of nodes and links respectively, so that $V=\\{1,2,\dots,N\\}$, and $E\subset
V\times V$ is the set of links. Then we can characterize the network’s
diameter as $D=\max\big{\\{}\mbox{dist}(j,k):j,k=1,2,\dots,N\big{\\}}$. Here
$\mbox{dist}(j,k)$ denotes the distance between the nodes $j$ and $k$.

The $d$-path adjacency matrix ${A}^{[d]}$ can then be written as,

${A}_{jk}^{[d]}=\begin{cases}1,&\text{if}\leavevmode\nobreak\
\mbox{dist}(j,k)=d,\\\ 0,&\text{otherwise}.\end{cases}$

For a clear perception, we choose a small exemplary network with $N=6$ nodes
with the diameter $D=3$, and in Fig. 16, we depict the associated $d$-path
$\leavevmode\nobreak\ (d=1,2,3)$ networks along with the adjacency matrices.
In Fig. 16(a), we display the original given network, i.e., the direct
$1$-path network, together with the associated adjacency matrix ${A}^{[1]}$ in
the lower panel. Similarly, the extracted $2$\- and $3$-path networks are
depicted in Figs 16(b) and 16(c) respectively, with the corresponding
adjacency matrices ${A}^{[2]}$ and ${A}^{[3]}$. Then casting each node by the
dynamics of Stuart-Landau oscillators, the dynamics of the $j$-th node in the
network subject to long-range communication can be described as
majhi2022dynamical ,

$\begin{array}[]{lcl}\dot{z}_{j}=(\alpha_{j}+i\omega-|z_{j}|^{2})z_{j}+\dfrac{1}{N}\sum\limits_{d=1}^{D}\sigma_{d}\sum\limits_{k=1}^{N}{A}^{[d]}_{jk}(z_{k}-z_{j}),j=1,2,\cdots,N,\end{array}$
(37)

In this context, $\sigma_{d}$ represents the interaction intensity between the
$j$-th and $k$-th nodes when the distance between them is $d$, where $D$
denotes the network’s diameter. Thus the strength of interaction between each
pair of nodes essentially depends on the distance between them. We also choose
$\alpha_{j}=a>0$ and $\alpha_{j}=-b<0$ respectively for the active and
inactive set of oscillators in the system.

Figure 16: A network featuring prominent long-distance connections: the
$d$-path networks and their corresponding adjacency matrices are displayed in
the upper and lower rows, respectively: (a) When $d=1$, it refers to the
original network provided, (b) When $d=2$, it represents the network with a
2-path derived from (a). (c) When $d=3$, it indicates the network with a
3-path derived from (a). Reprinted figure with permission from
Ref.majhi2022dynamical

Figure 17:  The order parameter $R$ for the system (37) with respect to the
simultaneous variations of (a) the inactivation ratio $p$ and the interaction
strength $\sigma$ where $\beta=2$ is fixed, (b) the rate of decay $\beta$ and
the inactivation ratio $p$ where $\sigma=20$ is fixed, (c) The critical
inactivation ratio $p_{c}$ as a function of $\beta$ and $\sigma$. The white
curves in (a) and (b) correspond to the the analytical $p_{c}$ (Eq. (42)) .
Reprinted figure with permission from Ref.majhi2022dynamical

We specifically focus on a power-law decay of the coupling strength in
relation to the distance between the nodes, which is best represented by
$\sigma_{d}=\dfrac{\sigma}{d^{\beta}}$. Here, the exponent $\beta$ governs the
rate of decay in the power law. We consider Erdos-Renyi random network
architecture ($G(N,q)$ graph erdHos1960evolution ) as the underlying network
with network size $N=200$ and the connection probability $q=0.05$, and we
choose $a=1$, $b=1$ with $\omega=3$. Figure 17(a) portrays the variation in
the order parameter $R$ as a function of the interaction strength
$\sigma\in[0,30]$ and the inactivation ratio $p$, in which the decay rate
$\beta=2$ is kept fixed. The phase diagram explains how $R$ transits from non-
null to null value and hence aging takes place in the entire range of
$\sigma\in[0,30]$ as the inactivation ratio $p$ increases. The increasing
values of both the interaction strength $\sigma$ and the fraction $p$ lead to
faster aging and hence the network becomes more dynamically vulnerable. On
this phase diagram acquired through numerical means, we depict the critical
value $p_{c}$, which is obtained through analytical methods (discussed later),
and aligns with the numerical findings. This analytical curve thus separates
the aging region from the oscillatory regime. In Fig. 17(b), order parameter
$R$ is plotted with respect to the simultaneous variation of $p$ and the decay
rate $\beta\in[0,4]$, with a fixed coupling strength $\sigma=20$. In contrast
to the effect of the interaction strength $\sigma$ on the dynamical
robustness, increasing values of $\beta$ inhibits aging transition and hence
dynamical robustness increases. The analytically obtained $p_{c}$ values are
plotted over the phase diagram, as before, that fits into the numerical
outcome. Ultimately, in Fig. 17(c), we illustrate the variations in the
critical inactivation ratio $p_{c}$ across the parameter plane
$(\sigma,\beta)$. The phase diagram describes how the robustness of the
network alters due to the simultaneity of the two crucial parameters $\sigma$
and $\beta$. The critical ratio $p_{c}$ remains unity for sufficiently small
interaction strength $\sigma$, irrespective of the value of $\beta$. However,
as $\sigma$ increases aging takes place in the networked system so that
$p_{c}$ decreases. On the contrary, $p_{c}$ increases for increasing values of
$\beta$. This indicates that the system is dynamically vulnerable, whenever
the coupling strength $\sigma$ is high and the long-range exponent $\beta$ is
low.

To establish a comprehensive analytical framework for evaluating the dynamic
robustness of the networked system under consideration, we employ a degree-
weighted mean field approximation. As a result, the system (37) can be
approximated as

$\begin{array}[]{lcl}\dot{z_{j}}=(\alpha+i\omega-|z_{j}|^{2})z_{j}+\dfrac{1}{N}\sum\limits_{d=1}^{D}\sigma_{d}{k_{j}}^{[d]}\Big{[}(1-p)M_{A}(t)+pM_{I}(t)-z_{j}\Big{]}.\\\
\end{array}$ (38)

Here,

$\begin{array}[]{lcl}M_{A}(t)=\dfrac{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{A}}{k_{j}}^{[d]}z_{j}(t)}{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{A}}{k_{j}}^{[d]}},\par\leavevmode\nobreak\ \leavevmode\nobreak\
M_{I}(t)=\dfrac{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{I}}{k_{j}}^{[d]}z_{j}(t)}{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{I}}{k_{j}}^{[d]}}\par\end{array}$ (39)

are the mean-fields with weighted degrees corresponding to the active and
inactive sets of dynamical units, respectively.

where ${k_{j}}^{[d]}(j=1,2,...,N)$ refers to the degree of the $j$-th node
related to the $d$-path network. Assuming the state variables in the form
$z_{j}(t)=r_{j}(t)e^{i(\omega t+\theta)}$, Eq. (38) can be expressed as

$\small{\begin{array}[]{lcl}\dot{r}_{j}=\Big{(}\alpha_{j}-\dfrac{1}{N}\sum\limits_{d=1}^{D}\sigma_{d}{k_{j}}^{[d]}-r_{j}^{2}\Big{)}r_{j}+\dfrac{1}{N}\sum\limits_{d=1}^{D}\sigma_{d}{k_{j}}^{[d]}\Big{[}(1-p)R_{A}(t)+pR_{I}(t)\Big{]},\\\
\end{array}}$ (40)

in which

$\begin{array}[]{lcl}R_{A}(t)=\dfrac{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{A}}{k_{j}}^{[d]}r_{j}(t)}{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{A}}{k_{j}}^{[d]}},\leavevmode\nobreak\ \leavevmode\nobreak\
R_{I}(t)=\dfrac{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{I}}{k_{j}}^{[d]}r_{j}(t)}{\sum\limits_{d=1}^{D}\sigma_{d}\sum_{j\in
S_{I}}{k_{j}}^{[d]}}.\end{array}$ (41)

Then, pursuing a similar procedure as above in Sec. II.2.2, we arrive at the
critical inactivation ratio

$\begin{array}[]{lcl}p_{c}=\dfrac{H(\sigma,a)-1}{H(\sigma,a)-H(\sigma,-b)},\end{array}$
(42)

where $H$ arises from the following equation,

$\begin{array}[]{lcl}H(\sigma,\alpha)=\dfrac{\sum\limits_{d=1}^{D}\sigma_{d}\sum\limits_{j=1}^{N}\Big{[}\dfrac{{k_{j}}^{[d]}L_{j}}{L_{j}-N\alpha_{j}}\Big{]}}{N^{2}\sum\limits_{d=1}^{D}\sigma_{d}s^{[d]}},\end{array}$
(43)

in which
$L_{j}=\sum\limits_{d=1}^{D}\sigma_{d}{k_{j}}^{[d]};\leavevmode\nobreak\
j=1,2,...,N$.

##  III Dynamical robustness of quantum oscillators

So far, we have considered classical Stuart-Landau oscillators to understand
different routes leading to aging transition. We now pose the question: “Can
the aging transition occur in the quantum domain, and if it does, how would it
manifest?” Motivated by this question, Bandyopadhyay et al.
bandyopadhyay2023aging  investigate the occurrence of aging transitions in
quantum systems, emphasizing the differences compared to classical models
Daido_2004_PRL . They examine a globally coupled quantum Stuart-Landau
oscillators, which can be either active or inactive. The classification of
nodes as active or inactive in the quantum context is determined by the
characteristics of the dissipators in the quantum master equation. The quantum
master equation for $N$ globally coupled quantum Stuart-Landau oscillators
with diffusive coupling is expressed as ishibashi2017oscillation ,

$\displaystyle\dot{\rho}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}\left(-i[H,\rho]+\underbrace{G_{j}\mathcal{D}[O_{j}](\rho)}+\kappa\mathcal{D}[{a_{j}}^{2}](\rho)\right)+\frac{V}{N}\sum_{j=1}^{N}\sum_{j^{\prime}=1}^{N}{}^{{}^{\prime}}\mathcal{D}[a_{j}-a_{j^{\prime}}](\rho),$
(44)

where $H=\omega{a_{j}}^{\dagger}a_{j}$. $a_{j}$ and ${a_{j}}^{{\dagger}}$
represent bosonic annihilation and creation operators of the $j$-th
oscillator, respectively. The Lindblad dissipator $\mathcal{D}[\hat{L}]$ have
the form
$\mathcal{D}[\hat{L}](\rho)=\hat{L}\rho\hat{L}^{\dagger}-\frac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\rho\\}$,
where $\hat{L}$ is an operator (we set $\hbar=1$, without any loss of
generality ). The operator $O_{j}$ of the second term (shown under the brace)
is introduced to incorporate the concept of active and inactive elements. As
system approaches the classical limit ($G_{j}>\kappa$), the quantum master
equation becomes equivalent to the classical Stuart-Landau equation by the
relation: $\langle\dot{a}\rangle=\mbox{Tr}(\dot{\rho}a)$. Here
$\sum_{j^{\prime}}{}^{{}^{\prime}}$ indicates that the sum does not include
the condition $j^{\prime}=j$. The concept of active and inactive elements can
be introduced into the quantum master equation through the properties of the
operator associated with the coefficient $G_{j}$ in the Lindblad dissipator as
follows:

$O_{j}=\begin{cases}{a_{j}}^{\dagger}&\text{for active
oscillators,}\\\\[10.0pt] {a_{j}}&\text{for inactive oscillators.}\end{cases}$

$O_{j}={a_{j}}^{\dagger}$ is the dissipator $G_{j}\mathcal{D}[O_{j}](\rho)$ in
Eq. (44) describes a single boson, gain with a rate of $G_{j}$ and it lead to
stable limit cycle for the jth oscillator to make it quantum active element.
When $O_{j}={a_{j}}$, the system experiences a single boson loss at a rate of
$G_{j}$, which show the non-oscillatory or inactive quantum behaviour of jth
oscillator. Similar to the case of the classical system Daido_2004_PRL , the
whole network is divided into two groups one group consists of $N_{a}$ active
elements and the other consists of $N_{i}$ inactive elements and calculate
fraction of inactive node $p=\frac{N_{i}}{N}$ in the network.

In the uncoupled state ($V=0$), the phase space representation of the Wigner
function for active and inactive elements is depicted in Fig.18(a) and (b),
respectively. The ring-shaped Wigner function carmichael2013statistical
indicates a quantum limit cycle at $G_{j}=4,O_{j}=a_{j}$ and a probability
blob at $G_{j}=2,O_{j}=a_{j}{\dagger}$ represents a non-self-oscillatory
element. When dealing with a large number of oscillators, the density matrix
of the many-body system can be approximately factorized as
$\rho\approx\otimes_{j=1}^{N}\rho_{j}$. This approach aligns with the mean-
field approximation, simplifying the master equation Eq. 44 into individual
master equations for each oscillator, which then interact with the mean-field
as follows ishibashi2017oscillation :

$\displaystyle\dot{\rho_{j}}$
$\displaystyle=-i[\omega{a_{j}}^{\dagger}a_{j},\rho_{j}]+G_{j}\mathcal{D}[O_{j}](\rho_{j})+\kappa\mathcal{D}[{a_{j}}^{2}](\rho_{j})$
$\displaystyle+\frac{2V(N-1)}{N}\mathcal{D}[a_{j}](\rho)+V\left(A[{a_{j}}^{\dagger},\rho_{j}]-A^{*}[a_{j},\rho_{j}]\right),$
(45)

where $A$ and $A^{*}$ are defined as follows:
$A=\frac{1}{N}\sum_{j^{\prime}=1}^{{}^{\prime}N}\langle
a_{j^{\prime}}\rangle_{j}$ and
$A^{=}\frac{1}{N}\sum_{j^{\prime}=1}^{{}^{\prime}N}\langle
a_{j^{\prime}}^{\dagger}\rangle_{j}$.

Figure 18:  (a,b) Distribution of Wigner function in phase space for uncoupled
(a) active oscillator ($G_{j}=4$ and $O_{j}={a_{j}}^{\dagger}$) and (b)
inactive oscillator ($G_{j}=2$ and $O_{j}=a$) at fixed parameters are
$\omega=2$ and $\kappa=0.2$. Reprinted figure with permission from Ref.
bandyopadhyay2023aging

.

Eq. III is numerically solved by self-consistent method using QuTiP
johansson2012qutip . We take $G_{j}=4$ for the active elements
$[j\in\\{1,2,...,N_{a}\\}]$, and $G_{j}=2$ for the inactive elements
$[j\in\\{N_{a}+1,...,N\\}]$. In the network, we differentiate between the
oscillatory state and the oscillation-collapsed state by calculating the
average boson number per oscillator:
$Q=\frac{\bar{n}_{mf}(p)}{\bar{n}_{mf}(0)}$, where, $\bar{n}_{mf}(p)$ is the
mean boson number per oscillator for a particular $p$ value.

To examine how the average mean boson number changes as the value of $p$
increases, $Q$ is plotted as function of $p$ for different coupling strengths
V in Fig. 19(a). For $V\leq 2.73$, $Q$ decreases monotonically as $p$
increases. However, once V exceeds approximately 2.73, the rate at which Q
decreases exhibits two distinct phases: initially, Q declines sharply as p
increases, but beyond a critical threshold $p_{cq}$, the rate of decrease
becomes nearly linear. The point on the curve that marks the transition
between the steep decline and the inclined linear region is referred to as the
knee point (for $V=5$ star mark in the Fig.19(a)). These knee points $p_{cq}$
consider as aging transition threshold and its associated order parameter
denoted as $Q_{c}$. As p increases further, the curve progressively approaches
Q=0 because of the growing number of inactive elements in the network. When p
reaches one, all oscillators become inactive, resulting in $Q=0$, which is a
straightforward scenario. The dependence of $p_{cq}$ and $Q_{c}$ on the
coupling strength V is illustrated in Figs. 19(b) and (c), respectively. The
results shows that $p_{cq}$ initially increases as V increases. However, for
strong coupling strength, $p_{cq}$ reaches a plateau and exhibits little
further change. This behavior differs from the classical case, where $p_{c}$
generally decreases with increasing coupling strength Daido_2004_PRL . It’s
notable that $Q_{c}$ diminishes as the coupling strength increases, aligning
with the anticipated effect of stronger coupling promoting aging.

Figure 19: (a) Order parameter $Q$ is plotted against $p$ for various values
of coupling strength $V$. For $V=5$, the star mark indicates the knee point of
the curve. The corresponding critical values of $p$ and $Q$ are indicated as
$p_{cq}$ and $Q_{c}$, respectively. (b-c) The variation of $p_{cq}$ and
$Q_{c}$ with coupling strength $V$ at parameter $\omega=2$ and $\kappa=0.2$.
Reprinted figure with permission from Ref. bandyopadhyay2023aging

.

The key observations that distinguish the quantum aging transition from its
classical counterpart are significant. Unlike in classical systems, where
aging transition is marked by the complete collapse of the network, the
quantum aging transition is characterized by a rapid decrease in the average
mean boson number. Furthermore, the quantum aging process involves two
distinct phases: initially, up to a critical “knee” point, the order parameter
decreases rapidly; beyond this point, the rate of decrease slows down. During
this latter phase, all inactive elements populate the ground state through the
single boson loss process. However, active oscillators never fully reach the
ground state, as their relaxation is governed solely by the two-boson
absorption process. This unique scenario, which has no classical counterpart,
underscores the distinct nature of quantum aging processes.

##  IV Dynamical robustness of biological networks

So far, we confined ourselves to the Stuart-Landau model which is a normal
form of a nonlinear oscillating system near the Hopf bifurcation point.
However, the scenario of loss of collective oscillatory activity in diverse
networked systems due to the failure of their components exhibits numerous
practical applications in specific biological systems as well. It is
critically important to investigate how various conditions affect macroscopic
activity when certain microscopic units fail and lose their self-oscillatory
behavior. This is especially pertinent in biological systems, since collective
functions play an integral part frequently resulting from interactions among
oscillatory units that can gradually deteriorate under pathological
conditions. Notable attempts have been made to further explore dynamical
robustness in terms of aging transition in realistic biological networks
composed of active and inactive elements, such as neuronal ensembles, and
ecological systems. Motivated by these facts, in this section we highlight
important results obtained for basic models in ecological and neuronal
systems.

###  IV.1 Spatial metapopulation networks

In recent past, the aging transition phenomenon has been studied in the
context of metapopulation survivability by Kundu et al. Rakshit_2017_PRE .
Metapopulation dynamics, a concept used to describe the movement of spatially
separated populations of one species in spatial ecologylevins1969some ;
hanski1999metapopulation , has shed light on the long-term dynamics of
structured populations. These studies have revealed that population densities
of a particular species often undergo synchronized fluctuations across
extensive geographic regions ranta1995synchrony . In this framework, a patch
is typically represented as a system of differential equations that displays
oscillatory solutions. Spatially structured metapopulations can be conceived
as a network composed of interconnected oscillators. In this context, nodes
correspond to viable habitat patches, and the links connecting these nodes
signify functional pathways. This conceptual framework enables the examination
of the ecological network’s dynamical robustness, particularly in the context
of predator-prey patches. The mathematical representation of the dynamics
within a single patch is as follows,

$\displaystyle\dot{x}=f(x,y)=\frac{1}{\epsilon}[x(1-x)(x-\theta)-xy],$
$\displaystyle\dot{y}=g(x,y)=xy-dy.$

In this context, the variables $x$ and $y$ represent the normalized prey and
predator population densities, respectively. The parameter $\epsilon\in(0,1]$
signifies the time scale separation between the prey and predator populations,
$\theta\in(0,1)$ represents the Allee threshold, and $d$ stands for the
natural mortality rate of the predator population. The nontrivial fixed point
$(d,(1-d)(d-\theta))$ exists under the condition $\theta<d<1$. This fixed
point is stable when $d>\frac{1+\theta}{2}$, and a supercritical Hopf
bifurcation occurs at $d=\frac{1+\theta}{2}$. When $d\leq\frac{1+\theta}{2}$,
the coexistence of oscillation (stable limit cycle) and a stable extinction
state (0,0) arises based on the initial population density. However, by
further reducing the predator mortality rate, species extinction occurs
through a boundary crisis of the limit cycle attractor.

The prey-predator model involving $N$ patches is mathematically described by
the following equation,

$\dot{\mathbf{X}_{i}}=\mathbf{F(X}_{i})+\mathbf{M}\sum_{j=1}^{N}A_{ij}(\mathbf{X}_{j}-\mathbf{X}_{i}),$
(46)

In this framework, $\mathbf{X}i=(x_{i},y_{i})^{T}$ represents the state
vector, and $\mathbf{F(X}i)=(f(x_{i},y_{i}),g(x_{i},y_{i}))^{T}$ describes the
inherent dynamics of the $i$-th patch. The second term signifies diffusive
coupling, illustrating interactions among species across different patches.
Here, $\mathbf{M}=(\frac{m}{deg(i)},\frac{m}{deg(i)})^{T}$ represents the
dispersal matrix, where $m$ denotes the dispersal rate between patches and the
term $deg(i)$ signifies the number of patches (degree) connected to the $i$-th
patch. Here $A_{ij}$ stands for the adjacency matrix. In this study, an active
patch indicates stable limit cycle oscillations in both populations, with
$d=0.5$ set. Conversely, an inactive patch implies species extinction, with
$d=0.3$.

To investigate the aging transition, we adopt the mathematical framework
introduced by Daido_2004_PRL . This approach involves considering a scenario
where, in a fraction $p$ of patches, species go extinct due to the absence of
dispersal between these patches. The order parameter $R$, which quantifies the
level of dynamical activity in the network, is defined as follows,

$R=\frac{1}{2}(R_{x}+R_{y}),$ (47)

where $R_{x}=\frac{1}{N}\sum_{i=1}^{N}(\langle x_{i,max}\rangle_{t}-\langle
x_{i,min}\rangle_{t})$, $R_{y}=\frac{1}{N}\sum_{i=1}^{N}(\langle
y_{i,max}\rangle_{t}-\langle y_{i,min}\rangle_{t})$, with $\langle...\rangle$
representing the long-time average. $R=0$ indicates the presence of stable
steady states. To distinguish between trivial (extinction state) and non-
trivial steady states, the concept of $\Delta=\Theta(\mathbf{X}_{i}-\delta)$
is introduced, where $\delta$ is a predefined threshold, and $\Theta(x)$ is
the Heaviside step function. Non-zero values of the order parameter $R$
indicate the continued existence of the metapopulation throughout the network,
whereas $R=0$ indicates extinction of this metapopulation.

Figure 20: The normalized order parameter $D=\frac{R(p)}{R(0)}$ for the system
(46) is plotted against the inactivation ratio $p$ in the globally coupled
network for various dispersal values of $m$. The critical ratio $p_{c}$ at
which $D$ reaches zero increases as the dispersal rate decreases. Reprinted
figure with permission from Ref.Rakshit_2017_PRE  Figure 21: A comparison of
the critical inactivation ratio ($p_{c}$) values for the system (46) in
relation to the dispersal rate ($m$) is presented for three distinct network
topologies: global, small-world, and scale-free. The metapopulation comprises
a total of $N=500$ patches, and for the scale-free network, a random
inactivation procedure has been employed. Reprinted figure with permission
from Ref.Rakshit_2017_PRE

We now examine three distinct types of dispersal networks: global, small
world, and scale-free.
i) In the case of an all-to-all network, as illustrated in Fig. 20, the
normalized order parameter $D=\frac{R(p)}{R(0)}$ is plotted against $p$ for
various dispersal rates $m$. The figure reveals that as the dispersal rate
decreases, the dynamical robustness increases until it reaches a value of
unity at a critical threshold of $m=0.03$, below which $p_{c}$ remains at
unity. This implies that a lower dispersal rate supports metapopulation
survivability. Additionally, we explore complex dispersal topologies,
specifically small-world and scale-free networks, among patches, aiming to
understand population revival in inactive patches through dispersal.
ii) The small-world dispersal topology exhibits aging transition behavior that
is qualitatively akin to the global (all-to-all) network.
iii) In case of a scale-free dispersal network we employ three distinct
inactivation strategies for the patches, namely random, targeted hub (highest
degree node), and targeted low-degree nodes. Analyzing the variation of $D$ in
relation to $p$ across the three different dispersal topologies, it becomes
apparent that in the case of random inactivation, the likelihood of
metapopulation persistence is higher compared to targeted inactivation.
Concurrently, it is evident that for targeted inactivation, particularly of
high-degree nodes, the critical value $p_{c}$ is notably lower. Across all
conceivable scenarios, there exists a substantial abundance of species within
the metapopulation until the inactivation ratio $p$ reaches the critical value
$p_{c}$, at which point sudden and explosive extinction occurs. In conclusion,
we undertake a comparison of the aging transition across all three dispersal
topologies. Figure 21 illustrates the variation of the critical inactivation
ratio $p_{c}$ with the dispersal rate $m$, clearly indicating that a small-
world network exhibits the highest ecological robustness. Conversely, in the
case of global dispersal, the likelihood of metapopulation extinction is
notably higher compared to the more intricate complex dispersal networks.

As an extension to this study, multilayer framework being capable of offering
an inherent structural framework to model diverse ecological systems, the
authors in the Ref kundu2021persistence  investigated the persistence in
multilayer ecological network consisting of harvested patches. They considered
small-world dispersal topologies in the layers for modeling the communications
between the prey-predator patches. The significant effect on the global
persistence of species caused by asymmetric intralayer and interlayer
dispersal strengths, along with the unique network topologies within the
layers, is examined in detail.

###  IV.2 Neuronal networks

We here present some significant results for studying neuronal ensembles. In
contrast to the above discussions of dynamical robustness of networked systems
in which the self-oscillatory dynamics of the inactive units is lost via
inverse Hopf bifurcation, we start our discussion by focusing on the
robustness of neuronal populations where the inactive neuronal systems lose
their dynamism through a saddle-node bifurcation on the invariant circle
(SNIC) pazo2006universal . To be precise, here the inactive units exhibit the
regime of class-I excitability. Analogous to the earlier approach, we split
the entire population of size $N$ into two sets $S_{E}$ and $S_{A}$,
comprising of $pN$ excitable and $(N-Np)$ acitive systems, respectively. The
dynamical evolution of the networked system reads as,

$\begin{array}[]{lcl}\dot{\bf x}_{j}={{\bf F}_{j}}({{\bf
x}_{j}})+\dfrac{K}{N}\sum\limits_{k=1}^{N}({\bf x}_{k}-{\bf
x}_{j});\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
j=1,2,\cdots,N,\end{array}$ (48)

in which ${\bf F}_{j}={\bf F}_{A(E)}$ whenever $j\in S_{A(E)}$. We,
specifically, choose the paradigmatic Morris-Lecar models for our analysis so
that system (48) becomes,

$\begin{array}[]{lcl}C\dot{\bf
V}_{j}=g_{L}(-V_{j}-V_{L})-w_{j}g_{K}(V_{j}+V_{K})-g_{Ca}m_{\infty}(V_{j})(V_{j}-V_{Ca})-\phi_{j}(V_{j}-0.2)+\dfrac{K}{N}\sum\limits_{k=1}^{N}(V_{k}-V_{j}),\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\dot{w_{j}}=\lambda(V_{j})\\{w_{\infty}(V_{j})-w_{j}\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ j=1,2,\cdots,N,\end{array}$ (49)

where $m_{\infty}(V_{j})=[1+\mbox{tanh}\\{(V_{j}-v_{1})/v_{2}\\}]/2$,
$\lambda(V_{j})=\lambda_{0}[1+\mbox{cosh}\\{(V_{j}-v_{3})/v_{4}\\}]$ and
$w_{\infty}(V_{j})=[1+\mbox{tanh}\\{(V_{j}-v_{3})/v_{4}\\}]/2$. The parameters
are chosen as $g_{L}=0.5$, $V_{L}=0.4$, $g_{K}=2$, $V_{K}=0.7$,
$g_{Ca}=V_{Ca}=C=1$. $\lambda_{0}=0.33$ and
$v_{1,2,3,4}=(-0.01,0.15,0.10,0.145)$. Moreover, $\phi_{j}=\phi_{A}$ assumes
the value $\phi_{A}>\phi_{*}\sim 0.076$ for an active (self-oscillatory)
individual dynamics whereas $\phi_{j}=\phi_{E}$ follows $\phi_{E}<\phi_{*}$
for an excitable cell. We next define the mean ensemble’s frequency as a
measurement of the global oscillation in the ensemble as follows,

$\begin{array}[]{lcl}\Omega=\dfrac{1}{N}\sum\limits_{k=1}^{N}\Omega_{k},\end{array}$
(50)

with its normalized value $R=\frac{\Omega(p)}{\Omega(0)}$, in which
$\Omega_{k}$ is the mean frequency of the $k$-th dynamical unit. For a network
of $N=500$ Morris-Lecar units, Fig. 22(a) displays the normalized average
frequency of the ensemble (along with the individual frequency $\Omega_{j}$ in
the right panels) as a function of the increasing inactivation ratio $p$, for
different values of the interaction strength $K$. It is reasonably discernible
that for low $K$, the transition to the global quiescence state (whenever)
takes place only when all the units are in the inactive (excitable) regime,
i.e., for $p_{c}=1$. This scenario alters when we consider a coupling strength
higher than $K_{c}\sim 0.144$ for which aging occurs even when not all the
elements are excitable. For high $K$, smooth profiles of $R$ and hence of
aging are observed, whereas intermediate values of $K$ leads to step-like
profiles of $R$. The corresponding step-like profiles of the individual
frequencies $\Omega_{j}$ are clearly visible in Figs. 22(c) and 22(d) for
$K=0.144$ and $K=0.2$, respectively.

Figure 22:  Normalized value of the ensemble’s average frequency $R$ for the
system (49) with respect to the inactivation ratio $p$, for several values of
the interaction strength $K$. The right panels depict the individual
frequencies $\Omega_{j}$ for (b) $K=0.05$, (c) $K=0.144$, (d) $K=0.2$ and (e)
$K=0.3$. Reprinted figure with permission from Ref.pazo2006universal

Figure 23:  The order parameter $R$ for the system (51) with respect to the
simultaneous variation in the electrical coupling strength $K_{el}$ and the
inactivation ratio $p$, for three different values of the chemical synaptic
strengths (a) $K_{ch}=0$, (b) $K_{ch}=0.5$ and (c) $K_{ch}=1.5$. Reprinted
figure with permission from Ref.kundu2019chemical

In contrast to the above study based upon the assumption of the presence of
only linear diffusive electrical coupling, we next examine the robustness of
another neuronal population which consists of neuronal systems interacting
through both diffusive gap junctional and non-linear chemical synaptic
communication kundu2019chemical . We particularly emphasize on the multi-layer
framework of the neuronal ensemble and demonstrate that the chemical synapses
acting through the inter-layer connections are sufficiently potential in
recovering the global oscillation and hence the dynamical rhythmicity of the
network. We now choose Hindmarsh-Rose neuron model in order to cast the nodes
in both layers organized in a framework of bi-layer multiplex network. We
presume that the neurons communicate via electrical coupling within each layer
while the neurons across the layers are connected through chemical synapses.
The intra-layer connectivities are considered to be of small-world topology,
which is evident in case of brain networks. The mathematical description of
the entire multi-layer network can then be given by the following equations
as,

$\begin{array}[]{lcl}\dot{x}_{i,k}=a{x^{2}}_{i,k}-{x^{3}}_{i,k}-y_{i,k}-z_{i,k}+\dfrac{K_{el}}{N}\sum\limits_{j=1,\leavevmode\nobreak\
j\neq i}^{N}A_{ij}(x_{j,k}-x_{i,k})+K_{ch}(v_{s}-x_{i,k})\Gamma(x_{i,l}),\\\
\dot{y}_{i,k}=(a+\alpha){x^{2}}_{i,k}-y_{i,k},\\\
\dot{z}_{i,k}=c(bx_{i,k}-z_{i,k}+e),\end{array}$ (51)

where $N$ is the number of neurons present in each layer with $i=1,2,\dots,N$,
$k,l=1,2$ and $k\neq l$. The parameters $K_{ch}$ and $K_{el}$ account for the
chemical and electrical synaptic strengths, respectively. Finally,
${(A_{ij})}_{N\times N}$ represents the adjacency matrix associated to a
small-world architecture considered for each layer. Besides, the variables
${x}_{i,k}$, ${y}_{i,k}$ and ${z}_{i,k}$ correspond to the membrane
potentials, ion transport across the membrane via the fast and slow channels,
respectively. The synapses are assumed to be excitatory with the reversal
potential $v_{s}>{x}_{i,k}(t)$ for all $t$. Further, the chemical synaptic
function is of the form $\Gamma(x)=\frac{1}{1+e^{-10(x+0.25)}}$. With the
parameters $a=2.2$, $c=0.005$, $e=5$ and $\alpha=1.6$, isolated neuronal
systems display a plateau bursting for $b=9$ and exhibit a stable equipoint
regime whenever $b=4$. Thus, the active units correspond to $b=b_{A}=9$,
whereas the inactive systems are associated to $b=b_{I}=4$. With a similar
set-up of dividing the whole ensemble into two groups of active ($(1-p)N$
number of neuronal systems) and inactive dynamical units ($pN$ number of
neuronal systems), we define the following order parameter for each layer
$\bar{R}_{k}=\sqrt{\langle({\bf X}_{c_{k}}-\langle{\bf
X}_{c_{k}}\rangle)^{2}\rangle}$, where ${\bf
X}_{c_{k}}=\frac{1}{N}\sum\limits_{j=1}{N}({x}_{j,k},{y}_{j,k},{z}_{j,k})$ is
the centroid of the $k$-th layer ($k=1,2$) and $\langle\cdots\rangle$ stands
for long-time average. The order parameter for the entire network then is
defined as $\bar{R}(p)=\frac{1}{2}[{\bar{R}}_{1}(p)+{\bar{R}}_{2}(p)]$ with
the normalized value $R=\bar{R}(p)/\bar{R}(0)$.

Assuming $N=200$ neuronal systems in each layer, and a small-world network
topology (with $p_{sw}=0.05$ and average-degree $\langle k\rangle=50$) in each
layer, we plot the order parameter $R$ as a function of the inactivation ratio
$p$ and the gap-junctional strength $K_{el}$ for different values of the
chemical synaptic strength $K_{ch}$, in Fig. 23. The black regions in each of
the plots correspond to $R=0$ reflecting the aging transition, that is, when
the entire neuronal ensemble loses its dynamism. We start with the no-
multiplexing case (i.e., with $K_{ch}=0$) in Fig. 23(a), and observe that
increasing $p$ leads to an aging transition, depending on the strength
$K_{el}$ of the electrical coupling. Higher the gap-junctional strength, the
earlier the global oscillation of the network vanishes and hence aging
transition takes place. However, as we introduce multiplexing in the ensemble
through the non-zero chemical synaptic strength $K_{ch}=0.5$ (cf. Fig. 23(b)),
we witness a significant improvement in the dynamical robustness of the
networked system. This is realized in the form of a shrinked black region
associated with the aging transition. We further increase the chemical
synaptic strength to $K_{ch}=1.5$ and depict a similar phase diagram in Fig.
23(c). With this higher $K_{ch}$, we encounter narrower black region
demonstrating a further enhancement in the robustness of the system. This is
how the chemical synapses are capable of enhancing the rhythmicity of the
multiplexed neuronal ensemble.

In addition to these studies, the dynamical robustness of neuronal networks is
further examined in terms of the phenomenon of aging transitions in an
Erdős–Rényi network of interacting Rulkov neurons based on network
connectivity, connection strength, and the ratio of inactive neurons
biswas2022ageing . Both noise-free and stochastic networks, with additive
noise affecting coupling strength, are investigated. Both smooth and explosive
aging transitions are witnessed in both noise-free and stochastic networks.
Although, noise is found to mitigate the impact of inactive neurons and reduce
the occurrence of explosive transition in the networked system.

Research in this area includes the illustrious study presented by Barać at al.
barac2023determinants  which explores dynamical robustness in terms of the
collective failures in networks of interacting heterogeneous excitable systems
using the FitzHugh-Nagumo neuronal model. These networks are assumed to
exhibit essential characteristics like broad-scale degree distribution, small-
world feature and high modularity. The proportion of inactive excitable units,
the interaction strength between them, and their proximity to the bifurcation
point all influence the network failure resulting in collective aging
transition. It is further demonstrated that intermediate coupling strengths
prolong global network activity when high-degree nodes are inactivated first.
Additionally, the most effective strategy for inducing collective failure
adepends non-monotonically on coupling strength and the distance from the
bifurcation point to the oscillatory dynamics of the units.

Furthermore, the Ref. liu2022analysis  studied the dynamical robustness of a
multilayer neuronal network with electrical intra-layer interaction and non-
synaptic ephaptic coupling between the layers. Ephaptic coupling arises due to
electromagnetic induction caused by extracellular electric fields faber2018two
, and is a form of non-synaptic interaction among neurons that plays crucial
role in neuronal communication. It is worth mentioning that the inter-layer
ephaptic interaction enhances the dynamical robustness of both individual
layers and the entire network, contrasting with electrical coupling, which
tends to weaken it. The network dynamics with such inter-layer ephaptic
coupling is also validated using an analog circuit built on Multisim.

##  V Enhancement and maintenance of dynamical robustness

In today’s interconnected world, where complex systems govern various aspects
of our lives, ensuring their stability and robustness has become paramount.
Enhancing dynamical robustness in the form of resurrecting oscillatory
activity is a crucial endeavor that safeguards the reliability and safety of
critical system, ranging from power grids and transportation networks to
financial systems and healthcare infrastructure. The oscillatory behavior in
neurons is pivotal for processing neural information and coordinating
processes related to cognitive functions and memory varela2001brainweb ;
schmidt2018network . Neurons thus demonstrate a pronounced inclination to
engage in rhythmic activity both at the individual and collective level
stiefel2016neurons ; singer2018neuronal . Therefore, the interruption in the
oscillatory behavior of neurons can directly impact essential neural
processes. Proper functioning of cardiac and respiratory systems
jalife1998self , as well as physiological processes like cell necrosis within
organs gurtner2007progress  relies on oscillatory dynamics. Power-grid
networks necessitate stable, synchronized rhythmic activity as a requirement
motter2013spontaneous ; gambuzza2017analysis . On the other hand, the crisis
of extinction of species is primarily attributed to factors such as climate
change or excessive utilization of natural resources, significantly impacting
the surrounding ecosystems on a broad scale garcia2018effect ;
may2019stability . In ecological networks, the extinction of patches within
the meta-population can thus result in significant alterations to its overall
sustainability gilarranz2012spatial ; allesina2012stability ;
kundu2017survivability ; kundu2021persistence . The significance and worth of
investigations on the enhancement of robustness is thus evident in both
natural and artificial contexts.

###  V.1 Controlled diffusion method

In the last few years, there have been a number of significant attempts in
reaching an enhancement of the dynamical robustness of complex networked
systems. In order to retrieve dynamism from the state of aging in damaged
dynamical networks, we first describe one of the most efficient and simple
procedure based upon controlled diffusion liu2016enhancing . We unravel that a
tiny deviation from the usual diffusive coupling among the dynamical units of
a network can increase the dynamical robustness of networks quite
comprehensively. Instead of easier aging, strong interaction under such
controlled diffusion can support robustness of the concerned system. Let us