A reaction-diffusion predator-prey model with pursuit, evasion, and nonlocal sensing

Abstract: In this paper, we propose and analyze a reaction-diffusion model for predator-prey interaction, featuring both prey and predator taxis mediated by nonlocal sensing. Both predator and prey densities are governed by parabolic equations. The prey and predator detect each other indirectly by means of odor or visibility fields, modeled by elliptic equations. We provide uniform estimates in Lebesgue spaces which lead to boundedness and the global well-posedness for the system. Numerical experiments are presented and discussed, allowing us to showcase the dynamical properties of the solutions.


Introduction
The mathematical modeling of predator-prey interactions has a long and rich history.The basic dynamics are given by the system of Lotka-Volterra (here in nondimensional form) where u(t) ≥ 0 represents the predator and w(t) ≥ 0 the prey at time t ≥ 0. According to system (1.1), the predator population decreases exponentially in the absence of prey, while the prey follows a logistic growth law.Interactions between predators and prey are modeled by a mass-action law benefitting the predator and depleting the prey.The system's main feature is the global asymptotic stability of its only nontrivial steady state (u, w) = ( α−1 α , 1 α ), which, for α > 1, expresses a balance between predation, prey reproduction, and predators' natural death rate.
When spatial movement is expected to influence the dynamics, it is natural to consider predator and prey densities u(t, x) and w(t, x) depending on a spatial variable x in some physical domain Ω ⊂ R 2 .Then, one can introduce spatial diffusion and advection to model foraging movement, the spreading of the population in a territory, and/or movement in a preferred direction.
Recently, indirect prey-and predator-taxis have been introduced as a mechanism allowing pursuit and evasion [3,4,10,15].This supposes that the advection velocities are mediated by some indirect signal, which may be an odor, a chemical, a field of visual detection, or seen as a potential.In this spirit, following [3,4,10,15], we consider the predator-prey system with pursuit, evasion and non-local sensing (already written in a non-dimensional form) for t > 0, x in a bounded, open Ω ⊂ R 2 , supplemented with boundary conditions and initial data u 0 (x), w 0 (x).Here, u(t, x) is the predator density, w(t, x) is the prey density, p(t, x) is the odor produced by the prey, q(t, x) is the odor produced by the predator, and α, β, D w , D p , D q , δ p , δ q are non-negative non-dimensional constants -see the Appendix for the physical meaning of the constants and details on the non-dimensionalisation procedure.System (1.2) states that the predator is attracted to the odor p of the prey w, which solves a (steady-state) diffusion equation with source proportional to w, while the prey is repelled by the odor q produced by the predator.Notice that the equations for the odors of the prey and predator are elliptic, rather than parabolic.This is justified in cases where the diffusion of the odor happens in a much faster time scale than the movement of individuals, which is reasonable on a variety of ecological settings.
Note also that we refer to p and q as "odors" but these quantities do not necessarily model chemical odors.They may be more generally interpreted as potentials representing the chance of an animal being detected at a distance by, e.g., visual means.
We remark on the Neumann boundary conditions (1.3).They model the fact that in the time evolution of system (1.2) there is no flow of individuals across the border of the physical domain.Depending on the particular application, this may be a natural assumption, modeling for instance an enclosed area.However, different boundary conditions could be envisaged (e.g., Dirichlet or mixed-type) reflecting distinct natural settings.Here we limit ourselves to the analisys with (1.3).Still, in the numerical experiments presented below, we show an example with Dirichlet boundary conditions.

Main contributions and outline of the paper
In light of the recent mathematical results in [4,15], our main contributions are the introduction of the predator-prey Lotka-Volterra dynamics and the numerical simulations.As we will see, the interaction terms on the right-hand side make the analysis more involved, in particular with respect to the derivation of L ∞ estimates.Indeed, while it is shown in the above-mentioned works that the system with no population or interaction dynamics does not originate finite time blow-up of the solutions, the predator equation dynamics introduces a quadratic term uw so it is not obvious that the boundedness property remains valid.We shall see that the attractive-repulsive nature of the advection terms continues to ensure boundedness of the solutions.Moreover, the property of instantaneous boundedness of the solution even for unbounded initial data, observed already in [4], remains valid in this setting.
An outline of the paper follows.In Section 2, we present our main well-posedness result.Next, in Section 3 we derive our main a priori estimates in L γ spaces, for γ ∈ [1, ∞].This will allow us, in Section 4, to construct strong and then weak solutions to the system (1.2), completing the proof of our well-posedness result.Finally, in Section 5 we detail an implicit-explicit two-step finite volume method for the approximation of system (1.2) and present some numerical experiments.

