Pattern Formation in a Predator-Prey Model with Allee Effect and Hyperbolic Mortality on Multiplex Networks

1. Introduction

As early as around 1925, Lotka and Volterra proposed a mathematical model to describe the population evolution of snow rabbits (prey) and lynx (predators), known as the Lotka-Volterra model [1]. Since then, researchers have conducted in-depth research and improvement on the model based on real-life backgrounds [2,3,4,5,6,7,8,9,10,11,12]. What interests us is that in 2012, Nagano and Maeda introduced diffusion terms into the classical Rosenzweig-MacArthur model and analyzed the reaction-diffusion Rosenzweig-MacArthur model [3,13]. They provided a phase diagram of the model and studied its lattice formation, revealing the existence of various spatiotemporal patterns, including well-known spiral and lattice patterns. In 2014, Zhang et al. pointed out that the model did not have the possibility of generating Turing pattern based on this, so they attempted to consider more complex inducing factors such as hyperbolic mortality [4] or delay [14]. They also displayed many beautiful spatial patterns, including hexagonal spot patterns, stripe patterns, micro spirals patterns, and target patterns. Liu and his colleagues have long been engaged in the study of the Allee effect on population model dynamics [5,6,7,8,11,12]. Therefore, in 2019, they attempted to introduce the predator-prey model with hyperbolic mortality proposed by Zhang et al. into the Allee effect. The study found that the Allee effect would increase the independence of pattern, in other words, it would make the spatial distribution of the population more concentrated, which is beneficial for the continuation of the population [7]. In 2010, the Turing pattern on large-scale random networks was first studied by Nakao and Mikhailov, and their findings demonstrated that early linear instability caused network nodes to naturally separate into wealthy and poor groups. In the following nonlinear stage, the freshly forming Turing pattern undergoes yet another significant reshaping. Several steady-state coexistence and hysteresis effects were observed [2]. Inspired by this, the work on Turing pattern diagrams on complex networks has been widely carried out, and these works can be subdivided into others from the perspective of theory and application background. In theory, it mainly focuses on the development of different types of network pattern formation theory and pattern dynamics [15,16,17,18,19,20,21,22,23,24,25]. In application, it mainly focuses on different research backgrounds, such as infectious disease models, population models, and rumor propagation models [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49]. Here, we choose the population model as the research background. In particular, we will describe in detail the work related to Turing pattern based on predator-prey models in complex networks. In 2012, Fernandes and collaborators generalized the work of Nakao and Mikhailov to food web models, and their results highlighted that differences in abundance distribution among patches could be dynamically generated through a self-organized Turing pattern [26]. Subsequently, in 2019 and 2021, Chang et al. considered the small delay and large delay factors in the Leslie-Gower model, respectively. The aim was to explore the evolution of rich spatiotemporal patterns caused by delay and further confirmed the impact of delay on Turing pattern [27,32]. Meanwhile, Liu et al. tested the Leslie-Gower model in LA4 ($2\mathrm{D}$ lattices with average degrees $\langle k\rangle =4$) network, LA12 network ($2\mathrm{D}$ lattices with average degrees $\langle k\rangle =12$), LA23 network ($2\mathrm{D}$ lattices with average degrees $\langle k\rangle =23$), ER random network [50] and BA (Barabási-Albert) scale-free network [51] in 2020, respectively. And found that the average degree of the network plays an important role in the pattern formation [31]. Influenced by the above works, Ye and Zhou constructed a predator-prey model with Allee effect and hyperbolic mortality on ER random network, and discovered the existence of the spatiotemporal pattern affected by the initial value [52]. Their model expression is as follows

