# Impact of Diffusion on synchronization pattern of epidemics in nonidentical
metapopulation networks

Anika Roy*
International Institute of Information Technology
Hyderabad, India - 500032
<EMAIL_ADDRESS>& Ujjwal Shekhar*
International Institute of Information Technology
Hyderabad, India - 500032
<EMAIL_ADDRESS>& Aditi Bose
International Institute of Information Technology
Hyderabad, India - 500032
<EMAIL_ADDRESS>& Subrata Ghosh
International Institute of Information Technology
Hyderabad, India - 500082
<EMAIL_ADDRESS>& Santosh Nannuru
International Institute of Information Technology
Hyderabad, India - 500032
<EMAIL_ADDRESS>& Syamal Kumar Dana
Department of Mathematics, Jadavpur University
Kolkata 700032, India
Division of Dynamics, Lodz University of Technology
Stefanowskiego 1/15, 90-924 Lodz, Poland
<EMAIL_ADDRESS>& Chittaranjan Hens
Centre for CNS and Bioinformatics
International Institute of Information Technology
Hyderabad, India - 500032
<EMAIL_ADDRESS>
###### Abstract

In a prior study, a novel deterministic compartmental model known as the
SEIHRK model was introduced, shedding light on the pivotal role of test kits
as an intervention strategy for mitigating epidemics. Particularly in
heterogeneous networks, it was empirically demonstrated that strategically
distributing a limited number of test kits among nodes with higher degrees
substantially diminishes the outbreak size. The network’s dynamics were
explored under varying values of infection rate. In this research, we expand
upon these findings to investigate the influence of migration on infection
dynamics within distinct communities of the network. Notably, we observe that
nodes equipped with test kits and those without tend to segregate into two
separate clusters when coupling strength is low, but beyond a critical
threshold coupling coefficient, they coalesce into a unified cluster. Building
on this clustering phenomenon, we develop a reduced equation model and
rigorously validate its accuracy through comprehensive simulations. We show
that this property is observed in both complete and random graphs.

11footnotetext: U.S. and A.R. have contributed equally to the work.

_Keywords_ Cluster Synchronisation $\cdot$ SEIHRK model $\cdot$ Coupling
Coefficient

##  1 Introduction

Over the past two decades, there has been a notable increase in various
infectious diseases within human populations. Examples include the Severe
Acute Respiratory Syndrome (SARS) outbreak in 2003 [1, 2], the swine flu
pandemic in 2009 [3, 4, 5], and the ongoing COVID-19 pandemic [6, 7, 8]. They
strain healthcare systems [9, 10], causing a surge in patients and potentially
high mortality rates [11, 12]. Thus, suitable control (pharmaceutical and
nonpharmaceutical) strategies are required to mitigate the severity of the
outbreaks [13, 14, 15]. Against this backdrop, to acquire a better
understanding of the disease transmission process and to assess the efficacy
of intervention methods, researchers use compartmental mathematical models
ranging from stochastic (Markovian or non-Markovian) to deterministic
frameworks [16, 13, 14, 17, 18]. In particular, computational epidemiological
modelling has proven to be a significant tool for predicting COVID virus onset
and spread [19, 20]. Note that, compartmental models, such as the SIR
(Susceptible-Infectious-Recovered) model, are useful for understanding disease
transmission in well-mixed populations. However, real-world interactions
frequently occur in complex networks in which individuals are linked in a more
detailed structure, and they can be effectively mapped with the real-world
scenarios of disease transmission from a small to a large population scale.
Thus, both the modeling techniques and the structural properties of complex
networks are essential for developing effective strategies to control the
spread of epidemics. This interdisciplinary approach allows researchers to
consider the nuances of real-world interactions and devise more realistic and
impactful public health interventions [21, 22, 23, 24, 4, 25, 26]. For
instance, complex networks play a crucial role in determining the optimal
vaccination strategy [15, 27]. Also, several structural properties of complex
networks are used as guidelines for shaping up the outbreak distribution [28]
and may provide suitable paths of disease extinction [29, 30]. Both networks
and underlying disease dynamics are extensively used for predicting the
epidemic propagation patterns [24, 26], and for identifying super-spreading
events [31, 32] among metacommunities. The general dynamics that capture the
states of compartments in a network of metacummunities can be written as
$\bf{{\dot{x}_{i}}={F(x_{i})}+\epsilon\sum_{j}A_{ij}G(x_{i},x_{j})}$, where
the vector ${\bf x_{i}}$ represents the current state of each compartment of
the $i$ th community, and $\epsilon$ determines the strength of the diffusion
(${\bf G(x_{i},x_{j})=x_{j}-x_{i}}$). Here, we investigate the strength of
migration (diffusion constant $\epsilon$) in a network of metacommunies.

The concept of migration highlights the significance of travel in the
transmission of infections. With advancements in aviation and transportation
networks, viruses can now spread more quickly than before [33, 34]. From a
theoretical perspective, these reaction-diffusion network dynamics are highly
beneficial for comprehending the transmission of epidemics via human migration
and predicting the arrival time [4] and velocity of the epidemic wave front
[25]. In this context, the central focus of our research is to examine how
epidemics spread and synchronize within complex networks, with the aim of
revealing their spreading patterns. In nonlinear oscillators having diffusive
coupling, it is already established that after a critical diffusion, the
coupled systems may follow the same trajectory, i.e update themselves
synchronously [35, 36, 37]. This phenomenon is common in a wide range of
coupled systems, from coordinated flight of birds [38] to rhythmic flashes of
fireflies [39], The synchronized beats of the heart [40], firing of neurons in
the brain [41], and organized migration of animals [42] are other examples
where the exchange of signals between units plays a crucial role in shaping
the synchronization dynamics. However, the role of diffusion in the
synchronization pattern of epidemics has not been explored rigorously. Also,
in the presence of non-identical meta communities (assume few communities have
a lower infection rate), no one ever studied clustered synchronization of
disease dynamics, except a recent study indicates that the interplay between
two diseases, seasonal influenza and COVID-19, can affect their behaviors.
Such coupling leads to the synchronization of both infections, with a more
significant impact observed on the dynamics of influenza [43].