Main Results
The main results of this paper are concerned with the well-posedness and Lebesgue integrability of weak solutions of the system (1.2).We follow a strategy similar to [4,16], making use of fine a priori estimates in Lebesgue spaces, and a De Giorgi level-set method [17] to obtain boundedness of the solutions.Still, the character of the present system introduces several changes in the analysis with respect to the results in [4,16].
The system (1.2) is, mathematically, of chemotaxis type [18,19].As is well known, such systems may exhibit blow-up of solutions in finite time, see for example the review [20].Therefore, it is not obvious at first glance whether solutions might also possess blow-up behavior.As we shall see in the following analysis, the indirect nature of the sensing, as well as the attraction-repulsion behavior, prevent the densities from becoming infinite in finite time.
We say that the quadruple (u, w, p, q) is a weak solution of the system (1.2) if it satisfies: = Ω w 0 (x)ξ(0, x) dx and Note that while the biologically relevant regime is when α > 1, the value of α > 0 does not change the results, so we present most of our analysis with α, as well as all the remaining constants in system (1.2), set to 1.
We will suppose throughout the paper that the initial data (u 0 , w 0 ) is non-negative and has finite mass and that Ω ⊂ R 2 is a bounded domain of class C 2 .The main results regarding the system (1.2) are collected here: Theorem 2.1.Let the initial data u 0 , w 0 be in L α (Ω) for some α > 2.Then, the system (1.2) has a unique non-negative weak solution.The following estimates are satisfied by the solutions for any In particular, we have the L ∞ -integrability property for some C independent of t > 0.

Analysis of the system (1.2)
In this section we provide a priori estimates which will be used in the well-posedness results.To establish the existence of a weak solution, we must first prove that strong, or classical, solutions exist.We say that (u, w, p, q) is a classical solution of the system (1.2) if it satisfies: 1. (u, w, p, q) ∈ C([0, T ]; L 2 (Ω)) and each of the terms in the system (1.2) are well defined functions in L 2 ((0, T ) × Ω), 2. The equations on (1.2) are satisfied almost everywhere, and 3.The initial data (u, w)| t=0 = (u 0 , w 0 ) and boundary conditions (1.3) are satisfied almost everywhere.
An essential feature of solutions of (1.2) is the following mass estimate: Proposition 3.1 (Mass estimate).Let (u, w, p, q) be sufficiently smooth non-negative solutions of the system (1.2) with the boundary conditions (1.3).Then there exists a constant M depending on u 0 1 , w 0 1 , α, β and |Ω|, but not on t, such that for all t > 0, Proof.Integrating the first and second equations of (1.2) and using the Neumann boundary conditions (1.3), we find From β(w − w 2 ) ≤ (β+1) The conclusion of the proposition readily follows.

