# Controlling species densities in structurally perturbed intransitive cycles
with higher-order interactions

Sourin Chatterjee∗ Department of Mathematics and Statistics, Indian Institute
of Science Education and Research, Kolkata, West Bengal 741246, India
Sayantan Nag Chowdhury∗<EMAIL_ADDRESS>Physics and Applied
Mathematics Unit, Indian Statistical Institute, 203 B. T. Road,
Kolkata-700108, India  Technology Innovation Hub (TIH), IDEAS (Institute of
Data Engineering Analytics and Science Foundation), @ Indian Statistical
Institute, 203 B. T. Road, Kolkata-700108, India  Dibakar Ghosh  Physics and
Applied Mathematics Unit, Indian Statistical Institute, 203 B. T. Road,
Kolkata-700108, India  Chittaranjan Hens<EMAIL_ADDRESS>Physics
and Applied Mathematics Unit, Indian Statistical Institute, 203 B. T. Road,
Kolkata-700108, India  International Institute of Information Technology,
Gachibowli, Hyderabad-500032, India

###### Abstract

111∗ Both the authors have contributed equally to the work.

The persistence of biodiversity of species is a challenging proposition in
ecological communities in the face of Darwinian selection. The present article
investigates beyond the pairwise competitive interactions and provides a novel
perspective for understanding the influence of higher-order interactions on
the evolution of social phenotypes. Our simple model yields a prosperous
outlook to demonstrate the impact of perturbations on intransitive competitive
higher-order interactions. Using a mathematical technique, we show how alone
the perturbed interaction network can quickly determine the coexistence
equilibrium of competing species instead of solving a large system of ordinary
differential equations. It is possible to split the system into multiple
feasible cluster states depending on the number of perturbations. Our analysis
also reveals the ratio between the unperturbed and perturbed species is
inversely proportional to the amount of employed perturbation. Our results
suggest that nonlinear dynamical systems and interaction topologies can be
interplayed to comprehend species’ coexistence under adverse conditions.
Particularly our findings signify that less competition between two species
increases their abundance and outperforms others.

> Recent trends in ecological research suggest the pairwise interaction among
> competitors may not readily infer the mechanism of tremendous biodiversity
> of species. Due to this difficulty, ecologists often resort to the structure
> of diverse competitive higher-order networks. Recently, Grilli et al. Grilli
> _et al._ (2017) shed light on the stability of biodiversity in higher-order
> networks, while pairwise competition allows only the neutral cycling around
> the stationary points. We study the influence of structural perturbations on
> the intransitive cycle of higher-order species interactions and attempt to
> capture the unexplained complexities of competitive ecosystems. Our
> analytical treatment provides a straightforward method to determine the
> globally stable stationary point without solving a large number of ordinary
> differential equations with higher-order terms. A tiny structural
> perturbation among the competitors in the higher-order structure leads to
> the cluster synchrony, where the unperturbed model settles into the complete
> synchronized death state. Our study suggests that two species can easily
> enhance their respective densities in comparison to others by reducing the
> pairwise winning probability from the stronger species. Our simple
> mathematical framework with the presence of intransitive competition and
> higher-order interactions advances our understanding of biodiversity
> maintenance among competitive communities.

##  I Introduction

Understanding how species interact, compete, and try to push each other to
extinction while preserving bio-diversity is a long-standing problem in
ecology. As a result, several theoretical and statistical approaches to
replicating biodiversity have been developed Lotka (1920); Volterra (1926);
May and Leonard (1975); Nag Chowdhury  _et al._ (2021a); Hubbell (2011); May
(1974); Nag Chowdhury  _et al._ (2021b); May (2019); Chesson (2000); Nag
Chowdhury  _et al._ (2020); Gao  _et al._ (2016); Hens  _et al._ (2019); Roy
_et al._ (2022).

One of the classical approaches is Niche theory, which predicts that in the
same niche, when compared to a strong competitor, the weak one is removed due
to the competitive exclusion principle, but large species-coexistence is
maintained due to habitat heterogeneity and appropriate adaptation and
exploitation of niche for each species Chase and Leibold (2009); Soberón and
Nakamura (2009); Tilman (2004). In contrast, in neutral theory, all species
are ecologically equivalent for a given trophic level, i.e., all have an equal
probability of colonizing in open space Bell (2000); Brokaw and Busing (2000),
and bio-diversity emerges as a stochastic balance between speciation and
extinction Hubbell (2011). According to neutral theory, there will be
overlapping niches. However, both classical theories have limitations as well
as accomplishments Wennekes  _et al._ (2012); Chisholm and Pacala (2010);
Gravel  _et al._ (2006); McGill (2003). For example, ’neutral drift’ is
incompatible with the stability of forest diversity Clark and McLachlan
(2003). Also, the paradox in the plankton world, where the number of
coexisting species exceeds the number of limiting resources Huisman and
Weissing (1999), cannot be explained by competition trade-off models.

Figure 1: Schematic representation of the impact of perturbation: The red
circle denotes species, and its size indicates the density of that species.
The two-sided arrow represents the interaction between two species, and its
width contemplates the winning probability of the species towards which it is
pointed. The self-loop depicts the intra-species interaction with a $50\%$
winning probability. In (a), every species has the same density. In (b), we
increase the winning likelihood of species 1 (Stronger species in the
interaction) against species 2, the density of both species decreases (size of
the circles are decreased compared to species 3,4, and 5). In (c), we increase
the winning probability of species 2 (Weaker species in the interaction)
against species 1, the density of both species increases (size of the circles
are now enhanced compared to others).

Motivated by the shortcomings of classical theories, people are looking for
alternative frameworks, particularly those based on game theory. One of them
is the widely studied intransitive cyclic competition He  _et al._ (2011);
Szabó and Fath (2007); Berr  _et al._ (2009); Sinervo and Lively (1996); Kerr
_et al._ (2002); Laird and Schamp (2006); Hibbing  _et al._ (2010); Gilpin
(1975), in which bio-diversity is maintained within a community despite strong
competition. Even in the presence of hierarchy, this type of formalism allows
a large number of species to converge, particularly into stable cycles or, in
some cases, into a stable equilibrium. This phenomenon is known as the
’stabilizing effect of intransitivities’ Allesina and Levine (2011); Gilpin
(1975). The rock-paper-scissors (RPS) type game Allesina and Levine (2011);
Gilpin (1975); Bhattacharyya  _et al._ (2020); Szolnoki  _et al._ (2014) is a
classic example of intransitivity: rock beats scissor, scissor beats paper,
but paper outperforms rock Bhattacharyya  _et al._ (2020); Szolnoki  _et al._
(2014); Reichenbach  _et al._ (2007); Islam  _et al._ (2022). Coexistence is
still possible when a large number of species participate in different
variations of such game theoretical models Park  _et al._ (2017); Allesina and
Levine (2011); Shi  _et al._ (2010). Empirical evidence suggests that
intransitive competition (RPS competition loop) appears to be common in
bacterial strains Kirkup and Riley (2004); Kerr  _et al._ (2002),
phytoplankton communities Soliveres  _et al._ (2015), Parasite-grass-forb
Cameron  _et al._ (2009) or in general consumer-resource problems Huisman and
Weissing (1999). However, such a game theoretic approach does not always
guarantee the stability of feasible equilibrium points Gilpin (1975); Allesina
and Levine (2011), but rather provides a stable cycle in the majority of
cases. Weak perturbation, on the other hand, destabilizes the desired
equilibrium points. However, theoretically one can show, beyond a pairwise
interaction, higher order intransitive mixing in ecological communities can
avoid such drawbacks by allowing the species to maintain stabilization Grilli
_et al._ (2017); Battiston  _et al._ (2020); Letten and Stouffer (2019);
Mayfield and Stouffer (2017); Levine  _et al._ (2017); Bairey  _et al._ (2016)
by settling into a robust equilibrium point. A species is affected by the
interaction of three or an even greater number of species in higher order
interactions Abrams (1983); Lambiotte  _et al._ (2019). It should be noted
that the role of intransitivity and intra as well as inter specific higher
order interactions are tested in density dependent cyclic dynamics Stouffer
_et al._ (2018), particularly in plant dynamics, using data from four co-
occurring annual plant species in SW Western Australia. In contrast, habitat
loss has the potential to reverse the effect of HOI Li  _et al._ (2020). Singh
et al.  Singh and Baruah (2021) demonstrated that the negative strength of
interspecific HOI determines the fate of the bio-diversity.

Here we revisit the classical higher order model developed by Grilli et al
Grilli  _et al._ (2017). The question here is, under what conditions can a
species or a group of species gain or lose abundance strength (in comparison
to others)? Our research provides a suitable and useful tweaking of the
intransitive matrix to provide the species’ target abundance states. The
abundance here denotes the stable equilibrium point. In particular, we focus
on the role of perturbation in the elements of the matrix $A$, which encodes
information about interactions between two species. If the system has $N$
species/seedlings in a forest, then $A$ is an $N\times N$ matrix, where
$A_{ij}$ describes the probability of winning for the $i$-th species over the
$j$-th species, i.e., the ability to fill the gap in a canopy among two
seedlings. Hence, $A_{ij}+A_{ji}=1$ $\forall$ $i,j$. Thus, the diagonal
entries of $A$ will be $0.5$. If all the elements of $A$ are equal, i.e.,
$A_{ij}=0.5$, it represents the neutral model, where as $0-1$ structure
provides a tournament matrix and coexistence is possible in the presence of
intransitive cycles of competitive dominance Allesina and Levine (2011);
Grilli  _et al._ (2017). In this paper, we look at the effect of perturbation
in a higher order interaction (HOI) where the intransitive cycle is captured
by a circulant matrix

$A=\rm{circ}\big{(}0.5,\alpha,1-\alpha,\cdots,\cdots,\cdots).$ (1)

Here the circulant matrix generated by a single element $\alpha\in[0,1]$. If
the square matrix $A$’s order $N$ is odd, the last entry will be $1-\alpha$.
However, if $N$ is an even number, this last entry will be $\alpha$. In the
presence of HOI, for $N$ species, this circulant matrix steers the whole
system to a globally stable equilibrium point $\frac{1}{N}$ (Sec. (III)),
where $N$ is odd. The choice of such a circulant matrix provides an utterly
ecological community with equal strength of interaction among two species,
similar up to the direction. Thus, initially, we do not yield any extra
privilege to any species in the unperturbed system. Furthermore, we try to
inspect the effect of feasible perturbation in a specific pair of species.