Here, to investigate the clustering of communities in terms of disease spread,
we employ a complex six-compartment meta-population epidemiological model that
includes the use of test kits, introduced in [13]. We explore the effect of
diffusion coupling ($\epsilon$) within the meta-population model and derive
the critical coupling strength ($\epsilon_{c}$) of cluster and global
synchronization. The rationale for selecting this model pertains to the
necessity of delineating two distinct categories of nodes, distinguishing
between those possessing diagnostic test kits and those lacking testing
equipment. Based on this, it is discerned that nodes furnished with test kits
and those devoid of test kits exhibit a tendency to segregate into discrete
clusters at low coupling strength. However, upon surpassing the critical
threshold in the coupling strength, all nodes amalgamate into a singular
global cluster.

Figure 1: An illustrative overview of our study depicts a meta-population
model wherein each community shares an identical population size.
Distinguished by colour, some communities possess access to test kits
(depicted in green), while others lack this resource (depicted in red). Within
each community, a subset of individuals are infected. In Case A, the network
exhibits low coupling strength (represented by dotted arrows) between
communities, while Case B portrays a scenario with increased coupling
strength. In Case B, the infection trends of the communities cluster, forming
distinct clusters under the influence of higher coupling strength.

Note that numerous synchronization patterns are conceivable, but the primary
focus of our investigation pertains to the examination of partial or cluster
synchronization [44, 36, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54]. In
ecological contexts, cluster synchronization is used [55, 56, 57, 58] to
investigate how the connectivity of a dispersal network and varying coupling
strengths affect the persistence of species in different patches or habitats.
In separate studies of bursting in neurons [59, 54, 60, 61], researchers
illustrate that a network of heterogeneous de-synchronized neurons, when
subject to diverse memory settings, can ultimately lead to bursting and the
formation of cluster synchronization beyond a certain coupling strength
threshold. By building upon these existing works and contributing new insights
into infection control strategies, our study aims to enhance our understanding
of disease dynamics and inform the development of more effective public health
interventions.

Our research specifically hones in on infection trends within interconnected
communities, taking into account human mobility patterns driven by diffusion-
like migration. Figure 1 illustrates our research motivation. As inter-
community migration strength increases, nodes segregate based on the presence
or absence of test kits. This clustering becomes more pronounced with higher
diffusion rates. Additionally, the migration term’s influence on the behavior
of $I_{\text{max}}$ (peak of the infection curve) becomes apparent as
diffusion strength changes. This influence stems from individuals in a
compartment diffusing from nodes with higher population to those with lower
population. Testkit communities, experiencing rapid hospitalization, have
fewer individuals in the infected compartment, resulting in a higher influx of
infected individuals from non-testkit to testkit communities. Thus, with
increasing diffusion strength, infection peaks for non-testkit nodes decrease
and testkit nodes increase as seen in the figure.

In this paper, we examine the process of cluster formation among nodes
possessing test kits and those lacking test kits. We describe the model and
quantities we work with in Sec. 2. The detailed analysis for the effect of
coupling strength in a complete graph is presented in Sec. 3. We analyse the
peak of infection ($I$max) for nodes with and without test kits, along with
their respective averages, as a function of coupling strength ($\epsilon)$,
for different values of initial fractions of nodes ($p$) equipped with test
kits. We also demonstrated the error of global, cluster (with test kits) and
cluster ( without test kits) synchronization in relation to coupling strength.
In addition to our inter-cluster observations, our analysis unveils intriguing
phenomena, including the emergence of two peaks in infection plots and an
unexpected dip in $I$max due to the inherent characteristics of our model. To
gain deeper insights, we have studied the phase plot of synchronization errors
(global and cluster-wise) in ($\mathcal{R}_{0}$, $\epsilon$) and
($p$,$\epsilon$) plane. Here, $\mathcal{R}_{0}$ represents the basic
reproduction number. Evidently, we construct a reduced equation [54] model and
employ it to compute the thresholds for critical coupling strength (Sec. 4).
Furthermore, we corroborated these findings by applying various network
topologies, including regular graphs(Sec. 3) and Erdős-Rényi network(Sec. 5).
Finally, conclusive remarks are made with the results and key takeaways in
Sec. 6.

##  2 A brief description of the model

We employ an epidemiological model, as detailed in a prior study [13]. This
model extends the standard S-E-I-R framework by introducing a ‘hospitalized’
compartment ($H$) to account for individuals requiring hospitalization due to
the disease. Additionally, it incorporates a novel state variable, $K$,
representing the quantity of available test kits. The motivation behind this
extension lies in recognizing the pivotal role of test kits in mitigating the
disease’s spread. The model comprises five distinct compartments: Susceptible
($S$), Exposed ($E$), Infected ($I$), Recovered ($R$), and Hospitalized ($H$),
each reflecting different stages of disease exposure and progression.
Specifically, individuals transition from the susceptible to the exposed state
through close contact with infected individuals, a process governed by the
infection transmission rate $\beta$. The latent period and recovery rate are
denoted by $\sigma$ and $\gamma$, respectively. Notably, the transition of
infected individuals to the hospitalized compartment depends on the
availability of test kits and is expressed as a function of $K$, denoted as
$\alpha(K)$. We posit that the rate of hospitalization is influenced by the
number of available test kits and is modelled as a linear function,
$\alpha(K)=\left(\alpha_{0}+g_{n}(K)\right)$. Here, $\alpha_{0}$ represents
the baseline hospitalization rate in the absence of any test kit. The
production of test kits is tied to the total infected population, with a
production rate $\zeta$, while the efficacy of these kits diminishes due to
production faults at a rate $\chi$. Furthermore, the distribution of test kits
across communities is determined by $g_{n}(K)=\alpha_{1}\times p_{n}K$, where
$\alpha_{1}$ signifies the effectiveness of the test kit and $p_{n}$ is
calculated as the fraction of test kits being allocated to a particular
community $n$ out of the $l$ communities with access to test-kits, here,
$p_{n}=\frac{d_{n}}{\sum_{i=1}^{l}d_{i}}$ while $d_{n}$ corresponds to the
degree of a specific patch, where $n\in{1,2,..,M}$. $\alpha_{1}$ is chosen
keeping $\alpha(K)$ bounded; e.g., the diagnosis rate remains less than or
equal to 1.

In the context of heterogeneous networks, the model’s dynamics are coupled,
resulting in a system of interconnected equations enumerated as follows:

$\displaystyle\frac{dS_{n}}{dt}$
$\displaystyle=-\beta\left(\frac{S_{n}I_{n}}{N_{n}}\right)+\frac{\epsilon}{d_{n}}\displaystyle\sum_{m=1}^{M}{A_{mn}(S_{m}-S_{n})}$
(1) $\displaystyle\frac{dE_{n}}{dt}$
$\displaystyle=\beta\left(\frac{S_{n}I_{n}}{N_{n}}\right)-\sigma
E_{n}+\frac{\epsilon}{d_{n}}\displaystyle\sum_{m=1}^{M}{A_{mn}(E_{m}-E_{n})}$
(2) $\displaystyle\frac{dI_{n}}{dt}$ $\displaystyle=\sigma
E_{n}-\left(\alpha_{0}+g_{n}(K)\right)I_{n}+\frac{\epsilon}{d_{n}}\displaystyle\sum_{m=1}^{M}{A_{mn}(I_{m}-I_{n})}$
(3) $\displaystyle\frac{dR_{n}}{dt}$ $\displaystyle=\gamma
H_{n}+\frac{\epsilon}{d_{n}}\displaystyle\sum_{m=1}^{M}{A_{mn}(R_{m}-R_{n})}$
(4) $\displaystyle\frac{dH_{n}}{dt}$
$\displaystyle=\left(\alpha_{0}+g_{n}(K)\right)I_{n}-\gamma H_{n}$ (5)
$\displaystyle\frac{dK}{dt}$
$\displaystyle=\zeta\displaystyle\sum_{m=1}^{M}{I_{m}}-\chi K.$ (6)

In this context, the unweighted adjacency matrix element $A_{mn}$ represents
the interconnectivity between $M$ distinct communities, and $\epsilon$ denotes
the strength of migration. The model accounts for migration dynamics through
the diffusive term $\sum_{m=1}^{M}{A_{mn}(X_{m}-X_{n})}$, which governs the
exchange of individuals between these communities. Two critical model
assumptions prevail: first, that the disease under examination does not result
in any fatalities, and second, that the timescale of birth and death processes
greatly exceeds that of the disease outbreak, rendering these demographic
processes negligible and thus disregarded within the model. Before delving
into the impact assessment of the coupling strength, it is also essential to
establish a standard mathematical formulation for quantifying clustering. The
following section outlines our formulation of the Synchronization Error, $E$.
We define synchronization errors at both global and cluster levels to
quantitatively assess clustering [52, 44, 50]. The Global Synchronization
Error is denoted as $E_{\mathrm{global}}$, while the cluster synchronization
error is denoted as $E_{\mathrm{cluster}}$, abbreviated as GSE and CSE
respectively.

This formulation of synchronization error employs a well-established
mathematical concept, asserting that when a series of values are closely
clustered, they are also close to their mean.

$E_{\mathrm{global}}=\Big{\langle}\sqrt{\langle(x_{i}-\bar{x})^{2}\rangle\stackengine{5pt}{}{\\!\scriptstyle\mathcal{N}}{U}{l}{F}{F}{L}}\Big{\rangle}_{t}$
(7)

$E_{\mathrm{cluster}}=\left\langle\sqrt{\Big{(}{\frac{1}{\mathcal{N}_{k}}\sum_{i\in
v_{k}}(x_{i}-\bar{x}_{k})^{2}\Big{)}}}\right\rangle_{t}$ (8)

Here, $x_{i}$ could represent any compartment of the community, however, in
our analysis, $x_{i}$ denotes the number of infected individuals in the
$i^{th}$ community, $\bar{x}$ denotes the mean of the entire dataset,
$\bar{x}_{k}$ is the mean of cluster $k$, and $\mathcal{N}_{k}$ is the size of
cluster $k$.

##  3 Impact of Coupling Strength on a Globally Connected Network

Our investigation begins by examining a compact network comprising 20 fully
interconnected communities, among which only 5 have access to test kits.
Through simulations employing the Ordinary Differential Equations (ODEs)
specified by the SEIHRK model over a span of 3000 time steps, we observe two
levels of synchronization in infection trends across communities. The first
level of synchronization manifests within the test-kit and non-test-kit
communities, demonstrating consistency in their infection trends. The second
level is a global synchronization of infection trends across all communities
at high coupling strengths. Given the inherent coupling of equations governing
all compartments, synchronization in one compartment naturally extends to the
others.

Figure 2: Figures A through I show snapshots of infection trends among
different communities with changing coupling strengths. There are $20$
communities, with $5$ having access to test kits. Each community has a
population of $10000$ people, and $E_{0}$ for each community is randomly
sampled between [1, 750]. The plot parameters used were
$\sigma=0.1,\alpha_{0}=0.01,\ \alpha_{1}=0.0001,\ \gamma=0.07,\ \zeta=0.02,\
\chi=0.1$ and $\beta=0.03$.

Figure 2 depicts the infection trends in each community as the coupling
coefficient magnitude ($\epsilon$) increases. Several observations emerge, and
we will explore and eliminate any other contributing parameters in subsequent
sections of our study. In Fig. 2 (A-C), we observe that communities exhibit a
tendency to draw closer together. The emergence of the second peak can be
explained by the varying dynamics between communities with and without test
kits. Non-test kit communities exhibit a higher susceptibility to infection,
which results in a rapid ascent of the first peak, compelling an increased
production of test kits. In response, the abundant availability of test kits
efficiently suppresses the initial peak within the test kit-equipped
communities.

However, as the first peak recedes within non-test kit communities, a critical
issue unfolds - the scarcity of available test kits. This dearth of resources
leads to the resurgence of infections within the test kit communities,
culminating in the appearance of the second peak. This second peak is not
observed in non-test kit communities as the number of susceptibles remaining
within the population is much less after the big first peak. Notably, as the
coupling coefficient ($\epsilon$) is increased, communities equipped with test
kits also converge, the second peak diminishes, and only the first peak
remains, reflecting the intricate interplay between community interactions and
the availability of test kits.

In Figure 2(D), we see that the second peak for the communities with test kits
has almost completely vanished. The infection trends for communities without
test kits are now closer. Figure 2(E) is when we can observe a synchronisation
of communities with and without test kits amongst themselves. In Fig. 2(F-H),
it can be observed that the two synchronised infection trends start to move
towards each other until finally, in Fig. 2(I), the two trends also
synchronise to form a global cluster.

###  3.1 Studying $E$ vs $\epsilon$ plots and corresponding clustering of
$I_{\text{max}}$