A priori estimates in L γ
In this section we prove that data (u 0 , w 0 ) ∈ L 1 (Ω) generate instantaneous L γ -integrability, with γ > 1, for classical non-negative solutions of (1.2).Proposition 3.2 (A priori estimates in L γ ).Let (u, w, p, q) be sufficiently smooth nonnegative solutions of the system (1.2) with boundary condition (1.3) and integrable initial data, and let t > 0 be arbitrary.Then, for any γ ∈ (1, ∞), > 0, we have the estimate Moreover, if u 0 , w 0 ∈ L γ (Ω), then actually for some C > 0 depending on M and L γ −norms of the data, but independent of t.
Proof.We will frequently use the Gagliardo-Nirenberg-Sobolev (GNS) inequality in two dimensions (see e.g.[21]), holding for any α ≥ 1, as well as the interpolation inequality We start by multiplying the second equation of (1.2) by w γ−1 and integrating, to find (discarding a nonpositive term) Now multiply the fourth equation of (1.2) by w γ and integrate by parts to get Also using we find, discarding nonpositive terms, Let us consider the terms on the right-hand side of the previous inequality.Take a small > 0 to be specified later.We use the following consequence of Young's inequality, qw γ ≤ w γ+1 + −γ q γ+1 , and also the inequality which is obtained from Young's inequality, the mass estimate (3.1), and the interpolation inequality (3.5).Therefore, for some constant C depending on γ, M, and , This way, we find and from the GNS inequality (3.4) and sufficiently small, Mathematical Biosciences and Engineering Volume 16, Issue 5, 5114-5145.
Again, from Ω w γ dx ≤ C + Ω w γ+1 dx, we get for some C depending on γ, M, the GNS constant, and the parameters of the system.
To deal with the last term on the right-hand side of (3.8), multiply the fourth equation of (1.2) by q γ−1 to get From the GNS inequality (3.4) we deduce that Choosing small, we find In view of (3.9), the estimate (3.8) becomes Note that the gain from w being the prey population, and thus having a repulsive behavior, is reflected in the lower power u γ .In contrast, performing very similar computations using the first and third equations of the system (1.2), so that instead of (3.6)only is valid, we find that for α ≥ 2 Adding (3.10) and (3.11) gives It is clear that to conveniently bound the terms on the right-hand side using the left-hand side, we should take α < γ < α + 1.Now, we make use of the interpolation inequalities Recalling the mass estimate (3.1), we get for C depending on γ, α, and M. Now, θ 2 (α + 1) < γ + 1 and θ 1 γ < α + 1, so using Young's inequality with a sufficiently small allows the terms on the right-hand side to be absorbed into the left-hand side.This gives for C depending on γ, α, and M. Next, use the inequality (3.5) to find γ−1 and convexity of the power function, we find, setting In view of the definition of η, one finds for C depending on γ, α, and M.This proves the estimate (3.2).The uniform estimate (3.3) follows from the afore-mentioned ODE comparison results in [16].This concludes the proof of Proposition 3.2.

L ∞ estimates
In this section we prove two boundedness results adopting De Giorgi's energy method (see [17] and [4,16,[22][23][24] for related applications of the method).Throughout this section, we will use the notation γ + to denote an arbitrary fixed number in (γ, +∞).
First, we consider initial data u 0 and w 0 only in L 1 (Ω), and obtain an estimate of the type for some constant C > 0 depending on M, but not on T > 0.Then, we will suppose that the initial data u 0 and w 0 are in L ∞ (Ω).Then we can upgrade the estimates above to where C now depends also on the L ∞ norms of initial data.
To further clarify the notation, let us spell out that, for instance, the estimate (3.12) precisely means that for any > 0, it holds with a constant C that may, however, depend on .

Initial data in L 1
The main result in this section is the following: Proposition 3.3.Let (u, w, p, q) be a sufficiently smooth non-negative solution of the system (1.2) with boundary conditions (1.3) and integrable, nonnegative initial data, and let T > 0 be arbitrary.Then, for all t ∈ (0, T ], we have the estimate where the constant C is independent of T > 0. Let us begin by recording here the following L ∞ estimates, which are a consequence of elliptic regularity for the last two equations of the system (1.2) and Proposition 3.2. and The first step in the proof of Proposition 3.3 is a boundedness result valid on each interval (t * , T ) with t * > 0.
Proof.Although the structure of the proof is the similar to corresponding results in [4,16], we show the details since the rather involved calculations depend heavily on the structure of the system.Consider (u, p) a non-negative, sufficiently smooth solution to the general problem where n is the outward unit normal vector to ∂Ω, f and w are known, and w satisfies Similarly, take a non-negative solution (w, q) to the problem g and u given and u satisfying We denote Multiplying the first equation of (3.15) by u λ , integrating in space and using (3.16), we obtain For the first term on the right-hand side note that using Young's inequality, we find Let M > 0 be a constant to be defined later, let t * > 0, and set where we use the notation u k = u λ k .Similarly, we define for w the energy functional with the same definitions of λ k (with N instead M) and t k , for some N > 0 to be chosen later.
With λ = λ k on (3.19), integrating over [s, t], we obtain We use this relation with Integrating with respect to s over [t k−1 , t k ], bearing in mind that t k − t k−1 = t * /2 k and only the first term on the right-hand side of this inequality depends on s, we obtain Now, we are going to introduce a series of estimates for I 1 , I 2 and I 3 .For this, we will use the Gagliardo-Nirenberg interpolation inequality in R n (see e.g.[25]) and a key estimate for 1 V k .Note that we are temporarily performing the analysis in n dimensions.We have which holds for any α ∈ [0, 1] and 1 ≤ p ≤ ∞.Choosing the parameter α ∈ [0, 1] such that αp = 2, it follows that p = 2 n+2 n and holds for all a ≥ 0.
• Estimate for I 1 First, note that for a ≥ 0 we have We choose a in (3.24) such that 2 + a = p = 2 n+2 n , so a = 4 n .Thus, 2 ) ds.