Figure 2: Oscillation quenching and global convergence: The system of
differential equations (4) with the interaction matrix $A$ given in Eq. (1)
converges to the stationary point
$\bigg{(}\dfrac{1}{3},\dfrac{1}{3},\dfrac{1}{3}\bigg{)}$. The initial
condition $(0.45,0.41,0.14)$ is chosen randomly within the
$(0,1)\times(0,1)\times(0,1)$ maintaining the constraint
$x_{1}(0)+x_{2}(0)+x_{3}(0)=1.0$. Here, $\alpha=0.75$.

We investigate the effect of localized perturbation on the identical
equilibrium points. Here local perturbation signifies that if $\alpha$ value
between two species are slightly changed, the rest of the interaction remains
unchanged/unperturbed. We have shown analytically that the density of both
perturbed species maintains a relation in the form of a rectangular hyperbola
with respect to the applicable perturbation’s strength for fixed $\alpha$.
Thus, the respective densities of the perturbed species are either increased
or decreased depending on the strength of the perturbation. Specifically, the
species density of each of those perturbed species increases if their internal
competition reduces, i.e., they try to reach neutral competition
($||A_{ij}-A_{ji}||\to 0$). In contrast, the density of both species decreases
if they increase the strength of competition. A schematic figure is drawn to
portray this scenario in Fig. (1) where in the subfigure (a), all species are
unperturbed, maintaining the original strength following the interaction
matrix $A$. The slight increment of winning probability of species 1 against
species 2, keeping the other species unperturbed, leads to a reduction in the
respective densities of the perturbed species (See Fig. (1) (b)). A reverse
scenario of enhancing the winning probability of species 2 against species 1
is contemplated in Fig. (1) (c), where interestingly, the density of both
species increases compared to the unperturbed species.

We have also found the two-clustered equilibrium points of perturbed (say $K$
species) and unperturbed ($N-K$) species analytically. As per our systematic
approach, the perturbed species’ density varies within the interval $[0.5,1]$.
Similarly, by applying $w$ number of perturbations to $2w$ species among
$N\geq(2w+1)$ species, we can split the unique equilibrium points into $(w+1)$
number of different clusters. Finally, we relax the perturbation scheme and
analytically calculate the inhomogeneous stable equilibrium of $N=3$ species
where all the elements in the upper triangular part of the circular matrix are
perturbed within a specific range. Thus, we are able to demonstrate our
findings in a general set-up, where the interaction matrix $A$ attains all
intermediate values along with the neutral and tournament matrix.
Interestingly, if more than two perturbations have been applied, then the
relative change in density of one perturbed group relative to another is
inversely proportional to its own perturbed strength. Note that perturbations
can be considered evolutions of traits in our system of interacting species.
In the presence of selection pressure, traits can be evolved in such a way
that the traits achieve a significant advantage during any competition. A
similar analogy is already well-established in predator-prey systems Abrams
(2000).

##  II Mathematical Model

We consider a dynamical model describing the temporal evolution of the density
$x_{i}(t)$ of the $i$-th species as follows Grilli  _et al._ (2017)

$\begin{split}{\dot{x}}_{i}=x_{i}\bigg{[}\sum_{j=1}^{N}\sum_{k=1}^{N}\big{(}2A_{ij}A_{ik}+A_{ij}A_{jk}+A_{ik}A_{kj}\big{)}x_{j}x_{k}-\eta_{i}\bigg{]},\\\
i=1,2,\dots,N.\end{split}$ (2)

Here, $\eta_{i}$ is the mortality rate of the $i$-th species. The death rate
$\eta_{i}=\eta=1.0$ is chosen to be identical for all the species. Since,
$x_{i}(t)$ represents the density of the $i$-th species, thus the sum of all
species density is one, i.e., $\sum_{i=1}^{N}x_{i}(t)=1$. Instead of the
pairwise interaction, we consider the interactions among the possible
triplets, where the two species compete among each other and the winner play
against the third one. The term within the bracket
$\big{(}2A_{ij}A_{ik}+A_{ij}A_{jk}+A_{ik}A_{kj}\big{)}$ portrays the winning
probability of $i$-th species over other two $j$ and $k$ species under all
possible circumstances. $A_{ij}A_{jk}$ describes the probability of winning of
the $j$-th species over the $k$-th species and, ultimately, the $i$-th species
dominate the $j$-th species. Similarly, the term $A_{ik}A_{kj}$ represents
$k$-th species win over the $j$-th species, and then, the $i$-th species beat
the $k$-th species. The probability of winning of $i$-th species over the
other two $j$ and $k$ species is encapsulated by $A_{ij}A_{ik}$. The factor
$2$ here depicts that $i$-th species can play first with $j$-th species, and
then with $k$-th species or vice-versa. Since, we choose $\eta_{i}=1$ and
$\sum_{i=1}^{N}x_{i}(t)=1$, thus Eq. (2) transforms to

${\dot{x}}_{i}=x_{i}\bigg{[}\sum_{j=1}^{N}\sum_{k=1}^{N}\big{(}2A_{ij}A_{ik}+A_{ij}A_{jk}+A_{ik}A_{kj}-1\big{)}x_{j}x_{k}\bigg{]}.$
(3)

Then, Eq. (3) reduces to the form,

${\dot{x}}_{i}=x_{i}\bigg{[}\sum_{j=1}^{N}\sum_{k=1}^{N}B_{ijk}x_{j}x_{k}\bigg{]},$
(4)

where
$B_{ijk}=2A_{ij}A_{ik}+A_{ij}A_{jk}+A_{ik}A_{kj}-1=2A_{ij}A_{ik}-A_{ji}A_{jk}-A_{ki}A_{kj}$,
since $A_{ij}+A_{ji}=1$ $\forall i,j=1,2,\dots,N$.

Note that there is a thin difference between our model and the RPS game. The
classical RPS game is basically a three-player zero-sum game where cyclic
dominance is observed. We consider a winning probability between two species
as opposed to the classical RPS, wherein in a game, one species beats another
with $100\%$ possibility. Moreover, our model describes the competition
between $N\geq 3$ species rather than three species. One may quickly obtain a
suitable variant of the RPS game by choosing suitable $\alpha$ and $N$.

Figure 3: Effect of localized perturbation: The species densities vary for
different choices of $\beta$ for fixed $\alpha=0.75$. (a) Initially, for
$\beta=0$, the unperturbed system leads to the coexistence of all species with
equal densities. However, the scenario changes with the choices of non-zero
$\beta$. (b) For our specific choice of $\alpha>0.5$, $\beta>0$ reduces the
first and second species’ density (orange and blue respectively) and,
consequently, raises $x_{3}$ (purple). (c) $\beta<0$ improves (raises) $x_{1}$
and $x_{2}$ compared to the unperturbed case, and diminishes $x_{3}$. The
results are valid for any choice of initial conditions maintaining the
relation $x_{1}(0)+x_{2}(0)+x_{3}(0)=1$. For further information, please see
the main text Sec. 4.(IV.1)

.

##  III The dynamics of the replicator equations

To determine the optimal strategy for the matrix $A$ by adopting the methods
provided in Ref. Spaniel (2014), we assume $p_{i}\in[0,1]$ is the probability
of the $i$-th species with density $x_{i}$, $i=1,2,3,\cdots,N$. Let us further
assume, $N\geq 3$ is odd. The existence of the interior stationary point with
all non-zero components is unachievable for the even number of species Fisher
and Reeves (1995). The impossibility of maintenance of biodiversity with an
even number of species is emphasized on the Appendix (Dynamics of even
species). Let $\mathbf{p}=\big{(}p_{1},p_{2},p_{3},\cdots,p_{N}\big{)}^{T}$ be
the stochastic vector with $0\leq p_{i}\leq 1$. The expected values of each
species’s density for the row player are

$\displaystyle\begin{split}E(x_{1})=(A\mathbf{p})_{1}=0.5p_{1}+\alpha(p_{2}+p_{4}+\cdots+p_{N-1})\\\
+(1-\alpha)(p_{3}+p_{5}+\cdots+p_{N}),\\\
E(x_{2})=(A\mathbf{p})_{2}=(1-\alpha)(p_{1}+p_{4}+p_{6}+\cdots+p_{N-1})\\\
+0.5p_{2}+\alpha(p_{3}+p_{5}+\cdots+p_{N}),\\\ \vdots\hskip 142.26378pt\\\
E(x_{N})=(A\mathbf{p})_{N}=\alpha(p_{1}+p_{3}+\cdots+p_{N-2})\\\
+(1-\alpha)(p_{2}+p_{4}+\cdots+p_{N-1})+0.5p_{N},\end{split}$ (5)

along with $\sum_{i=1}^{N}p_{i}=1$. Here, $(A\mathbf{p})_{i}$ represents the
$i$-th row of the matrix $A\mathbf{p}$. We mainly assume $\alpha\neq 0.5$ to
avoid the occurrence of the rank-deficient matrix for which the matrix $A$
provides multiple optimal strategies. Solving Eq. (5), we obtain
$p_{i}=\dfrac{1}{N}$, $i=1,2,\dots,N$. Thus, the optimal strategy suggests all
species will coexist with equal densities. Substituting $x_{i}=\dfrac{1}{N}$
in Eq. (4) for odd $N$, we find
$\bigg{(}\dfrac{1}{N},\dfrac{1}{N},\cdots,\dfrac{1}{N}\bigg{)}$ is a
stationary point of the Eq. (4), as

$\displaystyle\sum_{j=1}^{N}\sum_{k=1}^{N}B_{ijk}=0,\hskip
14.22636pti=1,2,\dots,N.$ (6)

Thus, the equilibrium $(p_{1},p_{2},\cdots,p_{N})$ calculated through the
interaction matrix $A$ serves as a stationary point
$(x_{1}=p_{1},x_{2}=p_{2},\cdots,x_{N}=p_{N})$ of the system (4). For the
validation, we first choose the simpliest case of $N=3$. We find the system
converges to the stationary point
$\bigg{(}\dfrac{1}{3},\dfrac{1}{3},\dfrac{1}{3}\bigg{)}$ for all initial
conditions chosen from $(0,1)\times(0,1)\times(0,1)$. The initial conditions
are chosen maintaining the constraint $x_{1}(0)+x_{2}(0)+x_{3}(0)=1.0$. The
local stability analysis provides the following three eigen values