Figure 3: Plots A and C show the spread of $I$max values for different
communities decreasing with increasing $\epsilon$. Plots B and D show the
decrease in Synchronisation Error with increasing $\epsilon$. There are 100
communities each with a population of $10000$ people. The plot parameters used
were $\sigma=0.1,\ \alpha_{0}=0.01\,\alpha_{1}=0.0001,\ \gamma=0.07,\
\zeta=0.02,\ \chi=0.1$ and $\beta=0.03$. $p$ denotes the fraction of
communities having test kits. We can observe the existence of a bicluster in
Plots A and C and confirm that the GSE vanishes only after CSE for both test-
kit and testkit-free communities vanishes in plots B and D.

In Fig. 3(A, C) we report that $I$max for each community vs $\epsilon$ as the
background patch, and the mean of $I$max across all communities is indicated
by the solid lines through the patch. In Fig. 3(B, D), it is discernible that
the $\epsilon$ values associated with errors $E$ falling below a pre-defined
empirical threshold ($1.5\times 10^{-4}$) partition the graph into three
distinct regions. These regions are identified by colour in Fig. 3(B, D): the
Multi-cluster region(characterized by a white backdrop), the Bi-cluster
region(characterized by a light grey backdrop) when both test and non-testkit
cluster errors are below the threshold but not the global error and finally,
the Mono-cluster region(distinguished by a dark grey backdrop) where the
global error is below the threshold as well. With increasing $\epsilon$,
errors for test-kit and without test-kit node clusters vanish after a
threshold epsilon while the Global Synchronisation Error is still non-zero.
Upon further increasing $\epsilon$, we observe complete clustering, and the
Global Synchronisation Error (GSE) also drops to zero. We can see now that the
GSE vanishes only after CSE for both test-kit and testkit-free communities
vanishes establishing the existence of a bicluster region.

Figure 3(A, C) show that for smaller orders of the coupling coefficient
$\epsilon$, the $I$max has a greater spread around its mean value - indicated
by a patch around the line showing the mean value of $I$max across all
communities. We can now see that the peak of infection curves form two
clusters and then merge into one. Employing our formulation for
synchronization error, we will ascertain that synchronization extends beyond
the $I$max values alone, encompassing both the spread and the overall shape of
the curve. Upon looking more closely, we observe that the mean $I$max for the
non-testkit communities falls like a smooth-step function with increasing
$\epsilon$, but the testkit communities first monotonically decrease and then
after reaching a minimum start monotonically increasing, creating a small-dip
in the plot. We can explain this using Fig. 2, where we observe two peaks for
lower $\epsilon$ values in the test-kit communities, with the first peak being
lower. As $\epsilon$ values increase, we observe the second peak reduces and
the first peak increases in magnitude, leading to a particular epsilon value
where the $I$max switches from the second to the first peak. As the first peak
increases in magnitude and the second one decreases, $I$max initially
decreases until the first peak overtakes the second peak to become the
determinant of $I$max and continues increasing. Thus, we observe a dip at this
point of overtake.

From this point on, we investigate the impact of other parameters in our study
and confirm that the observed behaviour is due to $\epsilon$ alone.

###  3.2 Investigating Impact of the fraction of communities with test kits
and reproductive number on clustering

Let’s take the fraction of communities with test kits to be denoted by $p$. We
report the Synchronisation errors(Global and Cluster-wise) as a function of
$\mathcal{R}_{0}$ and $\epsilon$ in Fig.4(A-C). The role of $p$ and $\epsilon$
is shown in Fig.4(D-F). In our analysis, $\mathcal{R}_{0}$ represents the
basic reproduction number, calculated as the ratio of the transmission rate
($\beta$) and the constant $\alpha_{0}$. In both Figure 4(A-C) and (D-F), our
observations reveal a consistent trend wherein the Cluster-wise
Synchronisation Error (CSE) converges to zero before the Global
Synchronisation Error (GSE). This notable pattern persists across various
values of $\mathcal{R}_{0}$ and $p$, indicating a prevalent emergence of the
bicluster regime preceding the monocluster regime. Consequently, we
confidently assert that the manifestation of the bicluster regime is solely
attributed to the influence of the migration strength parameter, $\epsilon$.

Figure 4: Here we take 51 communities, each with a population of $10000$
people. The plot parameters used were $\sigma=0.1,\ \alpha_{0}=0.01,\
\alpha_{1}=0.0001,\ \gamma=0.07,\ \zeta=0.02,\ \chi=0.1$ and $\beta=0.03$. The
2D projection of Synchronisation errors (Global and Cluster-wise) in the
$\mathcal{R}_{0}$ vs $\epsilon$ plane are shown from (A-C) and in the $p$ vs
$\epsilon$ plane are shown from (D-F). The dashed line indicates the contour
line when the CSE vanishes, and the solid line indicates the same for GSE.
Additionally, the dotted line in plot (D) shows the same threshold obtained
via the reduced equation model.

##  4 Reduced Equation model

The observations suggest that test kit and non-test kit communities exhibit
similar trends beyond a specific threshold value of $\epsilon$. Both test-kit
and non-test-kit communities cluster, reducing the number of differential
equations from $5M+1$ to only $11$. This simplification arises from
consolidating the $M$ communities into just two, as beyond a certain
$\epsilon$ value, distinct trends converge to either test-kit or non-test-kit.

###  4.1 Derivation of reduced model equations

To do this, we consolidate all the variables with equations containing a
migration term (refer section 2.1) into a vector $\textbf{X}_{i}$
corresponding to the $i^{th}$ community. Let $y$ and $z$ be the number of
communities with and without test kits, respectively.

$\textbf{X}_{i}=\begin{bmatrix}S_{i}\\\ E_{i}\\\ I_{i}\\\ R_{i}\\\
\end{bmatrix}$

We will start by looking at all $i\in[1,y]$, that is, all the communities
where test kits were distributed. According to the equations in section 2.1,
we can write the derivative of $\textbf{X}_{i}$ as

$\displaystyle\dot{\textbf{X}}_{i}$
$\displaystyle=f(\textbf{X}_{i})+\frac{\epsilon}{N-1}\displaystyle\sum_{j=1}^{N}{(\textbf{X}_{j}-\textbf{X}_{i})}$
(9)
$\displaystyle=f(\textbf{X}_{i})+\frac{\epsilon}{N-1}{\biggl{(}(\textbf{X}_{1}-\textbf{X}_{i})+(\textbf{X}_{2}-\textbf{X}_{i})+(\textbf{X}_{3}-\textbf{X}_{i})+\dots+(\textbf{X}_{N}-\textbf{X}_{i})\biggr{)}}$
(10)

