Numerical solution of a spatio-temporal predator-prey model with infected prey

1.   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 predator-prey 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 ^[22] (formulated in a non-spatial setting in that paper). The movement of predators is based on the non-local 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 ^[11]. 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 $\boldsymbol{x} \in \Omega$ and time $t\in \mathcal{T}: = [0, T]$ on a bounded domain $\Omega \subset \mathbb{R}^2$. The final model for

is given by a convection-diffusion reaction system of the type

supplied with initial and boundary conditions, where the convective fluxes $\boldsymbol{F}^{\mathrm{c}}(\boldsymbol{u})$, the diffusion matrix $\boldsymbol{\mathcal{D}}$, and the vector of reaction terms $\boldsymbol{s}(\boldsymbol{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 $\boldsymbol{\mathcal{D}} \Delta \boldsymbol{u}$. As in ^[11] 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 well-identified sub-populations or "patches" (cf., e.g., ^{[1,2,3,12,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 ^[22] 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 ^[22] 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 ^[23]: "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 $\nabla \cdot \boldsymbol{F}^{\mathrm{c}} (\boldsymbol{u})$. Related advection terms, for which the essential functional dependence for one compartment is $\boldsymbol{F}^{\mathrm{c}} = \boldsymbol{F}^{\mathrm{c}} (\boldsymbol{x}, u) = \boldsymbol{b}(\boldsymbol{x}) u$ with a given velocity $\boldsymbol{b}(\boldsymbol{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 $\boldsymbol{D} \Delta \boldsymbol{u}$ combined with an explicit one for the convective and reactive terms, $\boldsymbol{F}^{\mathrm{c}}(\boldsymbol{u})$ and $\boldsymbol{s} (\boldsymbol{u})$, respectively. To introduce the main idea, we assume that

represents a method-of-lines semi-discretization of (1.1), where $\Phi^*(\boldsymbol{v})$ and $\Phi(\boldsymbol{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 $y$-direction. Then the stability restriction on the time step $\Delta t$ that explicit schemes impose when applied to (1.2) is very severe ($\Delta t$ must be proportional to $h^2$), due to the presence of $\Phi(\boldsymbol{v})$. The implicit treatment of both $\Phi^*(\boldsymbol{v})$ and $\Phi(\boldsymbol{v})$ would remove any stability restriction on $\Delta t$. However, the upwind nonlinear discretization of the convective terms in $\Phi^*(\boldsymbol{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 $\Phi(\boldsymbol{v})$ and explicitly with $\Phi^*(\boldsymbol{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,7,8,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.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.

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 ^[22]), 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, $\beta$ is the force of infection or the rate of transmission, $\mu$ 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, $\delta$ 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 $\beta 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 righthand 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 [22,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 ^[23]. 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 susceptible 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,23,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 ^[22] as follows: if $I(0) = 0$ and $Y(0) \geq 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(\boldsymbol{x}, 0) > 0$ for $\boldsymbol{x} \in \Omega$.

Theorem 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 $R: = A^2\beta^2-4d^3\beta-8d^2\beta\mu-4d\beta\mu^2$. These equilibria are feasible if $R>0$, which occurs if $\beta> 4d(d+\mu)^{2}/A^{2}$.

$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 \boldsymbol{J}^* (E_4)>0$, and since the entries of $\boldsymbol{J}^* (E_4)$ are real, this means that at least one eigenvalue of $\boldsymbol{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 $\boldsymbol{J}^*$ as

where $e_1 = - {\rm{tr}}(\boldsymbol{J}^*)$, $e_2 = \sum^3_{i = 1} M_i$, where $M_i$ is the determinant of the matrix obtained by eliminating row $i$ and column $i$ from $\boldsymbol{J}^*$, and $e_3 = -\det \boldsymbol{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

This criterion is now applied to $\boldsymbol{J}^* (E_5)$. From (2.5) for $k = 5$ we get

so we always have $e_3>0$ for $\boldsymbol{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$, $\mu$, and $d$ are those proposed in ^[45] while the choices of $m$, $h$, and $\delta$ are assumed but are comparable with those utilized in ^[22]. 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 $\beta$, and the search rate $c$ of the predators towards infected prey. The stability of $E_5$ will be analyzed under variations of $\beta$ and $c$ only.

To determine suitable choices of $\beta$ and $c$, we assume that $e_1 = e_1(\beta, c)$, $e_2 = e_2(\beta, c)$ and $e_3 = e_3(\beta, 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 evaluation to deduce that $e_1$, $e_2$, and a coefficient function $\gamma$ that arises further below, are positive on a relevant $(\beta, 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 $(\beta, 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 [23,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 $(\beta, c) \in \mathcal{R}: = [30, 50] \times [0, 20]$, we also get $e_2 >0$ wherever $E_4$ and $E_5$ are feasible. (We limit the further discussion to $(\beta, c) \in \mathcal{R}$, and understand that $R_1$ to $R_4$ are subregions of $\mathcal{R}$.)

Let us show that across the curve $c \mapsto \beta_{\mathrm{H}} (c)$, which is marked in blue in Figure 1, defined implicitly by $e_1 (\beta _{\mathrm{H}} (c), c) e_2 (\beta_{\mathrm{H}} (c), c) - e_3 (\beta_{\mathrm{H}} (c), c) = 0$, a Hopf bifurcation occurs, i.e., we fix $c$ as a parameter and then consider $\beta$ as a bifurcation parameter. Clearly, the characteristic equation (2.6) takes the form $(\lambda^2 +e_2)(\lambda + e_1) = 0$ with the three roots $\smash{\lambda_1 = \mathrm{i} e_2^{1/2}}$, $\smash{\lambda_2 = -\mathrm{i} e_2^{1/2}}$ and $\smash{\lambda_3 = -e_1}$. Thus, in a neighborhood of $(\beta_{\mathrm{H}} (c), c)$, (2.6) cannot have real positive roots; rather, the roots are of the form $\lambda_1 = u + \mathrm{i} v$, $\lambda_2 = u - \mathrm{i} v$ and $\lambda_3 = - \phi$ with $\phi>0$, and where $u, v \in \mathbb{R}$ and $\phi$ depend on $(\beta_{\mathrm{H}}(c), c)$. It remains to verify the transversality condition, i.e.,

To this end, we substitute $\lambda = \lambda_1 = u + {\rm{i}}v$ into (2.6), equate real and imaginary parts, and calculate the derivative (see ^[23]) to obtain the following equations, where $u = u(\beta, c)$, $u' = \partial u(\beta, c) / \partial \beta$, and analogously for $v$, $e_1$, $e_2$, and $e_3$ and quantities $\rho = \rho (\beta, c)$, $\sigma$, $\tau$ and $\nu$ defined below. This yields

where

At $\beta = \beta_{\mathrm{H}} (c)$, $u(\beta_{\mathrm{H}} (c)) = 0$ and $v(\beta_{\mathrm{H}} (c)) \neq 0$. Since $e_1 >0$ in a neighborhood of $(\beta_{\mathrm{H}}(c), c)$, it follows that $\sigma \neq 0$. Applying Cramer's rule to (2.10) we get

We wish to show that $\sigma \nu + \rho \tau \neq 0$. Evaluating (2.11) at $\beta = \beta_{\mathrm{H}}(c)$, we obtain at that state $\sigma \nu + \rho \tau = 2 e_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' = \gamma (\beta_{\mathrm{H}} (c), c) \partial (\beta I_5^2)/ \partial \beta$, where we define

For the parameters (2.9), we obtain

as we infer from the rightmost plot of Figure 2. (Indeed, $\gamma$ is positive on a slightly larger subset of $\mathcal{R}$. However, positivity on $R_3 \cup 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

which proves the transversality condition. Thus, a Hopf bifurcation occurs across $c \mapsto \beta_{\mathrm{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 $\boldsymbol{J}^* (E_5)$ zero, we must have $\det \boldsymbol{J}^* (E_5) = 0$, but here we know that $\det \boldsymbol{J}^* (E_5) < 0$ wherever $E_5$ is feasible. In particular, it cannot happen that $\boldsymbol{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 \boldsymbol{J}^* (E_5) < 0$, which has been established symbolically and is therefore valid for any choice of parameters (instead of (2.9)) and $(\beta, c)$ for which $E_5$ is feasible.

2.3.   Spatio-temporal predator-prey model

The ODE model (2.1) exhibits a rich solution behaviour, which motivates us to analyze a spatial variant that in turn generalizes the spatio-temporal hyperbolic-parabolic predator-prey model studied by Colombo and Rossi ^[16]. Following ^[11,16], we then replace (2.1c) by a non-local spatio-temporal PDE that contains a term describing the spatial movement of predators towards infected prey, and equip (2.1a) and (2.1b) with diffusion terms. The resulting spatio-temporal model is given by

where $\nabla \cdot$ denotes the (spatial) divergence operator. Here $x$ and $y$ are the space variables, $\Delta = \partial^2/\partial x^2+\partial^2/\partial 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.

2.4.   Convective fluxes and diffusion matrix

The flux $\boldsymbol{F}_Y$ appearing in the left-hand side of (2.15c) is assumed to have the form $\boldsymbol{F}_Y(I, Y) = \kappa_{Y}Y \boldsymbol{V} (I)$, where $\kappa_Y \geq 0$ is constant and

is the non-local unscaled velocity function ^[16]. Here $\eta$ denotes a radial convolution kernel with radius $\varepsilon$, i.e., $\eta$ is a piecewise smooth function such that

i.e., for any function $w$ defined on $\Omega \times \mathcal{T}$ and $\boldsymbol{x} \in \Omega$ such that $B_{\varepsilon} (\boldsymbol{x}) : = \{ \boldsymbol{y} \in \mathbb{R}^2 : \|\boldsymbol{y}- \boldsymbol{x} \| < \varepsilon \} \subset \Omega$, we have

(with slight modifications for points $\boldsymbol{x}$ with ${\rm{dist}} (\boldsymbol{x}, \partial \Omega) < \varepsilon$.) Since $\nabla (w \ast \eta) = w \ast \nabla \eta$, the velocity $\boldsymbol{V}(w)$ indeed depends (non-locally) on $w$ and not on its gradient.

Summarizing, we obtain $\smash{\boldsymbol{F}^{\mathrm{c}}(\boldsymbol{u}) = \bigl(0, 0, F_{Y}(I, Y))^{\mathrm{T}}}$ and $\boldsymbol{\mathcal{D}} = {\rm{diiag}} (D_{S}, D_{I}, 0)$. The vector $\boldsymbol{s}(\boldsymbol{u}) = (s_1(\boldsymbol{u}), s_2(\boldsymbol{u}), s_3(\boldsymbol{u}))^{\mathrm{T}}$ is given by the right-hand sides of (2.15). The system (2.15) is considered on $\Omega \times \mathcal{T}$ along with the initial condition

where $\boldsymbol{u}_0$ is a given function, and zero-flux boundary conditions

where $\boldsymbol{n}$ is the unit exterior normal vector to the boundary $\partial \Omega$ of $\Omega$.

3.   Numerical method

3.1.   Spatial discretization of local convection and diffusion terms

We take $\Omega = [0, L]\times [0, L]$ and use a Cartesian grid with nodes $(x_i, y_j)$, $i, j = 1, \dots, M$, with $x_i = y_i = (i-1/2)h$, $h = L/M$. We discretize $\nabla \cdot \boldsymbol{F}^{\mathrm{c}}(\boldsymbol{u})$ by WENO finite differences and $\Delta \boldsymbol{u}$ by the standard second-order scheme with a five-point stencil to get a spatial semi-discretization of (1.1) for a ${\rm{3}} \times M \times M$-matrix $\smash{\boldsymbol{v} (t)}$ of unknown approximations $\smash{v_{\ell, i, j} (t)\approx u_{\ell}(x_i, y_j, t)}$ for $i, j = 1, \dots, M$ and $\ell = 1, \dots, 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 $\smash{\nabla_h \cdot \boldsymbol{\tilde{F}}^{\mathrm{c}}(\boldsymbol{v})}$ is the discretization of $\nabla \cdot \boldsymbol{F}^{\mathrm{c}}(\boldsymbol{u})$, to be defined in Section 3.3, and

denotes the discretization of the diffusion terms. Here $\boldsymbol{v}_{\ell}$ denotes the $M\times M$ submatrix given by $(\boldsymbol{v}_{\ell})_{i, j} = v_{\ell, i, j}$ and $\boldsymbol{\Delta}_h$ is the approximation of the standard two-dimensional Laplacian operator with Neumann boundary conditions. Furthermore, $\boldsymbol{S}(\boldsymbol{v})$ is the ${\rm{3}} \times M \times M$-matrix with components

with corresponding submatrices $\boldsymbol{S}_{\ell}(\boldsymbol{v})$, given by $\smash{\boldsymbol{S}_{\ell}(\boldsymbol{v})_{i, j} = s_\ell(v_{\ell, i, j})}$.

3.2.   Discretization of the convolutions

The following identity arises from (2.16) if we take into account that $\nabla (w \ast \eta) = w \ast \nabla \eta$:

The corresponding convolutions $w\ast \partial{\eta}/\partial{x}$ and $w \ast \partial{\eta}/ \partial{y}$ are calculated approximately on the discrete grid via a composite Newton-Cotes quadrature formula, such as the composite Simpson rule.

Since $B_{\varepsilon}(0)\subseteq [-rh, rh]^2$, $r = \lceil \varepsilon/h \rceil < M$, and according to boundary conditions, $w$ is extended to the exterior of $\Omega$ by reflection, e.g. by setting $w(-x, y) = w(x, y)$, $(x, y)\in\Omega$, we obtain

where $\alpha_{p}$ and $\alpha_{q}$ are the coefficients in the quadrature rule (e.g., for the composite Simpson rule, $\boldsymbol{\alpha} = (1, 4, 2, 4, \dots, 2, 4, 1)$). Consequently, the approximation (3.4) for $\boldsymbol{W} = (w_{i, j}) \in\mathbb R^{M\times M}$, $w_{i, j}\approx w(x_i, y_j)$, is given by

where we define

The discrete approximation of $\boldsymbol{V}(w)$ in (3.3) is then given by

Since $r\approx \varepsilon /h = \varepsilon M/L$, the computational cost of this discrete convolution is $M^2(2r+1)^2\approx 4 \varepsilon^2M^4/L^2$, which can be very high for large $M$. This cost can be substantially reduced to $\mathcal{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 $\boldsymbol{W} = (w_{i, j}) \in\mathbb R^{M\times M}$ a matrix $\smash{\boldsymbol{\tilde{W}}} = (\tilde{w}_{i, j}) \in\mathbb{R}^{2M\times 2M}$ such that

and use the notation $[i]'_{2M} = \mod (i-1, 2M)+1$, i.e., $[i]'_{2M} = i+2kM$, with $k$ being the integer such that $1\leq [i]'_{2M}\leq 2M$. It is then readily checked that

Therefore (3.5) for $i, j = 1, \dots, M$ can be rewritten as

The convolution on the right-hand side of (3.6) can be performed by FFTs applied to the $(2M)\times (2M)$ matrix $\smash{\boldsymbol{\tilde{W}}}$. To save further computational costs, the FFT of the kernel $\beta_{p, q}$ is performed only once, son each convolution entails two two-dimensional FFT of $(2M)\times (2M)$ matrices and a product of $4M^2$ numbers, with an overall computational cost of $\mathcal{O}(M^2\log M)$.

3.3.   Discretization of the convective term

The convective flux for the prey compartments ($S$ and $I$, corresponding to $\ell\in\{1, 2\}$) is zero and for the predator (species $Y$, $\ell = 3$) is given by $\smash{\boldsymbol{F}^{\mathrm{c}}_{3} (\boldsymbol{u}) = u_{3} (\kappa_{3} \boldsymbol{V} (u_2))}$. To discretize its divergence $\nabla\cdot \boldsymbol{F}^{\mathrm{c}}_{3} (\boldsymbol{u})$ for the approximation $\boldsymbol{v}$, we first approximate the convolution terms as expounded in Section 3.2 to obtain

We introduce the notation ${\smash (f_{i, j}^{x}, f_{i, j}^{y}): = \boldsymbol{\tilde{F}}^{\mathrm{c}}_{\ell} (\boldsymbol{v})_{i, j}}$, where we have dropped the index $\ell$ in $\smash{f_{i, j}^{x}}$ etc. to obtain clearer expressions. Our purpose is to use a fifth-order WENO finite difference discretization ^[26,34,42] of $\nabla\cdot \boldsymbol{F}^{\mathrm{c}}_{\ell} (\boldsymbol{u})$ for which

for suitable numerical fluxes $\smash{\hat f^{x}_{i+1/2, j}, \hat 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, \pm}$ is given by

with an analogous formula for $\smash{f^{y, \pm}_{i, j}}$. If $\mathcal{R}^{\pm}$ 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 $\Phi(\boldsymbol{v})$ is handled by an implicit $s$-stage diagonally implicit (DIRK) scheme with coefficients $\boldsymbol{A}\in\mathbb{R}^{s\times s}$, $\boldsymbol{b}, \boldsymbol{c}\in \mathbb{R}^s$, in the common notation, where $\boldsymbol{A} = (a_{ij})$ with $a_{ij} = 0$ for $j>i$. For the term $\Phi^*(\boldsymbol{v})$ we employ an $s$-stage explicit scheme with coefficients $\smash{\boldsymbol{\hat{A}}\in\mathbb{R}^{s\times s}}$, $\smash{\boldsymbol{\hat{b}}}$, $\boldsymbol{\hat{c}}\in\mathbb{R}^{s}$ and $\smash{\boldsymbol{\hat{A}} = (\hat{a}_{ij})}$ with $\hat{a}_{ij} = 0$ for $j\geq 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 IMEXRK scheme gives rise to the following algorithm (see ^[38]).

Algorithm 3.1.

    Input: approximate solution vector $\boldsymbol{v}^n$ of (1.2) for $t = t^n$

    for $p = 1, \dots, s$

        if $p = 1$ then

            ${\displaystyle \boldsymbol{\hat{v}}^{(1)} \leftarrow \boldsymbol{v}^n, \quad \bar{\boldsymbol{v}}^{(1)} \leftarrow \boldsymbol{v}^n }$

        else compute $\boldsymbol{\hat{v}}^{(p)}$ and $\bar{\boldsymbol{v}}^{(p)}$ from the known values $\boldsymbol{K}_1, \dots, \boldsymbol{K}_{p-1}$:

            ${\displaystyle \boldsymbol{\hat{v}}^{(p)} \leftarrow \boldsymbol{v}^n + \Delta t\sum_{j = 1}^{p-1}\hat{a}_{pj} \boldsymbol{K}_j}, \quad {\displaystyle \bar{\boldsymbol{v}}^{(p)} \leftarrow \boldsymbol{v}^n + \Delta t\sum_{j = 1}^{p-1} a_{pj} \boldsymbol{K}_j}$

        endif

        solve for $\boldsymbol{K}_p$ the linear system

                    ${\boldsymbol{K}_p} = \Phi^*\bigl(\boldsymbol{\hat{v}}^{(p)}\bigr) + \Phi\bigl(\bar{\boldsymbol{v}}^{(p)} + \Delta t a_{pp} \boldsymbol{K}_p \bigr)\quad \quad \quad \quad \quad \quad (3.8)$

    endfor

    ${\displaystyle \boldsymbol{v}^{n+1} \leftarrow \boldsymbol{v}^n+\Delta t \sum_{j = 1}^sb_j\boldsymbol{K}_j}$

    Output: approximate solution vector $\boldsymbol{v}^{n+1}$ of (1.2) for $t = t^{n+1} = t^n+\Delta t$.

To solve the linear equation (3.8) for $\boldsymbol{K}_{p}$, in view of (3.7) we rewrite it as

where $\boldsymbol{I}$ denotes the identity operator for ${\rm{3}} \times M \times M$ matrices. From the definition of the matrix $\boldsymbol{\mathcal{B}}$ in (3.2) and from the definition of $\Phi^*$ in (3.7), if we equate the $\ell$ submatrices along the first dimension of both sides of (3.9) we get

where $\boldsymbol{I}_{M\times M}$ is the $M\times M$ identity operator and, e.g., $(\boldsymbol{K}_{p})_{\ell}$ is the $\ell$ submatrix of $\boldsymbol{K}_p$ along the first dimension, i.e., $((\boldsymbol{K}_{p})_{\ell})_{i, j} = (\boldsymbol{K}_{p})_{\ell, i, j}$. If $\mu_{\ell} = 0$ (for $\ell = 3$), then

otherwise (3.10) is solved by Fast Cosine Transforms (due to boundary conditions), which entails a nearly optimal computational cost of $\mathcal{O}(M^2\log M)$.

4.   Numerical results

We wish to compare the ODE model (2.1) with the spatio-temporal model (2.15) in different scenarios.

4.1.   Solution to the ODE model

For this model we consider parameter values (2.9) along with

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 1 (besides $E_1 = (1, 0, 0)$) where the predator is absent:

as well as the two points

From Section 2.2 we know that all eigenvalues of $\boldsymbol{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.

4.2.   Solution to the spatio-temporal model

In Examples 1 to 10 we select the spatial domain $\Omega = [0, L]\times[0, L]$ with $L = 100$. For the spatial discretization we use $M = 400$, such that $\Delta x = L/M = 0.25$. For the time discretization we employ Algorithm 3.1. The parameters $\varepsilon = 4$ and $\kappa_Y = 1$ are chosen in the non-local 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 Ⅰ, where the values of $\beta$ and $\mu$ 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 Ⅱ, 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 Figure 4. The simulations are carried out until the non-local system (2.15) attains a stable non-homogeneous steady state.

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 \in \{ S, I, Y\}$ and time instants $t^n = n \Delta t$ the following quantity:

which represents the approximate total number in $\Omega$ of individuals of compartment $X$ at time $t^n$. We recall that in the PDE model, the unknowns $X \in \mathcal{C}$ are densities, and so an integration over $\Omega$ 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 Ⅰ 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 $\beta$ and $\mu$ 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 ^[45]. This distribution will remain constant over time. Moreover, Figure 6 shows the integral quantities $\mathcal{I} (S)$, $\mathcal{I} (I)$, and wherever appropriate $\mathcal{I} (Y)$, for Examples 1 to 6. For Example 1, and consistently with the nearly stationary spatial pattern, we observe (in Figures 6(a) and (b)) that $\mathcal{I} (S)$ and $\mathcal{I} (I)$ become practically stationary after very short time, and assume values close to $L^2 \times 0.2483389 \approx 2483$ and $L^2 \times 0.2684504 \approx 2684$ that would be obtained if the solution were equal to the stable equilibrium $E_3$ (see 4.2) on the whole domain. In this sense, in this example the temporal behavior of the model is consistent with (2.1).

Example 2: Case Ⅰ with a "triangular" initial distribution of the predator.

Here we introduce predators through the "triangular" initial condition of Figure 4(b). All other parameters are the same as in Example 1. In Figure 7 we display the dynamics of model (2.15). The convective term in (2.15c) defined in terms of the non-local velocity function (2.16) allows the movement of the predators to places where infective prey are present. At simulated time $t = 30$ we observe that the distribution of prey varies in comparison with Example 1 in places where predators are present. At $t = 90$ we observe that in places where the predators are absent, preys are distributed in forming spots as in Example 1, however, where predators are present, these patterns take the form of stripes. At $t = 200$, predators are present in the whole domain, and this has motivated that the distribution of the prey in the domain varies drastically from its natural distribution without predator (Example 1). This distribution will remain constant over time. The corresponding plots of $\mathcal{I} (S)$ and $\mathcal{I} (I)$ (Figure 6(c)) and $\mathcal{I} (Y)$ indicate a growth of the predator population over a very long period of time. Shortly before $t = 200$ these quantities become constant. It is worth noting that while $\mathcal{I} (S)$ and $\mathcal{I} (I)$ tend to values that are close to $L^2$ times the first and second component of $E_5$, namely 2986.8 and 2364.4, respectively, the limit value of $\mathcal{I} (Y)$ exceeds $L^2$ times the $Y$-component of $E_5$, that is 23644.3, by more than a factor of two.

Example 3: Case Ⅰ with a random initial distribution of the predator.

Here we consider a random initial condition for predators with the same parameters of Example 1. The initial condition is shown in Figure 4(c). In Figure 8 we display the dynamics of the model with this initial condition. At simulated time $t = 30$ we observe that stripes and some spots patterns emerge mixed in the distribution of each species. After the stripe patterns form, they grow steadily with time until they reach a certain arm length, and the spatial patterns become distinct at $t = 200$. It is observed that these stripe patterns are not similar to those obtained in Example 1, instead these patterns in form of stripes or filaments are very similar to those observed in Example 2. Furthermore, Figures 6(e) and (f) indicate that $\mathcal{I} (S)$, $\mathcal{I} (I)$ and $\mathcal{I} (Y)$ quickly attain steady states, with final values similar to those of Example 2.

Example 4: Case Ⅱ in the absence of the predator.

Here, we consider absence of predators like in Example 1 and we use here and in Examples 5 and 6 the parameter values

In this case the parameters $\beta$ and $\mu$ together with the diffusion constants are not in the Turing space. The steady states according to Theorem 1 are now, besides $E_1 = (1, 0, 0)$,

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 Ⅱ 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 appropriate 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 Ⅱ with a random initial distribution of the predator.

This case is the analogue of Example 2, namely we impose the initial condition 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 patterns become slightly distinct at $t = 200$ since at this time there are strips and spots. It is observed that these patterns differ from those obtained in Example 5, where only stripe-like patterns are formed. The observations corresponding to the evolution of $\mathcal{I}(S)$, $\mathcal{I}(I)$ and $\mathcal{I}(Y)$ are analogous to those of Example 3.

Example 7

Here and in Examples 8 and 9 we use the values

The equilibrium points are $E_1 = (1, 0, 0)$, $E_2$ and $E_3$ given by (4.5), and

The integral quantities for this case and Examples 8 to 10 are plotted in Figure 13.

As in Example 4, the parameters $\beta$ and $\mu$ together with the diffusion constants are not in the Turing space. Considering the initial condition of Figure 4(a), in this case no patterns are formed and the system rapidly tends to a constant equilibrium 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 $\mathcal{I}(S)$, $\mathcal{I} (I)$ and $\mathcal{I} (Y)$ (see Figure 13(c) and (d)) are similar to those for Example 5.

We now consider an analogue of Examples 3 and 6 corresponding to the initial configuration of Figure 4(c). In Figure 15 we display the dynamics of the model with this initial condition. At $t = 50$ we observe that some stripes and spots patterns emerge mixed in the distribution of each species. The pattern distribution is maintained over time after its formation. It is observed that these patterns differ from those obtained in Examples 5, 6, and 8. As in Examples 3 and 6, the quantities $\mathcal{I} (S)$, $\mathcal{I} (I)$ and $\mathcal{I} (Y)$ very quickly attain constant values (see Figures 13(e) and (f)).

Example 10: Effect of varying parameter $c$.

In this example, a constant value is considered for $\beta = 48$ and in combination with the three different values $c = 4$, $c = 6$ and $c = 10$ (all pairs $(48, c)$ belong, again, to $R_4$, see Figure 1). For these three pairs a variation is obtained in the formation of patterns for the different values of $c$. The initial datum is given by Figure 4(c). In Figure 16 the corresponding numerical results are shown. It turns out that with increasing $c$ the species tend to form more filament-like rather than spot-like spatial structures. Moreover, there is a marked difference in behavior of $\mathcal{I}(S)$, $\mathcal{I} (I)$ and $\mathcal{I} (Y)$, as can be inferred from Figures 13(g) to (l). Note that the steady-state value of $\mathcal{I} (I)$ decreases with increasing $c$, as does the steady-state value of $\mathcal{I}(Y)$.

4.3.   Convergence test and effect of the choice of $\varepsilon$

In Example 11, we investigate the sensitivity of (2.15) to the variation of discretization parameter $\Delta 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 Figure 4(c)) defined for a $50 \times 50$ discretization, which is also utilized for finer discretizations with $M = 100$, $200$, 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 $\Delta 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 $\varepsilon$, that is the radius of the convolution kernel of the non-local 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 $\varepsilon$ considered, while the prey distribution depends appreciably on $\varepsilon$. As $\varepsilon$ decreases, the regions in which predators concentrate become smaller and the densities become higher. Moreover, Figures 19(b), (d), (f), and (h) indicate a decrease of the steady state value of $\mathcal{I}$ as $\varepsilon$ is decreased.

5.   Conclusions

We have shown that a relatively simple three-species model can produce complex spatio-temporal patterns of predators and prey. As is stated in Section 1.1, the treatment combines the spatio-temporal epidemic model of ^[45], the non-spatial three-species predator-prey model analyzed in ^[22], the nonlocal velocity function introduced in ^[16] and the numerical method developed in ^[11]. 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 $\Omega$ 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 non-local velocity function and where the diffusion caused by the Laplacian in the prey equation counterbalances the first-order non-local differential operator in the predator equation ^[16]. A similar tendency was observed in some scenarios of the eight-compartment epidemiological model studied in ^[11], 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 $\varepsilon$. 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.^[23] that generalizes (2.1) and ^[22] 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,14,15,23,22,25,30,40]}. Concerning the numerical scheme, further analysis could compare results with higher-order discretizations of the diffusion operator and time-stepping scheme.

Acknowledgements

RB is supported by Fondecyt project 1170473; and Centro CRHIAM Proyecto Conicyt/Fondap/15130015. PM is supported by Spanish MINECO project MTM2017-83942-P and Conicyt (Chile), project PAI-MEC, folio 80150006. LMV is supported by Fondecyt project 1181511. RB and LMV are also supported by BASAL pro-ject CMM, Universidad de Chile and Centro de Investigación en Ingeniería Matemá-tica (CI$^{\mathrm{2}}$MA), Universidad de Concepción and the INRIA Associated Team "Efficient numerical schemes for non-local transport phenomena (NOLOCO; 2018--2020)". EG is supported by CONICYT scholarship. GC acknowledges financial support from grants NSF-IIS RAPID award #1518939, NSF grant 1318788 Ⅲ: Small: Data Management for Real-Time Data Driven Epidemic simulation, and Conicyt (Chile), project MEC80170119.

Conflict of interest

All authors declare no conflicts of interest in this paper.