$\displaystyle\lambda_{1,2}=-0.0625\pm 0.4315i\text{and}\hskip
2.84544pt\lambda_{3}=-0.0207,$ (7)

and hence, it suggests the stationary point
$\bigg{(}\dfrac{1}{3},\dfrac{1}{3},\dfrac{1}{3}\bigg{)}$ is a stable focus-
node. In fact, we have performed a global stability analysis by choosing the
following Lyapunov function,

$\displaystyle V(x_{i})=-\dfrac{1}{N}\sum_{i=1}^{N}\log(Nx_{i}),$ (8)

where $N\geq 3$ is an odd integer. Since $\log$ is a concave function, thus by
using Jensen’s inequality Jensen (1906); Dekking  _et al._ (2005); Ruel and
Ayres (1999), we find

$\displaystyle\sum_{i=1}^{N}\frac{1}{N}\log(Nx_{i})\leq\log\bigg{(}\sum_{i=1}^{N}\frac{1}{N}(Nx_{i})\bigg{)}=\log\sum_{i=1}^{N}x_{i}=\log(1)=0.$
(9)

Thus, clearly $V(x_{i})\geq 0$, for $x_{i}\in(0,1)$ and

$V\bigg{(}\dfrac{1}{N}\bigg{)}=-\dfrac{1}{N}\sum_{i=1}^{N}\log\bigg{(}\frac{N}{N}\bigg{)}=0.$
(10)

Let $\xi_{j}$ be a non-zero perturbation such that

$\xi_{j}=x_{j}-\dfrac{1}{N},j=1,2,\cdots,N.$ (11)

Then, $\sum_{j=1}^{N}\xi_{j}=0$. Now,

$\displaystyle\begin{split}\frac{dV}{dt}=-\frac{2}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}\big{(}A_{ij}A_{ik}+A_{ij}A_{jk}\big{)}x_{j}x_{k}+1\\\\[0.0pt]
=-2\sum_{j=1}^{N}\sum_{k=1}^{N}\big{(}\sum_{i=1}^{N}\frac{1}{N}A_{ij}\big{)}A_{jk}x_{j}x_{k}\\\
-\frac{2}{N}\sum_{i=1}^{N}\big{(}\sum_{j=1}^{N}A_{ij}x_{j}\big{)}\big{(}\sum_{k=1}^{N}A_{ik}x_{k}\big{)}+1\\\\[0.0pt]
\end{split}$

Since
$\sum_{i=1}^{N}\frac{1}{N}A_{ij}=\dfrac{\frac{N-1}{2}+0.5}{N}=\frac{1}{2}$ and
$\sum_{j=1}^{N}\sum_{k=1}^{N}A_{jk}x_{j}x_{k}=\frac{1}{2}$, we have

$\displaystyle\begin{split}\frac{dV}{dt}=-\frac{2}{N}\sum_{i=1}^{N}\big{(}\sum_{j=1}^{N}A_{ij}x_{j}\big{)}^{2}+\frac{1}{2}\\\
=-\frac{2}{N}\sum_{i=1}^{N}\bigg{(}\sum_{j=1}^{N}A_{ij}\bigg{(}\frac{1}{N}+\xi_{j}\bigg{)}\bigg{)}^{2}+\frac{1}{2}\\\\[0.0pt]
=-\frac{2}{N}\sum_{i=1}^{N}\bigg{(}\sum_{j=1}^{N}A_{ij}\xi_{j}\bigg{)}^{2}-\frac{2}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}A_{ij}\xi_{j}\\\
=-\frac{2}{N}\sum_{i=1}^{N}\bigg{(}\sum_{j=1}^{N}A_{ij}\xi_{j}\bigg{)}^{2}-\sum_{j=1}^{N}\xi_{j}\\\\[0.0pt]
=-\frac{2}{N}\sum_{i=1}^{N}\bigg{(}\sum_{j=1}^{N}A_{ij}\xi_{j}\bigg{)}^{2}<0,\end{split}$

unless, $A$ is a rank deficient matrix. Thus, the stationary point
$\bigg{(}\dfrac{1}{N},\dfrac{1}{N},\cdots,\dfrac{1}{N}\bigg{)}$ is globally
stable with odd $N$. Therefore, the system (4) may possess several stationary
points; nevertheless, the system always converges to the globally stable
stationary point predicted by the equilibrium calculated through the
interaction matrix $A$ given in (1) for odd $N\in\mathbb{N}$.

##  IV Consequences of perturbation among the interacting species

###  IV.1 Perturbation in single species-pair $(N=3)$

Now, we want to understand the impact of a perturbation at a single entry of
the interaction matrix $A$. Does it alter the globally stable stationary
point? If yes, can we predict it before investigating the dynamics through
numerical simulations? Again, we consider the simplified choice to investigate
this and begin with $N=3$. The perturbed matrix looks like

$\tilde{A}=\begin{bmatrix}0.5&\alpha+\beta&1-\alpha\\\
1-\alpha-\beta&0.5&\alpha\\\ \alpha&1-\alpha&0.5\end{bmatrix}.$ (12)

Still this new interaction matrix satisfies $\tilde{A}_{ij}+\tilde{A}_{ji}=1$,
as each entry here too contemplates the probability of winning among the
species. Note that we just perturb a single pair of species’ interaction,
which is $A_{12}$ and consequently $A_{21}$ without loss of any generality in
this case. Clearly, $\beta$ may be positive or negative, however we choose the
specific interval which can mimick the physically meaningful interaction
matrix. The choice of $\beta$ is based on the following two criteria,

  * •

One needs to choose $\beta$ in such a way that each entry of the interaction
matrix lies within $[0,1]$ and indicates a non-negative real number
representing a probability of winning the pairwise competition.

  * •

Applying $\beta$ does not change the role of stronger and weaker species; it
only increases or reduces the interaction between two competitors. This is to
make sure that every species coexist after applying perturbation. Otherwise,
before perturbation, every species was stronger and weaker in the same number
of interactions; reversing the role will form a disparity and eventually lead
to the destruction of all species’ coexistence.

Figure 4: Impact of a single perturbation on the first-second species
interaction: (a, d) A common convergence for all the species in an identical
density is observed for the unperturbed case with $\beta=0$. (b, e) For
$\beta=0.1$, the external perturbation reduces their respective densities of
first and second species. The emergence of two clusters is represented in
subfigure (e) after the initial transient. (c, f) For $\beta=-0.1$, the
internal perturbation enhances the densities of two perturbed species. The
unperturbed three species maintain a coherent state. All these observations
remain valid for any permissible choice of initial conditions. Here,
$\alpha=0.75$ and $N=5$.

Here $\beta=0$ leads to $A=\tilde{A}$. Irrespective of the sign of $\beta$,
one can determine the optimal strategy by equating the expected values of all
three species corresponding to the row. We have the following three equations

$\displaystyle E(x_{1})=(1-\alpha)+(\alpha-0.5)p_{1}+(2\alpha+\beta-1)p_{2},$
(13) $\displaystyle E(x_{2})=\alpha+(1-2\alpha-\beta)p_{1}+(0.5-\alpha)p_{2},$
$\displaystyle E(x_{3})=0.5+(\alpha-0.5)p_{1}+(0.5-\alpha)p_{2},$

where $p_{1}$, $p_{2}$, and $p_{3}=1-p_{1}-p_{2}$ are the probabilities of
species survivability. Solving the Eq. (13), we obtain

$\displaystyle p_{1}=\frac{\alpha-0.5}{3(\alpha-0.5)+\beta},$ (14)
$\displaystyle p_{2}=\frac{\alpha-0.5}{3(\alpha-0.5)+\beta},$ $\displaystyle
p_{3}=\frac{(\alpha-0.5)+\beta}{3(\alpha-0.5)+\beta}.$

Figure 5: The impact of localized perturbation on a single pair of
interaction: We have perturbed the first and second species’ density by
$\beta\neq 0$. Due to our choice of $\alpha=0.75>0.50$, $\beta$ can be varied
within $[-0.25,0.25]$. However, we ignore the case for $\beta=-0.25$ to avoid
the rank deficient matrix $\tilde{A}$. If we choose $\beta<-0.25$, the
interaction between the first and the second species gets reversed, so we
disregard these values of $\beta$. For $\beta>0.25$, $\alpha$ exceeds the
value $1$. (a) External perturbation ($\beta>0$). Here the unperturbed species
is marked by blue diamonds and first two perturbed species are marked by red
circles. Here the abundance of unperturbed species increases slowly and the
abundances of perturbed species are decreased slowly. (b) Internal
perturbation ($\beta<0$). Here the abundance of unperturbed species (blue
diamonds) is dropped rapidly. The analytical calculation fits perfectly with
the numerically simulated values shown by the markers.

Initially, we start with $\beta>0$ for $\alpha>0.5$ and $\beta<0$ for
$\alpha<0.5$. The valid range for these choices of $\beta$ is $[-\kappa,0)$ or
$(0,-\kappa]$, where, $\kappa={\rm min}\\{\alpha,1-\alpha\\}$. For these
choices of $\beta$, we have
$||A_{ij}-A_{ji}||<||\tilde{A}_{ij}-\tilde{A}_{ji}||$. We refer to this state
as the external perturbation because this choice of $\beta$ enhances the
surviving probability of the $i$-th species (which is the first species in the
case of the matrix (12)) when interacted with the $j$-th species (the second
species for our particular choice in matrix (12)). To scrutinize the impact of
$\beta$, we choose a particular value of $\beta=0.1$ and $\alpha=0.75$. The
numerical integration of Eq. (4) with random permissible initial conditions
provides $x_{1}=0.29412,x_{2}=0.29412,x_{3}=0.41176$ (See Fig. (3) (b),
orange, blue and purple respectively). It is a stable focus-node as the eigen
values of the corresponding Jacobian are

$\lambda_{1,2}=-0.0788\pm 0.4792i\hskip 5.69046pt\text{and}\hskip
5.69046pt\lambda_{3}=-0.0226.$ (15)

Again, we need to choose the initial conditions maintaining the constraint
$x_{1}(0)+x_{2}(0)+x_{3}(0)=1$. Interestingly, our numerically simulated
values, in this case too, are well fitted with our analytically predicted
values $\eqref{eq12}$. The densities of the first (orange) and second (blue)
species reduce compared to the unperturbed ($\beta=0$) case (See Fig. (3)
(a-b)). This scenario suggests external perturbation is helpful in gaining
species density for the unperturbed species as $x_{3}>x_{2}=x_{1}$. Apart from
the analytical understanding, we realize this finding is absolutely
recognizable through the physical interpretation. An external perturbation in
the first-second interaction increases the probability of the survivability of
the first species against the second species. This, in turn, decreases the
density of the second species, and as the second species dominates the third
species, the density of the third species goes up. On the other hand, as the
third species dominate the first species, increasing $x_{3}$ implies firm
control over the first species, lowering its density.