Note that
• Estimate for I 2 We have that Now, choosing a = p in (3.24), we obtain where we need q > 1.Thus, where we used (3.23) with the relation (3.22).We choose α ∈ (0, 1) satisfying αp q = 2 to find where q is such that 1 < q < n n−2 .Thus, where in this last step we proceeded as in the estimate for I 1 .Therefore • Estimates for I 3 Proceeding as we did for the estimate of I 2 , using that u k ≤ u1 V k , we have where α, p and q are the same of the estimate for I 2 .Since sup (3.28) Now we are going to restrict our reasoning to the case n = 2, where we have α = q q+1 and the advantage that we can take 1 < q < ∞, since 1 < q < n n−2 = ∞.Using Proposition 3.2, we have Likewise, and particularizing to f + = uw, Therefore, Now, substituting these results in (3.28), we obtain (3.31) We are going to prove that there exists a constant a ∈ (0, 1) depending only on q such that U k ≤ a k U 0 for all k ∈ N. First, set V k = a k U 0 .Then applying the recurrence relation defined by the right-hand side of (3.31) to V k gives (3.32) Now we choose a such that max{2 3 a, 2 2(q+1) q a 1 q } < 1.So, the last line of (3.32) is bounded by In other words, V k is a supersolution of the recurrence relation defined by (3.31).By a comparison principle, we have U k ≤ a k U 0 −→ k→+∞ 0. Thus, we obtain Mathematical Biosciences and Engineering Volume 16, Issue 5, 5114-5145.
Proof of Proposition 3.3.In Lemma 3.4 we found M, N > 0 such that 0 ≤ u(t, x) ≤ M and 0 ≤ w(t, x) ≤ N almost everywhere on (t * , T )×Ω.From (3.33) and (3.34) we see that M and N depend directly on U 0 and W 0 , respectively.First, we are going to obtain appropriate estimates for these terms, from which the L ∞ −bound for w and u in the time interval (t * , 1) follows.After this, we will extend the estimate for general large intervals.Proceeding as we did in (3.19) for the particular case f + = uw, we have where, in the last step, we used which follows from the first equation of (3.15).Observe that Substituting this in (3.35), we obtain Integrating over [s, t], we find Observe that from this inequality we conclude . Now, via (3.33), using the estimates made previously, we get for where we remember that 0 < t * < 1.For W 0 , a similar computation with (3.17) and particularizing g = w leads to where we recall that 0 < t * < 1.These estimates are valid whenever 0 < t * ≤ t ≤ T = 1, but the same reasoning can be applied for any T > 0 and it would provide a bound depending on T .In order to justify that the same L ∞ -estimate can be made uniform with respect to T , we proceeding by extending this estimate in a similar way as was done in [16]: let t 1 ∈ (t * , T − t * ) and note that the shifted functions w t 1 (t, x) = w(t + t 1 , x) and u t 1 (t, x) = u(t + t 1 , x) are still solutions of the same problem with initial data w t 1 (0, x) = w(t 1 , x) and u t 1 (0, x) = u(t 1 , x), and the appropriate right-hand side.Since the constant C doesn't change due to the Proposition 3.1, we pick t 1 ∈ (0, T ) and repeat the same arguments to w t 1 and u t 1 , which leads the same L ∞ -bounds (3.37)We suppose now we have initial data u 0 and w 0 in L α (Ω), for some α > 2. We will slightly modify the previous analysis in order to obtain better versions of the estimates in Propositions 3.2 and 3.3, as well as a uniform estimate in L ∞ in the case of bounded initial data.Proposition 3.5.Let u 0 and w 0 be initial data in L α (Ω), for some α ∈ (2, ∞].The estimates (3.2) of Propositions 3.2 and 3.3 can be upgraded by adding to the constant C the dependence of L α -norms on the initial data: for any γ > α, we have and, in particular,