Here, $f(\textbf{X}_{i})$ denotes the non-migration terms for each equations
for compartments present in $X_{i}$

$f(\textbf{X}_{i})=\begin{bmatrix}-\beta\left(\frac{S_{i}I_{i}}{N_{i}}\right)\\\
\beta\left(\frac{S_{i}I_{i}}{N_{i}}\right)-\sigma E_{i}\\\ \sigma
E_{i}-\left(\alpha_{0}+g_{i}(K)\right)I_{i}\\\ \gamma H_{i}\end{bmatrix}$

After entering the Bi-Cluster regime, communities corresponding to $i\in[1,y]$
and communities corresponding to $j\in[y+1,N]$ must have clustered amongst
themselves. This implies that for all $i\in[1,y]$

$\displaystyle\textbf{X}_{1}=\textbf{X}_{2}=\textbf{X}_{3}=\dots=\textbf{X}_{y}\implies\textbf{X}_{i}=\textbf{X}_{y}$
(11)

Similarly for all $j\in[y+1,N]$,

$\displaystyle\textbf{X}_{y+1}=\textbf{X}_{y+2}=\dots=\textbf{X}_{N}\implies\textbf{X}_{j}=\textbf{Y}_{z}$
(12)

Many terms from the summation collapse together using the above equalities,

$\displaystyle\dot{\textbf{X}}_{y}$
$\displaystyle=f(\textbf{X}_{y})+\frac{\epsilon}{N-1}\biggl{(}{\displaystyle\sum_{j=1}^{y}{(\textbf{X}_{j}-\textbf{X}_{y})}}+\displaystyle\sum_{j=y+1}^{N}{(\textbf{X}_{j}-\textbf{X}_{y})}\biggr{)}$
(13)
$\displaystyle=f(\textbf{X}_{y})+\frac{\epsilon}{N-1}\biggl{(}\displaystyle\sum_{j=y+1}^{N}{(\textbf{X}_{j}-\textbf{X}_{y})}\biggr{)}$
(14)

Since $y+z=N$,

$\displaystyle\dot{\textbf{X}}_{y}$
$\displaystyle=f(\textbf{X}_{y})+\frac{\epsilon}{N-1}\biggl{(}{\biggl{(}\displaystyle\sum_{j=y+1}^{N}{\textbf{X}_{j}}\biggr{)}}-z\times\textbf{X}_{y}\biggr{)}$
(15)
$\displaystyle=f(\textbf{X}_{y})+\frac{\epsilon}{N-1}\biggl{(}{\biggl{(}\displaystyle\sum_{j=y+1}^{N}{\textbf{Y}_{z}}\biggr{)}}-z\times\textbf{X}_{y}\biggr{)}$
(16)
$\displaystyle=f(\textbf{X}_{y})+\frac{\epsilon}{N-1}\biggl{(}z\times\textbf{Y}_{z}-z\times\textbf{X}_{y}\biggr{)}$
(17)
$\displaystyle=f(\textbf{X}_{y})+\frac{z\epsilon}{N-1}\biggl{(}\textbf{Y}_{z}-\textbf{X}_{y}\biggr{)}$
(18)

Without loss of generality, we can deduce a similar equation for non-test kit
communities too,

$\displaystyle\dot{\textbf{Y}}_{z}=f(\textbf{Y}_{z})+\frac{y\epsilon}{N-1}\biggl{(}\textbf{X}_{y}-\textbf{Y}_{z}\biggr{)}$
(19)

Along with these 8 equations, we have $2$ for $H_{i}$ and $H_{j}$ and $1$ more
for $K$,

$\displaystyle\dot{H}_{y}=\alpha I_{y}+\gamma H_{y}$ (20)
$\displaystyle\dot{H}_{z}=\alpha I_{z}+\gamma H_{z}$ (21)

Here, $H_{y}$ refers to the Hospitalization compartment for test kit
communities and $H_{z}$ refers to the Hospitalization compartment for non-test
kit communities. Similarly for $K$,

$\displaystyle\dot{K}=\zeta(yI_{y}+zI_{z})-\chi K$ (22)

Thus, we reduce $5M+1$ differential equations to 11 differential equations.

On plotting the infection values obtained from with and without test-kit
equations, they follow very closely after the threshold $\epsilon$ when
clustering begins as we expected.

Figure 5: Plotting $I$max obtained from simulation and reduced equation model
together to see areas of synchronisation. There are $100$ communities with a
population of $10,000$ people each. 0.2 and 0.8 fractions of communities get
test kits. Initial conditions were randomised as a uniform distribution
between 1 and 1000. The plot parameters used were $\sigma=0.1,\
\alpha_{0}=0.01,\ \alpha_{1}=0.0001,\ \gamma=0.07,\ \zeta=0.02,\ \chi=0.1$ and
$\beta=0.03$.

Figure 5 shows that the reduced model poorly approximates the value of
$\epsilon$ before entering the Bi-cluster regime. Part of the possible reasons
behind this poor approximation at low coupling strength comes from the assumed
synchronisation before the clusters actually form. However, the trend of the
reduced equation model closely follows the real model at higher values of
$\epsilon$ due to the formation of clusters.

##  5 Impact of the coupling coefficient on an Erdos-Renyi Network

In this investigation, we experimentally validate our conceptual framework on
an Erdos-Renyi Network to ascertain the presence of a bicluster regime.
Consistent with the observations in the preceding section, the clustering
dynamics in this network exhibit the formation of bi-clusters, culminating in
complete synchronisation. A nuanced distinction arises when compared to the
fully connected network; specifically, the two clusters within the formed
bicluster display subtle noise. This variation introduces an additional layer
of complexity in the clustering behaviour, highlighting the network topology’s
influence on the observed synchronisation patterns.

Figure 6: There are $100$ communities with a population of $10000$ people,
each. Test kits were distributed among $40\%$ of the communities. Initial
conditions were randomised as uniform distribution between 1 and 1000. The
plot parameters used were $\sigma=0.1,\ \alpha_{0}=0.01,\ \alpha_{1}=0.0001,\
\gamma=0.07,\ \zeta=0.02,\ \chi=0.1$ and $\beta=0.03$. We observe a
synchronisation of infection trends with an increasing coupling coefficient.
This aligns with what we see in Figure 2, where clusters show some noise near
the peaks. Despite this, a clear synchronization pattern remains evident.

