NUMERICAL SOLUTION OF A SPATIO-TEMPORAL PREDATOR-PREY MODEL WITH INFECTED PREY

A spatio-temporal eco-epidemiological model is formulated by combining an available non-spatial model for predator-prey dynamics with infected prey [D. Greenhalgh and M. Haque, Math. Meth. Appl. Sci., 30 (2007), 911– 929] with a spatio-temporal susceptible-infective (SI)-type epidemic model of pattern formation due to diffusion [G.-Q. Sun, Nonlinear Dynamics, 69 (2012), 1097–1104]. It is assumed that predators exclusively eat infected prey, in agreement with the hypothesis that the infection weakens the prey, making it available for predation otherwise we assume that the predator has essentially no access to healthy prey of the same species. Furthermore, the movement of predators is described by a non-local convolution of the density of infected prey as proposed in [R.M. Colombo and E. Rossi, Commun. Math. Sci., 13 (2015), 369–400]. The resulting convection-diffusion-reaction system of three partial differential equations for the densities of susceptible and infected prey and predators is solved by an efficient method that combines weighted essentially non-oscillatory (WENO) reconstructions and an implicit-explicit RungeKutta (IMEX-RK) method for time stepping. Numerical examples illustrate the formation of spatial patterns involving all three species. Future directions of research are suggested.