.40)
If u 0 , w 0 ∈ L ∞ (Ω), then we have the uniform estimate where the constant C > 0 is independent of T > 0.
Proof.First of all, by Proposition 3.2, there exist a constant A > 0 such that We change the definition t k = (1 − 1/2 k )t * , and observe now that t 0 = 0. Note first that (3.30) becomes In this way, (3.31) becomes with a constant C depending on M, A and T , but independent of t * and k.We are going to proceed as we did after (3.31), but now with (3.42), relying on the fact that here we have where the constant C depends on M, A and T .Let us now see that there exist a ∈ (0, 1), M > 0, so that U k ≤ a k U 0 , for all k.So, taking the right-hand side of (3.42) applied to V k := a k U 0 , we find Then, taking a so that 2 3 a < 1, Choosing M > 0 so that we get that V k is a supersolution of the recurrence defined by (3.42), and so Thus, we have 0 ≤ u(t, x) ≤ M, whence it follows that We get the same estimate for w proceeding exactly in the same way for W k , which leads to

Construction of classical and weak solutions
The a priori estimates of the previous sections will now allow us to prove the global well-posedness of the system (1.2).The first step is to prove existence and uniqueness of classical solutions with smooth initial data.For this, we will use the Banach fixed-point theorem.A stability result for such solutions will be obtained and the existence of a weak solution follows as a consequence.

