On the nonstandard numerical discretization of SIR epidemic model with a saturated incidence rate and vaccination

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 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.

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:

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

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

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.

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 \lbrace 1, \dots, n \rbrace$, we denote by $(x_i, y_i)$ the densities of preys and predators respectively. We assume that each patch $i \in \lbrace 1, \dots, n \rbrace$ can be connected to other patches and we denote by $\mathcal{N}_i \subset \lbrace 1, \dots, n \rbrace$ 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:

for $1 \leq i \leq n$, with $\sigma_1 \geq 0$ and $\sigma_2 \geq 0$. Note that the index $i$ ranges over the set $\lbrace 1, \dots, n \rbrace$, whereas the subscript $j$, in the two sums which determine the couplings of the network, ranges of the set $\mathcal{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 $\mathcal{N}_i = \lbrace 1, \dots, n \rbrace \setminus \lbrace i \rbrace$; 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 = \big((x_1, y_1), \dots, (x_n, y_n)\big)^\top \in \mathbb{R}^{2n}$. For each $i \in \lbrace 1, \dots, n \rbrace$, we denote

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

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

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

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

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

with $1 \leq i \leq n$ and

where $\left\vert{ \mathcal{N}_i}\right\vert$ 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 \leq i \leq n}$ satisfy $x_i(t) \geq 0$, $y_i(t) \geq 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 \lbrace 1, \dots, n \rbrace$, positive coefficients $a_i$, $b_i$ can be found such that

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

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

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

Next, we write

and analogously

We obtain

Applying Lemma 1 leads to

which implies that the region $\Theta$ defined by

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 \lbrace 1, \dots, n \rbrace$ 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

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 $\mathcal{N}_i = \lbrace 1, \dots, n \rbrace \setminus \lbrace i \rbrace$ for $1 \leq i \leq n$, where $\mathcal{N}_i$ denotes the finite set of patches which are connected to patch $i$. For all $i, j \in \lbrace 1, \dots, n \rbrace$, we introduce the energy functions $u_{i, j}$ defined along the trajectories of the complex network by

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

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

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

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 \leq r \leq 6$, such that the functions $f_1$ and $f_2$ defined in (2.5) satisfy

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

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

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

for all $i \in \lbrace 1, \dots, n \rbrace$.

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

which leads to

where $d^+ = \max\limits_{1 \leq i \leq n} d_i$. Similarly, we have

Now we write

from which it follows that

where $r^+ = \max\limits_{1 \leq i \leq n} r_i$. Finally, we compute

For the first term in the brackets, we write

For the second term in the brackets, we have

which leads to

where $\alpha^+ = \max\limits_{1\leq i \leq n} \alpha_i$, $c^+ = \max\limits_{1\leq i \leq n} c_i$ and $\alpha^- = \min\limits_{1\leq i \leq n} \alpha_i$.

Gathering inequalities (3.5) and (3.6) leads to

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

with positive constants $k_r$, $1 \leq r \leq 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 $\big(x_i(t), y_i(t)\big)_{1 \leq i \leq n}$ of the complex network (2.7), starting from an initial condition $\big(x_{i, 0}, y_{i, 0}\big)_{1 \leq i \leq n} \in \Theta$. In order to lighten our notations, we omit the time dependence:

Now, we observe that

By virtue of Lemma 2, we obtain

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

thus we have

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

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

which completes the proof of Theorem 2.

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 \lbrace 1, \dots, n \rbrace$, the energy function $u_{i, j}$ is bounded by the solution $w$ of the Bernoulli equation

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

which directly yields the following corollary.

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

If $\lambda_i = \lambda_j$ for all $i, j \in \lbrace 1, \dots, n \rbrace$, 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 \lbrace \sigma_1, \sigma_2 \rbrace$ 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 \lbrace 1, \dots, n \rbrace^2$. In that case, the synchronization state $\lbrace(x_i, y_i) = (x_j, y_j)\rbrace$ 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

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

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 $\lbrace(x_i, y_i) = (x_j, y_j)\rbrace$ can be reached, which justifies the expression near-synchronization. This remark can be formalized by the following definition.

Definition 2. Let $i, j \in \lbrace 1, \dots, n \rbrace$ 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

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

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$).

The initial conditions are defined by

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.

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

for each $i \in \lbrace 1, 2, \dots, 10 \rbrace$; 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 show 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.

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.

Appendix

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

Computation code for a two-patches network

Computation code for a ten-patches network