For $\beta<0$ for $\alpha>0.5$ and $\beta>0$ for $\alpha<0.5$, we can define
an internal perturbation for which we have
$||A_{ij}-A_{ji}||>||\tilde{A}_{ij}-\tilde{A}_{ji}||$. The choice of $\beta$
is restricted within the range $(-\rho,0)$ or $(0,\rho)$, where
$\rho=\frac{1}{2}|1-2\alpha|$. We call this as internal perturbation as it is
reducing the competition between them leading them to the equal cooperation.
We have encapsulated all the requisite information regarding the
classification of internal and external perturbation in Table (1). However, we
discard the case of $\alpha=0.5$, which creates a rank deficient matrix that
yields multistable system (4). We choose $\beta=-0.1$ for $\alpha=0.75$, we
find the densities of the first and second species increase (orange and blue
respectively) and the third species’ density (purple) decreases (See Fig. (3)
(c)). An internal perturbation in the first-second species interaction
decreases the surviving probability of the first species against the second
species. This, in turn, increases the density of the second species, and as
the second species dominates the third species, the density of the third
species goes down. On the other hand, as the third species dominate the first
species, a decrease in $x_{3}$ implies weak control over the first species,
raising its density. Moreover, the simulated result in Fig. (3) (c) yields
$x_{1}=x_{2}=0.38462,x_{3}=0.23077$, which fits with our derived equilibrium
(14) of $\tilde{A}$ for $\alpha=0.75$ and $\beta=-0.1$. This stationary point
is a focus-node as per the local stability analysis with

$\lambda_{1,2}=-0.0422\pm 0.3594i\hskip 5.69046pt\text{and}\hskip
5.69046pt\lambda_{3}=-0.0166.$ (16)

Table 1: Distinguishing between internal and external perturbation depending on the choice of $\alpha$ and the permissible range of $\beta$ applied on the first and second species’ interaction $\alpha<0.5$ | $\beta>0$ | Internal perturbation
---|---|---
$\beta<0$ | External perturbation
$\alpha=0.5$ | $\beta\neq 0$ | External perturbation
$\alpha>0.5$ | $\beta>0$ | External perturbation
$\beta<0$ | Internal perturbation

Besides, its stability can also be confirmed by the similar global stability
analysis provided in Sec. (III). To prove global stability, one needs to
choose the positive-definite Lyapunov function Rao  _et al._ (2021)

$V(x_{i})=-\sum_{i=1}^{N}{x}^{*}_{i}\log{\frac{x_{i}}{{x}^{*}_{i}}},$ (17)

where ${x}^{*}_{i}\in(0,1)$ is the $i$-th component of the stationary point of
the system (4) with odd $N\in\mathbb{N}$, and proceeds like our previous
analysis to show the time derivative of this function is globally negative
definite.

Figure 6: Repercussion of two perturbations: For $\alpha=0.75$, we represent
four different cases. (a) $\beta_{1}=\beta_{2}=0.1$ provides an improved
density for the fifth species, while the rest of the densities decline
compared to $\beta_{1}=\beta_{2}=0$. (b) Through the consideration of
$-\beta_{1}=\beta_{2}=0.1$, we introduce the interplay of internal and
external perturbations on the $N=5$ species. (c) Another different combination
of external and internal perturbation with $\beta_{1}=-\beta_{2}=0.1$
demonstrate the variation of species densities, where $x_{1}=x_{2}$ and
$x_{3}=x_{4}$. (d) $\beta_{1}=\beta_{2}=-0.1$ displays the influence of two
internal perturbations on the first two species. Simulations are done with
random permissible set of initial conditions. Initially for
$\beta_{1}=\beta_{2}=0$, all species densities
$x_{1}=x_{2}=x_{3}=x_{4}=x_{5}=0.2$ are equal. For further details, please see
the main text.

We confirm our claim of predicting the asymptotic values of each species’
densities using Eq. (14) without solving it (Eq. (4)) numerically. By keeping
fixed $\alpha=0.75$, we vary $\beta$ within the interval $[-0.2,0]$ for the
internal perturbation and within $[0,0.2]$ for the external perturbation. The
employed external perturbation ($\beta>0$) with $\alpha=0.75$ destroys the
coherence among each trajectory in the asymptotic state. Since we have applied
the perturbation only in the first and second species’ density, those two
trajectories split from the synchronous state. Satisfying the calculated
stationary points in (14) shown by solid lines in Fig. (5) (a), we find the
species density of the first and second species diminish (red circles) with an
increment of $\beta\in[0,0.2]$, and as a result of that, the density of third
species (blue diamonds) increases. Our analytical prediction through the
equilibrium (14) of $\tilde{A}$ suggests $x_{1}=x_{2}$ (red circles) increase
and $x_{3}$ (blue diamonds) diminishes for negative $\beta\in[-0.2,0]$ with
$\alpha=0.75$. Even the numerically simulated values shown by the marker in
Fig. (5) (b) perfectly match the analytical predictions. Thus a suitable
single perturbation impacting over a two-species interaction can control the
density of each species. Particularly, partially motivated by the analysis of
Connolly et al.  Connolly  _et al._ (2001) we can measure the relative steady
abundance between the perturbed and unperturbed species as follows (using Eq.
(14))

$\Delta_{R}=\frac{p_{\rm perb}}{p_{\rm
unperb}}=\frac{\alpha-0.5}{(\alpha-0.5)+\beta}.$ (18)

Figure 7: Effect of two perturbations on the first four out of five species:
The triadic interaction-induced system diversities, i.e., non-trivial steady
states, are explored using numerical methods (ODE45, Matlab R2021a). The
choice of $\beta_{1}\neq\beta_{2}$ allows three distinct groups, viz.
$x_{1}=x_{2}$, $x_{3}=x_{4}$, and $x_{5}$. We plot three different parameter
space $\beta_{1}-\beta_{2}$ with respect to (a) $x_{1}(=x_{2})$, (b)
$x_{3}(=x_{4})$, and (c) $x_{5}$, respectively for fixed $\alpha=0.75$. The
numerical simulations completely agree with our derived equilibrium (24) of
$\tilde{A}$ (See the matrix (22)). The system (4) is iterated for a
sufficiently long number of iterations to disregard the initial transient and
the final densities ($x_{i}$) are collected from the asymptotic behavior of
the long term stationary state. Furthermore, we repeat the process for $50$
times independently with random permissible initial conditions to establish
the generic behavior of the triple-wise competitive interaction and eliminate
the initial conditions dependent finding.

For a suitable $\alpha$, the ratio between the unperturbed and perturbed
species is inversely proportional to $\beta$, i.e., if $\beta$ increases in
negative direction this ratio (18) sharply increases where as it slowly decays
if $\beta$ increases in positive direction. This is also clear from the Fig.
(5) (a), where, in presence of positive $\beta$ the unperturbed species (red
diamonds) slowly grows compared to perturbed species. On the contrary, it
drops rapidly for negative $\beta$. Thus, negative perturbation, i.e., a less
competition helps the species to win strongly over the unperturbed species. In
the next section we have generalized the perturbed scheme for $N$ species.

###  IV.2 Perturbation in single species-pair (for any odd $N$)

Again, we only consider an odd number of species as a comprehensive analysis
is provided in Appendix (Dynamics of even species) to demonstrate the
impossibility of the existence of the interior stationary point with an even
number of species. Till now, although we have represented all the numerical
simulations only for $N=3$. Nevertheless, one can easily derive the
equilibrium of $\tilde{A}$ for any odd $N\geq 3$ as follows,

$\displaystyle
x_{1}=x_{2}=\frac{(\alpha-0.5)}{N(\alpha-0.5)+(N-2)\beta},\text{and}$ (19)
$\displaystyle
x_{i}=\frac{(\alpha-0.5)+\beta}{N(\alpha-0.5)+(N-2)\beta},\hskip
5.69046pt\text{when}\hskip 5.69046pti=3,4,\cdots,N.$

where the perturbation $\beta\neq 0$, chosen from the permissible range, is
applied only at ${\tilde{A}}_{12}$ and ${\tilde{A}}_{21}$. We calculate this
stationary point by following the similar analysis provided in Sec. (III). We
construct $N$ linear equations using the rows of $\tilde{A}$, and solving
those $N$ hyperplanes; we obtain Eq. (19). To check the analytical expression
of Eq. (19), we consider here 5 species intransitive cycle where the
perturbation $\beta\neq 0$ is employed at the interaction of the first and
second species. The corresponding matrix looks like

$\tilde{A}=\begin{bmatrix}[r]\phantom{-}0.5&\alpha+\beta&1-\alpha&\alpha&1-\alpha\\\
1-\alpha-\beta&0.5&\alpha&1-\alpha&\alpha\\\
\alpha&1-\alpha&0.5&\alpha&1-\alpha\\\ 1-\alpha&\alpha&1-\alpha&0.5&\alpha\\\
\alpha&1-\alpha&\alpha&1-\alpha&\phantom{-}0.5\end{bmatrix}.$ (20)