Introduction.
1.1.Scope.Mathematical models are powerful tools to characterize and gain insights into the spatial-temporal dynamics of complex systems.More specifically, mathematical models are frequently used to gain a qualitative understanding of the dynamics and control of natural and social dynamic phenomena [44].For example, systems of mathematical models are often simple and crude approximations that aim to capture key mechanisms underlying the dynamics of infectious disease transmission in human, animal, and plant populations.Such first approximation models can be subsequently refined to study the effects of additional factors and processes including spatial heterogeneity in density of hosts with different characteristics (e.g., fitness, susceptibility, infectivity) as well as the role of specific host behaviors (e.g, migration and diffusion patterns) in response to local host density or disease prevalence.
It is the purpose of this contribution to advance a spatio-temporal predatorprey model taking into account infected prey with a particular focus on its efficient numerical solution and dynamic characteristics.The resulting eco-epidemiological model combines several submodels that have been proposed in the recent literature.The population of prey is described by a spatio-temporal SIR-type epidemiological model similar to the one studied by Sun [45] while the predator-prey interaction follows the treatment by Greenhalgh and Haque [23] (formulated in a non-spatial setting in that paper).The movement of predators is based on the nonlocal velocity model advanced for the predator movement in a spatio-temporal predator-prey model by Colombo and Rossi [16].That velocity model was also used for the male compartments in a recently published model of hantavirus infection [12].In particular the implicit-explicit numerical technique for efficiently solving the governing partial differential equations (PDEs) is adapted to the new model presented herein.
In most general terms, the governing model considers a prey population with a susceptible and infective compartment, denoted by S and I, and a variable Y denoting the population of the predator.In the spatio-temporal context, these variables are understood as local densities that are functions of position x ∈ Ω and time t ∈ T := [0, T ] on a bounded domain Ω ⊂ R 2 .The final model for u := (u 1 , u 2 , u 3 ) T := (S, I, Y ) T , u = u(x, t) is given by a convection-diffusion reaction system of the type supplied with initial and boundary conditions, where the convective fluxes F c (u), the diffusion matrix D, and the vector of reaction terms s(u) are specified in later parts of the paper.
It is proposed to describe the movement of the predators in a particular way that depends non-locally on the density I of infected prey, while the movement of infected and susceptible prey, that is of members of the compartments I and S, is standard diffusion.This property calls for numerical methods that allow the efficient computation of numerical fluxes based on the non-local evaluation of data but avoid the severe time step restriction incurred by explicit time discretizations of the diffusion term D∆u.As in [12] the first goal can be achieved by a technique based on Fast Fourier Transforms (FFT) in combination with an implicit-explicit (IMEX) discretization to handle the second.The numerical method constructed in this way is applied to simulate several scenarios that allow comparing the effect of initial conditions and different parameter values.The numerical results exhibit a rich variety of behavior, including formation of patterns taking the shapes of spots, stripes, both from randomized perturbations of steady states and in a scenario that describes the invasion of a prey habitat by the predator.The structure of solutions obtained in the latter case are reminiscent of permanent wave fronts.
1.2.Related work.The spatial spread of infectious diseases is described mathematically in a number of monographs that include [18,35,36,41,48].The temporal evolution of diseases is also treated in [10,21].Our assumption that all epidemiological compartments are distributed over the whole spatial domain is opposed to the alternative metapopulation approach that describes spatial structure through wellidentified sub-populations or "patches" (cf., e.g., [1][2][3]11,46,47]).The description of spatial structure by explicitly specifying the mobilities between "patches" is typical for characterizing the behavior of humans who undertake directed travels, while a description through a convection-diffusion-reaction mechanism is more suitable for non-human infectious agents such as spores, insects, and bacteria that would disperse [18,36,37].A decisive advantage of the spatially continuous approach, namely its amenability to mathematical analysis, is emphasized in [49].
The introduction of [23] and Hadeler and Freedman [24] provide extensive references to real-world examples of three-species eco-epidemiological systems of sound prey, infected prey, and predators.Specifically, in [24] a predator-prey model is described in which the prey is infected by a parasite and the prey in turn infects the predator with that parasite.The models studied in [23] and herein are based on the general assumption that the infection weakens the prey and increases its susceptibility to predation.Furthermore, we mention that the eco-epidemiological system of the Salton Sea in the desert of Southern California, New Mexico, is an important eco-epidemiological system where birds (particularly pelicans) prey on fish (particularly tilapia).The mechanism of predation is explained in [22]: "The vibrio class of bacteria is very common in salt water fish, and those infected with vibrio may have salt water present in their tissue.As the fish struggle in their death process, they tend to rise to the surface of the sea for oxygen.When they do, they become very attractive to fish-eating birds, specifically the pelicans."Mathematical models of this system were proposed and analyzed by Chattopadhyay and collaborators in a series of papers [14,15,30,40].
A less common ingredient in mathematical epidemiology is the convective term ∇ • F c (u). Related advection terms, for which the essential functional dependence for one compartment is F c = F c (x, u) = b(x)u with a given velocity b(x), arise if the population under study is transported (as is the case with wind-borne infectious agents, plankton, etc.).In our model the main role of the convective term is similar to that of [36,Ch. 14] in that it imposes a preferred direction of movement of predators, where the global dependence of the direction is determined by convolution with population data within a certain horizon of the current position.This idea of non-local dependence of biological fluxes goes back to a non-local predator-prey model introduced and analyzed by Colombo and Rossi [16], and for which a numerical scheme was analyzed in [39].The scheme of [39] is based on the Lax-Friedrichs scheme for the ease of demonstrating convergence properties; we here employ higher-order weighted essentially non-oscillatory (WENO) reconstructions [26,34] to achieve high-order spatial accuracy.
Furthermore, IMEX Runge-Kutta (IMEX-RK) schemes play an important role for the efficient numerical solution of (1.1).Roughly speaking, an IMEX-RK method for the convection-diffusion-reaction equation (1.1) consists of a Runge-Kutta scheme with an implicit discretization of the diffusive term D∆u combined with an explicit one for the convective and reactive terms, F c (u) and s(u), respectively.To introduce the main idea, we assume that represents a method-of-lines semi-discretization of (1.1), where Φ * (v) and Φ(v) are spatial discretizations of the convective and reactive terms and of the diffusive term, respectively.Assume that the spatial mesh width is h > 0 in both the x-and ydirection.Then the stability restriction on the time step ∆t that explicit schemes impose when applied to (1.2) is very severe (∆t must be proportional to h 2 ), due to the presence of Φ(v).The implicit treatment of both Φ * (v) and Φ(v) would remove any stability restriction on ∆t.However, the upwind nonlinear discretization of the convective terms in Φ * (v) that is needed for stability makes its implicit treatment extremely involved.This situation becomes even more complicated due to the use of WENO reconstructions [26,34,42,51].However, numerical integrators that deal implicitly with Φ(v) and explicitly with Φ * (v) can be used with a time step restriction dictated by the convective-reactive term alone.These schemes have been profusely used in convection-diffusion problems and convection problems with stiff reaction terms [4, 6-9, 17, 19, 28, 38, 52].
1.3.Outline of the paper.The remainder of the paper is organized as follows.
In Section 2 the governing mathematical model is introduced.To this end, we first formulate (in Section 2.1) a non-spatial predator-prey model, defined by three nonlinear ordinary differential equations (ODEs), in which an epidemic spreads in the prey.This model gives rise to five equilibrium points whose respective stability is discussed in Section 2.2, along with a discussion of the Hopf bifurcation associated with the unique stable one of the five equilibria that includes all three species.(These stability results will motivate the choice of parameters of the simulated spatio-temporal scenarios.)Then, in Section 2.3, we obtain the full spatio-temporal model by equipping the ODE model with diffusive terms for the prey species and a non-local convection term for predators.These ingredients are specified in detail in Section 2.4.Section 3 is devoted to the description of the numerical scheme proposed to solve the spatio-temporal model developed in Section 2.3.Specifically, the (standard) spatial discretization of local convection and diffusion terms is introduced in Section 3. 2. Mathematical model.

