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

: With the rapid development of network science, Turing patterns on complex networks have attracted extensive attention from researchers. In this paper, we focus on spatial patterns in multiplex ER (Erdös-Rényi) random networks, taking the predator–prey model with Allee effect and hyperbolic mortality as an example. In theory, the threshold condition for generating Turing patterns is given using the Turing instability theory of multiplex networks. Numerically, we design relevant experiments to explore the impact of network topology on Turing patterns. The factors considered include model parameters, diffusion rate, average degree of the network, and differences in the average degree of different layers. The results indicate that the importance of diffusion rate and network average degree for Turing patterns is afﬁrmed on the single-layer network. For multiplex networks, the differentiation of average degrees in different layers controls the generation of Turing patterns, which are not affected by the diffusion rates of the two populations. More interestingly, we observe the switching of Turing patterns and spatiotemporal patterns. We believe that these ﬁndings contribute to a better understanding of self-organization on complex networks.


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]. Turing demonstrated in 1952 that variations in the diffusion constants of activator and inhibitor species might cause the uniform state to destabilize and result in the spontaneous formation of periodic spatial patterns [13]. 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,14]. 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 patterns based on this, so they attempted to consider more complex inducing factors such as hyperbolic mortality [4] or delay [15]. They also displayed many beautiful spatial patterns, including hexagonal spot patterns, stripe patterns, micro-spiral patterns, and target patterns. As an important ecological effect of population dynamics, the Allee effect reflects the correlation between population density and average individual fitness within a population. Specifically, 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 was first studied on large-scale random networks 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 patterns 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 [16][17][18][19][20][21][22][23][24][25][26]. In application, it mainly focuses on different research backgrounds, such as infectious disease models, population models, and rumor propagation models . Here, we choose the population model as the research background. In particular, we will describe in detail the work related to Turing patterns 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 [27]. 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 patterns [28,33]. Meanwhile, Liu et al. tested the Leslie-Gower model in the LA4 (2D lattices with average degrees k = 4) network, LA12 network (2D lattices with average degrees k = 12), LA23 network (2D lattices with average degrees k = 23), ER random network [52] and the BA (Barabási-Albert) scale-free network [53] in 2020, respectively. And found that the average degree of the network plays an important role in the pattern formation [32]. Influenced by the above works, Ye and Zhou constructed a predator-prey model with Allee effect and hyperbolic mortality on the ER random network, and discovered the existence of the spatiotemporal pattern affected by the initial value [54].
The above work is mostly based on single-layer networks. We note that generating Turing patterns on multiplex networks does not seem to require significant differences in diffusion rates between two populations, which was first mentioned in [18]. 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 [29]. They applied it to the predator-prey model in 2023 and also discovered these complex heterogeneous spatial patterns, which effectively confirmed their early findings [49]. To our knowledge, most of the research on patterns in multiplex networks currently focuses on Turing patterns, 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, our model can be written as follows The reaction term here can be expressed as 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 (u) and U (v) , respectively. There is an interactive relationship between two populations across the layers through the imaginary lines, which means the reaction term ( f (u i , v i ) and g(u i , v i )). 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 ij v j . 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 the U (u) layer and U (v) 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, . . . , N and i = j; otherwise, A ij = 0. We set the degree of node i, as shown by k i = ∑ N j=1 A ij . 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 δ ij is considered and δ ij illustrates the Kronecker δ, the diffusive flux of the predator population v to the node index i is written as ij v j and the diffusive flux of the prey population u to the node index i as ∑ N j=1 L (u) ij u j . The biological significance of other parameters is shown in Table 1. For hyperbolic mortality, it should be noted that when η = 0 and e = 0, it degenerates into a quadratic mortality rate. When η = 0 and e = 0, it degenerates into a linear mortality rate. When both η and e are not zero, this is a hyperbolic mortality rate. 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 α The ratio of intrinsic growth rate to predator birth rate A The weak Allee effect constant β The ratio of environmental capacity to half-saturation prey density e The coefficients of light attenuation by water η The coefficients of light attenuation by self-shading γ The death rate of the predator population L ij The elements in the network Laplacian matrix L 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 the 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.

Turing Instability on Multiplex Networks
In this section, we mainly focus on the critical conditions for the occurrence of Turing bifurcation in model (1). 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 (1). For the convenience of calculation, we set η = γ and e = 1 as in [4,7,54]. According to [54], model (1) may have four equilibria E 0 = (0, 0), Furthermore, we can obtain the following lemma.
, the stability conditions for the two positive equilibria in model (1) can be summarized as follows

•
The positive equilibrium E (1) * is always unstable, which is a saddle.
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 [55]. Refer to [2], the following condition 0 = Λ needs to be satisfied. The eigenvectors are orthonormalized as (2) * ) + (δu i , δv i ) and the linearizing form of model (1) can be written as where In particular, for a single-layer network, where Λ a detailed analysis can be seen in [54], which can be summarized as the following lemma. In addition, to have an intuitive feeling, we give the dispersion relation between Laplace eigenvalue (Λ α ) and model (1) eigenvalue (λ α ).