When $\beta=0$ with $\alpha=0.75$, the system (4) stabilizes to the stationary
point $(0.2,0.2,0.2,0.2,0.2)$ coinciding with the derived equilibrium of $A$
for odd $N$, as portrayed in Fig. (4) (a). Figure (4) (d) depicts the snapshot
at a particular time after the transients. However, this scenario alters for
non-zero $\beta$. For $\beta=0.1$ and $\alpha=0.75$, $x_{1}$ and $x_{2}$ split
from the coherent group and we find $x_{1}=x_{2}=0.16129<x_{i}=0.22581$, for
$i=3,4,5$. The asymptotic values of these two coherent groups fit with our
derived values $\big{(}\zeta_{1},\zeta_{1},\zeta_{2},\zeta_{2},\zeta_{2})$ at
Eq. (19), where $\zeta_{1}=\frac{(\alpha-0.5)}{5(\alpha-0.5)+3\beta}$ and
$\zeta_{2}=\frac{(\alpha-0.5)+\beta}{5(\alpha-0.5)+3\beta}$. We plot the
temporal evolution of each $x_{i}$ for $i=1,2,3,4,5$ in Fig. (4) (b), and the
snapshot at a particular time after the transient in Fig. (4) (e). Apart from
the external perturbation, we also examine the impact of internal perturbation
on the interaction among the first and second species in subfigures (c) and
(f) of Fig. (4). The choice of $\beta=-0.1$ leads to an increment of
$x_{1}=x_{2}=0.26316$ as predicted by Eq. (19) compared to the unperturbed
($\beta=0$) scenario for $\alpha=0.75$. The value of
$x_{3}=x_{4}=x_{5}=0.15789$ also satisfies our theoratical calculated value in
Eq. (19). So, applying the external perturbation at the first and second
species’ interaction will reduce the corresponding densities. Whereas the
internal perturbation always boosts the densities of the related species. The
behavior of relative steady abundance $\Delta_{R}$ also follows the same
trend. A relevant question we also raise here is what will happen if unequal
perturbations are used for each pair of the species? We will deal with this
concern in the following subsection. In the following subsection, we will
investigate the influence of perturbation in multiple species.

###  IV.3 Perturbation in multiple species-pair (for any odd $N$)

Still, now, we have discussed the results of a single perturbation on the
ensemble. However, one may introduce more than a single perturbation. We
assume that a single species cannot be perturbed twice to investigate this
scenario. By allowing $k$ perturbations, one can disturb the first $2k$
species. Thus, we strictly need $N>2k$. As $N$ is odd, we consider the greater
than sign. The theoretical calculation yields the following equations
predicting the asymptotic species densities of $N$ species as follows,

$\displaystyle
x_{i}=\frac{(\alpha-0.5)}{N(\alpha-0.5)+(N-2k)\beta},\text{when}\hskip
5.69046pti=1,2,\cdots,2k$ (21) $\displaystyle
x_{i}=\frac{(\alpha-0.5)+\beta}{N(\alpha-0.5)+(N-2k)\beta},\text{when}\hskip
5.69046pti=2k+1,2k+2,\cdots,N.$

We find this solution by solving $N$ equations as shown in detail in Sec.
(III). The denominator contains the term $(N-2k)\beta$. This term appears as
$k$ perturbations affect $2k$ species, and thus, the remaining $(N-2k)$
species will experience the influence of $\beta$. Thus perturbation in a
single pair or identical perturbation in multiple species-species interaction
split the system into two clusters steady states, and $\Delta_{R}$ follows the
same trend as discussed before. Equations (21) easily leads to the earlier
discussed results in Eq. (19) by setting $k=1$.

Thus our recipe provides a way to split the entire network into two clusters,
and abundance of perturbed species pair will be increased if the pair
decreases their competition. Against this backdrop one may wonder whether this
approach can increase the number of steady clusters and can control the
abundance of each pair in the network. The simple way is to introduce the
different amounts of perturbation without disturbing the same species twice.
For instance, let us consider the following interaction matrix for $N=5$

$\tilde{A}=\begin{bmatrix}0.5&\alpha+\beta_{1}&1-\alpha&\alpha&1-\alpha\\\
1-\alpha-\beta_{1}&0.5&\alpha&1-\alpha&\alpha\\\
\alpha&1-\alpha&0.5&\alpha+\beta_{2}&1-\alpha\\\
1-\alpha&\alpha&1-\alpha-\beta_{2}&0.5&\alpha\\\
\alpha&1-\alpha&\alpha&1-\alpha&0.5\end{bmatrix}.$ (22)

To calculate the equilibrium of this matrix (22), we equate the following five
expressions,

$\displaystyle
E(x_{1})=(A\mathbf{x})_{1}=0.5x_{1}+(\alpha+\beta_{1})x_{2}+(1-\alpha)x_{3}+\alpha
x_{4}+(1-\alpha)x_{5},$ (23) $\displaystyle
E(x_{2})=(A\mathbf{x})_{2}=(1-\alpha-\beta_{1})x_{1}+0.5x_{2}+\alpha
x_{3}+(1-\alpha)x_{4}+\alpha x_{5},$ $\displaystyle
E(x_{3})=(A\mathbf{x})_{3}=\alpha
x_{1}+(1-\alpha)x_{2}+0.5x_{3}+(\alpha+\beta_{2})x_{4}+(1-\alpha)x_{5},$
$\displaystyle E(x_{4})=(A\mathbf{x})_{4}=(1-\alpha)x_{1}+\alpha
x_{2}+(1-\alpha-\beta_{2})x_{3}+0.5x_{4}+\alpha x_{5},$ $\displaystyle
E(x_{5})=(A\mathbf{x})_{5}=\alpha x_{1}+(1-\alpha)x_{2}+\alpha
x_{3}+(1-\alpha)x_{4}+0.5x_{5},$

and derive a set of four distinct equations
$(A\mathbf{x})_{1}=(A\mathbf{x})_{2}$, $(A\mathbf{x})_{1}=(A\mathbf{x})_{3}$,
$(A\mathbf{x})_{1}=(A\mathbf{x})_{4}$, and
$(A\mathbf{x})_{1}=(A\mathbf{x})_{5}$ along with the constraint
$x_{1}+x_{2}+x_{3}+x_{4}+x_{5}=1$. Hence, we have five equations in five
variables $x_{1}$, $x_{2}$, $x_{3}$, $x_{4}$, and $x_{5}$. Solving these five
equations, we obtain

$\displaystyle x_{1}=x_{2}=\frac{d_{1}(d_{1}+\beta_{2})}{d_{2}},$ (24)
$\displaystyle x_{3}=x_{4}=\frac{d_{1}(d_{1}+\beta_{1})}{d_{2}},$
$\displaystyle x_{5}=\frac{(d_{1}+\beta_{2})(d_{1}+\beta_{1})}{d_{2}},$

where $d_{1}=\alpha-0.5$, and
$d_{2}=(\alpha-0.5)[5(\alpha-0.5)+3\beta_{1}+3\beta_{2}]+\beta_{1}\beta_{2}$.

Figure 8: Temporal evolution of three species under the triadic interaction
for diverse choices of perturbations: $N=3$ species are interacted as per the
interaction matrix $\tilde{A}$ given in Eq. (29) with fixed $\alpha=0.25$. We
consider eight different cases, viz. (a) $\beta_{i}>0$, $i=1,2,3$, (b)
$\beta_{1},\beta_{2}>0$, and $\beta_{3}<0$, (c) $\beta_{1}>0$, $\beta_{2}<0$,
and $\beta_{3}>0$, (d) $\beta_{1}>0$, $\beta_{2}<0$, and $\beta_{3}<0$, (e)
$\beta_{1}<0$, $\beta_{2}>0$, and $\beta_{3}>0$, (f) $\beta_{1}<0$,
$\beta_{2}>0$, and $\beta_{3}<0$, (g) $\beta_{1}<0$, $\beta_{2}<0$, and
$\beta_{3}>0$, and (h) $\beta_{1}<0$, $\beta_{2}<0$, and $\beta_{3}<0$. All
these possible combinations allow species coexistence with different densities
(a) $(0.08333,0.6667,0.25)$, (b) $(0.45,0.4,0.15)$, (c)
$(0.16667,0.33333,0.5)$, (d) $(0.64286,0.14286,0.21429)$, (e)
$(0.0625,0.5,0.4375)$, (f) $(0.375,0.33333,0.29167)$,(g) $(0.1,0.2,0.7)$, and
(h) $(0.5,0.1111,0.3889)$, respectively. Here, $|\beta_{1}|=0.1$,
$|\beta_{2}|=0.15$, and $|\beta_{3}|=0.2$. The blue, red, and yellow curves
contemplate temporal evolutions of the densities $x_{1}$, $x_{2}$, and
$x_{3}$, respectively.

This formula is equivalent to the relation (22) for $k=2$ and
$\beta_{1}=\beta_{2}$. Figure (6) portrays the impact of two different amounts
of perturbations on the first four species among $N=5$ species with
$\alpha=0.75$. When $\beta_{1}=\beta_{2}=0.1$, Fig. (6) (a) reflects the
decreament of first four species’ density due to the employed external
perturbations. This outcome attests to the essential signature of external
perturbation on the species. The asymptotic values of the
$x_{1}=x_{2}=x_{3}=x_{4}=0.18519$, and $x_{5}=0.25926$ fits absolutely
perfectly with the derived equilibrium (24) of $\tilde{A}$ (c.f. (22)) for
$\beta_{1}=\beta_{2}=0.1$. Interestingly, when $\beta_{2}=0.1>0$ and
$\beta_{1}=-0.1<0$, the first two species experience the internal
perturbation, and the third and fourth species feel the presence of the
external one. The fifth species is left undisturbed. A clear expectation from
our previous observation is that $x_{1}=x_{2}$ (orange and blue bar in Fig.
(6) (b)) will increase and $x_{3}=x_{4}$ (purple and brown) will decrease
compared to the unperturbed case (cyan). The steady state values of the three
clusters are $x_{1}=0.28926$, $x_{2}=0.28926$, $x_{3}=0.12397$,
$x_{4}=0.12397$, and $x_{5}=0.17355$. Similarly, here we can design the
relative density of first (or second) species with the last one as
$\Delta_{R}(1/2-5)=\frac{\alpha-0.5}{\alpha-0.5+\beta_{2}}$ and the relative
density of third (or forth) species with the last one as
$\Delta_{R}(3/4-5)=\frac{\alpha-0.5}{\alpha-0.5+\beta_{1}}$. Clearly, if we
increase the negative strength of $\beta_{1}(\beta_{2}=0.1)$, the
$\Delta_{R}(3/4-5)$ will increase, but $\Delta_{R}(1/2-5)$ will be fixed.
Also, if we increase the positive strength of $\beta_{2}(\beta_{1}=-0.1)$,
$\Delta_{R}(1/2-5)$ will decrease, however $\Delta_{R}(3/4-5)$ will be fixed.
Noticeably, the first two species always win in the competition due to the
internal perturbation.