$$\left\{\begin{array}{c}\frac{\mathrm{d}{u}_i}{\mathrm{d}t}=\overline{\alpha }{u}_i\left(1-u_i\right)\left(\frac{u_i}{u_i+A}\right)-\frac{\overline{\alpha }{u}_i{v}_i}{1+\beta u_i}+d_1\sum _{j=1}^N{L}_{ij}{u}_j,\hfill \\ \frac{\mathrm{d}{v}_i}{\mathrm{d}t}=\frac{\beta u_i{v}_i}{1+\beta u_i}-\frac{\gamma v_i^2}{e+\eta v_i}+d_2\sum _{j=1}^N{L}_{ij}{v}_j,\hfill \\ u_i\left(0\right)\ge 0,v_i\left(0\right)\ge 0.\hfill \end{array}\right.$$

(1)

The biological significance of each parameter and variable in model (1) is shown in Table 1 and [3,4,7,52].

Table 1. Notation used in model (1).

Parameter	Description
$u_i$	The prey population occupies every node index i in the network
$v_i$	The predator population occupies every node index i in the network
$d_1$	The self-diffusion rate of prey population
$d_2$	The self-diffusion rate of predator population
$\overline{\alpha }$	The ratio of intrinsic growth rate to predator birth rate
A	The weak Allee effect constant
$\beta$	The ratio of environmental capacity to half-saturation prey density
e	The coefficients of light attenuation by water
$\eta$	The coefficients of light attenuation by self-shading
$\gamma$	The death rate of the predator population
$L_{ij}$	The elements in the network Laplacian matrix $\mathbf{L}$

Table 1. Notation used in model (1).

Parameter	Description
$u_i$	The prey population occupies every node index i in the network
$v_i$	The predator population occupies every node index i in the network
$d_1$	The self-diffusion rate of prey population
$d_2$	The self-diffusion rate of predator population
$\overline{\alpha }$	The ratio of intrinsic growth rate to predator birth rate
A	The weak Allee effect constant
$\beta$	The ratio of environmental capacity to half-saturation prey density
e	The coefficients of light attenuation by water
$\eta$	The coefficients of light attenuation by self-shading
$\gamma$	The death rate of the predator population
$L_{ij}$	The elements in the network Laplacian matrix $\mathbf{L}$

The above work is mostly based on single-layer networks. We note that generating Turing pattern on multiplex networks does not seem to require significant differences in diffusion rates between two populations, which was first mentioned in [17]. Therefore, compared to single-layer networks, multiplex networks are more general. Soon after, in 2020, Gao et al. attempted to consider cross-diffusion mechanisms on multiplex networks and found many complex heterogeneous spatial patterns that had not yet been discovered in single-layer networks [28]. They applied it to the predator-prey model in 2023 and also discovered these complex heterogeneous spatial patterns, which effectively confirmed their early findings [48]. To our knowledge, most of the research on patterns in multiplex networks currently focuses on Turing pattern, while other types of spatiotemporal patterns are rarely mentioned. Therefore, this article takes the predator-prey model with Allee effect and hyperbolic mortality as an example to explore the possibility of spatial patterns on multiplex networks. Therefore, model (1) is rewritten as follows

$$\left\{\begin{array}{c}\frac{\mathrm{d}{u}_i}{\mathrm{d}t}=f\left(u_i,v_i\right)+d_1\sum _{j=1}^N{L}_{ij}^{\left(u\right)}{u}_j,\hfill \\ \frac{\mathrm{d}{v}_i}{\mathrm{d}t}=g\left(u_i,v_i\right)+d_2\sum _{j=1}^N{L}_{ij}^{\left(v\right)}{v}_j,\hfill \\ u_i\left(0\right)\ge 0,v_i\left(0\right)\ge 0.\hfill \end{array}\right.$$

(2)

The reaction term here can be expressed as

$$f(u_i,v_i)=\overline{\alpha }{u}_i(1-u_i)\left(\frac{u_i}{u_i+A}\right)-\frac{\overline{\alpha }{u}_i{v}_i}{1+\beta u_i},g(u_i,v_i)=\frac{\beta u_i{v}_i}{1+\beta u_i}-\frac{\gamma v_i^2}{e+\eta v_i}.$$

We divide the population on the multiplex ER random networks into two categories: prey population $u_i$ and predator population $v_i$. The schematic diagram is shown in Figure 1. The prey population $u_i$ and predator population $v_i$ occupy every node in layers $U^{\left(u\right)}$ and $U^{\left(v\right)}$, respectively. There is an interactive relationship between two populations across the layers through the imaginary lines which means the reaction term $\left(f\left(u_i,v_i\right)\right.$ and $\left.g\left(u_i,v_i\right)\right)$. Meanwhile, species on the network will also flow in or out of this node on their own layer through the solid lines which means the diffusion term $\left(d_1{\sum }_{j=1}^N{L}_{ij}^{\left(u\right)}{u}_j\right.$ and $\left.d_2{\sum }_{j=1}^N{L}_{ij}^{\left(v\right)}{v}_j\right)$.

Figure 1. Reaction-diffusion prey-predator model locates in a multiplex ER random network, which shows the prey population $u_i$ (green nodes) and predator population $v_i$ (orange nodes) occupy every node on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer, respectively. Among them, the imaginary lines represent the reaction term, and the solid lines represent the diffusion term.

Figure 1. Reaction-diffusion prey-predator model locates in a multiplex ER random network, which shows the prey population $u_i$ (green nodes) and predator population $v_i$ (orange nodes) occupy every node on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer, respectively. Among them, the imaginary lines represent the reaction term, and the solid lines represent the diffusion term.

Here, the total number of nodes in the network is N, and its topology is defined as a symmetric adjacency matrix whose elements satisfy $A_{ij}=1$ if the node i and j are connected where $i,j=1,\dots ,N$ and $i\ne j$, otherwise $A_{ij}=0$. We set the degree of node i as shown by $k_i={\sum }_{j=1}^N{A}_{ij}$. In addition, we sort network nodes index i in decreasing order of their degrees $k_i$. The aim is to ensure that $k_1\ge k_2\ge \cdots k_N$ is established [2,17]. The sum of the incoming fluxes from other connecting nodes index j to node index i provides the diffusion of species u and v at the node index i. Fick’s law states that the flux is inversely proportional to the difference in node concentrations. So, the network Laplacian matrix $L_{ij}=A_{ij}-k_i{\delta }_{ij}$ is considered and ${\delta }_{ij}$ illustrates the Kronecker product, the diffusive flux of the predator population v to the node index i is written as ${\sum }_{j=1}^N{L}_{ij}^{\left(v\right)}{v}_j$ and the diffusive flux of the prey population u to the node index i as ${\sum }_{j=1}^N{L}_{ij}^{\left(u\right)}{u}_j$. The biological significance of other parameters is shown in Table 1.

The main structure of this paper can be organized as follows. In Section 2, by utilizing the Turing instability theory on multiplex networks, the threshold conditions for Turing bifurcation were calculated. In Section 3, numerical experiments are carried out on Turing pattern in single-layer networks and multiplex networks, respectively. In Section 4, numerical experiments are carried out on spatiotemporal patterns in single-layer networks and multiplex networks, respectively. In Section 5, a brief conclusion is given.

2. Turing instability on multiplex networks

In this section, we mainly focus on the critical conditions for the occurrence of Turing bifurcation in model (2). Just as the prerequisite for Turing instability in continuous media is to ensure that the positive equilibrium is stable without diffusion. Therefore, on the networks, we need to first determine the stability of the positive equilibrium of the non-diffusion model (2). For the convenience of calculation, we set $\eta =\gamma$ and $e=1$ as in [4,7,52]. According to [52], model (2) may have four equilibria $E_0=(0,0)$, $E_1=(1,0)$, $E_{*}^{\left(1\right)}=(u_{*}^{\left(1\right)},v_{*}^{\left(1\right)})$, and $E_{*}^{\left(2\right)}=(u_{*}^{\left(2\right)},v_{*}^{\left(2\right)})$ when $\beta >\frac{\gamma }{1-\gamma },0<\gamma <1$, and $\frac{\gamma }{\beta }<A<\frac{\left({(\beta \gamma -\gamma -\beta )}^2+4\beta {\gamma }^2\right)}{4{\beta }^2\gamma }$. It should be emphasized that $u_{*}^{\left(1\right)}=\frac{\beta \gamma -\gamma -\beta -\sqrt{{(\beta \gamma -\gamma -\beta )}^2+4(\gamma -\beta A)\beta \gamma }}{2\beta \gamma }$, $u_{*}^{\left(2\right)}=\frac{\beta \gamma -\gamma -\beta +\sqrt{{(\beta \gamma -\gamma -\beta )}^2+4(\gamma -\beta A)\beta \gamma }}{2\beta \gamma }$, and $v_{*}^{\left(i\right)}=\frac{\beta }{\gamma }{u}_{*}^{\left(i\right)}$ with $i=1,2$. Furthermore, we can obtain the following lemma.

Lemma 1.

[52] When $\beta >\frac{\gamma }{1-\gamma },0<\gamma <1$, and $\frac{\gamma }{\beta }<A<\frac{\left({(\beta \gamma -\gamma -\beta )}^2+4\beta {\gamma }^2\right)}{4{\beta }^2\gamma }$, the stability conditions for the two positive equilibra in model (2) can be summarized as follows

• The positive equilibrium $E_{*}^{\left(1\right)}$ is always unstable which is a saddle.
• Let $H=\frac{v_{*}^{\left(2\right)}}{1+\beta u_{*}^{\left(2\right)}}-\frac{u_{*}^{\left(2\right)}{v}_{*}^{\left(2\right)}}{\left(1+\beta u_{*}^2\right)\left(1-u_{*}^{\left(2\right)}\right)}+\frac{A{v}_{*}^{\left(2\right)}}{\left(1+\beta u_{*}^{\left(2\right)}\right)\left(A+u_{*}^{\left(2\right)}\right)}-\frac{v_{*}^{\left(2\right)}}{{\left(1+\beta u_{*}^{\left(2\right)}\right)}^2}>0$, which can be obtained ${\overline{\alpha }}_H=\frac{\beta u_{*}^{\left(2\right)}}{H{\left(1+\beta u_{*}^{\left(2\right)}\right)}^2}$. So we have the following conclusion: if $\overline{\alpha }<{\overline{\alpha }}_H$, then $E_{*}^{\left(2\right)}$ is stable. If $\overline{\alpha }>{\overline{\alpha }}_H$, then $E_{*}^{\left(2\right)}$ is unstable. Particularly, when $\overline{\alpha }={\overline{\alpha }}_H$ Hopf bifurcation occurs since $\frac{\mathrm{d}Tr}{\mathrm{d}\overline{\alpha }}=H>0$ and

  $$\begin{array}{cc}\hfill Tr=& {\overline{\alpha }}_H\left(\frac{v_{*}^{\left(2\right)}}{1+\beta u_{*}^{\left(2\right)}}-\frac{u_{*}^{\left(2\right)}{v}_{*}^{\left(2\right)}}{\left(1+\beta u_{*}^{\left(2\right)}\right)\left(1-u_{*}^{\left(2\right)}\right)}+\frac{A{v}_{*}^{\left(2\right)}}{\left(1+\beta u_{*}^{\left(2\right)}\right)\left(A+u_{*}^{\left(2\right)}\right)}-\frac{v_{*}^{\left(2\right)}}{{\left(1+\beta u_{*}^{\left(2\right)}\right)}^2}\right)\hfill \\ \hfill & -\frac{\beta u_{*}^{\left(2\right)}}{{\left(1+\beta u_{*}^{\left(2\right)}\right)}^2}.\hfill \end{array}$$

When applied to classical continuous media, the non-uniform perturbation is typically decomposed as spatial Fourier modes that correspond to various wave numbers for plane waves. Inspired by this, Othmer and Scriven observed that the eigenvectors can indicate the influence of the plane wave and wave number on the network, which is ${\Phi }^{\left(\alpha \right)}=\left({\varphi }_1^{\left(\alpha \right)},\dots ,{\varphi }_N^{\left(\alpha \right)}\right)$ and eigenvalues ${\Lambda }_{\alpha }^{(u,v)}(\alpha =1,\dots ,N)$ of Laplace matrix, where ${\sum }_{j=1}^N{L}_{ij}^{(u,v)}{\varphi }_j^{\left(\alpha \right)}={\Lambda }_{\alpha }^{(u,v)}{\varphi }_i^{\left(\alpha \right)}$ [53]. Refer to [2], the following condition $0={\Lambda }_1^{(u,v)}\ge {\Lambda }_2^{(u,v)}\ge \cdots \ge {\Lambda }_N^{(u,v)}$ need to be satisfied. The eigenvectors are orthonormalized as ${\sum }_{i=1}^N{\varphi }_i^{\left(\alpha \right)}{\varphi }_i^{\left(\omega \right)}={\delta }_{\alpha ,\omega }$ and $\alpha ,\omega =1,\dots ,N$. The small perturbations $\delta u_i$ and $\delta v_i$ are introduced to the positive equilibrium like $\left(u_i,v_i\right)=(u_{*}^{\left(2\right)},v_{*}^{\left(2\right)})+\left(\delta u_i,\delta v_i\right)$ and the linearizing form of model (2) can be written as

$$\begin{array}{c}\frac{\mathrm{d}\delta u_i}{\mathrm{d}t}=a_{11}\delta u_i+a_{12}\delta v_i+d_1\sum _{j=1}^N{L}_{ij}^{\left(u\right)}\delta u_j,\\ \frac{\mathrm{d}\delta v_i}{\mathrm{d}t}=a_{21}\delta u_i+a_{22}\delta v_i+d_2\sum _{j=1}^N{L}_{ij}^{\left(v\right)}\delta v_j,\end{array}$$

(3)

where

$$a_{11}=\frac{\overline{\alpha }{v}_{*}^{\left(2\right)}}{1+\beta u_{*}^{\left(2\right)}}-\frac{\overline{\alpha }{u}_{*}^{\left(2\right)}{v}_{*}^{\left(2\right)}}{\left(1+\beta u_{*}^{\left(2\right)}\right)\left(1-u_{*}^{\left(2\right)}\right)}+\frac{\overline{\alpha }A{v}_{*}^{\left(2\right)}}{\left(1+\beta u_{*}^{\left(2\right)}\right)\left(A+u_{*}^{\left(2\right)}\right)}-\frac{\overline{\alpha }{v}_{*}^{\left(2\right)}}{{\left(1+\beta u_{*}^{\left(2\right)}\right)}^2},$$

$$a_{12}=\frac{-\overline{\alpha }{u}_{*}^{\left(2\right)}}{1+\beta u_{*}^{\left(2\right)}},a_{21}=\frac{{\beta }^2{v}_{*}^{\left(2\right)}}{{\left(1+\beta u_{*}^{\left(2\right)}\right)}^2},a_{22}=\frac{-\beta u_{*}^{\left(2\right)}}{{\left(1+\beta u_{*}^{\left(2\right)}\right)}^2}.$$

In particular, for a single-layer network, where ${\Lambda }_{\alpha }^{\left(u\right)}={\Lambda }_{\alpha }^{\left(v\right)}={\Lambda }_{\alpha }$. A detailed analysis can be seen in [52], which can be summarized as the following lemma. In addition, to have an intuitive feeling, we give the dispersion relation between Laplace eigenvalue (${\Lambda }_{\alpha }$) and model (2) eigenvalue (${\lambda }_{\alpha }$) in Figure 3 and Figure 5.

Lemma 2.

[52] There are two critical values $a_{11}+a_{22}=0$ and $a_{11}{a}_{22}-$$a_{12}{b}_{21}+a_{11}{d}_2{\Lambda }_{\alpha }+a_{22}{d}_1{\Lambda }_{\alpha }+d_1{d}_2{\Lambda }_{\alpha }^2=0$. They correspond to the critical value conditions for generating Hopf bifurcation and Turing bifurcation, respectively. And we can determine ${\Lambda }_{\alpha }^{\left(1\right)}$ and ${\Lambda }_{\alpha }^{\left(2\right)}$ as

$$\begin{array}{cc}\hfill {\Lambda }_{\alpha }^{\left(1\right)}& =\frac{-\left(a_{11}{d}_2+a_{22}{d}_1\right)-\sqrt{\Delta }}{2d_1{d}_2},\hfill \\ \hfill {\Lambda }_{\alpha }^{\left(2\right)}& =\frac{-\left(a_{11}{d}_2+a_{22}{d}_1\right)+\sqrt{\Delta }}{2d_1{d}_2}.\hfill \end{array}$$

Where $\Delta ={\left(a_{11}{d}_2+a_{22}{d}_1\right)}^2-4d_1{d}_2\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)$, Turing instability requires ${\Lambda }_{\alpha }^{\left(1\right)}<{\Lambda }_{\alpha }<{\Lambda }_{\alpha }^{\left(2\right)}$ to be maintained.

Next, we consider the more general case of multiplex networks, whose Turing-instability conditions are derived from [17]. We rewrite the linearized model (3) as follows

$$\frac{\mathrm{d}\mathrm{S}}{\mathrm{d}t}=(\mathrm{R}+\mathrm{D})S,$$

(4)

where $\mathrm{S}={\left(\delta u_1,\cdots ,\delta u_N,\delta v_1,\cdots ,\delta v_N\right)}^T$ is considered as the perturbation vector and

$$\mathrm{R}=\left(\begin{array}{cc}{a}_{11}{\mathrm{I}}_N\hfill & a_{12}{\mathrm{I}}_N\hfill \\ a_{21}{\mathrm{I}}_N\hfill & a_{22}{\mathrm{I}}_N\hfill \end{array}\right),\mathrm{D}=\left(\begin{array}{cc}{d}_1{\mathbf{L}}^{\left(u\right)}& 0\\ 0& d_2{\mathbf{L}}^{\left(v\right)}\end{array}\right).$$

Notice that ${\mathrm{I}}_N=diag(1,1,\dots ,1)$ be used as $N\times N$ identity matrix, ${\mathbf{L}}^{\left(u\right)}$ and ${\mathbf{L}}^{\left(v\right)}$ are the Laplacian matrix of $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer, respectively. Among them, $\mathrm{D}=\mathcal{A}-\mathcal{D}$ is established, here

$$\mathcal{A}=\left(\begin{array}{cc}{d}_1{\mathcal{A}}^{\left(u\right)}& 0\\ 0& d_2{\mathcal{A}}^{\left(v\right)}\end{array}\right),\mathcal{D}=\left(\begin{array}{cc}{d}_1{\mathcal{D}}^{\left(u\right)}& 0\\ 0& d_2{\mathcal{D}}^{\left(v\right)}\end{array}\right).$$

Specifically, ${\mathcal{A}}^{\left(u\right)}$ and ${\mathcal{A}}^{\left(v\right)}$ are the adjacency matrices of $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer, respectively. The ${\mathcal{D}}^{\left(u\right)}$ and ${\mathcal{D}}^{\left(v\right)}$ represent the degree matrices of $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer, respectively. For the matrices $\mathcal{A}$, which has elements with values of order $\mathcal{O}\left(d_1\right)$ or $\mathcal{O}\left(d_2\right)$. For the matrices $\mathcal{D}$, which has elements with values of order $\mathcal{O}\left(d_1{k}^{\left(u\right)}\right)$ or $\mathcal{O}\left(d_2{k}^{\left(v\right)}\right)$. Assuming $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are dense enough that $〈k^{\left(u\right)}〉\gg 1$ and $〈k^{\left(v\right)}〉\gg 1$. According to [17], it is easy to see that $\mathcal{A}$ can be neglected and $\mathrm{D}\approx -\mathcal{D}$. Therefore, Eq. (4) can be rewritten as

$$\frac{\mathrm{d}\mathrm{S}}{\mathrm{d}t}=\left(\begin{array}{cc}{a}_{11}{\mathrm{I}}_N-d_1{\mathcal{D}}^{\left(u\right)}& a_{12}{\mathrm{I}}_N\\ a_{21}{\mathrm{I}}_N& a_{22}{\mathrm{I}}_N-d_2{\mathcal{D}}^{\left(u\right)}\end{array}\right)\mathrm{S}.$$

(5)

It is not difficult to obtain

$$Det\left(\begin{array}{cc}{a}_{11}-d_1{k}^{\left(u\right)}-\lambda & a_{12}\\ a_{21}& a_{22}-d_2{k}^{\left(u\right)}-\lambda \end{array}\right)=0,$$

(6)

where $k^{\left(u\right)}$ and $k^{\left(v\right)}$ represent the degrees of nodes in $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer, respectively. In summary, the characteristic equation can be expressed as

$${\lambda }^2-Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)\lambda +Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)=0,$$

