# Complex network near-synchronization for non-identical predator-prey systems

Guillaume Cantin, Cristiana Silva

► To cite this version:

Guillaume Cantin, Cristiana Silva. Complex network near-synchronization for non-identical predator-prey systems. AIMS Mathematics, 2022, 7 (11), pp.19975-19997. 10.3934/math.20221093. hal-03778221

HAL Id: hal-03778221

https://hal.science/hal-03778221v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research article

Complex network near-synchronization for non-identical predator-prey
systems

Guillaume Cantin<fnref n="1" /> and Cristiana J. Silva<fnref n="2" />,<fnref n="3" />,*

**Abstract:** In this paper, we analyze the properties of a complex network of predator-prey systems, modeling the ecological dynamics of interacting species living in a fragmented environment. We consider non-identical instances of a Lotka-Volterra model with Holling type II functional response, which undergoes a Hopf bifurcation, and focus on the possible synchronization of distinct local behaviours. We prove an original result for the near-synchronization of non-identical systems, which shows how to and to what extent an extinction dynamic can be driven to a persistence equilibrium. Our theoretical statements are illustrated by appropriate numerical simulations.

**Keywords:** complex network; non-identical systems; synchronization; ecology

**Mathematics Subject Classification:** 34A34, 34C60, 92B05

# 1. Introduction

It is now widely admitted by the scientific community that the increase of the human activity over the last century exerts a high pressure on the equilibrium of ecological systems, which can account for a major defaunation and a fast loss of biodiversity, often qualified as the “sixth extinction” [7, 12]. One of the main causes of this biological crisis is the degradation of the ecological environments of wildlife, which takes various forms, in the forefront of which are deforestation and habitat fragmentation [8, 10, 13]. Facing the challenge of restoring biodiversity, while maintaining human activity at a reasonable level, a carefully considered solution consists in the implementation of ecological corridors between each component of the fragmented environment, so as to increase the connectivity of natural habitats and to avoid local extinction of several wildlife species. For example, the restoration of aquatic ecosystems

<footnote n="1">¹ Laboratoire des Sciences du Numérique de Nantes, Nantes Université, CNRS UMR 6004, France; guillaume.cantin@univ-nantes.fr</footnote>

<footnote n="2">² Instituto Universitário de Lisboa (ISCTE-IUL), Av. das Forças Armadas, 1649-026, Lisbon, Portugal</footnote>

<footnote n="3">³ Center for Research and Development in Mathematics and Applications (CIDMA), Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal</footnote>

<footnote>* Correspondence: Email: cjoaosilva@ua.pt.</footnote>
by maintaining hydrological connectivity with riparian corridors is studied in [15]; different complexity and connectivity patterns are investigated at multiple scales in [21].

In this paper, our aim is precisely to propose an original mathematical model, in order to study these ecological concerns. Therefore, we consider a network of ecological systems, which is intended to model the complex dynamics of trophic chains in a degraded area. The complex network structure reproduces the heterogeneous natural environment, which is perturbed by fragmentation, by coupling several patches on which interacting wild species are living. On each patch, the ecological inter-species dynamics are modeled by a Lotka-Volterra predator-prey model with Holling type II functional response, which is able to describe several biological dynamics, such as extinction, coexistence or ecological cycles [9, 17, 19] (see also e.g. [3, 18, 23, 24] for other Lotka-Volterra type models). Next, the complex network is constructed so that each patch can admit its own dynamic. In this way, the local components of the network can for instance exhibit an extinction equilibrium on some places, whereas other places can present cycles. Furthermore, migrations of biological individuals in space, between each component of the fragmented environment, are taken into account by coupling the patches of the network, as schematized in Figure 1, where the disks model the patches of the fragmented habitat, and the arrows model the ecological corridors between these patches.

With this mathematical model in hand, we analyze the effect of the couplings on the local dynamics of each patch. In particular, we investigate the possibility to modify a local dynamic of extinction of the species, by increasing the couplings with patches on which persistence, with or without oscillations, is ensured. More generally, we search for sufficient conditions of *synchronization* of the local dynamics, under a variation of the couplings, which roughly means that a unique local dynamic is imposed on each patch. On this point, we establish a novel theorem of *near-synchronization*, which guarantees that the complex network remains in a neighborhood of a synchronization state, provided the coupling strength is strong enough, even if the local behaviors are non-identical.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of a complex network of predator-prey systems. The disks model the patches of a fragmented environment, where the inter-species dynamics of Lotka-Volterra type occur. The arrows model the ecological corridors which can be implemented between these patches, so as to increase the migrations in space of the species between each patch.</figcaption></figure>

Our choice to model the dynamics of interacting species living in a fragmented environment by a complex network structure follows a line of recent works. Indeed, synchronization in complex networks has been studied in a great number of papers, for various real-world applications, among them coupled oscillators, networks of chemical reactions, neural networks or meta-populations models (see for instance [1, 2, 5, 14] and the references therein). Recently, complex networks of Lotka-Volterra
models have been studied in [6]; however, the sufficient conditions of synchronization which have been established in this paper, correspond only to the particular case of identical dynamics. This state of the art highlights the main contributions of the present paper: it is the first time, at our knowledge, that a near-synchronization result for complex networks of non-identical systems is established at a theoretical level.

Our paper is organized as follows. In the next Section, we show how to construct a complex network of predator-prey systems, stemming from a Lotka-Volterra model embedded with a Holling type II functional response. In Section 3, we establish our main theoretical result, with Theorem 2, which establishes sufficient conditions of near-synchronization in a network of non-identical systems. Finally, in Section 4, we illustrate our theoretical statements by relevant numerical simulations.

## 2. Setting of the complex network of predator-prey models

In this section, we present the construction of a complex network of Lotka-Volterra systems, which describes the dynamics of interacting species living in a fragmented environment.

### 2.1. Lotka-Volterra predator-prey model with Holling type II functional response

Let us consider a biological environment in which two species interact. We assume that the densities of the species are determined by a predator-prey model of Lotka-Volterra type, which can be written by:

$$ egin{cases} \dot{x} = rx(1-x) - \frac{cxy}{\alpha + x}, \\ \dot{y} = -dy + \frac{cxy}{\alpha + x}. \end{cases} \qquad (2.1) $$

Here, $x$ and $y$ denote the prey and predator density, respectively; $\dot{x}$ and $\dot{y}$ denote their derivatives with respect to the time variable $t$. The parameters $r, c, d$ and $\alpha$ are positive coefficients; $r$ is the birth rate of the preys, $d$ is the mortality rate of predators, and $c, \alpha$ determine the non-linear interaction between preys and predators. As mentioned in our introduction, the dynamics of the predator-prey system (2.1) have been widely studied (see for instance [11]). Depending on the values of the parameters $r, c, d, \alpha$, the solutions of system (2.1) can be attracted to a coexistence equilibrium, to an extinction equilibrium or to a limit cycle. The extinction equilibrium is denoted $E_0 = (0, 0)$. Under the parameter conditions $c > d$ and $\alpha < \frac{c-d}{d}$, system (2.1) admits a coexistence equilibrium $E_1$, which implies persistence of each species, given by

$$ E_1 = \left( \frac{\alpha d}{c-d}, \frac{r\alpha}{c-d} \left( 1 - \frac{\alpha d}{c-d} \right) \right). \qquad (2.2) $$

System (2.1) also admits the equilibrium $E_2 = (1, 0)$. Let us introduce the critical value $\alpha_0$ given by

$$ \alpha_0 = \frac{c-d}{c+d}. \qquad (2.3) $$