A similar but reverse scenario is found in Fig. (6) (c), where the external
perturbation in the first two species (orange and blue bar) with
$\beta_{1}=0.1>0$ reduces their density effectively. The internal perturbation
with $\beta_{2}=-0.1<0$ positively affects the third and fourth species
(purple and brown bar), and their densities enhance compared to the
undisturbed case (cyan). We find all the species settle down to a three-
cluster state with $x_{3}=x_{4}=0.28926>x_{5}=0.17355>x_{1}=x_{2}=0.12397$.
The influence of internal perturbations on the first four species is
investigated in Fig. (6) (d) by considering $\beta_{1}=\beta_{2}=-0.1$. The
system (4) converges to the stable stationary point
$x_{1}=x_{2}=x_{3}=x_{4}=0.21739,x_{5}=0.13043$, which matches excellently
with our analytical derivation (24). The impact of such negative perturbations
improves the species density and consequently, diminishes the unperturbed
species density $x_{5}$.

Thus, the previous analysis confirms the appearance of three distinct clusters
$x_{1}=x_{2}$, $x_{3}=x_{4}$ and $x_{5}$ for non-zero perturbations with
$\beta_{1}\neq\beta_{2}$. To further scrutinize, we plot three two dimensional
parameter spaces in Fig. (7) by varying $\beta_{1}$ and $\beta_{2}$ within
$[-0.2,0.2]$ for $N=5$. $\alpha$ is kept fixed here at $0.75$. The obtained
results are completely fitted with our derived species densities at Eq. (24).
For a fixed value of $\alpha=0.75$ and fixed value of $\beta_{1}$,
$x_{1}=x_{2}$ will magnify for increament of $\beta_{2}$ as revealed through
Fig. (7) (a). A similar increment of densities in $x_{3}=x_{4}$ is observed
with enhancement of $\beta_{1}$ in Fig. (7) (b) for fixed values of
$\alpha=0.75$ and $\beta_{2}$. Figure (7) (c) uncovers the reciprocity of
$\beta_{1}-\beta_{2}$ for $\alpha=0.75$ for the unperturbed fifth species.

In fact, one can quickly determine the effect of $k<N$ perturbations
$\beta_{1}$, $\beta_{2}$, $\cdots$, $\beta_{k}$ on first $2k$ out of $N$
species. Here, we assume $N$ is odd, and a single species cannot be perturbed
more than once. Under these circumstances, the species densities will be
derived analytically by calculating the equilibrium of $\tilde{A}$ as follows,

$\displaystyle
x_{1}=x_{2}=\frac{(\alpha-0.5)(\alpha-0.5+\beta_{2})(\alpha-0.5+\beta_{3})\cdots(\alpha-0.5+\beta_{k})}{\chi},$
(25) $\displaystyle
x_{3}=x_{4}=\frac{(\alpha-0.5)(\alpha-0.5+\beta_{1})(\alpha-0.5+\beta_{3})\cdots(\alpha-0.5+\beta_{k})}{\chi},$
$\displaystyle
x_{5}=x_{6}=\frac{(\alpha-0.5)(\alpha-0.5+\beta_{1})(\alpha-0.5+\beta_{2})\cdots(\alpha-0.5+\beta_{k})}{\chi},$
$\displaystyle\vdots$ $\displaystyle
x_{2k-1}=x_{2k}=\frac{(\alpha-0.5)(\alpha-0.5+\beta_{1})(\alpha-0.5+\beta_{2})\cdots(\alpha-0.5+\beta_{k-1})}{\chi}.$

Here,

$\displaystyle\chi=(\alpha-0.5)\cdots[(\alpha-0.5)[(\alpha-0.5)[N(\alpha-0.5)+(N-2)B_{1}]$
(26) $\displaystyle+(N-4)B_{2}]+(N-6)B_{3}]+(N-8)B_{4}]\cdots+(N-2k)B_{k},$

where,

$\displaystyle\begin{split}B_{1}=\sum_{i=1}^{k}\beta_{i},B_{2}=\sum_{i=1}^{k}\sum_{j\leq
i}\beta_{i}\beta_{j},\\\ B_{3}=\sum_{i=1}^{k}\sum_{j\leq i}\sum_{p\leq
j}\beta_{i}\beta_{j}\beta_{p},\cdots,\hskip
42.67912ptB_{k}=\beta_{1}\beta_{2}\beta_{3}\cdots\beta_{k}.\end{split}$ (27)

The densities of the remaining $(N-2k)$ unperturbed species are equal and it
can be calculated from the constraint $\sum_{i=1}^{N}x_{i}=1$. Thus, we have

$\displaystyle x_{\rm
unperturbed}=\dfrac{(\alpha-0.5+\beta_{1})(\alpha-0.5+\beta_{2})(\alpha-0.5+\beta_{3})\cdots(\alpha-0.5+\beta_{k})}{\chi}.$
(28)

More generally relative density of a species having $\beta_{l}$ perturbation
with respect to a unperturbed species can be written as $\Delta_{R}(\rm
Pert(\beta_{i})-Unpert)=\frac{\alpha-0.5}{\alpha-0.5+\beta_{i}}$. So, with the
increase of positive strength of $\beta_{k}$, this relative density will
decrease and reverse in the case of an increase in negative strength.

###  IV.4 Random perturbation in $\dfrac{N(N-1)}{2}$ species-pair (Case:
$N=3$)

Still, now we have applied the perturbations so that the same species is not
perturbed more than once. What if we disturb the same species more than once?
Does it affect the prediction of steady states through the equilibrium of the
interaction matrix? To investigate this particular case, we consider $N=3$
species and examine with the following interaction matrix,

$\tilde{A}=\begin{bmatrix}0.5&\gamma_{1}&\gamma_{2}\\\
1-\gamma_{1}&0.5&\gamma_{3}\\\ 1-\gamma_{2}&1-\gamma_{3}&0.5\end{bmatrix}.$
(29)

Figure 9: Non-existence of an even number of species: (a) $\alpha=0.75>0.5$,
and (b) $\alpha=0.25<0.5$ lead to the extinction of the fourth and first
species, respectively. Both cases allow the coexistence of the remaining
$(N-1)$ species with equal densities. Here, $N=4$.

where $\gamma_{1}=\alpha+\beta_{1}$, $\gamma_{2}=1-\alpha+\beta_{2}$, and
$\gamma_{3}=\alpha+\beta_{3}$. Each $\gamma_{i}$ lies within
$[0,1]-\\{0.5\\}$. We can simulate the equilibrium of this matrix similarly by
following the earlier process in Eq. (5), and we obtain

$\displaystyle p_{1}=\frac{\nu_{3}}{\nu},\hskip
5.69046ptp_{2}=\frac{\nu_{2}}{\nu},\hskip 5.69046pt\text{and}\hskip
5.69046ptp_{3}=\frac{\nu_{1}}{\nu},$ (30)

where $\nu_{1}=(\alpha-0.5)+\beta_{1}$, $\nu_{2}=(\alpha-0.5)-\beta_{2}$,
$\nu_{3}=(\alpha-0.5)+\beta_{3}$, and
$\nu=\nu_{1}+\nu_{2}+\nu_{3}=3(\alpha-0.5)+(\beta_{1}-\beta_{2}+\beta_{3})$.
Here we have fixed $\alpha$ at $0.25$, and we have considered one particular
case where $|\beta_{1}|=0.1<|\beta_{2}|=0.15<|\beta_{3}|=0.2$. We choose
random initial conditions from the permissible range and plot the variation of
their corresponding densities (Fig. (8)), including the transient. Depending
upon the condition that $|\beta_{1}|<|\beta_{2}|<|\beta_{3}|$, we observe that
the order of species density is same for $\beta_{1}<0$ and $\beta_{1}>0$ while
keeping $\beta_{2}$ and $\beta_{3}$ unchanged (see column-wise in Fig. (8) i.e
sub-figure "(a) and (e)" or "(b) and (f)" etc.). At first, we focus on the
Fig. (8) (a). Here all $\beta$ s are positive. Therefore, at $\alpha=0.25$,
using the Eq. (30) we can easily write $|\nu_{2}|>|\nu_{2}|>|\nu_{3}|$. This
condition leads us to the comparative relation among the species as:
$x_{1}({\rm blue})<x_{3}({\rm yellow})<x_{2}({\rm red})$. Noticeably, due to
the presence of all perturbation, the enhanced cooperation between the first
two species (as $\beta_{1}<0$) is not helping the first species to gain more
density, rather it loses its density (Fig. (8) (a)), however species wins in
the competition. Also in the case of $\beta_{1}<0$, $\beta_{2,3}>0$,
$|\nu_{2}|>|\nu_{1}|$ as $|\beta_{1}|<|\beta_{2}|$, and the condition
$|\nu_{1}|>|\nu_{3}|$ is obvious as the sign of $\beta_{1}$ is opposite ($<0$)
to $\beta_{3}$ (>0) which makes $\nu_{1}$ more negative. Thus, the same
relation among the species density appears in the Fig. (8)(e). Now let’s
consider the case where $\beta_{1}>0$ (or $\beta_{1}<0$), $\beta_{2}<0$, and
$\beta_{3}>0$. Using the (30) we can see that $|\nu_{1}|>|\nu_{2}|>|\nu_{3}|$
and hence we observe the species density having the order $x_{1}<x_{2}<x_{3}$.
Using the same logic along with (30) we can explain the scenarios for another
particular order of $\beta_{1}$, $\beta_{2}$ and $\beta_{3}$.

##  V Conclusion: Current status and future prospects

To understand the maintenance of biodiversity observed in a realistic
scenario, recently, an increasing trend has been observed to study ecological
systems in the light of higher-order interaction (HOI) Bairey  _et al._
(2016); Grilli  _et al._ (2017); Mayfield and Stouffer (2017); Levine  _et
al._ (2017); Singh and Baruah (2021); Majhi  _et al._ (2022). Letten et al.
Letten and Stouffer (2019) have not only shown how these models approximate
the ecological system with better accuracy but also provided motivation for
how this indirect effect arises from the phenomenological models. So this
paves the way for choosing HOI in our proposed model. It is still debated
whether pairwise interactions are distinct from HOI and how to define and
implement these interactions. Levine et al. Levine  _et al._ (2017) have
argued that both interactions operate over different timescales, making them
distinct, but sometimes both timescales are comparable Grilli  _et al._
(2017). Actually, looking at this through the eyes of a mechanical model of
resource competition, no distinction can be found Letten and Stouffer (2019).
Hence, we have considered only the HOI to describe the model.