###  5.1 Derivation of reduced model equations

In alignment with the derivation for a complete graph, we pursue a comparable
reduced model formulation. This requires the assumption that the degree in
each node is approximately equal to the average degree across all nodes,
denoted as $\langle d\rangle$.To do this, we again consolidate all the
variables with equations containing a migration term (refer section 2.1) into
a vector $\textbf{X}_{i}$ corresponding to the $i^{th}$ community. Let $y$ and
$z$ be the number of communities with and without test kits, respectively.

$\textbf{X}_{i}=\begin{bmatrix}S_{i}\\\ E_{i}\\\ I_{i}\\\ R_{i}\\\
\end{bmatrix}$

Our derivation starts by looking at all $i\in[1,y]$, that is, all the
communities where test kits were distributed. According to the equations in
section 2.1, we can write the derivative of $\textbf{X}_{i}$ as

$\displaystyle\dot{\textbf{X}}_{i}$
$\displaystyle=f(\textbf{X}_{i})+\frac{\epsilon}{N-1}\displaystyle\sum_{j=1}^{N}{(\textbf{X}_{j}-\textbf{X}_{i})}$
(23)
$\displaystyle=f(\textbf{X}_{i})+\frac{\epsilon}{d_{i}}{\biggl{(}(\textbf{X}_{1}-\textbf{X}_{i})+(\textbf{X}_{2}-\textbf{X}_{i})+(\textbf{X}_{3}-\textbf{X}_{i})+\dots+(\textbf{X}_{N}-\textbf{X}_{i})\biggr{)}}$
(24)

$\displaystyle\dot{\textbf{X}}_{i}$ $\displaystyle\approx
f(\textbf{X}_{i})+\frac{\epsilon}{\langle
d\rangle}{\biggl{(}(\textbf{X}_{1}-\textbf{X}_{i})+(\textbf{X}_{2}-\textbf{X}_{i})+(\textbf{X}_{3}-X_{i})+\dots+(\textbf{X}_{N}-\textbf{X}_{i})\biggr{)}}$
(25)

The main difference lies in how we approximate the degree of each node in the
network to its expected value, denoted as $\langle d\rangle$. This
approximation is based on the Poissonian model applied to the Erdos-Renyi
network. This choice reflects our aim to capture the network’s inherent
structure, emphasising the utility of the Poissonian model in representing the
expected degree distribution.

Figure 7: Plotting $I$max obtained from simulation and reduced equation model
together to see areas of synchronisation. There are $100$ communities with a
population of $10,000$ people each. 0.2 and 0.8 fractions of communities get
test kits. The plots here are for an Erdos-Renyi network. Initial conditions
were randomised as a uniform distribution between 1 and 1000. The plot
parameters used were $\sigma=0.1,\ \alpha_{0}=0.01,\ \alpha_{1}=0.0001,\
\gamma=0.07,\ \zeta=0.02,\ \chi=0.1$ and $\beta=0.03$.

Similar to our earlier derivation for the complete graph, we can condense all
differential equations into the previously discussed 11 equations. This
simplification provides a more concise representation of the system dynamics,
following the approach employed in our analysis of the complete graph.

$\displaystyle\dot{\textbf{X}}_{y}=f(\textbf{X}_{y})+\frac{z\epsilon}{\langle
d\rangle}\biggl{(}\textbf{Y}_{z}-\textbf{X}_{y}\biggr{)}$ (26)
$\displaystyle\dot{\textbf{Y}}_{z}=f(\textbf{Y}_{z})+\frac{y\epsilon}{\langle
d\rangle}\biggl{(}\textbf{X}_{y}-\textbf{Y}_{z}\biggr{)}$ (27)
$\displaystyle\dot{H}_{y}=\alpha I_{y}+\gamma H_{y}$ (28)
$\displaystyle\dot{H}_{z}=\alpha I_{z}+\gamma H_{z}$ (29)
$\displaystyle\dot{K}=\zeta(yI_{y}+zI_{z})-\chi K$ (30)

###  5.2 Comparing infection trends with reduced model prediction

The deviation of the behaviour of the reduced model in the bi-cluster regime
becomes increasingly apparent upon inspecting the plots. The reduced model
lags behind the true behaviour of the two clusters. As shown in the derivation
of the reduced model equation, this follows from approximating the degree
distribution to be exactly equal to its expected value.

Figure 8: There are $51$ communities with a population of $10000$ people,
each. 0.4 fraction of communities got test kits. Initial conditions were
randomised as a uniform distribution between 1 and 1000. This bears
resemblance with what was observed in Fig. 3.

As previously observed, the synchronisation error for individual clusters
diminishes before the Global Synchronisation Error reaches zero, aligning with
our anticipated pattern. Consequently, we assert the definitive presence of a
bi-cluster phenomenon in the Erdos-Renyi graph. In Figure 8, the bi-cluster
regime manifests over a more confined range of $\epsilon$ than noted in the
case of a complete graph.

##  6 Conclusion

The efficacy of the test-kit intervention strategy is evident, particularly
when deploying test-kits in nodes with high degrees. Our study demonstrates
that this approach, coupled with a substantial coupling strength, effectively
mitigates and synchronises infection trends within the population. We reveal a
critical range of coupling coefficients wherein infection trends bifurcate
into two distinct clusters based on the presence or absence of test-kit
distribution.

Furthermore, our investigation establishes that synchronisation is solely
contingent on coupling strength. We analyse this synchronisation across
various fractions of distributed test-kits and different transmission rates,
emphasising its dependency on coupling strength. The observed synchronisation
behaviour facilitates the reduction of numerous differential equations to a
mere 11, offering a simplified model applicable to both complete graphs and
Erdos-Renyi formulations. This approximation proves effective, especially as
the coupling coefficient enters the bi-cluster regime. Reducing the number of
differential equations streamlines the model and enhances its practicality in
the face of increasing complexity.

Future directions of work could involve carrying out the outlined analysis for
scale-free networks. Additionally, delving into the captivating landscape of
multilayer graphs, where nodes hold dual opinions, promises to offer a richer
understanding of how infections propagate across diverse layers of influence.
These intriguing directions not only expand the breadth of our research but
also hold the potential to unveil novel insights into the complex dynamics of
contagion within evolving network landscapes.