Construction of classical solutions
Let w 0 ∈ C ∞ c (Ω) and define the set Note that Υ with the metric induced by • Υ is a complete metric space.In this section we prove the following theorem: Theorem 4.1.Let u 0 and w 0 ∈ C ∞ c (Ω) be non-negative initial dada.Then, for all T > 0 the system (1.2) supplemented with the boundary condition (1.3) admits a unique non-negative classical solution.This solution satisfies the estimates obtained in Propositions 3.2, 3.3 and 3.5.
The main step to prove this theorem is the use of the following lemma, whose proof is standard and can be found in [16,Thm. 3.1].: Lemma 4.2.Let ψ be a smooth solution of Now let φ ∈ Υ and let p = p[φ] be the solution of Linear theory guarantees that there exists a unique solution p ∈ H 1 (Ω), and, since φ(t) ∈ L ∞ (Ω) and Ω is smooth, we have p(t) ∈ H 2 (Ω) with p(t) H 2 ≤ C φ(t) L 2 almost everywhere in time.Therefore we have p(t), ∇p(t) and ∆p(t) ∈ L ∞ (Ω).Now we associate u = u[φ] the solution of By linear theory, u is unique and such that u ∈ L 2 (0, T ; H 1 (Ω)) ∩ C([0, T ], L 2 (Ω)).Linear theory still guarantees that there exists some constant C > 0 depending on T and Ω such that u L ∞ (0,T ;H 1 (Ω)) ≤ C u 0 H 1 (Ω) .Now, we associate q = q[φ] the solution of The same arguments lead to the existence of a unique solution q ∈ H 1 (Ω), which satisfies where w[φ] is the only weak solution and such that w ∈ L 2 (0, T ; Proof.Applying Lemma 4.2 for w[φ], Note that all terms in the exponential can be bounded by C φ ∞ , with C > 0 depending on the data of the problem.Thus, for T > 0 small enough we have 0 ≤ w(t, x) ≤ 2 w 0 ∞ , which means w[φ] ∈ Υ. Proof.Let φ 1 , φ 2 ∈ Υ and define w = w 1 − w 2 , where w 1 and w 2 are the respective associated solutions, and so forth.It's easy to check that Multiplying by w and integrating in space, we obtain 1 2 where we used the estimates for w i , u i and φ i guaranteed by Lemma 4.2 and by the definition of Υ, respectively.Substituting and By Gronwall's lemma, it follows that for some constant K > 0 which may change from line to line.Likewise, for u = u 1 − u 2 , we get Note that Combining (4.5) and (4.6) with (4.4), we conclude that for all t ∈ [0, T ].Thus, for T > 0 small enough, Φ is a contraction.This is enough to prove Theorem 4.1.We take T > 0 being the smallest guaranteed by Lemma 4.3 and Lemma 4.4.By the fixed point theorem, the result follows for small time.Extension to [0, T ] is done in a standard way.
Let α > 2. We say that the system of equations (1.2) is stable on L ∞ (0, T ; L α (Ω)) when, given two pairs of initial data u 0,i , w 0,i ∈ L α (Ω), the respective classical solutions u i and w i admit C > 0 depending only on M and on L α −norms of the data such that Proposition 4.5.Let α > 2. The system (1.2) is stable (in the sense of (4.7)) on L ∞ (0, T ; L α (Ω)).
Proof.Let u = u 1 − u 2 and similarly for the other differences.The system for the differences reads Multiplying (4.8) by u and integrating in space, we get 1 2 It is easy to check using (3.40) that Hölder inequality guarantees that Observe that using 1 and (3.40), we find Then, Choosing q big enough, we have from (4.9) and the previous estimates that for some small .A similar result is obtained for w using the same arguments, where we get Defining Z(t) = u 2 2 + w 2 2 , we have Integration of this differential equation leads to the result.

Global well-posedness of weak solutions
Now we are ready to state and prove a well-posedness result for weak solutions: Theorem 4.6.Fix an arbitrary T > 0 and assume non-negative initial data u 0 , w 0 ∈ L α (Ω), for some α > 2.Then, there exists a unique non-negative weak solution for the system (1.2).This solution satisfies the estimates of Proposition 3.2 and 3.3-3.5.
Proof.Take a sequence of non-negative initial data u 0,k , w 0,k ∈ C ∞ c (Ω) with u 0,k → u 0 and w 0,k → w 0 both strongly in L α (Ω).Theorem 4.1 guarantees sequences u k , w k , p k and q k in C([0, T ], L 2 (Ω)) strong solutions of the system (1.2) for each pair of data u The following convergence properties hold: ii) p k → p and q k → q strongly in L ∞ (0, T ; H 1 (Ω)).
iii) u k u and w k w weakly in L 2 (0, T ; In fact, estimate (3.3) guarantees that u k and w k are uniformly bounded in L ∞ (0, T ; L α (Ω)) with respect to k ≥ 1.Therefore, Proposition 4.5 applied to the Cauchy differences u k − u l , w k − w l , k, l ∈ N guarantees that u k and w k are Cauchy sequences in L ∞ (0, T ; L 2 (Ω)), giving i).Now, let p = p k − p l with k, l ∈ N, and denote w = w k − w l .Multiplying by p the difference of the third equations of (1.2) and integrating in space, we get Similarly, with the fourth equation we compute Thus, using Proposition 4.5, we have p(t) H 1 + q(t) H 1 ≤ C u 0,k (t) − u 0,l (t) 2 + w 0,k (t) − w 0,l (t) 2 e C(M)T .(4.11)This means that p k and q k are Cauchy sequences in L ∞ (0, T ; H 1 (Ω)), so we get ii).Next, multiply the first and second lines of (1.2) by u k and w k , respectively, and integrate in space.Using (3.13), (3.14), we find 1 2 Integrating in time gives that ∇u k , ∇w k are in L 2 (0, T ; L 2 (Ω)) uniformly in k.From this we find iii).
For instance, if ϕ ∈ L 2 (0, T ; H 1 ), then we find (to take just the term involving ∇u k , and using that ∇u k is in Treating the other terms similarly, we get that ∂ t u k and ∂ t w k are bounded in L 2 (0, T ; [H 1 (Ω)] * ), so iv) holds.The last convergence follows from the previous ones.Thus, we can pass to the limit k → ∞ in (4.10) to find that (u, w, p, q) is a weak solution of (1.2).The condition (u(0), w(0)) = (u 0 , w 0 ) is satisfied by continuity at t = 0 which follows from the estimate on the time derivatives.Finally, using the approximating by classical solutions we can check that the stability result from Proposition 4.5 and estimate (4.11) hold for weak solutions.This can be used to prove uniqueness: using the notations in Proposition 4.5, if (u i , w i , p i , q i ), i = 1, 2 are two weak solutions associated to the same initial data, then (4.7) proves that (u 1 , w 1 ) = (u 2 , w 2 ) a.e. on [0, T ] × Ω.In turn, the estimate (4.11) (which, since it is a direct consequence of Proposition 4.5, also holds for weak solutions) gives uniqueness of p and q.