Furthermore, in the present article, HOIs are fully characterized by the
pairwise interaction without adding new parameters and making this model
empirically testable Grilli  _et al._ (2017). On the other hand, in ecological
systems, RPS dynamics are very well studied. It is essential from the
perspective of describing intransitive cycles, cyclic dominance of species,
and also describing the theory of species co-existence and diversityBuss and
Jackson (1979); Sinervo and Lively (1996); Paquin and Adams (1983); Kerr  _et
al._ (2002); Allesina and Levine (2011). Hence, we developed our model based
on generalized RPS game which is really in the center of interest for
ecologists. The present manuscript reflects how a network of interacting
species can lead to species co-existence and diversity. Also, how changing one
interaction can lead to different indirect effects in other species is
characterized here. In the presence of higher-order interaction, we
investigate the structurally perturbed intransitive cycle of $N$ species. We
have determined analytically how the globally stable coexisting equilibrium
points are arranged for a long-term evaluation. We have found that when two
species reduce their competition (winning probability is decreased slightly
from the stronger species), their density increases compared to other species.
The analytically calculated globally stable stationary points match perfectly
with the numerical derivations. We reveal how the feasible perturbation in the
interaction matrix splits the long-time behavior of the replicator equations
into a few synchronized clusters.  Finally, these results indicate the
interaction of how weighted network structures between species can determine
the population density in the ecological dynamics Rao  _et al._ (2021).
Internal perturbation is essentially the positive evolution of weaker species,
so the winning probability of stronger decreases. Similarly, external
perturbation is the positive evolution of stronger, so the winning likelihood
of stronger increases. We have shown how the interplay of this evolution leads
to diverse species density in a system. Biologists can further test this for
competing bacteria in laboratory conditions Friedman  _et al._ (2017).

Future work relies on relaxation of the assumption like species interacting
with a particular strength and constraint on the perturbation range and number
of perturbations can be given. This will allow us to arrive at the results in
a more general setting. Now, depending on the environmental changes, this
perturbation can exceed the range we have considered in this model. This will
lead to a reversal in dominance, ultimately implying the extinction of some
species, which will be more realistic in the large ecological system. One may
wish to look at the spread of the stationary points in a large system when the
perturbation is chosen from a truncated normal distribution with having a
fixed interaction matrix and finally in a more general case where the
interaction matrix is also a random matrix May (1972); Holling (1973);
Allesina and Tang (2012); Grilli  _et al._ (2017). This spread of equilibrium
points corresponds to different species having very diverse species densities.
One may further consider those fundamental concepts like resource consumption
and other vast range of physical parameters that can help better understand
how small perturbations, i.e., evolution or environmental changes, can drive a
species to extension. We conclude with the hope that our perceived theoretical
model offers a feasible step toward revealing the influence of higher-order
interactions on ecological communities, and our mathematical convenience may
be pertinent to much applied biological realism.

## Appendix

### Dynamics of even species

Nevertheless, when $N\geq 4$ is even, then it is impossible to calculate the
equilibrium of ${A}$ given in Eq. (1). For example, let us consider $N=4$.
Then the interaction matrix is given by

${A}=\begin{bmatrix}0.5&\alpha&1-\alpha&\alpha\\\
1-\alpha&0.5&\alpha&1-\alpha\\\ \alpha&1-\alpha&0.5&\alpha\\\
1-\alpha&\alpha&1-\alpha&0.5\end{bmatrix}.$ (31)

We omit the case of $\alpha=0.5$ to avoid multiple optimal strategies.
Proceeding as the previous way given in Eq. (5), we get

$\displaystyle E(x_{1})=(0.5-\alpha)p_{1}+(1-2\alpha)p_{3}+\alpha,$ (32)
$\displaystyle E(x_{2})=(\alpha-0.5)p_{2}+(2\alpha-1)p_{3}+(1-\alpha),$
$\displaystyle E(x_{3})=(1-2\alpha)p_{2}+(0.5-\alpha)p_{3}+\alpha,$
$\displaystyle
E(x_{4})=(0.5-\alpha)p_{1}+(\alpha-0.5)p_{2}+(0.5-\alpha)p_{3}+0.5,$

where we use $p_{1}+p_{2}+p_{3}+p_{4}=1$ to eliminate $p_{4}$. To evaluate the
equilibrium of the above matrix $A$, we equate each of these expected values

$E(x_{1})=E(x_{2})=E(x_{3})=E(x_{4}).$ (33)

For $E(x_{1})=E(x_{4})$, we get

$p_{2}+p_{3}=1,$ (34)

since $\alpha\neq 0.5$. Similarly, for $E(x_{2})=E(x_{3})$, we get

$p_{2}+p_{3}=\frac{2}{3}.$ (35)

Figure 10: Effect of localized perturbation on even number species: We choose
$N=4$ species and apply the perturbation of strength $\beta=0.15$ on the
interaction of first-second species. Depending on the choice $\alpha\neq 0.5$,
the system converges to different species densities. (a) When
$\alpha=0.25<0.5$, the first species goes extinct similar to the unperturbed
case ($\beta=0$), as shown in Fig. (9) (b). The system settles down to the
stationary point $\big{(}0,\frac{1}{3},\frac{1}{3},\frac{1}{3}\big{)}$. (b)
The $\alpha=0.75>0.5$ leads to the extinction of the fourth species in the
long run. The remaining three species will survive maintaining Eq. (14), and
the system (4) converges to $(0.2778,0.2778,0.4444,0)$. The result is
independent of the choice of initial conditions chosen from the permissible
domain.

This clearly indicates that Eq. (33) is a set of inconsistent equations.
Hence, few of these $p_{i}$’s must be zero. But, as soon as any one of these
$p_{i}$ is zero, we have a system of $(N-1)$ species, and since $N$ is even,
$(N-1)$ becomes odd. Thus, as per our previous analysis, the system settles
down to a stable equilibrium with one zero component and $(N-1)$ non-zero
component. Our numerical simulation suggests depending on the choice of
$\alpha\neq 0.5$, either the first or the $N$-th species will die when $N$ is
even. For $\alpha>0.5$, we have

$\displaystyle x_{i}=\dfrac{1}{N-1},\text{for}\hskip
5.69046pti=1,2,3,\cdots,N-1$ (36) $\displaystyle x_{N}=0.$

Similarly, we obtain for $\alpha<0.5$ the following stationary point,

$\displaystyle x_{1}=0,$ (37) $\displaystyle
x_{i}=\dfrac{1}{N-1},\text{for}\hskip 5.69046pti=2,3,\cdots,N$

Figure (9) demonstrates the extinction of one species and the remaining
$(N-1)=3$ species coexist with equal densities $\frac{1}{N-1}=\frac{1}{3}$
irrespective of the choice of $\alpha\neq 0.5$. Figure (9)(a) for
$\alpha=0.75$, the fourth species dies out in the long term, and in Fig. (9)
(b) for $\alpha=0.25$, the first species extincts. The choice of initial
conditions maintaining the constraint $x_{1}(0)+x_{2}(0)+x_{3}(0)=1$ does not
affect the final outcome. Thus, when $N$ is even, any one of the species
(either first or the $N$-th depending on the $\alpha$) die out, and the
prediction of the stationary point containing the remaining $(N-1)$ species
can be quickly done using the equilibrium of the interaction matrix $A$, as
$(N-1)$ is odd.

We consider the following approach to prove the impossibility of gaining a
feasible solution with $A$ of even order. Let us consider the matrix $C$ with
entries

$C_{ij}=\begin{cases}0\quad&\text{if}\,i=j\\\ 1\quad&\text{if}\,A_{ij}>0.5\\\
-1\quad&\text{if}\,A_{ij}<0.5\\\ \end{cases}$ (38)

The entries of the matrix $C$ can be interpreted in the way that if $i$-th
species win over the $j$-th, then we assign $C_{ij}=1$. That is in the case of
a directed graph; if the arc is pointed from $i$ to $j$, we have $C_{ij}=1$.
Similarly, if the arc is from $j$ to $i$, i.e., $i$ is dominated by $j$, then
we assign $C_{ji}=-1$. The diagonal entries are considered as zero since there
is no own competition. Let
$\mathbf{p}=\big{(}p_{1},p_{2},p_{3},\cdots,p_{N}\big{)}^{T}$ be the strategy
vector with $\sum_{i=1}^{N}p_{i}=1$ and $p_{i}\geq 0$ $\forall$ $i$. To
calculate the optimal strategy, i.e., the non-zero solution $\mathbf{p}$, we
must have

$C\mathbf{p}=\mathbf{0}\\\ \implies detC=0.$ (39)

Since $C$ is a skew symmetric matrix, if $N$ is even, then $detC=1$ is a non-
zero perfect square. Hence, there does not exist any non-zero solution
$\mathbf{p}$ of (39). On the other hand, the determinant of a skew symmetric
matrix of odd order equals zero. So, calculation of the non-zero $\mathbf{p}$
with $\sum_{i=1}^{N}p_{i}=1$ and $p_{i}\geq 0$ $\forall$ $i$ is possible, only
if $N\geq 3$ is odd.

### Perturbation in single species for even number of species

We now consider the impact of perturbation on an even number of species. For
simplicity, we choose $N=4$. When $\alpha=0.25<0.5$, the influence of
$\beta>0$ is neutralized as the first species die out. The remaining species
coexist with equal densities as reflected through Fig. (10) (a). However,
$\alpha>0.5$ and $\beta>0$ reduces $x_{1}$ and $x_{2}$ similar to the effect
of external perturbations on triadic interaction as shown in Fig. (5) (b). We
choose $\beta=0.15$ and $\alpha=0.75$ in Fig. (10) (b). The tripple-wise
interaction in this case leads to the cessation of fourth species’ density.
$x_{4}$ converges to zero in the long run and positive $\beta$ acts like an
external perturbation on the remaining three species for $\alpha>0.5$. The
asymptotic values of $x_{1},x_{2},$ and $x_{3}$ can be predicted using the
theoretically calculated equilibrium in Eq. (14) due to the extinction of the
fourth species. The survivability of the first three species now becomes
equivalent to studying the system (4) with $(N-1)=3$ species.

## AUTHORS’ CONTRIBUTIONS

Both the authors Sourin Chatterjee and Sayantan Nag Chowdhury have contributed
equally to this work.
Sourin Chatterjee: Conceptualization, Methodology, Software, Validation,
Formal analysis, Investigation Visualization, Writing - review & editing.
Sayantan Nag Chowdhury: Methodology, Validation, Formal analysis,
Investigation, Writing, Visualization. Dibakar Ghosh: Validation,
Visualization, Writing - review & editing. Chittaranjan Hens:
Conceptualization, Methodology, Software, Validation, Formal analysis,
Investigation, Writing, Visualization.

