Pattern dynamics of a delayed eco-epidemiological model with disease in the predator

The eco-epidemiology, combining interacting species with epidemiology, can describe some complex phenomena in real ecosystem. Most diseases contain the latent stage in the process of disease transmission. In this paper, a spatial eco-epidemiological model with delay and disease in the predator is studied. By mathematical analysis, the characteristic equations are derived, then we give the conditions of diffusion-driven equilibrium instability and delay-driven equilibrium instability, and find the ranges of existence of Turing patterns in parameter space. Moreover, numerical results indicate that a parameter variation has influences on time and spatially averaged densities of pattern reaching stationary states when other parameters are fixed. The obtained results may explain some mechanisms of phenomena existing in real ecosystem.


1.
Introduction. From the earliest Kermack-McKendric compartment model, epidemiological studies have made great progress, which provides some guidance for the prevention and control of disease [20,2,4,16,8,35,56,1]. The susceptible can be infected by direct contact with the infective, or by encountering with the freeliving pathogen in the environment. In real ecosystems, species does not live alone. There are mutualistic symbiosis, competition and predation relationships between interacting species [7,13,9]. Thus, it is more biologically meaningful to consider the interactions between species. Now, the population dynamics represented by Lotka-Volterra system has received great attentions [47,48,12,23,33,54,37,25]. Population dynamics provides some suggestions for the development and utilization of resources, conservation of biodiversity, and biological control. For example, A fishery's model with two types of fishermen was proposed, which evaluated the effects of their activities on the sustainability to fishing stocks at the Amazonian floodplain lakes [34].
Most species suffer from various infectious disease, and these diseases have significant effects on population size. Therefore, the combination of epidemiology and the interacting species has aroused more and more interest, which produce a new branch of eco-epidemiology [45,14,10,53]. For example, some scholars have studied how control the number of harmful insects by biologically-based technology. The combination of the prey-predator model and the epidemic model is one of the eco-epidemiology models. In recent decades, some scholars have considered that disease spreads only in the prey species [17,50,24], or only in the predator species [46,53,15], or in both the prey and the predator [19]. In this article, the classical Rosenzweig-MacArthur prey-predator model with Holling type II functional response is [31]: Hilker and Schmitz introduced a general SI epidemic model with standard incidence rate in Rosenzweig-MacArthur prey-predator model [18,6], then a prey-predator model with disease in the predator species was formulated: , Since most species live in spatial environments, they do not always stay in one place, but go to other places. In nature, a common phenomenon is that the predator searches for and catches the prey through movement, in turn the prey avoids being eaten by escape. So the diffusion of the species in space should be concerned. Previous papers considered the interactions between the species are independent of space position [22,52,17,46,18,6]. In fact, the distribution of the population in space is not homogeneous because of the interactions between the dynamics and the movement of the species. The spatial distribution of the species may form a variety of patterns, such as spotted patterns, strip patterns, and the coexistence of spotted and stripe patterns [30,5,43,27,38,36,42,41]. However, most of the previous studies discussed spatial pattern formations for the two-variable system. Sun et. al investigated the spatiotemporal complexity in a predator-prey system with Holling-Tanner form, and they demonstrated that migration had important effect on the pattern formation of the population [43]. Since phenomena of the interactions among multiple species exist widely. Although some researchers have discussed the spatial pattern formations of the three-variable systems [12,49,51,44]. In real ecosystem, predator can not immediately reproduce the next generation after eating prey [11,29]. In this paper, we therefore introduce time delay into a spatial eco-epidemiological model with disease in the predator.
The structure of this paper is as below. In section 2, a spatial eco-epidemiological model with delay and disease in the predator is formed. Moreover, we analyze the local stability of equilibria. In section 3, we derive the characteristic equation of the spatial eco-epidemiological model with delay and disease in the predator by the assumption of small time delay. The instability conditions of the positive constant stationary state are given. In section 4, numerical simulations exhibit patterns. Finally, we give some conclusions and discussion.
2. Mathematical modeling and analysis.