## References

  * [1] Roy M Anderson et al.  “Epidemiology, transmission dynamics and control of SARS: the 2002–2003 epidemic”  In _Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences_ 359.1447 The Royal Society, 2004, pp. 1091–1105
  * [2] Chris Dye and Nigel Gay  “Modeling the SARS epidemic”  In _Science_ 300.5627 American Association for the Advancement of Science, 2003, pp. 1884–1885
  * [3] Brian J Coburn, Bradley G Wagner and Sally Blower  “Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1)”  In _BMC medicine_ 7.1 BioMed Central, 2009, pp. 1–8
  * [4] Dirk Brockmann and Dirk Helbing  “The hidden geometry of complex, network-driven contagion phenomena”  In _science_ 342.6164 American Association for the Advancement of Science, 2013, pp. 1337–1342
  * [5] Denise J Jamieson et al.  “H1N1 2009 influenza virus infection during pregnancy in the USA”  In _The Lancet_ 374.9688 Elsevier, 2009, pp. 451–458
  * [6] Benjamin F Maier and Dirk Brockmann  “Effective containment explains subexponential growth in recent confirmed COVID-19 cases in China”  In _Science_ 368.6492 American Association for the Advancement of Science, 2020, pp. 742–746
  * [7] Marco Ciotti et al.  “The COVID-19 pandemic”  In _Critical reviews in clinical laboratory sciences_ 57.6 Taylor & Francis, 2020, pp. 365–388
  * [8] Duduzile Ndwandwe and Charles S Wiysonge  “COVID-19 vaccines”  In _Current opinion in immunology_ 71 Elsevier, 2021, pp. 111–116
  * [9] Brett CA Stekelenburg, Harald De Cauwer, Dennis G Barten and Luc J Mortelmans  “Attacks on health care workers in historical pandemics and COVID-19”  In _Disaster medicine and public health preparedness_ 17 Cambridge University Press, 2023, pp. e309
  * [10] Emanuele Preti et al.  “The psychological impact of epidemic and pandemic outbreaks on healthcare workers: rapid review of the evidence”  In _Current psychiatry reports_ 22 Springer, 2020, pp. 1–22
  * [11] Levent Akin and Mustafa Gökhan Gözel  “Understanding dynamics of pandemics”  In _Turkish journal of medical sciences_ 50.9, 2020, pp. 515–519
  * [12] Lone Simonsen et al.  “Pandemic versus epidemic influenza mortality: a pattern of changing age distribution”  In _Journal of infectious diseases_ 178.1 The University of Chicago Press, 1998, pp. 53–60
  * [13] Subrata Ghosh et al.  “Optimal test-kit-based intervention strategy of epidemic spreading in heterogeneous complex networks”  In _Chaos_ , 2021
  * [14] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem and Alessandro Vespignani  “Epidemic processes in complex networks”  In _Rev. Mod. Phys._ 87 American Physical Society, 2015, pp. 925–979  DOI: 10.1103/RevModPhys.87.925
  * [15] Zhen Wang et al.  “Statistical physics of vaccination”  In _Physics Reports_ 664 Elsevier, 2016, pp. 1–113
  * [16] Alex Arenas et al.  “Modeling the Spatiotemporal Epidemic Spreading of COVID-19 and the Impact of Mobility and Social Distancing Interventions”  In _Phys. Rev. X_ 10 American Physical Society, 2020, pp. 041055  DOI: 10.1103/PhysRevX.10.041055
  * [17] Linda J.. Allen  “An Introduction to Stochastic Epidemic Models”  In _Mathematical Epidemiology_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, pp. 81–130  DOI: 10.1007/978-3-540-78911-6_3
  * [18] Zhen Wang et al.  “Statistical physics of vaccination” Statistical physics of vaccination  In _Physics Reports_ 664, 2016, pp. 1–113  DOI: https://doi.org/10.1016/j.physrep.2016.10.006
  * [19] Alessandro Vespignani et al.  “Modelling COVID-19”  In _Nature Reviews Physics_ 2.6, 2020, pp. 279–281  DOI: 10.1038/s42254-020-0178-4
  * [20] Gabriel M Leung Joseph T Wu 1  “Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study”  In _Lancet (London, England)_ , 2019
  * [21] Romualdo Pastor-Satorras and Alessandro Vespignani  “Epidemic spreading in scale-free networks”  In _Physical review letters_ 86.14 APS, 2001, pp. 3200
  * [22] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem and Alessandro Vespignani  “Epidemic processes in complex networks”  In _Reviews of modern physics_ 87.3 APS, 2015, pp. 925
  * [23] Romualdo Pastor-Satorras and Alessandro Vespignani  “Epidemic dynamics in finite size scale-free networks”  In _Physical Review E_ 65.3 APS, 2002, pp. 035108
  * [24] Chittaranjan Hens et al.  “Spatiotemporal signal propagation in complex networks”  In _Nature Physics_ 15.4 Nature Publishing Group UK London, 2019, pp. 403–412
  * [25] Vitaly Belik, Theo Geisel and Dirk Brockmann  “Natural human mobility patterns and spatial spread of infectious diseases”  In _Physical Review X_ 1.1 APS, 2011, pp. 011001
  * [26] Flavio Iannelli et al.  “Effective distances for epidemics spreading on complex networks”  In _Physical Review E_ 95.1 APS, 2017, pp. 012313
  * [27] Zhen Wang, Yamir Moreno, Stefano Boccaletti and Matjaž Perc  “Vaccination and epidemics in networked populations—an introduction”  In _Chaos, Solitons & Fractals_ 103 Elsevier, 2017, pp. 177–183
  * [28] Jason Hindes, Michael Assaf and Ira B Schwartz  “Outbreak size distribution in stochastic epidemic models”  In _Physical Review Letters_ 128.7 APS, 2022, pp. 078301
  * [29] Petter Holme  “Extinction times of epidemic outbreaks in networks”  In _PloS one_ 8.12 Public Library of Science San Francisco, USA, 2013, pp. e84429
  * [30] Jason Hindes and Ira B Schwartz  “Epidemic extinction and control in heterogeneous networks”  In _Physical review letters_ 117.2 APS, 2016, pp. 028302
  * [31] Alex James, Jonathan W Pitchford and Michael J Plank  “An event-based model of superspreading in epidemics”  In _Proceedings of the Royal Society B: Biological Sciences_ 274.1610 The Royal Society London, 2007, pp. 741–747
  * [32] Bjarke Frost Nielsen, Lone Simonsen and Kim Sneppen  “COVID-19 superspreading suggests mitigation by social network modulation”  In _Physical Review Letters_ 126.11 APS, 2021, pp. 118301
  * [33] Muylaert Renata L. John Reju Sam and Hayman David T. S.  “High connectivity and human movement limits the impact of travel time on infectious disease transmission”  In _Proceedings of the National Academy of Sciences of the United States of America_ , 2006
  * [34] David J Rogers Andrew J Tatem 1  “Global traffic and disease vector dispersal”  In _J. R. Soc. Interface_ , 2024
  * [35] Stefano Boccaletti et al.  “The synchronization of chaotic systems”  In _Physics reports_ 366.1-2 Elsevier, 2002, pp. 1–101
  * [36] Alex Arenas et al.  “Synchronization in complex networks”  In _Physics reports_ 469.3 Elsevier, 2008, pp. 93–153
  * [37] Marcin Kapitaniak et al.  “Synchronization of clocks”  In _Physics Reports_ 517.1-2 Elsevier, 2012, pp. 1–69
  * [38] Máté Nagy et al.  “Synchronization, coordination and collective sensing during thermalling flight of freely migrating white storks”  In _Philosophical Transactions of the Royal Society B: Biological Sciences_ 373.1746 The Royal Society, 2018, pp. 20170011
  * [39] Gonzalo Marcelo Ramírez-Ávila, Jürgen Kurths, Stéphanie Depickere and Jean-Louis Deneubourg  “Modeling fireflies synchronization”  In _A mathematical modeling approach from nonlinear dynamics to complex systems_ Springer, 2019, pp. 131–156
  * [40] Plamen Ch Ivanov, Qianli DY Ma and Ronny P Bartsch  “Maternal–fetal heartbeat phase synchronization”  In _Proceedings of the National Academy of Sciences_ 106.33 National Acad Sciences, 2009, pp. 13641–13642
  * [41] John J Hopfield and Andreas V Herz  “Rapid local synchronization of action potentials: toward computation with coupled integrate-and-fire neurons.”  In _Proceedings of the National Academy of Sciences_ 92.15 National Acad Sciences, 1995, pp. 6655–6662
  * [42] Julien Cote et al.  “Behavioural synchronization of large-scale animal movements–disperse alone, but migrate together?”  In _Biological reviews_ 92.3 Wiley Online Library, 2017, pp. 1275–1296
  * [43] Jorge P. Rodríguez and Víctor M. Eguíluz  “Coupling between infectious diseases leads to synchronization of their dynamics”  In _Chaos_ , 2023
  * [44] Pitambar Khanra et al.  “Identifying symmetries and predicting cluster synchronization in complex networks”  In _Chaos, Solitons & Fractals_ 155 Elsevier, 2022, pp. 111703
  * [45] Nir Lahav et al.  “Topological synchronization of chaotic systems”  In _Scientific reports_ 12.1 Nature Publishing Group UK London, 2022, pp. 2508
  * [46] Pitambar Khanra, Prosenjit Kundu, Chittaranjan Hens and Pinaki Pal  “Explosive synchronization in phase-frustrated multiplex networks”  In _Physical Review E_ 98.5 APS, 2018, pp. 052315
  * [47] Prosenjit Kundu, Pitambar Khanra, Chittaranjan Hens and Pinaki Pal  “Transition to synchrony in degree-frequency correlated Sakaguchi-Kuramoto model”  In _Physical Review E_ 96.5 APS, 2017, pp. 052216
  * [48] Arindam Mishra, Suman Saha and Syamal K Dana  “Chimeras in globally coupled oscillators: A review”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 33.9 AIP Publishing, 2023
  * [49] Louis M Pecora et al.  “Cluster synchronization and isolated desynchronization in complex networks with symmetries”  In _Nature communications_ 5.1 Nature Publishing Group UK London, 2014, pp. 4079
  * [50] Francesco Sorrentino et al.  “Complete characterization of the stability of cluster synchronization in complex dynamical networks”  In _Science advances_ 2.4 American Association for the Advancement of Science, 2016, pp. e1501737
  * [51] Sarika Jalan and Aradhana Singh  “Cluster synchronization in multiplex networks”  In _Europhysics Letters_ 113.3 IOP Publishing, 2016, pp. 30002
  * [52] P. Khanra et al.  “Endowing networks with desired symmetries and modular behavior”  In _Phys. Rev. E_ 108 American Physical Society, 2023, pp. 054309  DOI: 10.1103/PhysRevE.108.054309
  * [53] B Rakshit, AJ Kartha and C Hens  “Predicting aging transition using Echo state network.”  In _Chaos (Woodbury, NY)_ 33.8, 2023, pp. 081102–081102
  * [54] Subrata Ghosh et al.  “Emergence of mixed mode oscillations in random networks of diverse excitable neurons: the role of neighbors and electrical coupling”  In _Frontiers in Computational Neuroscience_ 14 Frontiers Media SA, 2020, pp. 49
  * [55] Samali Ghosh et al.  “Dispersal induced persistent dynamics in interacting ecological patches under randomness of network topology”, 2023
  * [56] Alan Hastings Matthew D.  “Strong effect of dispersal network structure on ecological dynamics”  In _Nature_ , 2008
  * [57] Huawei Fan et al.  “Synchronization within synchronization: transients and intermittency in ecological networks”  In _National Science Review_ , 2020
  * [58] S Saha et al.  “Infection spreading and recovery in a square lattice.”  In _Physical review. E_ 102.5-1, 2020, pp. 052307–052307
  * [59] Sanjeev Kumar Sharma et al.  “Emergence of bursting in a network of memory dependent excitable and spiking leech-heart neurons”  In _J. R. Soc. Interface_ , 2020
  * [60] Suman Saha et al.  “Predicting bursting in a complete graph of mixed population through reservoir computing”  In _Phys Rev Research_ 2 American Physical Society, pp. 033338
  * [61] Chittaranjan Hens, Pinaki Pal and Syamal K Dana  “Bursting dynamics in a population of oscillatory and excitable Josephson junctions”  In _Physical Review E_ 92.2 APS, 2015, pp. 022915