2.1.
Dynamical system for a predator-prey model.Combining the SIR epidemic model of [45] with a second species corresponding to a predator (as in the model of [23]), and assuming that the predator only eats infected prey.That is, illness arising from infection debilitates the prey, making it available for predation otherwise we assume that the predator is not sufficiently fit to access healthy prey.Accordingly, we obtain the following model: where t is time, S(t) is the population of susceptible prey, I(t) is the population of infected prey, A is the recruitment rate of the prey population, d is the natural death rate of the prey population, β is the force of infection or the rate of transmission, µ is the disease-related death from the infected, Y (t) is the population of predators, c is the search rate of the predators towards infected prey, δ is the per capita growth rate of the predators, h is a constant related to the density dependent mortality of the predator population, and m > 0 is a constant.The incidence rate in (2.1a) and (2.1b) is chosen here as the nonlinear expression βSI 2 in agreement with Sun [45], who in turn appeals to the justification by Liu et al. [32,33].The last term in the right-hand side of (2.1b) is a ratio-dependent predation term; see, e.g., [29] for further justification of such expressions.Overall, the right-hand sides of (2.1a) and (2.1b) have been chosen such that they coincide with the reactive (algebraic, i.e., non-differential) terms of the governing model of [45] (Equations (1a) and (1b)) of that paper, so that results can in principle be compares.Moreover, the assumption of a constant recruitment rate A (rather than a per capita growth rate for the prey population) simplifies the subsequent stability analysis.Furthermore, note that (2.1c), that is the assumption of a logistical growth rate of the predator population with a carrying capacity proportional to the population size of prey, coincides with the predator equation [23, Eq. ( 1)], and that hY /I in that equation is the Leslie-Gower term [31] which measures the loss in the predator population due to the relative scarcity of the (infected) prey [22].An example of a real biological scenario for which (2.1) would be an adequate description is the "pelican-tilapia" eco-epidemiological model of the Salton Sea mentioned in Section 1.2, where S and I correspond to the number of suseptible and infected tilapia (fish) and Y to that of the predating pelicans (birds).In particular, Chattopadhyay and Bairagi [14] and Sarkar et al. [40] argue that in this scenario the predator only eats infected prey (although later works study variants of the model advanced in [14] in which the predator also eats susceptible (sound) prey, see [15,22,30]).
It is well known that the use of a ratio-dependent terms requires carefully handling situations of zero denominators.The right-hand sides of the ODE version (2.1) are potentially ill-defined if I(0) = 0.This issue is treated in [23] as follows: if I(0) = 0 and Y (0) ≥ 0, then (2.1c) is interpreted as implying that Y (t) = 0 for t > 0 and if I(0) = Y (0) = 0, then (2.1b) is interpreted as implying that I(t) = 0 for all t.With respect to the PDE version (2.15), we have not encountered any difficulty of division by zero in our numerical experiments since in all cases we choose the initial datum I(x, 0) > 0 for x ∈ Ω.
Theorem 2.1.The system (2.1) has the following equilibrium points: i) The equilibrium E 1 = (A/d, 0, 0) corresponding to the extinction of the epidemic.
ii) The following equilibria in absence of the predator: where we define iii) The following equilibria with presence of the predator: where I 4 < I 5 are the solutions of the quadratic equation They are feasible if if this condition is satisfied, then we have