Model formulation.
In real ecosystem, the growth of the species and the interactions between the species depend on the spatial environment. Hence, the spatial Rosenzweig-MacArthur prey-predator model with disease in the predator is considered. Since we are interested in the self-organization of patterns, then the nonzero initial condition and Neumann boundary conditions are chosen. Neumann boundary conditions denotes that the species lives in a closed space domain. It is well known that the predator takes some time to produce the next generation after eating the prey, that is, the growth rate of the predator depends on the number of the predator and of the prey at some previous time [11,29]. So we introduce time delay into the interaction term between the predator and the prey in the healthy predator equation. By incorporating time delay into the spatial Rosenzweig-MacArthur preypredator model with disease in the predator, the following model in two-dimensional space is formed: where P 1 (x, y, t), S 1 (x, y, t), I 1 (x, y, t) represent the numbers of the prey, the healthy predator and the infected predator, respectively. b and K are the intrinsic growth rate and the carrying capacity of logistic equation, respectively. The functional response of predators is considered as Holling type II with a maximum consumption rate a, the half-saturation constant H and the conversion efficiency e. The contact transmission rate is β, d denotes the natural mortality, and α is disease induced mortality. In order to describe the diffusion process, the usual Laplacian operator ∇ 2 = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 in two-dimensional space is given. d 1 , d 2 and d 3 are diffusion coefficients for different species. n is the outside normal vector of ∂Ω. φ 1 (x, y, t), φ 2 (x, y, t), φ 3 (x, y, t) denote the initial number of the prey, the healthy predator and the infected predator, respectively.
2.2. Local stability. Next, system (4) in the absence of diffusion and delay is corresponding to the following ordinary differential system: The system (5) has two equilibria, one is the predator without disease E 1 = (− hm m−1 , − hr(hm+m−1) , 0), which needs to be satisfied m < 1/(h + 1). The other is the predator with disease equilibrium E * = (P * , S * , I * ), where And the condition of the existence of E * is .
The jacobian matrix of the system (5) at E 1 is the corresponding characteristic equation is it is obvious that a eigenvalue is λ 11 = β 1 − m − µ, the remaining two eigenvalues λ 12 , λ 13 are determined by holds, then the predator without disease equilibrium E 1 is locally asymptotic stability.
The Jacobian of system (5) around the predator with disease equilibrium E * is given by where a 1 = −(a 11 + a 22 + a 33 ), a 2 = a 11 a 22 + a 11 a 33 − a 12 a 21 + a 22 a 33 − a 23 a 32 , a 3 = −a 11 a 22 a 33 + a 11 a 23 a 32 + a 12 a 21 a 33 − a 12 a 21 a 32 .
By Hurwitz criterion, the necessary and sufficient conditions, under which all the roots of the characteristic polynomial (7) have negative real parts, are given by which indicates that the equilibrium E * is locally asymptotic stability 3. Analysis of the delayed eco-epidemiological model with disease in the predator. In this section, we mainly analyze the stability of the predator with disease equilibrium E * = (P * , S * , I * ) in system (4). Based on the assumption of small time delay, P (x, y, [32,26].
The following system is given by: By using Taylor series, we expand the three equations of the system (9) in (P, S, I) and neglect the higher order nonlinear terms, the system (9) is transformed into: where g(P, S, I) = P h+P (S + I), g P = ∂g(P, S, I)/∂P , g S = ∂g(P, S, I)/∂S, g I = ∂g(P, S, I)/∂I. The system (10) can be written as Since the predator with disease equilibrium E * = (P * , S * , I * ) satisfies f 1 (P * , S * , I * ) = 0, f 2 (P * , S * , I * ) = 0 and f 3 (P * , S * , I * ) = 0. By expanding the three equation of the system (11) with Taylor series in E * and neglecting the higher order nonlinear terms, then the system (11) becomes where g P = a 21 , g S = g I = −a 12 . Let P = P * + δP , S = S * + δS and I = I * + δI, where δP , δS and δI denote small spatiotemporal perturbations at homogenous steady state E * = (P * , S * , I * ). Furthermore, one can derive: Under the Neumann boundary conditions, spatiotemporal perturbations δP , δS and δI can be expressed as the following forms: By inserting (14a)-(14c) into system (13), we can derive matrix equation with respect to eigenvalues: here k 2 x + k 2 y = k 2 . Next, the characteristic equation of system (4) is derived: where a 0 = −a 12 τ + 1 > 0 because of a 12 < 0, here Then we can get that: where  In the following, we define f (k, τ ) := a 1 (k, τ )a 2 (k, τ ) − a 3 (k, τ ).
Proof. a 1 (k, τ ) > 0 because of a 12 < 0. However, the signs of a 3 (k, τ ) and f (k, τ ) are indefinite. If at least one of a 3 (k, τ ) < 0 for some k = 0 and f (k, τ ) < 0 for some k = 0 is established, then the predator with disease equilibrium E * is unstable based on Hurwitz criterion.
Analogously, the other Turing bifurcation function T2 is: On the basis of Theorem 3.2 and Theorem 3.3, the conditions of diffusion induced Turing instability can be given by:   FIG. 2(a). In part 1, the positive constant stationary state E * is unstable because of Hopf bifurcation instability. In part 2, Turing instability and Hopf bifurcation instability lead to instability of E * . In part 3, diffusion induces the instability of E * , which is Turing instability. FIG. 2(b) is divided into four regions by Hopf bifurcation function H and Turing bifurcation function T1. Since Turing instability and Hopf bifurcation instability occur in region 1, then E * loses stability. E * becomes unstable because of Turing instability in region 2. E * is unconditional stable in region 3. E * loses stability only for Hopf bifurcation instability in region 4.
According to Theorem 3.3, if at least one of the conditions a 3 (k, τ ) < 0 for some values of k = 0 and a 1 (k, τ )a 2 (k, τ ) − a 3 (k, τ ) < 0 for some values of k = 0 holds, the parameter domains of Turing instability can be found. Thus, FIG.3  and FIG. 4 give coefficients of the dispersion relation of the characteristic equation  . 3 and FIG. 4, it is easy to see that a 1 (k, τ ) and a 1 (k, τ )a 2 (k, τ ) − a 3 (k, τ ) are always greater than 0 for all k 2 > 0, but there exist some k 2 = 0 satisfying a 3 (k, τ ) < 0.
Consequently, the conditions of delay induced instability are: 4. Numerical simulations. In order to verify theoretical results of the stability of the positive constant stationary state E * in system (4). Through using explicit Euler method, we have carried out a series of numerical simulations in two dimensional spaces. The reaction-diffusion system (4) is calculated in a discrete domain with N x × N y lattice units, and the space between the lattice points is constant dx. The time evolution is also discrete, that is, the time goes in steps of dt. In numerical simulations, we set N x = N y = 200, grid space dx = 1.0 and time step dt = 0.01. We mainly focus on the distribution of healthy predator in two dimension space domain. Therefore, the pattern formations of healthy predator are presented. If simulation results reach stationary states or pattern behaviors do not change with time, then simulations finish. In the following figures, the scale of the colorbar denotes the density of healthy predator. The density of healthy predator increases with the color of the colorbar changing from blue to red. In this section, diffusion-driven instability of the equilibrium E * arises Turing patterns which is symmetry-breaking spatial structures. According to conditions (25), the following two sets of parameters are given by:  6(b). This set of parameters shows small "black-eye" patterns and the corresponding spatially averaged population densities for different the values of r. At the initial time, the healthy predator uniformly distributes in space domain. As the time is being increased, the spatial distribution of the healthy predator changes. At last, the healthy predator gathers together, which forms small "black-eye" patterns. But there are some slight differences in the local space distributions of small "black-eye" patterns for (a) and (b). From figures of spatially averaged population densities, it can be found that the time and spatially averaged densities of pattern reaching stationary states are different. β 1 = 1.8, µ = 0.6, m = 0.7, r = 0.1, D 1 = 10, D 2 = 0.1, D 3 = 4, τ = 0.01. This set of parameters presents big "black-eye" patterns and the corresponding spatially averaged population densities for different the values of h in FIG. 7. At first, the healthy predator is randomly distributed in space. The spatial distribution of healthy predator changes with time evolution. Finally, the healthy predator is aggregated, which gives rise to big "black-eye" pattern. Figures of spatially averaged population densities demonstrate that the time and spatially averaged densities of pattern reaching stationary states are not same. Moreover, the local space distributions of big "black-eye" patterns are slightly different.   . Spatial patterns (top) and the corresponding spatially averaged population density (bottom). (a) Big "black-eye" pattern (h=0.1), (b) big "black-eye" pattern (h=0.14).

5.
Conclusions and discussion. The interacting species play important roles in real ecosystem. Furthermore, the disease widely exists in different species. Thus, the eco-epidemiology well combines population dynamics with epidemiology. In this paper, we consider a spatial eco-epidemiology with delay and disease in the predator species under the Neumann boundary condition, which is a three-variable reactiondiffusion system with delay. But a little attention had been paid to study the pattern dynamics of three-variable reaction-diffusion system with delay. Through analyzing the characteristic equation of system (4), diffusion and delay can lead to the results that the predator with disease equilibrium E * loses stability, respectively. However, diffusion-driven instability and delay-driven instability is two different mechanisms. Besides, Turing pattern region is obtained in parameter space r − h, that is, we can find the parameter ranges of pattern formations. From numerical simulations, it is found that a parameter variation affects time and spatially averaged densities of pattern reaching stationary states when other parameters are fixed.
In this paper, in the process of analysis of instability of the equilibrium E * , we derive the characteristic equation of system (2) by assuming small delay τ . But for general delay, the characteristic equation with e −λτ is derived, then it is difficult to analyze the conditions of the equilibrium losing stability. Hence, this work is worth investigating. Pattern phenomena widely exist in nature, such as, the patterns of animal skins, the regularly spatial structure distribution of the infected, mussel bed on intertidal flats and so on [3,21,28,39,38,27,40]. It is very meaningful to study the mechanism of the existence of various patterns in the natural world.