Lemma 2 ([54]
). There are two critical values a 11 + a 22 = 0 and a 11 a 22 − a 12 a 21 + a 11 They correspond to the critical value conditions for generating Hopf bifurcation and Turing bifurcation, respectively. And we can determine Λ Next, we consider the more general case of multiplex networks, whose Turinginstability conditions are derived from [18]. We rewrite the linearized model (2) as follows where S = (δu 1 , · · · , δu N , δv 1 , · · · , δv N ) T is considered as the perturbation vector and R = a 11 I N a 12 I N a 21 I N a 22 Notice that I N = diag(1, 1, . . . , 1) can be used as N × N identity matrix, L (u) and L (v) are the Laplacian matrix of the U (u) layer and U (v) layer, respectively. Among them, Specifically, A (u) and A (v) are the adjacency matrices of the U (u) layer and U (v) layer, respectively. The D (u) and D (v) represent the degree matrices of the U (u) layer and U (v) layer, respectively. For the matrices A, which has elements with values of order O(d 1 ) or O(d 2 ). For the matrices D, which has elements with values of order O d 1 k (u) or . Assuming the U (u) layer and U (v) layer are dense enough that k (u) 1 and k (v) 1. According to [18], it is easy to see that A can be neglected and D ≈ −D. Therefore, Equation (3) can be rewritten as It is not difficult to obtain Det where k (u) and k (v) represent the degrees of nodes in the U (u) layer and U (v) layer, respectively. In summary, the characteristic equation can be expressed as where Therefore, we have the following discussion: if Tr k (u) , k (v) < 0 and Det k (u) , k (v) < 0, Turing instability will occur. When E is stable in the non-diffusion model (1), obviously, the condition of Tr k (u) , k (v) < 0 holds because of k (u) ≥ 0 and k (v) ≥ 0. Consequently, Turing instability needs Det k (u) , k (v) < 0 to be established, that is to say, for for Similar to the idea of a single-layer network, we give the relationship between Tr k (u) , k (v) and Det k (u) , k (v) with the degrees of nodes in the U (u) layer and U (v) layer (k (u) , k (v) ), respectively (see Figure 2). We can find that when the control parameterᾱ = 1.2, there is the possibility of Tr k (u) , k (v) > 0. That is to say, when the diffusion item is considered, model (1) becomes unstable. Note that it is a non-Turing pattern at this time.

Turing Pattern on Networks
In this section, we discuss the evolution of Turing patterns 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 Λ α , 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 (0) = u rand(0, 1), and the remaining parameters are d 1 = 0.05, A = 0.01,ᾱ = 0.65, β = 6, γ = 0.5, e = 1, η = 0.5, where rand(0, 1) is a random perturbation. u

Diffusion Rate Induces Turing Pattern on Single-Layer Network
We first discuss the effect of predator diffusion rate d 2 on the Turing pattern. According to Lemma 2, we give the dispersion relation between the Laplace eigenvalue (Λ α ) and the model (1) eigenvalue (λ α ), 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 Λ α so that Re(λ α ) > 0. However, when we increase the predator diffusion rate (d 2 = 1), there are critical thresholds Λ (1) α ≈ −3.01 and Λ (2) α ≈ −0.546 such that Re(λ α ) = 0. Therefore, it is not difficult to see that when Λ α ∈ (Λ (1) α , Λ α ), Turing instability occurs where Re(λ α ) > 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 varying 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 a single-layer network, the diffusion rate plays a decisive role in the Turing pattern.

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 ( k = k (u) = k (v) ). As we gradually increase the value of the average degree of the network k , we find that, although the Turing instability region always exists, the eigenvalues Λ α (red dots) of the Laplacian matrix do not fall within this region. Therefore, when k = k (u) = k (v) = 15, the 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 varying 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

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 patterns 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, a Turing pattern will be generated [18]. 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 k (u) = 4 on the U (u) layer. We observe the changes in the spatial distribution pattern of species by changing the average degree k (v) on U (v) layer. The threshold conditions for generating Turing bifurcation can be obtained from Equations (7) and (8), 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 ( k (u) = k (v) = 4). So, we increase the average degree on U (v) layer, which leads to Turing instability. By comparing Figures 4d,e,f and 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.

Spatiotemporal Pattern on Networks
When we theoretically analyze the Turing instability on the network, we find that the premise of the 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 the ER random network? It should be pointed out that this spatial pattern type is not a Turing pattern but another spatiotemporal pattern. Therefore, the following experiments were carried out. The value of the parameter is the same as Section 3, except forᾱ.

Parameter Induces Spatiotemporal Pattern on Single-Layer Network
In the beginning, we hope to conduct our testing on a single-layer network and by changing the value of the control parameterᾱ. 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ᾱ increases, this phenomenon becomes more apparent. However, when looking at Figure 4a-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ᾱ over time, we find the existence of an obvious irregular periodic oscillation phenomenon. Similar spatiotemporal pattern phenomena have also received attention in some work [41][42][43].

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 k (u) = k (v) = 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 a Turing pattern generated. We want to know what changes occur if the generation conditions of Turing patterns and spatiotemporal patterns 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 (v) layer increased (from k (v) = 13 to k (v) = 20), the two spatial patterns switched, as shown in Figure 11. In other words, the spatiotemporal pattern transforms into a 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.

Conclusions
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 a Turing pattern on multiplex networks, which effectively validates the conclusion of [18,56]. More interestingly, we capture a switch between two types of spatial patterns. Specifically, when model (1) 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 a Turing pattern. Meanwhile, this phenomenon seems to be unique to multiplex networks and has not been found on single-layer networks. To our knowledge, the switching of this phenomenon has not been found in previous work. These findings once again emphasize the important impact of network topology on species spatial distribution.
In future work, we 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 [34][35][36]. 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 [26,48,56]. 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.