Delay-induced patterns in a predator–prey model on complex networks with diffusion

Reaction-diffusion (RD) systems with time delays have been commonly used in modeling biological systems and can significantly change the dynamics of these systems. For predator–prey model with modified Leslie–Gower and Holling-type III schemes governed by RD equations, instability induced by time delay can generate spiral waves. Considering that populations are usually organized as networks instead of being continuously distributed in space, it is essential to study the predator–prey model on complex networks. In this paper, we investigate instability induced by time delay for the corresponding network organized system and explore pattern formations on several different networks including deterministic networks and random networks. We firstly obtain instability condition via linear stability analysis and then the condition is applied to study pattern formations for the model in question. The simulation results show that wave patterns can be generated on different networks. However, wave patterns on random networks differ significantly from patterns on deterministic networks. Finally, we discuss the influences of network topology on wave patterns from the aspects of amplitude and period, and reveal the ecology significance implied by these results.


Introduction
Reaction-diffusion (RD) system with time delay is used to describe evolution system which is affected by the present as well as the past state. Its characteristic is that the influences of history on the present situation of system are fully taken into account. As a powerful tool for describing spatiotemporal dynamics, the delay RD systems have been widely applied in many fields such as ecology [1], epidemiology [2], chemistry [3], physics [4], evolutionary games [5], and so on. Taken ecology as an example, factors that introduce time delay may include age structure of the population (influencing the birth and death rates), maturation periods, incubation periods, food storage or ingestion delays, resource regeneration times [6].
In continuous space, a RD system often leads to nonlinear RD equations-a class of typical partial differential equations (PDEs)-with reaction terms corresponding to the local dynamics and diffusion terms describing the motion of matter in space. Models governed by RD equations can show rich dynamic behaviors involving traveling waves, self-organized spatial patterns or chaoses. Specifically, the introduction of time delay causes these models to exhibit much more complex dynamical behaviors than the counterpart without time delay. Therefore, RD models with time delays have received wide attention over the past years [7][8][9][10][11]. For example, the authors in [12,13] derived the critical threshold for the spatiotemporal Hopf bifurcation by considering the delay as bifurcation parameter and further analyzed the influence of delay on pattern formations. Ghosh et al [3,14] investigated RD systems with small time delays in which the authors reported that small delays may result in a Turing instability of such systems so that spatial patterns were formed. The inadequacy of the method in [3,14] is that only systems with small delays are valid. Subsequently, Zhang et al [15] proposed a method to deal with cases with larger time delays. With the help of this method, they studied Turing instability induced by large time delays and analyzed the effect of delay on this instability. Moreover, wave patterns of RD equations with time delays were studied in [8,16,17].
Many phenomena in nature can be modeled by continuous RD processes. However, there are some other RD processes occur in discrete spaces (random networks) rather than in continuous spaces. One typical example is that the species are usually organized in patches in a predator-prey system [18], instead of being continuously distributed in space. The patches can be represented as nodes of a network such that predator and prey interact locally in each patch and diffuse through connected nodes. Another example for discrete version of RD process is found in metapopulation epidemic models [19,20]. In this case, network nodes represent different locations such as cities or urban areas, individuals-who are divided into classes denoting their state with respect to the modeled disease such as infected, susceptible, immune and so on-moving between connected locations, and the reaction processes account for the possibility that individuals in the same location may get in contact and change their state according to the infection dynamics.
With the development of network science, the dynamical behaviors of RD systems defined on networks have received considerable attention. Early in 1971 Othmer et al [21] pointed out that Turing instability can occur in network organized systems and proposed a general mathematical framework for the analysis of such instability. Earlier research, however, were limited to regular lattices or small networks [22][23][24], which cannot account for the effect of network scale and topology on pattern formation. Recently, the key research results given by Nakao et al [25] changed this situation. They studied Turing patterns in large random networks, which reveal striking differences from the classical behavior and pave the way for novel discoveries in an area of widespread interest. Subsequently, further exploration has been done based on this seminal work. Fernandes et al [26] investigated Turing patterns in a predator-prey food webs with more than two species on networks. Asllani et al [27] studied the theory of pattern formation on directed networks. Perc et al [28] reviewed the cyclic dominance in evolutionary games focusing on pattern formation in which species contact in different networks. Kouvaris et al [29] extended a predator-prey model on multiplex networks.
Among these references mentioned above only RD systems on networks without time delay were considered. Here we investigate Turing patterns of a Leslie-Gower Holling-type III predator-prey model with time delay on networks. We have investigated the continuous version with delay in [30] and network version without delay. The results showed that these models have complex dynamical properties and instability induced by diffusion can generate stationary patterns. Moreover, for the model in continuous media, spiral waves could be resulted in if the instability is induced by time delay. Now that there exists wave patterns for continuous RD model when delay satisfies certain conditions, could the model on networks have similar phenomenon? If so, what forms of wave patterns will display in? What are the effects of network topology on wave patterns? Through the research we try to address these problems.
The paper is structured in the following way. In section 2, we introduce the Leslie-Gower Holling-type III predator-prey model with time delay on networks. In section 3, we analyze instability of the model induced by time delay with the help of linear stability method. Wave pattrerns and the effects of network topology on them are given in section 4 through simulations. Finally, we make some discussions in section 5.