2.2.
Behavior of the dynamical system.The point E 1 corresponds to the extinction of the epidemic.On the other hand, the detailed analysis in [45] shows that E 2 is unstable and E 3 is stable.In these states the predator is absent, but we are especially interested in equilibria in which the predator does participate.To handle the dynamics of the system around E 4 and E 5 , we first recall that at E k , k = 4 or k = 5, the stability matrix is given by A straightforward calculation reveals that For I 4 given by the first expression of (2.4), we obtain, after rearranging terms, Since the term in squared brackets is positive, it follows that det J * (E 4 ) > 0, and since the entries of J * (E 4 ) are real, this means that at least one eigenvalue of J * (E 4 ) has positive real part.Thus, the equilibrium E 4 is always unstable.
To analyze the stability of E 5 , we utilize the Routh-Hurwitz criterion.To this end, we write the characteristic equation of J * as det(J * − λI) = 0 ⇔ λ 3 + e 1 λ 2 + e 2 λ + e 3 = 0, ( where e 1 = − tr(J * ), e 2 = 3 i=1 M i , where M i is the determinant of the matrix obtained by eliminating row i and column i from J * , and e 3 = − det J * .Here, the Routh-Hurwitz conditions (cf., e.g., [20,Sect. 6.4]; see also [50,Th. 2]) mean that there are only roots with strictly negative real parts, and hence the corresponding state is stable, if and only if e 1 > 0, e 3 > 0, and e 1 e 2 > e 3 . (2.7) This criterion is now applied to J * (E 5 ).From (2.5) for k = 5 we get so we always have e 3 > 0 for J * (E 5 ) wherever E 5 is feasible.
The subsequent analysis of e 1 and e 2 , as well as that of the existence and possible type of a Hopf bifurcation around E 5 , will be simplified if we fix the parameters where the values of A, µ, and d are those proposed in [45] while the choices of m, h, and δ are assumed but are comparable with those utilized in [23].In order identify scenarios of interest for the simulation of the spatio-temporal model we will only vary the infectious force of the epidemic expressed by the parameter β, and the search rate c of the predators towards infected prey.The stability of E 5 will be analyzed under variations of β and c only.
To determine suitable choices of β and c, we assume that e 1 = e 1 (β, c), e 2 = e 2 (β, c) and e 3 = e 3 (β, c) and analyze satisfaction of the Routh-Hurwitz condition (2.7).In addition, and to make the discussion as short as possible, we utilize numerical ecaluation to deduce that e 1 , e 2 , and a coefficient function γ that arises further below, are positive on a relevant (β, c)-region in each case.Combined with the feasibility condition (2.3) and the fact that e 3 > 0 wherever (2.3) holds, this discussion gives rise to four distinct regions as subregions of the region of the (β, c)-plane where the feasibility condition (2.3), the first condition in (2.7), and all conditions of (2.7) are satisfied, see Figure 1.
Finally, we demonstrate that a Hopf bifurcation occurs across the curve that separates R 3 from R 4 .Our reasoning partly follows that of [22,Th. 6]. Figure 1 illustrates that at points where the stability changes (either from unstable to stable or from stable to unstable), it is the constraint e 1 e 2 > e 3 that changes rather than e 1 > 0, where we recall that e 3 > 0 is generically satisfied wherever E 4 and E 5 are feasible.Furthermore, for the parameters (2.9) and (β, c) ∈ R := [30,50] × [0, 20], we also get e 2 > 0 wherever E 4 and E 5 are feasible.(We limit the further discussion to (β, c) ∈ R, and understand that R 1 to R 4 are subregions of R.) Let us show that across the curve c → β H (c), which is marked in blue in Figure 1, defined implicitly by e 1 (β H (c), c)e 2 (β H (c), c) − e 3 (β H (c), c) = 0, a Hopf bifurcation occurs, i.e., we fix c as a parameter and then consider β as a bifurcation parameter.Clearly, the characteristic equation (2.6) takes the form (λ 2 + e 2 )(λ + e 1 ) = 0 with the three roots and λ 3 = −e 1 .Thus, in a neighborhood of (β H (c), c), (2.6) cannot have real positive roots; rather, the roots are of the form λ 1 = u + iv, λ 2 = u − iv and λ 3 = −φ with φ > 0, and where u, v ∈ R and φ depend on (β H (c), c).It remains to verify the transversality condition, i.e.,

∂ ∂β
Re λ 1 (β, c) To this end, we substitute λ = λ 1 = u + iv into (2.6), equate real and imaginary parts, and calculate the derivative (see [22]) to obtain the following equations, where u = u(β, c), u = ∂u(β, c)/∂β, and analogously for v, e 1 , e 2 , and e 3 and quantities ρ = ρ(β, c), σ, τ and ν defined below.This yields where (2.11) ), it follows that σ = 0. Applying Cramer's rule to (2.10) we get ∂u(β, c) ∂c . (2.12) We wish to show that σν + ρτ = 0. Evaluating (2.11) at β = β H (c), we obtain at that state σν + ρτ = 2e 2 (e 1 e 2 + e 2 e 1 − e 3 ).We already know that e 2 > 0 in the region of interest.Furthermore, so that e 1 e 2 + e 2 e 1 − e 3 = γ(β H (c), c)∂(βI 2 5 )/∂β, where we define For the parameters (2.9), we obtain as we infer from the rightmost plot of Figure 2. (Indeed, γ is positive on a slightly larger subset of R.However, positivity on R 3 ∪ R 4 is sufficient for the discussion of the bifurcation across the boundary between R 3 and R 4 .)Furthermore, we observe that Combining (2.13) and (2.14), we deduce from (2.12) that ∂u(β, c)/∂c| β=βH(c) < 0, which proves the transversality condition.Thus, a Hopf bifurcation occurs across c → β H (c).We have conducted this analysis for E 5 as we are interested in the dynamics around an equilibrium when all three species are present.To make one of the eigenvalues of J * (E 5 ) zero, we must have det J * (E 5 ) = 0, but here we know that det J * (E 5 ) < 0 wherever E 5 is feasible.In particular, it cannot happen that J * (E 5 ) has a complex conjugate pair of purely imaginary eigenvalues with the third one being zero.Thus, the system does not experience a supercritical bifurcation, and in particular (2.1) does not admit time periodic solutions.The latter property is a mere consequence of det J * (E 5 ) < 0, which has been established symbolically and is therefore valid for any choice of parameters (instead of (2.9)) and (β, c) for which E 5 is feasible.
where ∇• denotes the (spatial) divergence operator.Here x and y are the space variables, ∆ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the two-dimensional Laplace operator, and D S and D I are the diffusion coefficients, i.e. we assume standard linear diffusion for the prey compartments S and I.
The right-hand sides of (2.15) are identical to that of the non-spatial ODE model (2.1), i.e., this model is recovered if all divergence terms on the left-hand sides are set to zero and variables are considered to depend on t only, and the unknowns represent suitably scaled densities.
is the non-local unscaled velocity function [16].Here η denotes a radial convolution kernel with radius ε, i.e., η is a piecewise smooth function such that i.e., for any function w defined on Ω × T and w(y, t)η(x − y) dy (with slight modifications for points x with dist(x, ∂Ω) < ε.) Since ∇(w * η) = w * ∇η, the velocity V (w) indeed depends (non-locally) on w and not on its gradient.Summarizing, we obtain F c (u) = 0, 0, F Y (I, Y )) T and D = diag(D S , D I , 0).The vector s(u) = (s 1 (u), s 2 (u), s 3 (u)) T is given by the right-hand sides of (2.15).The system (2.15) is considered on Ω × T along with the initial condition u(x, 0) = u 0 (x), x ∈ Ω, (2.18) where u 0 is a given function, and zero-flux boundary conditions where n is the unit exterior normal vector to the boundary ∂Ω of Ω.
3.1.Spatial discretization of local convection and diffusion terms.We take Ω = [0, L] × [0, L] and use a Cartesian grid with nodes (x i , y j ), i, j = 1, . . ., M , with u) by WENO finite differences and ∆u by the standard second-order scheme with a five-point stencil to get a spatial semi-discretization of (1.1) for a 3 × M × M -matrix v(t) of unknown approximations v ,i,j (t) ≈ u (x i , y j , t) for i, j = 1, . . ., M and = 1, . . ., 3.This spatial semi-discretization is given by Applying a suitable IMEX-RK scheme to (3.1) yields the final fully-discrete scheme (see Section 3.4).Here , to be defined in Section 3.3, and denotes the discretization of the diffusion terms.Here v denotes the M × M submatrix given by (v ) i,j = v ,i,j and ∆ h is the approximation of the standard twodimensional Laplacian operator with Neumann boundary conditions.Furthermore, S(v) is the 3 × M × M -matrix with components S(v) ,i,j = s (v ,i,j ), i, j = 1, . . ., M, = 1, . . ., 3, with corresponding submatrices S (v), given by S (v) i,j = s (v ,i,j ).