It is well-known (see for instance [11], Chapter 3 or [4], Section 3.4.1) that system (2.1) undergoes a Hopf bifurcation at $\alpha = \alpha_0$. For $\alpha < \alpha_0$, a stable limit cycle bifurcates from the persistence equilibrium $E_1$. This Hopf bifurcation is illustrated in Figure 2. Therefore, for $\alpha$ small enough, system (2.1) presents oscillations, which are interpreted as healthy ecological cycles.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Hopf bifurcation occurring in the predator-prey model (2.1) when $\alpha$ crosses the critical value $\alpha_0 = \frac{c-d}{c+d}$. For $\alpha > \alpha_0$, the orbits converge to a stable focus. For $\alpha < \alpha_0$, the system admits a stable limit cycle.</figcaption></figure>

## 2.2. Complex network of predator-prey models for a fragmented environment

Next, we assume that the geographical habitat of the species is perturbed by the anthropic extension, so that it is fragmented in several patches. This fragmentation is likely to alter the equilibrium of the ecological system. In order to model such a fragmented environment, we construct a complex network of predator-prey models as follows.

First, let $n > 0$ denote the number of patches on the fragmented environment. On each patch $i \in \{1, \dots, n\}$, we denote by $(x_i, y_i)$ the densities of preys and predators respectively. We assume that each patch $i \in \{1, \dots, n\}$ can be connected to other patches and we denote by $N_i \subset \{1, \dots, n\}$ the set of patches which are connected to patch $i$. We assume that migrations of biological individuals can occur between two connected patches, at rates $\sigma_1$ for preys and $\sigma_2$ for predators. In this way, the dynamics of the fragmented environment are determined by the following complex network:

$$ \left\{ \begin{aligned} \dot{x}_i &= r_i x_i (1 - x_i) - \frac{c_i x_i y_i}{\alpha_i + x_i} - \sigma_1 \sum_{j \in N_i} (x_i - x_j), \\ \dot{y}_i &= -d_i y_i + \frac{c_i x_i y_i}{\alpha_i + x_i} - \sigma_2 \sum_{j \in N_i} (y_i - y_j), \end{aligned} \right. \qquad (2.4) $$

for $1 \le i \le n$, with $\sigma_1 \ge 0$ and $\sigma_2 \ge 0$. Note that the index $i$ ranges over the set $\{1, \dots, n\}$, whereas the subscript $j$, in the two sums which determine the couplings of the network, ranges of the set $N_i$ of the neighbors of $i$.

We emphasize that the parameters $r_i$, $c_i$, $d_i$, $\alpha_i$ can differ from one patch to another, which means that the ecological dynamics are non-identical within the fragmented environment. For instance, some patches could present limit cycles, whereas other patches could exhibit an extinction of both species. Note also that the couplings are symmetric, which means that if the species $x_i, y_i$ of patch $i$ can move towards some patch $j$, then the species $x_j, y_j$ of patch $j$ can conversely move towards patch $i$.

One remarkable case of fragmented environment is that of a complete graph topology, for which we have $N_i = \{1, \dots, n\} \setminus \{i\}$; this situation means that each patch is connected to all other patches. At the opposite, if the coupling parameters $\sigma_1, \sigma_2$ are equal to 0, then no migration of individuals occur in the
network.

Let us now introduce some notations. Let $X = ((x_1, y_1), \dots, (x_n, y_n))^\top \in \mathbb{R}^{2n}$. For each $i \in \{1, \dots, n\}$, we denote

$$
\begin{align}
\lambda_i &= (r_i, c_i, d_i, \alpha_i)^\top \in \mathbb{R}^4, \nonumber \\
f_1(x_i, y_i, \lambda_i) &= r_i x_i (1 - x_i) - \frac{c_i x_i y_i}{\alpha_i + x_i}, \nonumber \\
f_2(x_i, y_i, \lambda_i) &= -d_i y_i + \frac{c_i x_i y_i}{\alpha_i + x_i}, \tag{2.5} \\
g_1(x_i, X, \sigma_1) &= -\sigma_1 \sum_{j \in \mathcal{N}_i} (x_i - x_j), \nonumber \\
g_2(y_i, X, \sigma_2) &= -\sigma_2 \sum_{j \in \mathcal{N}_i} (y_i - y_j). \nonumber
\end{align}
$$

We also denote $\sigma = (\sigma_1, \sigma_2)^\top \in \mathbb{R}^2$ and

$$
\begin{align}
\Lambda &= (\lambda_1, \dots, \lambda_n)^\top \in \mathbb{R}^{4n}, \nonumber \\
F(X, \Lambda) &= (f_1(x_1, y_1, \lambda_1), f_2(x_1, y_1, \lambda_1), \dots, f_1(x_n, y_n, \lambda_n), f_2(x_n, y_n, \lambda_n))^\top \in \mathbb{R}^{2n}, \tag{2.6} \\
G(X, \sigma) &= (g_1(x_1, X, \sigma_1), g_2(y_1, X, \sigma_2), \dots, g_1(x_n, X, \sigma_1), g_2(y_n, X, \sigma_2))^\top \in \mathbb{R}^{2n}. \nonumber
\end{align}
$$

With these notations, the complex network (2.4) can be written under the following short form

$$
\dot{X} = F(X, \Lambda) + G(X, \sigma). \tag{2.7}
$$

Our first result guarantees that the complex network (2.7) admits global solutions, whose components
are non-negative and bounded.

**Theorem 1.** Let $X_0 \in (\mathbb{R}^+)^{2n}$. Then the complex network problem determined by equation (2.7) and $X(0) = X_0$ admits a unique global solution $X(t, X_0)$ defined on $[0, +\infty)$, whose components are non-negative.

Furthermore, the flow induced by equation (2.7) admits a positively invariant region $\Theta$ which is compact in $(\mathbb{R}^+)^{2n}$.

Before giving the proof of Theorem 1, we recall a Gronwall type inequality, whose proof can be
found in [22] for instance.

**Lemma 1.** Let $v$ denote a continuously differentiable function defined on $[0, T]$, with $T > 0$. Assume that $v$ satisfies

$$
\dot{v}(t) + \delta v(t) \leq \gamma,
$$

for all $t \in [0, T]$, with $\delta > 0$ and $\gamma > 0$. Then we have

$$
v(t) \leq \left[ v(0) - \frac{\gamma}{\delta} \right] e^{-\delta t} + \frac{\gamma}{\delta}, \quad t \in [0, T].
$$

We continue with the proof of Theorem 1.
*Proof of Theorem 1.* Let $X_0 \in (\mathbb{R}^+)^{2n}$. Standard results of the theory of differential equations (see for instance [16]) guaranty that the Cauchy problem determined by equation (2.7) and the initial condition $X(0) = X_0$ admits a unique local solution, which we denote $X(t, X_0)$, defined on $[0, T]$, with $T > 0$.

Let us now prove that the components of the local solution $X(t, X_0)$ are non-negative on $[0, T]$. First, we recall that the initial condition $X_0$ belongs to $(\mathbb{R}^+)^{2n}$. Next, we observe that the equations of system (2.4) can be rewritten

$$ \begin{cases} \dot{x}_i = x_i\phi_1(x_i, y_i) + \sigma_1 \sum_{j \in \mathcal{N}_i} x_j, \\ \dot{y}_i = y_i\phi_2(x, i, y_i) + \sigma_2 \sum_{j \in \mathcal{N}_i} y_j, \end{cases} $$

with $1 \le i \le n$ and

$$ \phi_1(x_i, y_i) = x_i(1 - x_i) - \frac{c_i y_i}{\alpha_i + x_i} - \sigma_1 |\mathcal{N}_i|, \quad \phi_2(x_i, y_i) = -d_i y_i + \frac{c_i y_i}{\alpha_i + x_i} - \sigma_2 |\mathcal{N}_i|, $$

where $|\mathcal{N}_i|$ denotes the cardinal of the finite set $\mathcal{N}_i$. Hence, by virtue of Proposition A.17 in [20], the components $(x_i, y_i)_{1 \le i \le n}$ satisfy $x_i(t) \ge 0, y_i(t) \ge 0$, for $t \in [0, T]$.

Finally, let us prove that the flow induced by equation (2.7) admits a positively invariant region $\Theta$, which is compact in $(\mathbb{R}^+)^{2n}$; as a consequence, every local solution $X(t, X_0)$ will be global in time. We first remark that for each $i \in \{1, \dots, n\}$, positive coefficients $a_i, b_i$ can be found such that

$$ r_i s(1-s) \le a_i - b_i s, $$

for all $s \in \mathbb{R}$. Afterwards, we introduce

$$ d_0 = \min_{1 \le i \le n} d_i, \quad a_0 = \sum_{i=1}^{n} a_i, \quad b_0 = \min_{1 \le i \le n} b_i, \quad c_0 = \min(b_0, d_0), \qquad (2.8) $$

and we denote by $P$ the total population of preys and predators in the complex network:

$$ P(t) = \sum_{i=1}^{n} [x_i(t) + y_i(t)], \quad t \in [0, T]. \qquad (2.9) $$

Since the couplings are symmetric, the sum over $i \in \{1, \dots, n\}$ of all equations of system (2.4) leads to

$$ \dot{P}(t) = \sum_{i=1}^{n} r_i x_i(t) [1 - x_i(t)] - \sum_{i=1}^{n} d_i y_i(t). $$

Next, we write

$$ - \sum_{i=1}^{n} d_i y_i(t) \le -d_0 \sum_{i=1}^{n} y_i(t), $$

and analogously

$$ \sum_{i=1}^{n} r_i x_i(t) [1 - x_i(t)] \le \sum_{i=1}^{n} [a_i - b_i x_i(t)] \le a_0 - b_0 \sum_{i=1}^{n} x_i(t). $$

We obtain

$$ \dot{P}(t) + c_0 P(t) \le a_0. $$
Applying Lemma 1 leads to

$$P(t) \leq \left[ P(0) - \frac{a_0}{c_0} \right] e^{-c_0 t} + \frac{a_0}{c_0},$$

which implies that the region $\Theta$ defined by

$$\Theta = \left\{ X = (x_i, y_i)_{1 \le i \le n} \in (\mathbb{R}^+)^{2n} ; \sum_{i=1}^{n} (x_i + y_i) \le \frac{a_0}{c_0} \right\} \quad (2.10)$$

is a positively invariant and compact region. The proof is complete. □

**Remark 1.** It is worth emphasizing that the bound $\frac{a_0}{c_0}$ of the positively invariant region $\Theta$ defined by (2.10) does not depend on the coupling parameters $\sigma_1, \sigma_2$. This uniform bound will be of great interest in Section 3 for studying the global dynamics of the complex network (2.7).

Now that the solutions of the complex network (2.7) are proved to be well defined and global in time, we give the definition of synchronization.

**Definition 1.** Let $i, j \in \{1, \dots, n\}$ such that $i \neq j$. We say that the patches $i$ and $j$ of the complex network (2.7) synchronize in $\Theta$ if, for any initial condition $X_0 \in \Theta$, the solution of (2.7) starting from $X_0$ satisfies

$$\lim_{t \to +\infty} \left( |x_i(t) - x_j(t)|^2 + |y_i(t) - y_j(t)|^2 \right) = 0.$$

We say that the complex network (2.7) synchronizes in $\Theta$ if every pair $(i, j)$ of patches synchronizes in $\Theta$.

**Remark 2.** Note that synchronization can occur in the complex network (2.7) without imposing any particular asymptotic dynamics; for example, the complex network could synchronize towards a global dynamic of extinction, towards a global dynamic of coexistence, or towards a global dynamic of limit cycles (it could even happen that a new dynamic emerges from the complex network structure). We will establish in Section 3 sufficient conditions of synchronization and discover that a complex network of non-identical instances of the predator-prey model (2.1) is likely not to fully synchronize.

### 3. Main result: near-synchronization of the complex network predator-prey model

In this section, we search for sufficient conditions of synchronization in the complex network (2.7). In a great number of papers, sufficient conditions of synchronization are established in the case of complex networks of identical systems. Here, we investigate the case of non-identical systems. We prove a novel theorem which shows that a complex network of non-identical systems can reach a near-synchronization state. For the sake of simplicity, we assume that the graph underlying the complex network (2.7) is complete, that is, each patch is connected to all other patches; equivalently, we have $N_i = \{1, \dots, n\} \setminus \{i\}$ for $1 \le i \le n$, where $N_i$ denotes the finite set of patches which are connected to patch $i$. For all $i, j \in \{1, \dots, n\}$, we introduce the energy functions $u_{i,j}$ defined along the trajectories of the complex network by

$$u_{i,j}(t) = \frac{1}{2} \left[ |x_i(t) - x_j(t)|^2 + |y_i(t) - y_j(t)|^2 \right], \quad (3.1)$$

and for $\lambda_i = (r_i, d_i, c_i, \alpha_i)$, $\lambda_j = (r_j, d_j, c_j, \alpha_j) \in \mathbb{R}^4$, we denote

$$\|\lambda_i - \lambda_j\|_{\infty} = \max \{ |r_i - r_j|, |d_i - d_j|, |c_i - c_j|, |\alpha_i - \alpha_j| \}.$$
The next theorem establishes an estimate of the energy functions $u_{i,j}$ defined by (3.1).

**Theorem 2.** There exist positive constants $\eta, \delta$ such that, for any initial condition $X_0 \in \Theta$, the solution of the complex network (2.7), starting from $X_0$, satisfies

$$ \dot{u}_{i,j}(t) \leq \eta \|\lambda_i - \lambda_j\|_{\infty} u_{i,j}^{1/2}(t) + [\delta - 2(n-1)\tilde{\sigma}] u_{i,j}(t), \quad t > 0, \quad (3.2) $$

where $\tilde{\sigma} = \min\{\sigma_1, \sigma_2\}$.

Furthermore, the constants $\eta$ and $\delta$ do not depend on the coupling parameters $\sigma_1, \sigma_2$.

We begin with a technical lemma.

**Lemma 2.** There exist positive constants $k_r, 1 \le r \le 6$, such that the functions $f_1$ and $f_2$ defined in (2.5) satisfy

$$ \begin{aligned} |f_1(x_i, y_i, \lambda_i) - f_1(x_j, y_j, \lambda_j)| &\le k_1 \|\lambda_i - \lambda_j\|_{\infty} + k_2 |x_i - x_j| + k_3 |y_i - y_j|, \\ |f_2(x_i, y_i, \lambda_i) - f_2(x_j, y_j, \lambda_j)| &\le k_4 \|\lambda_i - \lambda_j\|_{\infty} + k_5 |x_i - x_j| + k_6 |y_i - y_j|, \end{aligned} \quad (3.3) $$

for all $(x_i, y_i)_{1 \le i \le n}$, in the invariant region $\Theta$ defined by (2.10) and for all $i, j \in \{1, \dots, n\}$.

Furthermore, the constants $k_r, 1 \le r \le 6$, do not depend on the coupling parameters $\sigma_1, \sigma_2$.

**Proof of Lemma 2.** Since $(x_i, y_i)_{1 \le i \le n}$ belongs to the invariant region $\Theta$ defined by (2.10), there exists a positive constant $K$ such that

$$ |x_i| \le K, \quad |y_i| \le K, $$

for all $i \in \{1, \dots, n\}$.

In order to establish the estimates (3.3), we first write

$$ |d_i y_i - d_j y_j| \le |d_i y_i - d_i y_j| + |d_j y_i - d_j y_j|, $$

which leads to

$$ |d_i y_i - d_j y_j| \le d^+ |y_i - y_j| + K |d_i - d_j|, \quad (3.4) $$

where $d^+ = \max_{1 \le i \le n} d_i$. Similarly, we have

$$ \begin{aligned} |r_i x_i(1 - x_i) - r_j x_j(1 - x_j)| &\le |r_i x_i(1 - x_i) - r_i x_j(1 - x_j)| + |r_i x_j(1 - x_j) - r_j x_i(1 - x_i)| \\ &\le r_i |x_i(1 - x_i) - x_j(1 - x_j)| + K(1 + K) |r_i - r_j|. \end{aligned} $$

Now we write

$$ \begin{aligned} |x_i(1 - x_i) - x_j(1 - x_j)| &\le |x_i - x_j| + |x_i^2 - x_j^2| \\ &\le |x_i - x_j| + |x_i + x_j| \times |x_i - x_j| \\ &\le (1 + 2K) |x_i - x_j|, \end{aligned} $$

from which it follows that

$$ |r_i x_i(1 - x_i) - r_j x_j(1 - x_j)| \le r^+(1 + 2K) |x_i - x_j| + K(1 + K) |r_i - r_j|, \quad (3.5) $$
where $r^+ = \max_{1 \le i \le n} r_i$. Finally, we compute

$$ \left| \frac{c_i x_i y_i}{\alpha_i + x_i} - \frac{c_j x_j y_j}{\alpha_j + x_j} \right| \le \left| \frac{(\alpha_j + x_j) c_i x_i y_i - (\alpha_i + x_i) c_j x_j y_j}{(\alpha_i + x_i)(\alpha_j + x_j)} \right| \\ \le \frac{1}{\alpha_i \alpha_j} \left[ |\alpha_j c_i x_i y_i - \alpha_i c_j x_j y_j| + |c_i x_i x_j y_i - c_j x_i x_j y_j| \right]. $$

For the first term in the brackets, we write

$$ \begin{aligned} |\alpha_j c_i x_i y_i - \alpha_i c_j x_j y_j| & \le |\alpha_j c_i x_i y_i - \alpha_i c_j x_j y_j| + |\alpha_j c_i x_j y_j - \alpha_i c_j x_j y_j| \\ & \le \alpha_j c_i |x_i y_i - x_j y_j| + K^2 |\alpha_j c_i - \alpha_i c_j| \\ & \le \alpha_j c_i |x_i y_i - x_i y_j| + \alpha_j c_i |x_i y_j - x_j y_j| + K^2 |\alpha_j c_i - \alpha_i c_j| \\ & \le \alpha_j c_i K (|y_i - y_j| + |x_i - x_j|) + K^2 \alpha_j |c_i - c_j| + K^2 c_j |\alpha_i - \alpha_j|. \end{aligned} $$

For the second term in the brackets, we have

$$ \begin{aligned} |c_i x_i x_j y_i - c_j x_i x_j y_j| & \le |c_i x_i x_j y_i - c_i x_i x_j y_j| + |c_i x_i x_j y_j - c_j x_i x_j y_j| \\ & \le c_i K^2 |y_i - y_j| + K^3 |c_i - c_j|, \end{aligned} $$

which leads to

$$ \left| \frac{c_i x_i y_i}{\alpha_i + x_i} - \frac{c_j x_j y_j}{\alpha_j + x_j} \right| \le \frac{\alpha^+ c^+ K}{(\alpha^-)^2} |x_i - x_j| + \frac{\alpha^+ c^+ K + c^+ K^2}{(\alpha^-)^2} |y_i - y_j| \\ + \frac{K^2 \alpha^+ + K^3}{(\alpha^-)^2} |c_i - c_j| + \frac{c^+ K^2}{(\alpha^-)^2} |\alpha_i - \alpha_j|, \quad (3.6) $$

where $\alpha^+ = \max_{1 \le i \le n} \alpha_i$, $c^+ = \max_{1 \le i \le n} c_i$ and $\alpha^- = \min_{1 \le i \le n} \alpha_i$.

Gathering inequalities (3.5) and (3.6) leads to

$$ |f_1(x_i, y_i, \lambda_i) - f_1(x_j, y_j, \lambda_j)| \le k_1 \|\lambda_i - \lambda_j\|_{\infty} + k_2 |x_i - x_j| + k_3 |y_i - y_j|, $$

whereas gathering inequalities (3.4) and (3.6) leads to

$$ |f_2(x_i, y_i, \lambda_i) - f_2(x_j, y_j, \lambda_j)| \le k_4 \|\lambda_i - \lambda_j\|_{\infty} + k_5 |x_i - x_j| + k_6 |y_i - y_j|, $$

with positive constants $k_r$, $1 \le r \le 6$, which do not depend on $\sigma$ neither on the number $n$ of patches in the network. The proof of Lemma 2 is complete. □

We continue with the proof of Theorem 2.

**Proof of Theorem 2.** We compute the derivative of $u_{i,j}(t)$ along a trajectory $(x_i(t), y_i(t))_{1 \le i \le n}$ of the complex network (2.7), starting from an initial condition $(x_{i,0}, y_{i,0})_{1 \le i \le n} \in \Theta$. In order to lighten our notations, we omit the time dependence:

$$ \begin{aligned} \dot{u}_{i,j} &= (x_i - x_j)(\dot{x}_i - \dot{x}_j) + (y_i - y_j)(\dot{y}_i - \dot{y}_j) \\ &= (x_i - x_j)\biggl[f_1(x_i, y_i, \lambda_i) - \sigma_1 \sum_{k \ne i} (x_i - x_k) - f_1(x_j, y_j, \lambda_j) + \sigma_1 \sum_{k \ne j} (x_j - x_k)\biggr] \\ &\quad + (y_i - y_j)\biggl[f_2(x_i, y_i, \lambda_i) - \sigma_2 \sum_{k \ne i} (y_i - y_k) - f_2(x_j, y_j, \lambda_j) + \sigma_2 \sum_{k \ne j} (y_j - y_k)\biggr]. \end{aligned} $$
Now, we observe that

$$ \sigma_1 \sum_{k \neq i} (x_i - x_k) - \sigma_1 \sum_{k \neq j} (x_j - x_k) = \sigma_1 (n-1)(x_i - x_j), $$

$$ \sigma_2 \sum_{k \neq i} (y_i - y_k) - \sigma_2 \sum_{k \neq j} (y_j - y_k) = \sigma_2 (n-1)(y_i - y_j). $$

By virtue of Lemma 2, we obtain

$$ \begin{aligned} \dot{u}_{i,j} &\le |x_i - x_j| |k_1| \|\lambda_i - \lambda_j\|_{\infty} + k_2 |x_i - x_j| + k_3 |y_i - y_j| - \sigma_1(n-1)(x_i - x_j)^2 \\ &\quad + |y_i - y_j| |k_4| \|\lambda_i - \lambda_j\|_{\infty} + k_5 |x_i - x_j| + k_6 |y_i - y_j| - \sigma_2(n-1)(y_i - y_j)^2 \\ &\le \max\{k_1, k_4\} \|\lambda_i - \lambda_j\|_{\infty} \times [ |x_i - x_j| + |y_i - y_j| ] \\ &\quad + [k_2 - \sigma_1(n-1)] |x_i - x_j|^2 + [k_6 - \sigma_2(n-1)] |y_i - y_j|^2 \\ &\quad + (k_3 + k_5) |x_i - x_j| \times |y_i - y_j|. \end{aligned} $$

Next, we use the standard inequality $a + b \le \sqrt{2} \times \sqrt{a^2 + b^2}$, which is valid for all $a, b \in \mathbb{R}^+$, to write

$$ |x_i - x_j| + |y_i - y_j| \le \sqrt{2} \sqrt{|x_i - x_j|^2 + |y_i - y_j|^2} \le 2 \sqrt{u_{i,j}}, $$

thus we have

$$ \max\{k_1, k_4\} \|\lambda_i - \lambda_j\|_{\infty} \times [ |x_i - x_j| + |y_i - y_j| ] \le \eta \|\lambda_i - \lambda_j\|_{\infty} u_{i,j}^{1/2}, $$

with $\eta = 2 \times \max\{k_1, k_4\}$. In parallel, the Young inequality $a \times b \le \frac{a^2}{2} + \frac{b^2}{2}$ yields

$$ |x_i - x_j| \times |y_i - y_j| \le \frac{|x_i - x_j|^2}{2} + \frac{|y_i - y_j|^2}{2} \le u_{i,j}. $$

Finally, we set $\delta = 2 \max\{k_2, k_6\} + k_3 + k_5$ and obtain

$$ \dot{u}_{i,j} \le \eta \|\lambda_i - \lambda_j\|_{\infty} u_{i,j}^{1/2} + [\delta - 2(n-1)\tilde{\sigma}] u_{i,j}, $$

which completes the proof of Theorem 2. $\square$

Let us now discuss on some consequences of Theorem 2. We observe that estimate (3.2) is a differential inequality, which implies that for all $i, j \in \{1, \dots, n\}$, the energy function $u_{i,j}$ is bounded by the solution $w$ of the Bernoulli equation

$$ \dot{w} = \eta \|\lambda_i - \lambda_j\|_{\infty} w^{1/2} + [\delta - 2(n-1)\tilde{\sigma}] w. \quad (3.7) $$

If $\lambda_i = \lambda_j$, the latter differential equation can be simplified, so that estimate (3.2) becomes

$$ \dot{u}_{i,j}(t) \le [\delta - 2(n-1)\tilde{\sigma}] u_{i,j}(t), $$

which directly yields the following corollary.

**Corollary 1.** Assume that $\lambda_i = \lambda_j$ for some $i, j \in \{1, \dots, n\}$. Then the patches $i$ and $j$ synchronize if the following condition is fulfilled:

$$ (n-1)\tilde{\sigma} > \delta. \quad (3.8) $$
If $\lambda_i = \lambda_j$ for all $i, j \in \{1, \dots, n\}$, then obviously the whole network synchronizes under condition (3.8). Next, since the constant $\delta$ does not depend on the coupling parameters $\sigma_1, \sigma_2$, the sufficient condition (3.8) can easily be satisfied, provided the number $n$ of patches in the network is sufficiently large, or provided the minimum coupling strength $\tilde{\sigma} = \min\{\sigma_1, \sigma_2\}$ is sufficiently large.

From the ecological point of view, increasing the number $n$ of patches in the network would correspond to a worse fragmentation of the habitat, which is not a reasonable strategy for our purposes. However, increasing the minimum coupling strength $\tilde{\sigma}$ can be realized by providing wider ecological corridors.

The non trivial case of Theorem 2 corresponds to a complex network of non-identical patches, for which we have $\lambda_i \neq \lambda_j$ for at least one pair $(i, j) \in \{1, \dots, n\}^2$. In that case, the synchronization state $\{(x_i, y_i) = (x_j, y_j)\}$ is likely to present a soft loss of stability. Indeed, it is well-known that the solution $w$ of the Bernoulli equation (3.7) converges towards a positive limit given by

$$ \lim_{t \to +\infty} w(t) = \left( \frac{\eta \|\lambda_i - \lambda_j\|_{\infty}}{\delta - (n-1)\tilde{\sigma}} \right)^2, $$

provided $w(0) > 0$. We obtain the following corollary.

**Corollary 2.** *The energy function $u_{i,j}$ defined by (3.1), along with the solution of the complex network (2.7) starting, from $X_0 \in \Theta$, satisfies*

$$ 0 \le \limsup u_{i,j}(t) \le \left( \frac{\eta \|\lambda_i - \lambda_j\|_{\infty}}{\delta - (n-1)\tilde{\sigma}} \right)^2. \quad (3.9) $$

Finally, since the constants $\eta$ and $\delta$ do not depend on the coupling parameters $\sigma_1, \sigma_2$, then a sufficiently large value of the minimum coupling strength $\tilde{\sigma}$ ensures that any neighborhood of the synchronization state $\{(x_i, y_i) = (x_j, y_j)\}$ can be reached, which justifies the expression near-synchronization. This remark can be formalized by the following definition.

**Definition 2.** Let $i, j \in \{1, \dots, n\}$ such that $i \neq j$. We say that the patches $i$ and $j$ of the complex network (2.7) nearly synchronize in $\Theta$ with respect to $\tilde{\sigma}$ if, for any initial condition $X_0 \in \Theta$, and for any $\varepsilon > 0$, the solution of (2.7) starting from $X_0$ satisfies

$$ 0 \le \lim_{t \to +\infty} \left( |x_i(t) - x_j(t)|^2 + |y_i(t) - y_j(t)|^2 \right) < \varepsilon, $$

for $\tilde{\sigma}$ sufficiently large.

We say that the complex network (2.7) nearly synchronizes in $\Theta$ if every pair $(i, j)$ of patches nearly synchronizes in $\Theta$.

The discussion above can now be formulated as follows.

**Corollary 3.** *The complex network (2.7) nearly synchronizes in $\Theta$ with respect to the minimum coupling strength $\tilde{\sigma}$.*

The latter corollary means that a sufficiently large minimal coupling strength $\tilde{\sigma}$ ensures that the dynamics of each patch of the complex network (2.4) will be almost identical. Furthermore, by virtue of (3.9), if the minimal coupling strength $\tilde{\sigma}$ increases, or if the number of patches $n$ increases, then the near-synchronization state is getting closer to a synchronization state. In other words, $\tilde{\sigma}$ and $n$ can be viewed as parameters to control the level of near-synchronization in the network.
# 4. Numerical simulations

In this section, we provide two examples that illustrate the near-synchronization pattern of a complex network, under the effect of the couplings. We consider networks with two and ten patches, with different coupling strengths, so as to highlight various emergent dynamics. For the sake of simplicity, we assume that $\sigma_1 = \sigma_2$ and write $\sigma$ instead of $\sigma_1, \sigma_2$. Our complete computation codes are provided in the Appendix.

## 4.1. Near-synchronization in a two-patches network

We start by showing the effect of the coupling in a simple two-patches network. Although it may appear simple, the case of a two-patches network is fundamental, since it models a natural habitat which is divided by a single obstacle. Notably, the construction of a single road crossing a forest ecosystem exerts a high perturbation on the biodiversity hosted by this environment. Similarly, the establishment of a single river barrier or of a single maritime dyke profoundly modifies the behavior of a water hosted ecosystem.

Therefore, we consider the system given by

$$ egin{cases} 
\dot{x}_1 = r_1 x_1 (1 - x_1) - \frac{c_1 x_1 y_1}{\alpha_1 + x_1} - \sigma(x_1 - x_2), \\ 
\dot{y}_1 = -d_1 y_1 + \frac{c_1 x_1 y_1}{\alpha_1 + x_1} - \sigma(y_1 - y_2), \\ 
\dot{x}_2 = r_2 x_2 (1 - x_2) - \frac{c_2 x_2 y_2}{\alpha_2 + x_2} - \sigma(x_2 - x_1), \\ 
\dot{y}_2 = -d_2 y_2 + \frac{c_2 x_2 y_2}{\alpha_2 + x_2} - \sigma(y_2 - y_1), 
\end{cases} 
\quad (4.1) $$

with the parameter values from Table 1 and $\sigma > 0$. On patch 1, the parameters are chosen so that the system converges to a coexistence steady state ($\alpha = 0.5$) or to a limit cycle ($\alpha = 0.05$): the birth rate $r_1$ of preys and the death rate $d_1$ of predators are of the same order; in both cases, extinction is avoided. On patch 2, the parameters are chosen so that the system converges to the extinction state (it suffices to decrease the birth rate of preys: $r_1 = 0.2$).

<table><caption>**Table 1.** Parameter values for the two patches network given by system (4.1).</caption><thead><tr><th colspan="2">Patch 1</th><th colspan="2">Patch 2</th></tr><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>r<sub>1</sub></td><td>1</td><td>r<sub>2</sub></td><td>0.2</td></tr><tr><td>d<sub>1</sub></td><td>1</td><td>d<sub>2</sub></td><td>1</td></tr><tr><td>c<sub>1</sub></td><td>2</td><td>c<sub>2</sub></td><td>3</td></tr><tr><td>α<sub>1</sub></td><td>0.5 or 0.05</td><td>α<sub>2</sub></td><td>0.05</td></tr></tbody></table>

The initial conditions are defined by

$$ (x_1(0), y_1(0)) = (0.3, 0.2), \quad (x_2(0), y_2(0)) = (0.5, 0.3). $$

We show in Figure 3 the corresponding orbits of the two-patches complex network (4.1), in absence of coupling ($\sigma = 0$), and in Figure 4 the orbits for a positive coupling strength ($\sigma = 0.3$); two cases are
investigated: in the first case, we have $\alpha_1 = 0.5$, whereas in the second case, we have $\alpha_1 = 0.05$. Let us discuss the numerical results in regard of Theorem 2.

For $\sigma = 0$ and $\alpha_1 = 0.5$ (Figure 3(a)), the orbit $(x_1, y_1)$ of the first patch is attracted to a persistence equilibrium, while the orbit $(x_2, y_2)$ of the second patch is attracted to the extinction equilibrium. For $\sigma = 0$ and $\alpha_1 = 0.05$ (Figure 3(b)), the orbit of the second patch is the same, but the orbit $(x_1, y_1)$ of the first patch is now attracted to a limit cycle. We can predict, by virtue of Theorem 2, that the two-patches network will present a near-synchronization dynamic if the coupling strength $\sigma$ is sufficiently large. However, we aim to describe the common dynamic which is reached under this near-synchronization pattern. Indeed, for $\sigma = 0.3$, Figure 4 shows that the near-synchronization can hide various emergent dynamics that might occur. For $\sigma = 0.3$ and $\alpha_1 = 0.5$ (Figure 4(a)), the two patches nearly synchronize towards a persistence equilibrium, whereas for $\sigma = 0.3$ and $\alpha_1 = 0.05$ (Figure 4(b)), the two patches nearly synchronize towards a limit cycle. In both cases, the dynamics are obviously modified: the level of persistence is decreased on patch 1 for $\alpha_1 = 0.5$ (Figure 4(a)), as well as the amplitude of the oscillations on patch 1 for $\alpha_1 = 0.05$ (Figure 4(b)); however, the extinction on patch 2 is avoided. Overall, these results that show the good health of the ecosystem (persistence or ecological oscillations) can be recovered by the setting of a strong connection between the patches of the fragmented habitat.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Orbits of the two-patches network (4.1) in absence of coupling. (a) For $\sigma = 0$ and $\alpha_1 = 0.5$, the orbit $(x_1, y_1)$ of the first patch is attracted to a persistence equilibrium, while the orbit $(x_2, y_2)$ of the second patch is attracted to the extinction equilibrium. (b) For $\sigma = 0$ and $\alpha_1 = 0.05$, the orbit of the second patch is the same, but the orbit $(x_1, y_1)$ of the first patch is now attracted to a limit cycle.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>(b) $\sigma = 0.3$, $\alpha_1 = 0.05$</figcaption></figure>

<figure><img src="image_7.png" /></figure>

**Figure 4.** Orbits of the two-patches network (4.1) with a positive coupling strength. (a) For $\sigma = 0.3$ and $\alpha_1 = 0.5$, the two patches nearly synchronize towards a persistence equilibrium. (b) For $\sigma = 0.3$ and $\alpha_1 = 0.05$, the two patches nearly synchronize towards a limit cycle corresponding to ecological oscillations.

## 4.2. Synchronization of oscillatory behaviors in small networks

Finally, we aim to experiment an increase of the number $n$ of patches in the network. Therefore, we consider a ten-patches network with a complete graph topology; we choose randomly a set of parameters for each patch, such that

$$0 \le r_i \le 1, \quad 0 \le d_i \le 1, \quad 0 \le c_i \le 3, \quad 0.01 \le \alpha_i \le 0.41,$$

for each $i \in \{1, 2, \dots, 10\}$; the initial conditions are also chosen randomly in the interval $[0.1, 0.6]$. In this way, the patches of the complex network exhibit various local dynamics in absence of coupling. The corresponding orbits are shown in Figure 5; the corresponding time series are also provided in Figure 6. For instance, the orbit $(x_3, y_3)$ of patch 3 is attracted to a limit cycle; the orbit $(x_4, y_4)$ of patch 4 converges to the equilibrium $E_2(1, 0)$ (persistence of the preys and extinction of the predators); the orbit $(x_6, y_6)$ of patch 6 converges to the extinction equilibrium.

We depict in Figure 7 the orbits of the same ten-patches complex network, with a relatively large coupling strength ($\sigma = 1$); the corresponding time series are again provided in Figure 8. According to Theorem 2, we observe that the orbits nearly synchronize. The coupling strength is large enough to ensure that the near-synchronization is very closed to a synchronization state. Furthermore, the near-synchronization is characterized by a limit cycle; the amplitude of that limit cycle can slightly vary from one patch to another (for instance, $y_7$ is less than 0.7, while $y_8$ reaches 0.8).
From the ecological point of view, these results show again that healthy biological oscillations can be recovered by the setting of numerous and efficient connections between the patches of a fragmented environment.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Local dynamics in a ten-patches network in absence of coupling ($\sigma = 0$). The orbits of patches 1, 2, 3, 5, 9, 10 are attracted to a limit cycle; the orbits of patches 4, 7 converge to the equilibrium $E_2(1, 0)$ (persistence of the preys and extinction of the predators); the orbit of patch 6 converges to the extinction equilibrium.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Time series of a ten-patches network in absence of coupling ($\sigma = 0$). Depending on the values of the parameters, the solutions are attracted to a limit cycle, to the equilibrium $E_2(1, 0)$ (persistence of the preys and extinction of the predators) or to the extinction equilibrium.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Near-synchronization in a ten-patches network with a relatively large coupling strength ($\sigma = 1$). The non-identical dynamics are nearly synchronized towards a limit cycle.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Time series of a ten-patches network with a relatively large coupling strength ($\sigma = 1$), showing the near-synchronization of the local dynamics.</figcaption></figure>

**5. Conclusions and future work**

In this paper, we proposed a complex network to model a heterogeneous geographical habitat of species which is perturbed by an anthropic extension, being fragmented in several patches, where the fragmentation is likely to alter the equilibrium of the ecological system. The complex network was constructed by coupling several patches on which interacting wild species are living and where, for each patch, the ecological inter-species dynamics were modeled by a Lotka-Volterra predator-prey model with Holling type II functional response. An important feature of the complex network is that each patch can admit its own dynamic and migrations of biological individuals in space, between each component of the fragmented environment, are taken into account by coupling the patches of the network. We proved new sufficient conditions for the *near-synchronization* of the complex network, which guarantees that the complex network remains in a neighborhood of a synchronization state, provided the coupling strength is strong enough, even if the local behaviors are non-identical. This result allows us to modify
the local dynamic of extinction of the species, by increasing the couplings with patches on which persistence, with or without oscillations, is ensured.

When proving the sufficient conditions of synchronization, we discovered that a complex network of non-identical instances of the predator-prey model (2.1) is likely not to fully synchronize. This feature motivates the setting of an optimal control problem, so as to exert a command on the dynamics of the complex network (2.7) and to reach a synchronization state, even in the case of non-identical patches. As future work, we intend to apply optimal control theory to remedy the default of synchronization, where the coupling strengths will be represented by external controls acting on the dynamics of the network.

Another interesting research perspective corresponds to the fact that the dispersion of the species from one patch to another could be modeled at a refined level by adding time delays in the coupling terms; in particular, it will be natural to investigate the impact of time delays on the near synchronization dynamic that we have established in the present paper.

## Acknowledgments

This work was partially supported by Portuguese funds through CIDMA, The Center for Research and Development in Mathematics and Applications of University of Aveiro, and the Portuguese Foundation for Science and Technology (FCT–Fundação para a Ciência e a Tecnologia), within project UIDB/04106/2020. Silva was also supported by the FCT Researcher Program CEEC Individual 2018 with reference CEECIND/00564/2018.

## Conflict of interest

All authors declare no conflicts of interest in this paper.

## References

1. A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, C. Zhou, Synchronization in complex networks, *Physics reports*, **469** (2008), 93–153. https://doi.org/10.1016/j.physrep.2008.09.002

2. M. Barahona, L. M. Pecora, Synchronization in small-world systems, *Phys. rev. lett.*, **89** (2002), 054101. https://doi.org/10.1103/PhysRevLett.89.054101

3. D. Barman, J. Roy, S. Alam, Modelling hiding behaviour in a predator-prey system by both integer order and fractional order derivatives, *Ecol. Inform.*, **67** (2022), 101483. https://doi.org/10.1016/j.ecoinf.2021.101483

4. A. D. Bazykin, *Nonlinear dynamics of interacting populations*, World Scientific, 1998. https://doi.org/10.1142/2284

5. I. Belykh, M. Hasler, M. Lauret, H. Nijmeijer, Synchronization and graph topology, *Int. J. Bifurcat. Chaos*, **15** (2005), 3423–3433. https://doi.org/10.1142/S0218127405014143

6. G. Cantin, M. Aziz-Alaoui, Dimension estimate of attractors for complex networks of reaction-diffusion systems applied to an ecological model, *Commun. Pur. Appl. Anal.*, **20** (2021), 623. https://doi.org/10.3934/cpaa.2020283
7. R. Dirzo, H. S. Young, M. Galetti, G. Ceballos, N. J. Isaac, B. Collen, Defaunation in the Anthropocene, *Science*, **345** (2014), 401–406. https://doi.org/10.1126/science.1251817

8. N. M. Haddad, L. A. Brudvig, J. Clobert, K. F. Davies, A. Gonzalez, R. D. Holt, et al., Habitat fragmentation and its lasting impact on earth's ecosystems, *Sci. adv.*, **1** (2015), e1500052. https://doi.org/10.1126/sciadv.1500052

9. C. S. Holling, The functional response of predators to prey density and its role in mimicry and population regulation, *The Memoirs of the Entomological Society of Canada*, **97** (1965), 5–60. https://doi.org/10.4039/entm9745fv

10. C. N. Johnson, A. Balmford, B. W. Brook, J. C. Buettel, M. Galetti, L. Guangchun, et al., Biodiversity losses and conservation responses in the Anthropocene, *Science*, **356** (2017), 270–275. https://doi.org/10.1126/science.aam9317

11. Y. A. Kuznetsov, I. A. Kuznetsov, Y. Kuznetsov, *Elements of applied bifurcation theory*, volume 112. Springer, 1998.

12. R. E. Leakey, R. Lewin, The sixth extinction: patterns of life and the future of humankind, *J. Leisure Res.*, **29** (1997), 476. https://doi.org/10.1080/00222216.1997.11949812

13. B. J. McGill, M. Dornelas, N. J. Gotelli, A. E. Magurran, Fifteen forms of biodiversity trend in the Anthropocene, *Trends ecol. evol.*, **30** (2015), 104–113. https://doi.org/10.1016/j.tree.2014.11.006

14. A. Miranville, G. Cantin, M. Aziz-Alaoui, Bifurcations and synchronization in networks of unstable reaction-diffusion systems, *J. Nonlinear Sci.*, **31** (2021), 1–34. https://doi.org/10.1007/s00332-021-09701-9

15. R. J. Naiman, H. Decamps, M. Pollock, The role of riparian corridors in maintaining regional biodiversity, *Ecol. appl.*, **3** (1993), 209–212. https://doi.org/10.2307/1941822

16. L. Perko, *Differential equations and dynamical systems*, volume 7, Springer Science & Business Media, 2013.

17. L. A. Real, The kinetics of functional response, *The American Naturalist*, **111** (1977), 289–300. https://doi.org/10.1086/283161

18. J. Roy, D. Barman, S. Alam, Role of fear in a predator-prey system with ratio-dependent functional response in deterministic and stochastic environment, *Biosystems*, **197** (2020), 104176. https://doi.org/10.1016/j.biosystems.2020.104176

19. G. T. Skalski, J. F. Gilliam, Functional responses with predator interference: viable alternatives to the holling type ii model, *Ecology*, **82** (2001), 3083–3092. https://doi.org/10.1890/0012-9658(2001)082[3083:FRWPIV]2.0.CO;2

20. H. L. Smith, H. R. Thieme, *Dynamical systems and population persistence*, volume 118, American Mathematical Soc., 2011.

21. L. R. Tambosi, A. C. Martensen, M. C. Ribeiro, J. P. Metzger, A framework to optimize biodiversity restoration efforts based on habitat amount and landscape connectivity, *Restor. ecol.*, **22** (2014), 169–177. https://doi.org/10.1111/rec.12049

22. A. Yagi, *Abstract parabolic evolution equations and their applications*, Springer Science & Business Media, 2009. https://doi.org/10.1007/978-3-642-04631-5_4
23. R. Yang, C. Nie, D. Jin, Spatiotemporal dynamics induced by nonlocal competition in a diffusive predator-prey system with habitat complexity, *Nonlinear Dynam.*, (2022), 1–22.
https://doi.org/10.21203/rs.3.rs-1141642/v1

24. R. Yang, F. Wang, D. Jin, Spatially inhomogeneous bifurcating periodic solutions induced by nonlocal competition in a predator–prey system with additional food, *Math. Method. Appl. Sci.*, (2022). https://doi.org/10.1002/mma.8349

Appendix

In this appendix, we provide the computation codes of our numerical simulations. These codes are written in Python3 and require the scientific libraries ` matplotlib`, ` numpy` and ` scipy`.

Computation code for a two-patches network

#!/usr/bin/env python3

# Scientific libraries
from matplotlib import plot as plt
import numpy as np
from scipy.integrate import odeint
from random import random

# Parameters
r1 = 1; d1 = 1; c1 = 2; alpha1 = 0.05 # or 0.5
sigma = 0.3 # or 0
r2 = 0.2; d2 = 1; c2 = 3; alpha2 = 0.05

def lotka(X, t):
    x1, y1, x2, y2 = X
    dx1 = r1*x1*(1-x1) - c1*x1*y1/(alpha1+x1) - sigma*(x1-x2)
    dy1 = -d1*y1 + c1*x1*y1/(alpha1+x1) - sigma*(y1-y2)
    dx2 = r2*x2*(1-x2) - c2*x2*y2/(alpha2+x2) + sigma*(x1-x2)
    dy2 = -d2*y2 + c2*x2*y2/(alpha2+x2) + sigma*(y1-y2)
    return [dx1, dy1, dx2, dy2]

# Phase portrait
time = np.arange(0, 50, 0.01)
plt.figure()
x10 = 0.3; y10 = 0.2; x20 = 0.5; y20 = 0.3
orbit = odeint(lotka, [x10, y10, x20, y20], time)
x1, y1, x2, y2 = orbit.T
plt.plot(x1, y1, 'b', lw=0.5)
plt.plot(x2, y2, 'r', lw=1)
plt.show()
Computation code for a ten-patches network

#!/usr/bin/env python3

# Scientific libraries

from matplotlib import pyplot as plt

import numpy as np

from scipy.integrate import odeint

from random import random

# Parameters

r1 = random(); d1 = random(); c1 = 3*random(); alpha1 = 0.01 + 0.4*random()

r2 = random(); d2 = random(); c2 = 3*random(); alpha2 = 0.01 + 0.4*random()

r3 = random(); d3 = random(); c3 = 3*random(); alpha3 = 0.01 + 0.4*random()

r4 = random(); d4 = random(); c4 = 3*random(); alpha4 = 0.01 + 0.4*random()

r5 = random(); d5 = random(); c5 = 3*random(); alpha5 = 0.01 + 0.4*random()

r6 = random(); d6 = random(); c6 = 3*random(); alpha6 = 0.01 + 0.4*random()

r7 = random(); d7 = random(); c7 = 3*random(); alpha7 = 0.01 + 0.4*random()

r8 = random(); d8 = random(); c8 = 3*random(); alpha8 = 0.01 + 0.4*random()

r9 = random(); d9 = random(); c9 = 3*random(); alpha9 = 0.01 + 0.4*random()

r10 = random(); d10 = random(); c10 = 3*random(); alpha10 = 0.01 + 0.4*random()

def lotka(X, t):

x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, y6, x7, y7, x8, y8, x9, y9, x10, y10 = X

dx1 = r1*x1*(1-x1) - c1*x1*y1/(alpha1+x1) - sigma*(9*x1-x2-x3-x4-x5-x6-x7-x8-x9-x10)

dy1 = -d1*y1 + c1*x1*y1/(alpha1+x1) - sigma*(9*y1-y2-y3-y4-y5-y6-y7-y8-y9-y10)

dx2 = r2*x2*(1-x2) - c2*x2*y2/(alpha2+x2) - sigma*(9*x2-x1-x3-x4-x5-x6-x7-x8-x9-x10)

dy2 = -d2*y2 + c2*x2*y2/(alpha2+x2) - sigma*(9*y2-y1-y3-y4-y5-y6-y7-y8-y9-y10)

dx3 = r3*x3*(1-x3) - c3*x3*y3/(alpha3+x3) - sigma*(9*x3-x1-x2-x4-x5-x6-x7-x8-x9-x10)

dy3 = -d3*y3 + c3*x3*y3/(alpha3+x3) - sigma*(9*y3-y2-y1-y4-y5-y6-y7-y8-y9-y10)

dx4 = r4*x4*(1-x4) - c4*x4*y4/(alpha4+x4) - sigma*(9*x4-x1-x2-x3-x5-x6-x7-x8-x9-x10)

dy4 = -d4*y4 + c4*x4*y4/(alpha4+x4) - sigma*(9*y4-y2-y3-y1-y5-y6-y7-y8-y9-y10)

dx5 = r5*x5*(1-x5) - c5*x5*y5/(alpha5+x5) - sigma*(9*x5-x2-x3-x4-x1-x6-x7-x8-x9-x10)

dy5 = -d5*y5 + c5*x5*y5/(alpha5+x5) - sigma*(9*y5-y2-y3-y4-y1-y6-y7-y8-y9-y10)

dx6 = r6*x6*(1-x6) - c6*x6*y6/(alpha6+x6) - sigma*(9*x6-x2-x3-x4-x5-x1-x7-x8-x9-x10)

dy6 = -d6*y6 + c6*x6*y6/(alpha6+x6) - sigma*(9*y6-y2-y3-y4-y5-y1-y7-y8-y9-y10)

dx7 = r7*x7*(1-x7) - c7*x7*y7/(alpha7+x7) - sigma*(9*x7-x2-x3-x4-x5-x6-x1-x8-x9-x10)

dy7 = -d7*y7 + c7*x7*y7/(alpha7+x7) - sigma*(9*y7-y2-y3-y4-y5-y6-y1-y8-y9-y10)

dx8 = r8*x8*(1-x8) - c8*x8*y8/(alpha8+x8) - sigma*(9*x8-x2-x3-x4-x5-x6-x7-x1-x9-x10)

dy8 = -d8*y8 + c8*x8*y8/(alpha8+x8) - sigma*(9*y8-y2-y3-y4-y5-y6-y7-y1-y9-y10)

dx9 = r9*x9*(1-x9) - c9*x9*y9/(alpha9+x9) - sigma*(9*x9-x2-x3-x4-x5-x6-x7-x8-x1-x10)

dy9 = -d9*y9 + c9*x9*y9/(alpha9+x9) - sigma*(9*y9-y2-y3-y4-y5-y6-y7-y8-y1-y10)

dx10 = r10*x10*(1-x10) - c10*x10*y10/(alpha10+x10) - sigma*(9*x10-x2-x3-x4-x5-x6-x7-x8-x9-x1)

dy10 = -d10*y10 + c10*x10*y10/(alpha10+x10) - sigma*(9*y10-y2-y3-y4-y5-y6-y7-y8-y9-y1)
return [dx1, dy1, dx2, dy2, dx3, dy3, dx4, dy4, dx5, dy5,
dx6, dy6, dx7, dy7, dx8, dy8, dx9, dy9, dx10, dy10]

# Phase portrait

time = np.arange(0, 200, 0.01)
plt.figure()
X0 = [0.1 + 0.5*random() for k in range(20)]
sigma = 0 # or 1
orbit = odeint(lotka, X0, time)
x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, y6, x7, y7, x8, y8, x9, y9, x10, y10 = orbit.T
plt.plot(x1, y1)
plt.plot(x2, y2)
plt.plot(x3, y3)
plt.plot(x4, y4)
plt.plot(x5, y5)
plt.plot(x6, y6)
plt.plot(x7, y7)
plt.plot(x8, y8)
plt.plot(x9, y9)
plt.plot(x10, y10)
plt.show()

<figure><img src="image_12.png" /></figure>

AIMS Press

©2022 the Author(s), licensee AIMS Press. This
is an open access article distributed under the
terms of the Creative Commons Attribution License
(http://creativecommons.org/licenses/by/4.0)