(7)

where

$$\begin{array}{cc}& Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)=a_{11}+a_{22}-(d_1{k}^{\left(u\right)}+d_2{k}^{\left(v\right)}),\hfill \\ & Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)=a_{11}{a}_{22}-a_{12}{a}_{21}-(a_{11}{d}_2{k}^{\left(v\right)}+a_{22}{d}_1{k}^{\left(u\right)})+d_1{d}_2{k}^{\left(u\right)}{k}^{\left(v\right)}.\hfill \end{array}$$

Therefore, we have the following discussion, if $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)<0$ and $Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)<0$, Turing instability will occur. When $E_{*}^{\left(2\right)}=\left(u_{*}^{\left(2\right)},v_{*}^{\left(2\right)}\right)$ is stable in non-diffusion model (2). Obviously, the condition of $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)<0$ holds because of $k^{\left(u\right)}\ge 0$ and $k^{\left(v\right)}\ge 0$. Consequently, Turing instability need $Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)<0$ is established, that is to say

$$k^{\left(u\right)}>\frac{\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)-a_{11}{d}_2{k}^{\left(u\right)}}{a_{22}{d}_1-d_1{d}_2{k}^{\left(v\right)}}$$

(8)

for $a_{22}{d}_1-d_1{d}_2{k}^{\left(v\right)}>0$ or