3.2.
Discretization of the convolutions.The following identity arises from (2.16) if we take into account that ∇(w * η) = w * ∇η: The corresponding convolutions w * ∂η/∂x and w * ∂η/∂y are calculated approximately on the discrete grid via a composite Newton-Cotes quadrature formula, such as the composite Simpson rule.Since B ε (0) ⊆ [−rh, rh] 2 , r = ε/h < M , and according to boundary conditions, w is extended to the exterior of Ω by reflection, e.g. by setting w(−x, y) = w(x, y), (x, y) ∈ Ω, we obtain (w * χ)(x i , y j ) ≈ r p=−r r q=−r β p,q w(x i−p , y j−q ), β p,q = h 2 α p α q χ(x p , y q ), (3.4) where α p and α q are the coefficients in the quadrature rule (e.g., for the composite Simpson rule, α = (1, 4, 2, 4, . . ., 2, 4, 1)).Consequently, the approximation (3.4) for W = (w i,j ) ∈ R M ×M , w i,j ≈ w(x i , y j ), is given by where we define The discrete approximation of V (w) in (3.3) is then given by Since r ≈ ε/h = εM/L, the computational cost of this discrete convolution is M 2 (2r + 1) 2 ≈ 4ε 2 M 4 /L 2 , which can be very high for large M .This cost can be substantially reduced to O(M 2 log M ) by performing a convolution with periodic data by Fast Fourier Transforms (FFTs) (see [27,43]).To achieve this goal, we define from W = (w i,j Therefore (3.5) for i, j = 1, . . ., M can be rewritten as The convolution on the right-hand side of (3.6) can be performed by FFTs applied to the (2M ) × (2M ) matrix W .To save further computational costs, the FFT of the kernel β p,q is performed only once, son each convolution entails two twodimensional FFT of (2M ) × (2M ) matrices and a product of 4M 2 numbers, with an overall computational cost of O(M 2 log M ).
3.3.Discretization of the convective term.The convective flux for the prey compartments (S and I, corresponding to ∈ {1, 2}) is zero and for the predator (species Y , = 3) is given by for the approximation v, we first approximate the convolution terms as expounded in Section 3.2 to obtain We introduce the notation (f x i,j , f y i,j ) := F c (v) i,j , where we have dropped the index in f x i,j etc. to obtain clearer expressions.Our purpose is to use a fifth-order WENO finite difference discretization [26,34,42] h for suitable numerical fluxes f x i+1/2,j , f y i,j+1/2 obtained by WENO reconstructions of split fluxes.For the numerical flux in the x-direction, the Lax-Friedrichs-type flux splitting f x,± is given by with an analogous formula for f y,± i,j .If R ± denotes fifth-order WENO upwind biased reconstructions and we use matlab-type notation for submatrices, then 3.4.Implicit-explicit Runge-Kutta schemes.To specify the IMEX-RK integrators for (3.1), we rewrite (3.1) as (1.2),where The diffusive part Φ(v) is handled by an implicit s-stage diagonally implicit (DIRK) scheme with coefficients A ∈ R s×s , b, c ∈ R s , in the common ¡ notation, where A = (a ij ) with a ij = 0 for j > i.For the term Φ * (v) we employ an s-stage explicit scheme with coefficients Â ∈ R s×s , b, ĉ ∈ R s and Â = (â ij ) with âij = 0 for j ≥ i.In our simulations, we limit ourselves to the second-order IMEX-RK scheme H-DIRK2(2,2,2) that corresponds to the pair of Butcher arrays Alternative choices are provided and discussed in [5,6,38].If applied to the equation (1.2), the IMEX-RK scheme gives rise to the following algorithm (see [38]).
Algorithm 3.1.Input: approximate solution vector v n of (1.2) for t = t n for p = 1, . . ., s a pj K j endif solve for K p the linear system Output: approximate solution vector v n+1 of (1.2) for t = t n+1 = t n + ∆t.
To solve the linear equation (3.8) for K p , in view of (3.7) we rewrite it as where I denotes the identity operator for 3 × M × M matrices.From the definition of the matrix B in (3.2) and from the definition of Φ * in (3.7), if we equate the submatrices along the first dimension of both sides of (3.9) we get where I M ×M is the M × M identity operator and, e.g., (K p ) is the submatrix of K p along the first dimension, i.e., ((K p ) ) i,j = (K p ) ,i,j .If µ = 0 (for = 3), then  in accordance with the region of stability shown in Figure 1 (these parameters will be varied within the numerical solution of spatio-temporal models).For the parameters (2.9) and (4.1), we obtain the following equilibria of Theorem 2.1 (besides E 1 = (1, 0, 0)) where the predator is absent: From Section 2.2 we know that all eigenvalues of J * (E 5 ) have negative real part, and that the equilibrium E 4 is unstable.As a consequence, whenever predators and prey are present, the solution orbits of (2.1) approach the stable equilibrium E 5 , as is illustrated in Figure 3.   Furthermore, we wish to compare the numerical results with the predictions obtained using the the non-spatial ODE model.To this end we determine for each compartment X ∈ {S, I, Y } and time instants t n = n∆t the following quantity: which represents the approximate total number in Ω of individuals of compartment X at time t n .We recall that in the PDE model, the unknowns X ∈ C are densities, and so an integration over Ω is necessary to make results comparable with those of a model that predicts the total number of individuals in each compartment (as does the dynamical system (2.1)).
Example 1: Case I in the absence of the predator.Here, only the equations (2.15a), (2.15b) are considered, with the parameters given by (2.9) and (4.1).The values of β and µ in combination with D S = 6 and D I = 1 have been chosen such that they lie in the Turing region, see [45,Example 6], which means that the formation of a permanent spatial pattern by the standard Turing mechanism in reaction-diffusion systems (see, e..g., [36]) is expected.In our simulation the initial conditions for S and I are those of Figure 4 (a).In Figure 5 we display the numerical solution at three different times.We observe the formation of a pattern of spots similar to that    values In this case the parameters β and µ together with the diffusion constants are not in the Turing space.The steady states according to Theorem 2.1 are now, besides E 4 = (0.7510048, 0.0848976, 0.8489759), E 5 = (0.2489953, 0.2560630, 2.5606299).( Considering the initial conditions of Figure 4 (a), we obtain in this case that no patterns are formed and the system quickly arrives at a constant equilibrium (uniform distribution for S and I in whole domain), which is maintained over time, as can be observed in Figure 9.This figure, as well as Figures 6 (g) and (h), indicate that the global equilibrium is E 3 given by (4.5).Example 5: Case II with a "triangular" initial distribution of the predator.In this case we consider the same scenario as in Example 2, namely the triangular initial condition of Figure 4 (b), but utilize the parameters (4.4) instead of (4.1).In this case, as in Example 2, predators are heading towards the infected prey, and as is shown in Figure 10, the whole domain becomes successively filled with a pattern formed of stripes.Note that in regions occupied by the prey before the arrival of predators the solution is constant (see Example 4), and in particular does not exhibit any formed pattern (in contrast to Example 2).In a sense this scenario illustrates how the model (2.15), under the appropiate choice of parameters, predicts the formation of spatial patterns among the prey upon arrival of the predators.Furthermore, we observe that a "front" of predators is moving into the domain initially occupied by prey only, which seems to move at a lower velocity than in Example 2. It is worth noting that the final distribution of prey corresponds to total S and I populations close to L 2 times the corresponding entries of E 5 given by (4.6), while as in previous cases the final population of predators is much higher than L 2 times the Y -component of E 5 .
Example 6: Case II with a random initial distribution of the predator.This case is the analogue of Example 2, namely we impose the initial conditiom of Figure 4 (c).
In Figure 11 we display the dynamics of the model with this initial condition.We observe that at simulated time t = 50, some stripes and spots emerge as patterns in the distribution of each species.Comparing the successive solutions for t = 50, t = 100 and t = 200, we observe that once the stripes form, they grow steadily in time incorporating "spots", until they reach a certain arm length, and the spatial which is maintained in time, see Figure 12.For this case we observe a behavior similar to that of Example 4.
Example 8.This case is based on the same initial "triangular" scenario (Figure 4 (b)) as Examples 2 and 5.The numerical solution (Figure 14) shows that the predators are heading towards the infected prey, and for t = 50 and t = 150 we observe that the distribution of the prey varies in comparison with the uniform distribution of Example 7 in those places were the predators are present.The final spatial configuration is marked by spots and stripes for the prey as well as for the predators, and remains in time.However, the structure of patterns that are formed differs from the results of Example 5.It is interesting to note that the stripes visible at t = 600 are still roughly aligned with the original "front" of predators separating the triangular region from the rest of the domain.The conclusions concerning the evolution of I(S), I(I) and I(Y ) (see Figure 13     and 400, so the initial condition is exactly the same in all cases.We compute the solution at time t = 15 until we obtain a steady state.The results are displayed in Figure 17.It is observed that the resolution in the solution is improved as ∆x is reduced, and that the steady-state solution is apparently the same (modulo, of course, sharpness of resolution) for all values of M .These results do not provide a rigorous convergence proof but indicate that the numerical simulator is reliable.Finally, in Example 12 we utilize the same scenario as in Example 11, and now fix the spatial discretization M = 400 but vary the parameter ε, that is the radius of the convolution kernel of the nonlocal velocity function (see (2.17)).The corresponding results are given in Figures 18 and 19.We observe that the results for the prey compartments S and I are practically the same for all values of ε condered, while the prey distribution depends appreciably on ε.As ε decreases, the regions in which predators concentrate become smaller and the densities become higher.Moreover, Figures 19 (b Section 1.1, the treatment combines the spatio-temporal epidemic model of [45], the non-spatial three-species predator-prey model analyzed in [23], the non-local velocity function introduced in [16] and the numerical method developed in [12].The numerical method is not fully analyzed herein, and certainly leaves potential for further improvements.For instance, its accuracy deserves further investigation particularly to establish that its order of convergence corresponds to its design order.Further, we use a second-order spatial discretization of the Laplacian and a second-order time-stepping, which limit the high-order accuracy of the fifth-order WENO spatial semidiscretization.Nevertheless, on the basis of our (limited) evidence of robustness of the scheme (provided by Examples 11 and 12 along with the fact that numerical solutions of Examples 1, 4 and 7 are consistent with those of [45], where, however, the domain Ω has slightly different dimensions), the numerical results allow us to draw some conclusions and conjectures about solutions of (2.15).Let us first point out that long-time stable patterned configurations obtained in our numerical experiments are formed by stripes.This contrasts with the results in [16] obtained for a predator-prey model with one single prey species only, where the system tends to an asymptotic state with regular arrangements of spots of high predator concentrations, in agreement with the radial structure introduced by the nonlocal velocity function and where the diffusion caused by the Laplacian in the prey equation counterbalances the first order nonlocal differential operator in the predator equation [16].A similar tendency was observed in some scenarios of the eight-compartment epidemiological model studied in [12], where male individuals are supposed to move non-locally in response to the sampled density of their female counterparts.Furthermore, it is calling to attention that in our Example 12, the dominant mechanism of pattern formation seems to be the interaction between susceptible and infected prey, since the solution (see Figure 18) for the prey compartments S and Y is practically independent of ε.That the effect of predator location on the formation of the prey pattern should be weak, at least for the parameters (2.9) and (4.4), can probably be explained by the relative smallness of the predation term on the right-hand side of (2.15b).Finally, it would be interesting to investigate whether the invasion of predators into a domain initially occupied by prey solely, as in Examples 2, 5 and 8, can possibly be modeled by simpler, spatially one-dimensional version of (2.15).Clearly, a more in-depth mathematical analysis to describe the mechanisms that compel these and other phenomena is necessary.With respect to further extensions, we mention that the spatio-temporal approach could also be applied to the model by Greenhalgh et al. [22] that generalizes (2.1) and [23] in the sense that the assumption that the predator eats infected prey only is replaced by the assumption that the susceptible and infected populations are both exposed to the predator but to varying degrees.Another model variant whose predictions should be compared with the findings of the present work arises from replacing the constant recruitment rate A of susceptible prey (see (2.1a) or (2.15a)) by a per capita growth rate, for instance by a logistic growth term as is done in [13-15, 22, 23, 25, 30, 40].Concerning the numerical scheme, further analysis could compare results with higher-order discretizations of the diffusion operator and time-stepping scheme.

2. 4 .
Convective fluxes and diffusion matrix.The flux F Y appearing in the left-hand side of (2.15c) is assumed to have the form F Y (I, Y ) = κ Y Y V (I), where κ Y ≥ 0 is constant and

4 .
10) is solved by Fast Cosine Transforms (due to boundary conditions), which entails a nearly optimal computational cost of O(M 2 log M ).Numerical results.We wish to compare the ODE model (2.1) with the spatiotemporal model (2.15) in different scenarios.

4. 2 .
Solution to the spatio-temporal model.In Examples 1 to 10 we select the spatial domain Ω = [0, L]×[0, L] with L = 100.For the spatial discretization we use M = 400, such that ∆x = L/M = 0.25.For the time discretization we employ Algorithm 3.1.The parameters ε = 4 and κ Y = 1 are chosen in the nonlocal term, and the diffusion constants are chosen as D S = 6 and D I = 1, corresponding to the assumption that infected individuals of the prey population exhibit a lower degree of mobility than their susceptible counterparts.To observe the effect of the non-local movement of predators in the spatio-temporal dynamics of the prey, we consider two cases: Case I, where the values of β and µ in combination with D S and D I have been chosen such that they lie in the Turing region (see[45, Example 6]), and Case II, where this combination is not in the Turing space.The initial condition for S and I is always a spatially distributed random perturbation of the respective values 0.3 and 0.2.For both cases, three scenarios display three different initial conditions of the distributions of the predators: (a) absence of predators, (b) a "triangular" initial condition for predators to describe how predators invade a region that is initially only occupied by the prey, and (c) a random initial distribution of the predators.These three scenarios are illustrated in Figure4.The simulations are carried out until the nonlocal system (2.15) attains a stable non-homogeneous steady state.

Figure 5 . 1 I
Figure 5. Example 1: numerical solution of (2.15) in the absence of predators with parameters (2.9) and (4.1) at three different times.The initial datum given in Figure 4 (a) to the formation of a pattern in the domain.

Figure 7 .Figure 8 .
Figure 7. Example 2: numerical solution of (2.15) with parameters (2.9) and (4.1) at three different times.The initial datum given in Figure 4 (b) leads to the formation of patterns which are different from those obtained in Figure 6.
Figure 8. Example 3: numerical solution of (2.15) with parameters (2.9) and (4.1) at three different times.The initial datum given in Figure 4 (c) leads to the formation of patterns which are different from those obtained in Figures 6 and 7.

Figure 9 .
Figure 9. Example 4: numerical solution of (2.15) with parameters (2.9) and (4.4) at three different times.The initial datum given in Figure 4 (a) leads to a constant equilibrium in the whole domain.

Figure 10 .Example 7 .Figure 11 .
Figure 10.Example 5: numerical solution of (2.15) with parameters (2.9) and (4.4) at three different times.The initial datum given in Figure 4 (b) leads to the formation of patterns that contrast those observed in Figure 9.
(c) and (d)) are similar to those for Example 5.

Figure 12 .Example 9 .
Figure 12.Example 7: numerical solution of (2.15) with parameters (2.9) and (4.7) at three different times.The initial datum given in Figure 4 (a) leads to a constant equilibrium in the whole domain.

4. 3 .
Convergence test and effect of the choice of ε.In Example 11, we investigate the sensitivity of (2.15) to the variation of discretization parameter ∆x = L/M .We take L = 25 and the model parameters (2.9), (4.4) as in Examples 4 to 6.The initial datum is a random distribution (similar to Figure4 (c)) defined for a 50 × 50 discretization, which is also utilized for finer discretizations with M = 100, 200,

Figure 14 . 8 :
Figure 14.8: numerical solution of (2.15) with parameters (2.9) and (4.7) at three different times.The initial datum given in Figure 4 (b) leads to the formation of patterns that contrast those observed in Figure 12.These patterns are different from those shown in Figure 10.

Figure 15 .
Figure 15.Example 9: numerical solution of (2.15) with parameters (2.9) and (4.7) at three different times.The initial datum given in Figure 4 (c) leads to the formation of patterns that contrast those observed in Figure 12 and Figures 10 and 14.

Figure 16 .
Figure 16.Example 10: variation of the parameter c for β = 48, showing the numerical solution at t = 50 for (left) c = 4, (middle) c = 6 and (right) c = 10.All the remaining parameters are as in (2.9).The initial datum is shown in Figure 4 (c).

Figure 17 .
Figure 17.Example 11: variation of the spatial discretization ∆x = L/M of the numerical solution at t = 15 on a domain of side length L = 25.

Figure 18 .
Figure 18.Example 12: effect of variation of convolution radius ε on the numerical solution at t = 15 on a domain of side length L = 25.
1.The spatial discretization of the convolution defining the non-local predator velocities is described in Section 3.2.These velocities, in turn, are utilized in the definition of the final discretization of the convective flux by means of a fifth-order WENO discretization, as is detailed in Section 3.3.The time discretization of the complete model by an IMEX-RK scheme is outlined in Section 3.4.In Section 4, which is at the core of this paper, numerical solutions to the spatio-temporal model are presented.To motivate the choice of parameters, we focus on two equilibrium points of the three-equation ODE model of Section 2.1, as is briefly discussed in Section 4.1.In Section 4.2 we present a total number of 10 examples of numerical solutions of the spatio-temporal model that give rise to patterns.In Section 4.3 we present Example 11 that has been designed to demonstrate in a purely visual way the convergence of the numerical scheme to a definite solution upon refinement of the computational grid, while Example 12 illustrates the influence of the choice of the radius of the convolution kernel of the non-local prey velocity.Finally, some conclusions and possible directions of future work are collected in Section 5.