## Competing interests.

We declare we have no competing interests.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the anonymous referees for their insightful
suggestions that helped considerably improve the manuscript. We are indebted
to Jacopo Grilli of The Abdus Salam International Centre for Theoretical
Physics for valuable discussions. Sayantan Nag Chowdhury would like to
acknowledge the CSIR (Project No. 09/093(0194)/2020-EMR-I) for financial
assistance. Sayantan Nag Chowdhury also wants to thank the Department of
Science and Technology, Govt. of India, for the financial support through
Grant No. NMICPS/006/MD/2020-21 during the end of this work. Chittaranjan Hens
is financially supported by the INSPIRE-Faculty grant (Code: IFA17-PH193).

## DATA AVAILABILITY

The data supporting this study’s findings are openly available on Github at
web .

## References

  * Grilli  _et al._ (2017) J. Grilli, G. Barabás, M. J. Michalska-Smith, and S. Allesina, Nature 548, 210 (2017).
  * Lotka (1920) A. J. Lotka, Proceedings of the National Academy of Sciences 6, 410 (1920).
  * Volterra (1926) V. Volterra, Animal Ecology , 409 (1926).
  * May and Leonard (1975) R. M. May and W. J. Leonard, SIAM Journal on Applied Mathematics 29, 243 (1975).
  * Nag Chowdhury  _et al._ (2021a) S. Nag Chowdhury, S. Kundu, M. Perc, and D. Ghosh, Proceedings of the Royal Society A 477, 20210397 (2021a).
  * Hubbell (2011) S. P. Hubbell,  _The unified neutral theory of biodiversity and biogeography (MPB-32)_ (Princeton University Press, 2011).
  * May (1974) R. M. May, Science 186, 645 (1974).
  * Nag Chowdhury  _et al._ (2021b) S. Nag Chowdhury, S. Kundu, J. Banerjee, M. Perc, and D. Ghosh, Journal of Theoretical Biology 518, 110606 (2021b).
  * May (2019) R. M. May,  _Stability and complexity in model ecosystems_ (Princeton university press, 2019).
  * Chesson (2000) P. Chesson, Annual review of Ecology and Systematics 31, 343 (2000).
  * Nag Chowdhury  _et al._ (2020) S. Nag Chowdhury, S. Kundu, M. Duh, M. Perc, and D. Ghosh, Entropy 22, 485 (2020).
  * Gao  _et al._ (2016) J. Gao, B. Barzel, and A.-L. Barabási, Nature 530, 307 (2016).
  * Hens  _et al._ (2019) C. Hens, U. Harush, S. Haber, R. Cohen, and B. Barzel, Nature Physics 15, 403 (2019).
  * Roy  _et al._ (2022) S. Roy, S. Nag Chowdhury, P. C. Mali, M. Perc, and D. Ghosh, PLOS ONE 17, 1 (2022).
  * Chase and Leibold (2009) J. M. Chase and M. A. Leibold,  _Ecological Niches_ (University of Chicago Press, 2009).
  * Soberón and Nakamura (2009) J. Soberón and M. Nakamura, Proceedings of the National Academy of Sciences 106, 19644 (2009).
  * Tilman (2004) D. Tilman, Proceedings of the National Academy of Sciences 101, 10854 (2004).
  * Bell (2000) G. Bell, The American Naturalist 155, 606 (2000).
  * Brokaw and Busing (2000) N. Brokaw and R. T. Busing, Trends in Ecology & Evolution 15, 183 (2000).
  * Wennekes  _et al._ (2012) P. L. Wennekes, J. Rosindell, and R. S. Etienne, Acta Biotheoretica 60, 257 (2012).
  * Chisholm and Pacala (2010) R. A. Chisholm and S. W. Pacala, Proceedings of the National Academy of Sciences 107, 15821 (2010).
  * Gravel  _et al._ (2006) D. Gravel, C. D. Canham, M. Beaudet, and C. Messier, Ecology Letters 9, 399 (2006).
  * McGill (2003) B. J. McGill, Nature 422, 881 (2003).
  * Clark and McLachlan (2003) J. S. Clark and J. S. McLachlan, Nature 423, 635 (2003).
  * Huisman and Weissing (1999) J. Huisman and F. J. Weissing, Nature 402, 407 (1999).
  * He  _et al._ (2011) Q. He, M. Mobilia, and U. C. Täuber, The European Physical Journal B 82, 97 (2011).
  * Szabó and Fath (2007) G. Szabó and G. Fath, Physics Reports 446, 97 (2007).
  * Berr  _et al._ (2009) M. Berr, T. Reichenbach, M. Schottenloher, and E. Frey, Physical Review Letters 102, 048102 (2009).
  * Sinervo and Lively (1996) B. Sinervo and C. M. Lively, Nature 380, 240 (1996).
  * Kerr  _et al._ (2002) B. Kerr, M. A. Riley, M. W. Feldman, and B. J. Bohannan, Nature 418, 171 (2002).
  * Laird and Schamp (2006) R. A. Laird and B. S. Schamp, The American Naturalist 168, 182 (2006).
  * Hibbing  _et al._ (2010) M. E. Hibbing, C. Fuqua, M. R. Parsek, and S. B. Peterson, Nature Reviews Microbiology 8, 15 (2010).
  * Gilpin (1975) M. E. Gilpin, The American Naturalist 109, 51 (1975).
  * Allesina and Levine (2011) S. Allesina and J. M. Levine, Proceedings of the National Academy of Sciences 108, 5638 (2011).
  * Bhattacharyya  _et al._ (2020) S. Bhattacharyya, P. Sinha, R. De, and C. Hens, Physical Review E 102, 012220 (2020).
  * Szolnoki  _et al._ (2014) A. Szolnoki, M. Mobilia, L.-L. Jiang, B. Szczesny, A. M. Rucklidge, and M. Perc, Journal of the Royal Society Interface 11, 20140735 (2014).
  * Reichenbach  _et al._ (2007) T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 (2007).
  * Islam  _et al._ (2022) S. Islam, A. Mondal, M. Mobilia, S. Bhattacharyya, and C. Hens, Physical Review E 105, 014215 (2022).
  * Park  _et al._ (2017) J. Park, Y. Do, B. Jang, and Y.-C. Lai, Scientific Reports 7, 1 (2017).
  * Shi  _et al._ (2010) H. Shi, W.-X. Wang, R. Yang, and Y.-C. Lai, Physical Review E 81, 030901 (2010).
  * Kirkup and Riley (2004) B. C. Kirkup and M. A. Riley, Nature 428, 412 (2004).
  * Soliveres  _et al._ (2015) S. Soliveres, F. T. Maestre, W. Ulrich, P. Manning, S. Boch, M. A. Bowker, D. Prati, M. Delgado-Baquerizo, J. L. Quero, I. Schöning, _et al._ , Ecology Letters 18, 790 (2015).
  * Cameron  _et al._ (2009) D. D. Cameron, A. White, and J. Antonovics, Journal of Ecology 97, 1311 (2009).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Physics Reports 874, 1 (2020).
  * Letten and Stouffer (2019) A. D. Letten and D. B. Stouffer, Ecology Letters 22, 423 (2019).
  * Mayfield and Stouffer (2017) M. M. Mayfield and D. B. Stouffer, Nature ecology & evolution 1, 1 (2017).
  * Levine  _et al._ (2017) J. M. Levine, J. Bascompte, P. B. Adler, and S. Allesina, Nature 546, 56 (2017).
  * Bairey  _et al._ (2016) E. Bairey, E. D. Kelsic, and R. Kishony, Nature Communications 7, 1 (2016).
  * Abrams (1983) P. A. Abrams, The American Naturalist 121, 887 (1983).
  * Lambiotte  _et al._ (2019) R. Lambiotte, M. Rosvall, and I. Scholtes, Nature physics 15, 313 (2019).
  * Stouffer  _et al._ (2018) D. B. Stouffer, C. E. Wainwright, T. Flanagan, and M. M. Mayfield, Journal of Ecology 106, 838 (2018).
  * Li  _et al._ (2020) Y. Li, D. Bearup, and J. Liao, Proceedings of the Royal Society B 287, 20201571 (2020).
  * Singh and Baruah (2021) P. Singh and G. Baruah, Theoretical Ecology 14, 71 (2021).
  * Abrams (2000) P. A. Abrams, Annual Review of Ecology and Systematics 31, 79 (2000).
  * Spaniel (2014) W. Spaniel,  _Game theory 101: the complete textbook_ (CreateSpace, 2014).
  * Fisher and Reeves (1995) D. C. Fisher and R. B. Reeves, Linear Algebra and its Applications 217, 83 (1995).
  * Jensen (1906) J. L. W. V. Jensen, Acta Mathematica 30, 175 (1906).
  * Dekking  _et al._ (2005) F. M. Dekking, C. Kraaikamp, H. P. Lopuhaä, and L. E. Meester,  _A Modern Introduction to Probability and Statistics: Understanding why and how_, Vol. 488 (Springer, 2005).
  * Ruel and Ayres (1999) J. J. Ruel and M. P. Ayres, Trends in Ecology & Evolution 14, 361 (1999).
  * Rao  _et al._ (2021) S. Rao, N. Muyinda, and B. De Baets, Scientific Reports 11 (2021), https://doi.org/10.1038/s41598-021-93438-8.
  * Connolly  _et al._ (2001) J. Connolly, P. Wayne, and F. A. Bazzaz, The American Naturalist 157, 107 (2001).
  * Majhi  _et al._ (2022) S. Majhi, M. Perc, and D. Ghosh, Journal of the Royal Society Interface 19, 20220043 (2022).
  * Buss and Jackson (1979) L. Buss and J. Jackson, The American Naturalist 113, 223 (1979).
  * Paquin and Adams (1983) C. E. Paquin and J. Adams, Nature 306, 368 (1983).
  * Friedman  _et al._ (2017) J. Friedman, L. M. Higgins, and J. Gore, Nature Ecology & Evolution 1 (2017), https://doi.org/10.1038/s41559-017-0109.
  * May (1972) R. M. May, Nature 238, 413 (1972).
  * Holling (1973) C. S. Holling, Annual Review of Ecology and Systematics 4, 1 (1973).
  * Allesina and Tang (2012) S. Allesina and S. Tang, Nature 483, 205 (2012).
  * (69) https://github.com/CRHENS/HOI-cyclic-interaction.