$$k^{\left(u\right)}<\frac{\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)-a_{11}{d}_2{k}^{\left(u\right)}}{a_{22}{d}_1-d_1{d}_2{k}^{\left(v\right)}}$$

(9)

for $a_{22}{d}_1-d_1{d}_2{k}^{\left(v\right)}<0$. Similar to the idea of a single-layer network, we give the relationship between $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ and $Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ with the degrees of nodes in $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer $(k^{\left(u\right)},k^{\left(v\right)})$, respectively (see Figure 2). We can find that when the control parameter $\overline{\alpha }=1.2$, there is the possibility of $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)>0$. That is to say, when the diffusion item is considered, model (2) becomes unstable. Note that it is a non-Turing pattern at this time.

Figure 2. When $\overline{\alpha }=0.65$, panels (a) and (b) represent the relationship between $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ and $Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ with the degrees of nodes in $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer $(k^{\left(u\right)},k^{\left(v\right)})$, respectively. When $\overline{\alpha }=1.2$, panels (a) and (b) represent the relationship between $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ and $Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ with the degrees of nodes in $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer $(k^{\left(u\right)},k^{\left(v\right)})$, respectively. The rest of the parameters are $d_1=0.05,d_2=0.05,A=0.01,\beta =6,\gamma =0.5,e=1$, and $\eta =0.5$.