Numerical experiments
In order to show some of the relevant features of the system, we provide in this section the details of a numerical simulation of (1.2).Our goal is to present an implicit-explicit finite volume scheme and showcase some numerical results exhibiting the system's main features, namely, evasive behavior of the prey, and chasing by the predator.We consider the system (1.2) in a rectangular domain Ω = [0, L x ] × [0, L y ] , where we introduce a cartesian mesh consisting of the cells I i, j := [x i− 1 2 , x i+ 1  2 ] × [y j− 1 2 , y j+ 1  2 ], which for the sake of simplicity, are assumed square with uniform size, so |I i, j | := h 2 for all i and j.Consider a step size ∆t > 0 to discretize the time interval (0, T ).Let N > 0 the smallest integer such that N∆t ≤ T and set t n := n∆t for n ∈ {0, N}.The cell average of a quantity v at time t is defined by and define v n i, j := v i, j (t n ).Note that in this section we use x = (x, y) to denote the spatial variable.Let f k (u, w), k = 1, 2, be the reactive terms in the right-hand side of the first two equations in (1.2).Then, the terms 1 h 2 are approximated by f k,i, j := f k (u i, j , w i, j ), k = 1, 2. The Laplacian on a Cartesian grid is discretized via The numerical fluxes in the x-and y-directions are respectively and and in a similar way for F 2 i+ 1 2 , j (q) and F 2 i, j+ 1 2 (q).Finally we incorporate a first-order Euler time integration for the u and w components.The diffusive terms are treated in an implicit form and an explicit form is used for the convective and reactive terms.The initial data are approximated by their cell averages, To advance the numerical solution from t n to t n+1 = t n + ∆t, we use the following finite volume scheme: given u n = (u n i, j ) and w n = (w n i, j ) for all cells I i, j at time t = t n , the unknown values u n+1 and w n+1 are determined by the following two steps implicit-explicit scheme: Step 1 solve for p = (p i, j ) and q = (q i, j ) (5.3b) Step 2 solve for u n+1 = (u n+1 i, j ) and w n+1 = (w n+1 i, j ) u n+1 i, j − ∆t∆ i, j u n+1 = u n i, j + ∆t f n 1,i, j (5.4a) where we have used the notation ∆ h = (∆ i, j ) to indicate the matrix of the discrete Laplacian operator, and I is the identity matrix.
Theorem 5.1.Suposse that f 1 (u, w) = u(αw − 1) and f 2 (u, w) = βw(1 − u − w) Then the solutions p i, j , q i, j and u n+1 i, j , w n+1 i, j of the finite volume scheme (5.3a)-(5.3b)and (5.4a)-(5.4b)respectively, are nonnegatives for all i, j provided u n i, j , w n i, j are nonnegative for all i, j, and the following CFL-like condition is satisfied: where a = max i, j Proof.In (5.3a)-(5.3b),we have a linear system of algebraic equations for p i, j and q i, j which need to be solve in each time step t n .However, observe that the matrix of these linear systems are diagonally dominant, which guarantee the existence of solution and the positivity of p i, j and q i, j .Each system of equations (5.4a)-(5.4b)can be seen as a linear system for u n+1 i, j and w n+1 i, j respectively, where the right side is positive according with the CFL condition (5.5) (see [26][27][28]) which guarantees the positivity of u n+1 i, j and w n+1 i, j .
Each system of linear algebraic equations for p i, j , q i, j and u n+1 i, j , w n+1 i, j can be solved by using an accurate and efficient linear algebraic solver.