Model with delay on networks
Holling-type III response exists universally in population dynamics. A predator-prey model incorporates the Leslie-Gower functional response as well as Holling-type III functional response can be written as with initial condition u(0)=u 0 >0, v(0)=v 0 >0. This model describes a prey population u which serves as food for a predator with population v. The parameters α, β, γ are assumed to be only positive values. It should be noted that the first equation of system (2.1) is standard while the second equation is not because that the righthand side in the second equation contains a Leslie-Gower term. The Leslie-Gower formulation is based on the assumption that reduction in a predator population has a reciprocal relationship with per capita availability of its preferred food. Indeed, Leslie introduced a predator-prey model where the carrying capacity of the predator environment is proportional to the number of prey. He stresses the fact that there are upper limits to the rates of increase of both prey u and predator v, which are not recognized in the Lotka-Volterra model (see, e.g. [31]).
Including time delay in model (2.1) is a more realistic approach to the understanding of predator-prey dynamics. We introduce a single discrete delay τ>0 in the negative feedback of the predator's density, then obtain the following system , .
Taking into account the inhomogeneous distribution of the predator and its prey in different spatial locations within a fixed bounded domain  W Ì 2 at any given time, and the natural tendency of each species to diffuse to areas of smaller population concentration, we arrive at a continuous RD system defined on = W´¥ ( Here the positive constants d and σ d denote diffusion coefficients, D = ¶ ¶ + ¶ ¶ x y 2 2 2 2 is the Laplacian operator, n is the outward unit normal vector of the boundary ∂Ω, u 0 , v 0 are continuous positive functions. Time delay plays an important role in dynamical system (2.3), where it has been recognized to contribute critically to the stable or unstable outcome of prey density due to predation (see [30]). Now we extend model (2.3) to the network analogue version where prey and predator occupy discrete nodes of a network and are diffusively transported over links connecting them. The links represent diffusive connections between species from one habitat to another (see, e.g. [25]). In continuous media the Laplacian operator account for the diffusion of populations in space, while this operator is replaced by Laplacian matrix of a network in discrete media. In fact, Laplacian operator is the limit case of Laplacian matrices of some especial networks, we illustrate it in appendix. For a given undirected network with N nodes, the Laplacian matrix holds. Equations describing network organized prey-predator system are thus given by , i=1, 2, L, N represent populations of prey and predator at time t in the habitat i, respectively.
In the present paper we study the instability induced by time delay for model (2.4) and investigate corresponding patterns on some networks with different topologies.

Delay induced instability
In this section we give the stability analysis of network organized systems with time delay. The analysis method is performed in close analogy to the method proposed in [3,14] applied to continuous model except some details. Assuming τ to be small we replace Expanding in Taylor series and neglecting the high-order nonlinearities, equation (3.1) becomes be the positive equilibrium of system (2.1) defined by the positive solution of We now consider small spatiotemporal perturbations to the fixed point * * By expanding the reaction terms around this equilibrium in a Taylor series up to first order, it follows from the first equation of (2.4) and equation Note that the Laplacian of a network L is a real, symmetric and negative semi-definite matrix, its eigenvalues are real, non-positive. The eigenvalue μ s and its corresponding eigenvector Φ s are determined by LΦ s =μ s Φ s with we obtain, for each mode s, the following matrix equation for linear growth rates λ s : Thus λ s are given by the roots of characteristic polynomial l t m l t m - It is well known that the stability condition for system In the absence of time delay, following Turing's idea [25,32], the equilibrium was originally stable, but become unstable under the influence of diffusion so that biological patterns formed. The former is equivalent to b g The latter requires that B(0, μ s )<0 for some mode s. We now consider the first condition in (3.4) is not satisfied to explore instability induced by time delay for system (2.4). Let the homogeneous equilibrium for the system without delay be stable, i.e. condition (3.5) is satisfied. Note that  m 0 s , the inequality A(τ, μ s )>0 holds if In summary, (3.5)-(3.6) are sufficient conditions for instability induced by time delay. Under these two conditions Turing patterns can be observed.

Wave patterns on networks
In this section, we perform extensive numerical simulations on three different types of network with N=10 000 nodes to investigate pattern formations of system (2.4). Each type of network include several cases with different average degrees as follows.
(1) Two-dimensional (2D) lattices: N nodes are equally arranged in a lattice with M columns, M rows and the distances between columns (rows) are h=1. In the first case, two nodes are connected if the distance between them equals to 1. We name this case 'LA4' because that its average degree is á ñ = k 4. In the second and third cases, two nodes are connected if the distance between them is less than or equal to 2 and 3 respectively, which leads to á ñ = k 12 and á ñ = k 23, we name them 'LA12' and 'LA23' respectively. See  Here randN (0, 1) denotes random number of standard normal distribution. The initial condition (I) means that the density of prey and predator in all nodes is a random perturbation around * * ( ) u v , . The condition (II) we assume that the density is a small perturbation around * * ( ) u v , except the nodes on some lines. For 2D lattices, if one node in m column n row then i=(m−1)M+n with m, n=1, L, M, N=M 2 . For ER and BA networks, we firstly arrange the nodes equally in a M×M lattice and order them according to decreasing degree from left to right, bottom to top, then m and n are column and row indices of node i respectively.

Wave patterns on 2D lattices
We study simulation results on 2D lattices in this subsection. Figures 2(a)-(c) show that periodic spiral waves are generated starting from initial data (I) on 2D lattices. The spiral waves appear as polycentric microspirals owing to the influence of initial condition (I) which is a small random perturbation around the positive equilibrium for all nodes. Comparing figure 2(a)-(c), it can be easily noticed that the spatial frequency of microspirals on LA12 is smaller than that on LA4 and larger than that on LA23. Specifically, the microspirals on LA4 have more centers than those on LA12 and LA23. In figures 2(d)-(f), single-arm spiral waves are generated on 2D lattices with initial data (II). Here again we see that the spatial frequency of spiral wave on LA12 is smaller than spatial frequency on LA4 while larger than that on LA23. The most essential reason for this phenomenon is that the nodes on LA12 and LA23 communicate more with other nodes. We show average densities of all nodes for prey on 2D lattices with initial values (I) and (II) in figures 3(a)-(c) and (d)-(f) respectively. It is apparent that average densities display cyclic variation over time and the oscillations of average density on LA12 and LA23 are more regular.
Changes in the average density of prey have followed a cyclical wave pattern. The two main characteristics of wave pattern are amplitude (i.e. peak and trough) and period. We next explore the influences of network topology on wave pattern from these two characteristics. Numerical results in figure 4 (see black lines) illustrate that average densities on 2D lattices maintaining a steady state when Î [ ] t 300, 1000 . This stabilization displays in two aspects: Firstly, with the increases of average degree of 2D lattice, the differences of maximum values, minimum values and amplitudes of average density are very small. In fact, the maximums range from 0.29 to 0.33 and the minimums from 0.26 to 0.29 on LA4, LA12 and LA23 (see figure 4(a), (c)). Secondly, the periods of average density are almost constant (approximately 8, see figures 4(b), (d)) as the average degree changes.

Wave patterns on random networks
In this subsection we investigate simulation results on random networks. Under the present parameters setting, the prey densities, similar to results on 2D lattices, display periodic variation on ER and BA networks. Figures 5  and 6 show prey densities at different times with initial condition (I) on ER4, ER10 and BA4, BA10 networks respectively; figures 9 and 10 show the corresponding results with initial condition (II). The prey distributions  on network nodes are roughly separated into one group, while only a limited fraction of nodes do not belong to this group (see panels (a)-(c) in figures 5, 6 and 9, 10). Along with the increases of average degree of random network, the number of nodes which dissociate from the group becomes fewer and fewer (see panels (d)-(f) in figures 5, 6 and 9, 10). In other words, the differences of prey density among most nodes are very small at a definite time, especially on BA10 network, all nodes have the same prey density (see panels (d)-(f) in figures 6 and 10). Figures 7 and 11 show average densities of all nodes for prey on ER =  i i 2 , 2, , 7 networks with initial conditions (I) and (II) respectively, from which periodic wave patterns can be observed. It is intuitively plausible that, on ER4, ER6, ER8, ER10 and ER12, the waves with larger average degrees have higher peaks, lower troughs as well as longer periods. Indeed, red lines in figure 4 confirm these observation results. Take simulation results on ER networks with initial condition (II) for example, the maximum values of prey average density on ER4, ER6, ER8, ER10, ER12 are 0. 6780, 0.8329, 0.9062, 0.9658, 1.0000, the minimum values are 0.0883, 0.0320,  0.0192, 0.0164, 0.0157, the periods are 9.72, 10.58, 11.61, 13.45, 152. It is worth emphasizing that the period has  a rapid increase when the maximum value increases to 1, thereafter the prey-predator system reaches a steady state: the maximum is stabilized at 1, the minimum at 0.016 and the period at 152. Moreover, all nodes have an identical density at any given time once the system enters into the steady state, which implies that figures 7, 11(e), (f) are, in fact, density evolution of each node. Figures 8 and 12 show prey average densities on BA2i, i=2, L,7 networks with initial conditions (I) and (II) respectively. One can observe that changing initial condition or network type do not alter the steady state. The differences are that the amplitudes and periods of average density on BA networks approach steady state more quickly (see blue lines in figure 4) which can also be confirmed via numerical results on ER10 network and BA10 network (see panels (d) in figures 7, 8 and 11, 12). These results demonstrate that average degree of random network has a significant effect on wave patterns, and the effect becomes more prominent on BA networks.
Through extensive simulation experiments, we found that the prey-predator system (2.4) can achieve stability both on 2D lattices and on random networks with some certain average degrees. However, there are fundamental differences between these two steady states. For steady state on 2D lattices, the amplitudes and  periods of wave patterns are small and the prey can keep a certain number. While for steady state on random networks, the wave patterns have long periods, the prey density in a period is kept to peak value 1 most of the time, and near 0 in a very short time. In addition, average degree has only a small effect on the steady state for 2D lattices, but has significant effects for random networks. Trace its root, nodes of 2D lattices have the same degree and share the same important place in the RD process. However, node degrees of random networks especially for BA networks can vary widely which leads to these nodes have different status and play different roles in the RD process.
From the viewpoint of ecology, the steady state on random networks is a double-edged sword: on the one hand, population number maintains a high level for much of the time is benefit to protect species diversity. On the other hand, biological invasion may lead to species extinction at the time when prey density near 0. Thus, understanding network topology of population distribution is significant for us to maintain the ecological balance.

Discussion
We have studied instability induced by time delay for a predator-prey model with modified Leslie-Gower and Holling-type III schemes on complex networks. Instability conditions were obtained via linear stability analysis of network organized systems which are performed similarly with continuous model. Extensive simulations show that wave patterns can be generated for the present model if time delay satisfies instability conditions. Because of considering network organized systems, wave patterns are influenced by network topology inevitably. We explored the influences mainly from two aspects, i.e. wave amplitude (the size of peak and trough) and wave period using deterministic networks-2Dlattices and random networks-ER networks and BA networks.
Spiral waves for prey densities can be generated on 2D lattices and spatial frequencies of waves decrease with the increases of average degree. Even so, prey average densities on LA4, LA12 and LA23 are still in a steady state.  This is shown in two aspects: the amplitudes of average density are very small (the peaks and troughs of average density equal approximately to 0.3) and the periods change a little (equal to about 8). In this sense, 2D lattices have little effect on wave patterns. In contrast, numerical results on ER and BA networks particularly the latter suggest that the effects of these two networks on wave patterns are huge. In the beginning, the amplitudes and periods increase with the increases of average degree of ER/BA networks. Then, there is a sudden-change of period (increase from 14 on ER10 and 12 on BA8 to 152 on ER12 and BA10) when the maximum value of average density up to 1. After that, the prey average density enters into another steady state. Compared with the steady state on 2D lattices, the former has a higher peak (equals to 1), lower trough (equals to about 0.016) as well as longer period (equals to about 152). The characteristic of results on ER and BA networks has two sides from the point of view of ecosystem equilibrium and environmental protection. In particular, protection efforts should be enhanced when population density at the wave trough to prevent population extinction.
In the process of instability conditions derivation we assume that the time delay is small so that a Taylor expansion could be employed to reduce system with time delay to those without time delay. Then Turing's idea can be followed. In practice, however, time delays in biological systems may be large. One of our future work is to study the influences of large time delay and network topology (such as community structure, clustering coefficient and average degree) on pattern formations. Besides, metapopulation epidemic models can be described by RD process on networks. Study of disease spreading using network organized system is also our further work.
Appendix. From continuous model to discrete model Generally, the process of pattern forms in space for continuous model is an asymptotic behavior with large scale and long time, therefore setting the mesh size h=1 is acceptable when b?1 (see figure 1(a)). According to Taylor expansion, one has ¶ ¶ = -+ + + - Denoting u mn is the finite difference approximation to u(x m , y n ) and using (A.1) with h=1, the spatially semidiscrete scheme for prey equation in ( Thereby u 01 =u 10 =u 11 and (A.3) can be rewritten as A . 4 11 11 11 11 12 21 Other boundary nodes are handled similarly. For node (m, n) we denote i=(m−1)M+n with m, n=1, L, M, then the spatially semidiscrete scheme for prey equation is given as follows, Here the matrix L=(l ij ) N×N is discrete Laplacian operator, its elements satisfy