Figure 2. When $\overline{\alpha }=0.65$, panels (a) and (b) represent the relationship between $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ and $Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ with the degrees of nodes in $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer $(k^{\left(u\right)},k^{\left(v\right)})$, respectively. When $\overline{\alpha }=1.2$, panels (a) and (b) represent the relationship between $Tr\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ and $Det\left(k^{\left(u\right)},k^{\left(v\right)}\right)$ with the degrees of nodes in $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer $(k^{\left(u\right)},k^{\left(v\right)})$, respectively. The rest of the parameters are $d_1=0.05,d_2=0.05,A=0.01,\beta =6,\gamma =0.5,e=1$, and $\eta =0.5$.

3. Turing pattern on networks

In this section, we discuss the evolution of Turing pattern on single-layer ER random networks and multiplex ER random networks. In particular, we take the predator diffusion rate $d_2$, the average degree of the network ${\Lambda }_{\alpha }$, and the ratio of average degrees on different networks as the research objects, respectively. The purpose is to study the factors affecting the Turing pattern on the ER random network and discuss how these factors affect the generation of the pattern. Before starting the numerical experiment, we need to emphasize some basic settings, where the total number of nodes in the network is $N=1600$, the initial value is $u_i\left(0\right)=u_{*}^{\left(2\right)}+0.1\times rand(0,1),v_i\left(0\right)=v_{*}^{\left(2\right)}+0.1\times rand(0,1)$, and the remaining parameters are $d_1=0.05,A=0.01,\overline{\alpha }=0.65,\beta =6,\gamma =0.5,e=1,\eta =0.5$, where $rand(0,1)$ is a random perturbation. $u_{*}^{\left(2\right)}=u_{*}$ and $v_{*}^{\left(2\right)}=v_{*}$ always holds in all figures.