Test 1: chasing and evasion
In this numerical test, shown in Figure 1, we suppress the terms on the right-hand side of (1.2), to ignore the population dynamics and emphasize the effect of the pursuit and evasion.We can see that the predator starts to chase the prey even though at first any direct contact with it would be very small (due to diffusion only).The prey takes evasive action immediately.Note that by choosing large δ u , δ w in this example, we see from Tables 1 and 2 that this may be interpreted as saying that the chemical sensitivity of the predator and prey are large compared to their diffusion rates.Therefore, we expect that the movement observed in Figure 1 is due to the attraction and repulsion terms and not so much to the diffusion.The numerical parameters are a 400 by 400 spatial cell grid, and a time step of 0.01.Mathematical Biosciences and Engineering Volume 16, Issue 5, 5114-5145.

Test 2: full dynamics
In this test, shown in Figure 2, we set some generic parameters in the system (1.2) in order to observe the full behavior.We can see now the predator-prey interaction taking place, as the densities of the two species fluctuate more widely in relation to the previous example, due to the predator's population growth from predation.After some time, the solution seems to exhibit wave-like interaction patterns with decreasing amplitudes, stabilizing around the values predicted by the equilibrium point ( α−1 α , 1 α ) = (0.9, 0.1).Although here we show the solution only until t = 10, computation up to larger times, not shown here, confirm this behavior.The numerical parameters are a 400 by 400 spatial cell grid, and a time step of 0.01.

Test 3: Dirichlet boundary conditions
This test, which we show in Figure 3, implements a variant of the proposed method for system (1.2) with Dirichlet boundary conditions instead of the conditions (1.3).The parameters are the same as in Test 1.We can see that although the general dynamics remains similar, the boundary behavior affects the solution, in particular making the maximum density smaller.This in natural, since a Dirichlet condition corresponds, physically, to an absorption, or death, at the boundary.
One can also see, especially on the prey column, the formation of steep gradients along the border (boundary layers).This, again, is natural, since with the parameters of the simulation, the evasion dynamics of the prey should be dominant.Since evasion corresponds in (1.2) to an advection term, the formation of boundary layers when considering Dirichlet conditions is to be expected.

Conclusions and outlook
In this work, we have studied analytically and numerically a system of parabolic-elliptic PDEs modeling predator-prey dynamics including prey-and predator-taxis, indirect detection by means of a "potential" or odor, diffusion, and Lotka-Volterra dynamics.We have seen that the system is well posed and established boundedness properties of the solution in Lebesgue spaces.Although these purely mathematical properties may not be of immediate biological relevancy, the reasonings and results establish that the model does not lead to non-biological phenomena such as blow-up in finite timethis being by no means a given property of systems of parabolic equations.
The Lotka-Volterra dynamics with logistic term featured in the right-hand side of our system was chosen since it is the simplest model not leading to unrealistic population explosions.Therefore, possible extensions of our work include the question of whether the results still hold for more realistic population dynamics, namely Holling-type, or after the inclusion of Allee effects, prey handling time, etc.
One other obvious extension to the system (1.2) would be to consider a fully parabolic system where the odor diffusion time scale is of the same order than the population diffusion time scale.In that case, the analysis is more involved.We plan to address this question in future work.
Table 1.Physical parameters in system (6.1).Here, denotes length, t denotes time, bio denotes some measure of predator or prey quantity or mass, and odor denotes some measure of "odor".Prey odor diffusion rate relative to predator diffusion rate D q = α q /α u Predator odor diffusion rate relative to predator diffusion rate δ w = β u K w δ w /(α u β)

Parameter
Normalized prey odor production rate δ u = β w γδ u /(α u βδ) Normalized prey odor production rate δ p = δ p /β Prey odor degradation rate relative to predator death rate δ q = δ q /β Predator odor degradation rate relative to predator death rate

Figure 3 .
Figure 3. Numerical solution of system (1.2) with only pursuit and evasion and Dirichlet boundary conditions.Left column is the predator and right column is the prey.Shown times are, from top to bottom, t = 0, 1, 4, 8. Parameters in this simulation are the same as in Test 1.

Table 2 .
Dimensionless parameters in system (1.2).Here, denotes length, t denotes time, bio denotes some measure of predator or prey quantity or mass, and odor denotes some measure of "odor".
w = α w /α u Prey diffusion rate relative to predator diffusion rate α = αK w /β Predator efficiency relative to death rate β = γ/β Prey growth rate relative to predator death rate D p = α p /α u