3.1. Diffusion rate induces Turing pattern on single-layer network

We first discuss the effect of predator diffusion rate $d_2$ on Turing pattern. According to Lemma 2, we give the dispersion relation between Laplace eigenvalue (${\Lambda }_{\alpha }$) and model (2) eigenvalue (${\lambda }_{\alpha }$) as shown in Figure 3. It can be seen that when the two populations have the same diffusion rate ($d_1=d_2=0.05$), there is no ${\Lambda }_{\alpha }$ so that $Re\left({\lambda }_{\alpha }\right)>0$. However, when we increase the predator diffusion rate ($d_2=1$), there are critical thresholds ${\Lambda }_{\alpha }^{\left(1\right)}\approx -3.01$ and ${\Lambda }_{\alpha }^{\left(2\right)}\approx -0.546$ such that $Re\left({\lambda }_{\alpha }\right)=0$. Therefore, it is not difficult to see that when ${\Lambda }_{\alpha }\in ({\Lambda }_{\alpha }^{\left(1\right)},{\Lambda }_{\alpha }^{\left(2\right)})$, Turing instability occurs where $Re\left({\lambda }_{\alpha }\right)>0$. In addition, we also provide the density (red dots) of predator population ($v_i$) with different predator diffusion rates, the predator population density ($v_i$) on the ER random network varies with time under different predator diffusion rates, and the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different predator diffusion rates over time, respectively (see Figure 4). The results indicate that on single-layer network, the diffusion rate plays a decisive role in Turing pattern.

Figure 3. Dispersion relation between Laplace eigenvalue (${\Lambda }_{\alpha }$) and model (2) eigenvalue (${\lambda }_{\alpha }$) with different diffusion rate. The predator diffusion rate is $d_2=0.05$ for (a). The predator diffusion rate is $d_2=1$ for (b). Where the average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =4$. The rest of the parameters are $d_1=0.05,A=0.01,\overline{\alpha }=0.65,\beta =6,\gamma =0.5,e=1$, and $\eta =0.5$.

Figure 3. Dispersion relation between Laplace eigenvalue (${\Lambda }_{\alpha }$) and model (2) eigenvalue (${\lambda }_{\alpha }$) with different diffusion rate. The predator diffusion rate is $d_2=0.05$ for (a). The predator diffusion rate is $d_2=1$ for (b). Where the average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =4$. The rest of the parameters are $d_1=0.05,A=0.01,\overline{\alpha }=0.65,\beta =6,\gamma =0.5,e=1$, and $\eta =0.5$.

Figure 4. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different predator diffusion rates. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different predator diffusion rates. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different predator diffusion rates over time. The predator diffusion rate is $d_2=0.05$ for (a), (b), and (c). The predator diffusion rate is $d_2=1$ for (d), (e), and (f). Where the average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =4$.

Figure 4. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different predator diffusion rates. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different predator diffusion rates. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different predator diffusion rates over time. The predator diffusion rate is $d_2=0.05$ for (a), (b), and (c). The predator diffusion rate is $d_2=1$ for (d), (e), and (f). Where the average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =4$.

3.2. Network average degree induces Turing pattern on single-layer network

The above experiments have well verified the importance of diffusion rate. In this group of experiments, we investigate how the spatial distribution of species changes when the diffusion rate meets the Turing instability condition and changes the average degree of the network ($\langle k\rangle =\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle$). As we gradually increase the value of the average degree of the network $\langle k\rangle$, we find that although the Turing instability region always exists, the eigenvalues ${\Lambda }_{\alpha }$ (red dots) of the Laplacian matrix do not fall within this region. Therefore, when $\langle k\rangle =\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =15$, Turing pattern is not captured in Figure 5. Similarly, in this set of experiments, we also provide the density (red dots) of predator population ($v_i$) with different network average degrees, predator population density ($v_i$) on the ER random network varies with time under different network average degrees, and curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different average degrees over time (see Figure 6). The results indicate that as the average degree of the network increases, the spatial distribution of species gradually concentrates at equilibrium $E_{*}^{\left(2\right)}=(u_{*}^{\left(2\right)},v_{*}^{\left(2\right)})$ until the Turing pattern disappears.

Figure 5. Dispersion relation between Laplace eigenvalue (${\Lambda }_{\alpha }$) and model (2) eigenvalue (${\lambda }_{\alpha }$) with different average degree. The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =10$ for (a). The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =15$ for (b). The rest of the parameters are $d_1=0.05,d_2=0.05,A=0.01,\overline{\alpha }=0.65,\beta =6,\gamma =0.5,e=1$, and $\eta =0.5$.

Figure 5. Dispersion relation between Laplace eigenvalue (${\Lambda }_{\alpha }$) and model (2) eigenvalue (${\lambda }_{\alpha }$) with different average degree. The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =10$ for (a). The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =15$ for (b). The rest of the parameters are $d_1=0.05,d_2=0.05,A=0.01,\overline{\alpha }=0.65,\beta =6,\gamma =0.5,e=1$, and $\eta =0.5$.

Figure 6. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different network average degrees. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different network average degrees. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different average degrees over time. The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =10$ for (a), (b), and (c). The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =15$ for (d), (e), and (f).

Figure 6. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different network average degrees. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different network average degrees. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different average degrees over time. The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =10$ for (a), (b), and (c). The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =15$ for (d), (e), and (f).

3.3. Ratio of network average degree induces Turing pattern on multiplex networks

It is well known that diffusivity plays a decisive role in the formation of Turing pattern on continuous media, but is it the same on discrete media or complex networks? Kouvaris et al. provide an answer to this in their work, stating that even if two populations have the same diffusion rate on multiplex networks, Turing pattern will be generated [17]. Therefore, we hope to verify whether this conclusion is equally applicable to our model. Before starting the experiment, we set the diffusion rate of the two populations to $d_1=d_2=0.05$ and determine the average degree $\langle k^{\left(u\right)}\rangle =4$ on $U^{\left(u\right)}$ layer. Observing the changes in the spatial distribution pattern of species by changing the average degree $\langle k^{\left(v\right)}\rangle$ on $U^{\left(v\right)}$ layer. The threshold conditions for generating Turing bifurcation can be obtained from Eq. (8) and (9), as shown in Figure 7. It is not difficult to find that the Turing instability condition does not hold when the average degree of the two layers network is the same ($\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =4$). So, we increase the average degree on $U^{\left(v\right)}$ layer, which led to Turing instability. By comparing Figure 4((d), (e), and (f)) and Figure 8, it can be intuitively seen that the Turing pattern has evolved from scratch. In addition, we also find that the spatial distribution of species changes with the increase of average degree differentiation in different layers.

Figure 7. (a), (b), and (c) represent degree combination for pairs of nodes in plane $k^{\left(v\right)}-k^{\left(u\right)}$ together with the curve of Turing bifurcation (black) on multiplex ER random networks. The average degree on $U^{\left(v\right)}$ layer is $\langle k^{\left(v\right)}\rangle =4$, $\langle k^{\left(v\right)}\rangle =20$, and $\langle k^{\left(v\right)}\rangle =80$ for (a), (b), and (c). Where have same average degree on $U^{\left(u\right)}$ layer $\langle k^{\left(u\right)}\rangle =4$.

Figure 7. (a), (b), and (c) represent degree combination for pairs of nodes in plane $k^{\left(v\right)}-k^{\left(u\right)}$ together with the curve of Turing bifurcation (black) on multiplex ER random networks. The average degree on $U^{\left(v\right)}$ layer is $\langle k^{\left(v\right)}\rangle =4$, $\langle k^{\left(v\right)}\rangle =20$, and $\langle k^{\left(v\right)}\rangle =80$ for (a), (b), and (c). Where have same average degree on $U^{\left(u\right)}$ layer $\langle k^{\left(u\right)}\rangle =4$.

Figure 8. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different ratios of network average degree. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different ratios of network average degree. (c) and (f) represent the curves of the maximum, minimum, and average values for the predator population density ($v_i$) in all nodes on the network over time. The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =20$ for (a), (b), and (c). The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =80$ for (d), (e), and (f).

Figure 8. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different ratios of network average degree. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different ratios of network average degree. (c) and (f) represent the curves of the maximum, minimum, and average values for the predator population density ($v_i$) in all nodes on the network over time. The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =20$ for (a), (b), and (c). The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =80$ for (d), (e), and (f).

4. Spatiotemporal pattern on networks

When we theoretically analyze the Turing instability on the network, we find that the premise of Turing pattern always requires that the positive equilibrium in the non-diffusion model is always stable. If the positive equilibrium in the non-diffusion model is not stable at this time, what kind of spatial distribution state will the model present? On continuous media, Nagano and Maeda discovered the existence of spiral and lattice patterns [3]. What will happen on ER random network? It should be pointed out that this spatial pattern type is not Turing pattern but another spatiotemporal pattern. So, we conduct the following numerical experiments. The value of the parameter is the same as Section 3, except for $\alpha$.

4.1. Parameter induces spatiotemporal pattern on single-layer network

In the beginning, we hope to conduct our testing on single-layer network and by changing the value of control parameter $\overline{\alpha }$. It turned out that even if the diffusion rates of two populations are the same, we could still discover the existence of spatial heterogeneity patterns in Figure 9. And as $\overline{\alpha }$ increases, this phenomenon becomes more apparent. However, when looking at Figure 4 (a), (b), and (c), no such spatiotemporal patterns are observed. So we speculate that this is a spatiotemporal pattern phenomenon caused by parameters. By observing the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different parameters $\overline{\alpha }$ over time. We find the existence of an obvious irregular periodic oscillation phenomenon. Similar spatiotemporal pattern phenomena have also received attention in some work [40,41,42].

Figure 9. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different parameters $\alpha$. (b) and (e) represent the density of the predator population ($v_i$) on the ER random network varies with time under different parameters $\overline{\alpha }$. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different parameters $\overline{\alpha }$ over time. The parameter is $\overline{\alpha }=0.85$ for (a), (b), and (c). The parameter is $\overline{\alpha }=1.2$ for (d), (e), and (f). Where the average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =4$ and same diffusion rate is $d_1=d_2=0.05$.

Figure 9. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different parameters $\alpha$. (b) and (e) represent the density of the predator population ($v_i$) on the ER random network varies with time under different parameters $\overline{\alpha }$. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different parameters $\overline{\alpha }$ over time. The parameter is $\overline{\alpha }=0.85$ for (a), (b), and (c). The parameter is $\overline{\alpha }=1.2$ for (d), (e), and (f). Where the average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =4$ and same diffusion rate is $d_1=d_2=0.05$.

4.2. Network average degree induces spatiotemporal pattern on single-layer network and multiplex networks

Then, we consider whether this spatiotemporal pattern is still affected by the average degree of single-layer network and multiplex networks. Therefore, we provide a comparative experiment of spatiotemporal patterns under different average degrees (see Figure 10). When the average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =2$, spatial heterogeneity is particularly evident. If the diffusion rates of two populations on the single-layer network are the same, there will not be Turing pattern generated. We want to know what changes occur if the generation conditions of Turing pattern and spatiotemporal pattern are generated simultaneously. Due to this issue, we conduct our experiment on multiplex networks. To our surprise, we find that as the average value on the controlled $U^{\left(v\right)}$ layer increased (from $\langle k^{\left(v\right)}\rangle =13$ to $\langle k^{\left(v\right)}\rangle =20$), the two spatial patterns switched as shown in Figure 11. In other words, the spatiotemporal pattern transforms into Turing pattern, and even the oscillations in time disappear and become stable. These results indicate that network topology also has a significant impact on the generation of spatiotemporal patterns.

Figure 10. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different network average degrees. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different network average degrees. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different average degrees over time. The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =2$ for (a), (b), and (c). The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =5$ for (d), (e), and (f). Where have the same diffusion rate $d_1=d_2=0.05$ and $\overline{\alpha }=1.2$.

Figure 10. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different network average degrees. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different network average degrees. (c) and (f) represent the curves of the maximum, minimum, and average values of the predator population density ($v_i$) in all nodes on the network at different average degrees over time. The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =2$ for (a), (b), and (c). The average degree is $\langle k^{\left(u\right)}\rangle =\langle k^{\left(v\right)}\rangle =5$ for (d), (e), and (f). Where have the same diffusion rate $d_1=d_2=0.05$ and $\overline{\alpha }=1.2$.

Figure 11. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different ratios of network average degree. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different ratios of network average degree. (c) and (f) represent the curves of the maximum, minimum, and average values for the predator population density ($v_i$) in all nodes on the network over time. The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =13$ for (a), (b), and (c). The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =20$ for (d), (e), and (f). Where have the same diffusion rate $d_1=d_2=0.05$ and $\overline{\alpha }=1.2$.

Figure 11. (a) and (d) represent the density (red dots) of predator population ($v_i$) with different ratios of network average degree. (b) and (e) represent the predator population density ($v_i$) on the ER random network varies with time under different ratios of network average degree. (c) and (f) represent the curves of the maximum, minimum, and average values for the predator population density ($v_i$) in all nodes on the network over time. The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =13$ for (a), (b), and (c). The average degree on $U^{\left(u\right)}$ layer and $U^{\left(v\right)}$ layer are $\langle k^{\left(u\right)}\rangle =4$ and $\langle k^{\left(v\right)}\rangle =20$ for (d), (e), and (f). Where have the same diffusion rate $d_1=d_2=0.05$ and $\overline{\alpha }=1.2$.

5. Conclusion

Taken together, this paper considers the spatial patterns of a predator-prey model with weak Allee effect and hyperbolic mortality on multiplex networks. The main content of our research is to theoretically analyze the critical conditions for generating Turing patterns on multiplex networks. Moreover, we also conduct numerical experiments on several factors that affect spatial patterns, such as diffusion rate, the average degree of network, and the ratio of average degree. Our results indicate that the diffusion rate of two populations is not a necessary condition for generating Turing pattern on multiplex networks, which effectively validates the conclusion of [17,54]. More interestingly, we also found a switch between two types of spatial patterns. Specifically, when model (2) exhibits spatiotemporal patterns under specific parameters, as the difference in the average degree of the different layers on the network increases, the pattern transforms into Turing pattern. Meanwhile, this phenomenon seems to be unique to multiplex networks and has not been found on single-layer networks. These findings once again emphasize the important impact of network topology on species spatial distribution.

In future work, we will continue our work from two aspects: (1) The practical significance of spatiotemporal patterns in biological systems often represents the spatial distribution of species. In recent years, there have been some control methods for Turing patterns in reaction-diffusion systems on networks, which can effectively present the spatial distribution of species in a predetermined pattern [33,34,35]. Therefore, we hope to apply this theory to our model. (2) The higher-order structure of networks has gradually become a new research hotspot in the field of network science. Considering that higher-order interactions affect the topological properties of the network, which can interfere with the formation of spatiotemporal patterns [25,47,54]. Therefore, we also hope to introduce higher-order networks into our model, to understand the potential role of higher-order structures in species spatial distribution.