Analysis of a predator-prey model with disease in the predator species.

Piotr Radziński (Warsaw) Urszula Foryś (Warsaw)

Analysis of a predator-prey model with disease in the predator species.

Abstract In the paper we analyse a diffusive predator-prey model with disease in predator species proposed byQiao et al.(2014). In the original article there appears a mistake in the procedure of the model undimensionalisation. We make a correction in this procedure and show that some changes in the model analysis are necessary to obtain results similar to those presented byQiao et al.(2014)

We propose corrected conditions for global stability of one of the existing equilibria – disease free steady state and endemic state in the case without diffusion as well as in the model with diffusion. On the basis of the corrected analysis we present new stability results.

2010 Mathematics Subject Classification: Primary: 92D40; Secondary: 92D30..

Key words and phrases: predator-prey model, eco-epidemiology, local and global stability.

1. Introduction For many years mathematical modelling in ecology and epidemiology appeared as two different branches of applied mathematics. Re- cently, a new field called eco-epidemiology has started to develop and models describing interactions between preys and predators in which at least one of the species is divided into subpopulations of healthy and infected individuals have been proposed and analysed in the literature. In this paper we focus on one of such models with a disease in the predator species; cf. [5]. The model proposed by Qiao et al. [5] is based on the ideas presented in [1], [2] and [3] and includes simple spatial effects reflected by diffusion. During our stud- ies we have found that the model was incorrectly undimensionalised which resulted in sticking two parameters together. Therefore, the whole analysis presented in [5] needs some corrections. In this paper we show that the results obtained by Qiao et al. could be transferred into the corrected version of the undimensional model. However, assumptions need to be adjusted.

2. Model description Qiao et al. [5] considered a model describing interactions between a population of preys and a population of predators in

which some disease spreads. This population is therefore divided into two sub- populations of healthy and infected individuals. The dynamics of preys is de- scribed by a classic Lotka-Volterra model with carrying capacity for preys; c.f.

e.g. [4] for more details on such type of models. Preys are eaten only by healthy predators. The underlying dynamics of healthy predators is governed by May model ([3]) with the standard bilinear infection term. Infected predators ap- pear after infection of healthy individuals by infected ones, while they die due to natural reasons and due to overcrowding. Qiao et al. [5] included also simple spacial effects due to diffusion.

Eventually, the model we are interested in is expressed in the following way:

∂u

∂t = d1∆u + ru

 1 − u

K



− kuv, x ∈ Ω,

∂v

∂t = d2∆v + sv

 1 −hv

u



− βvw, x ∈ Ω,

∂w

∂t = d3∆w + βvw − dw − µw², x ∈ Ω,

∂u

∂η = ∂v

∂η = ∂w

∂η = 0, x ∈ ∂Ω,

u(x, 0) = u0(x) ≥ 0, v(x, 0) = v0(x) ≥ 0, w(x, 0) = w0(x) ≥ 0, x ∈ Ω,

(1)

where u, v and w represent densities of preys, susceptible predators and in- fected predators, respectively. Parameters r and s are responsible for repro- duction of preys and predators, K denotes carrying capacity for preys, capac- ity of predators is dependent on the prey population with parameter h, while k is hunting coefficient. The rate of infection is denoted by β, the linear death rate is reflected by d, while µ is responsible for density-dependent death rate.

Parameters d₁, d₂ and d₃ are diffusion coefficients. Ω is an open domain in Rⁿ with a smooth boundary ∂Ω, and _∂η^∂ is the outward derivative normal to

∂Ω.

To undimensionalise Eqs. (1) we use the following change of variables:

u 7−→ Ku, v 7−→ rK

shv, w 7−→ r

µw, t 7−→ t r, and we rename the parameters:

kK

sh 7−→ k, s

r 7−→ s, β

µ 7−→ β₁, βK

sh 7−→ β₂, d

r 7−→ d, d1

r 7−→ d₁, d2

r 7−→ d₂, d3

r 7−→ d₃,

obtaining the following undimensional model:

∂u

∂t = d₁∆u + u(1 − u) − kuv, x ∈ Ω,

∂v

∂t = d₂∆v + v s − v

u

− β₁vw, x ∈ Ω,

∂w

∂t = d3∆w + β2vw − dw − w², x ∈ Ω,

∂u

∂η = ∂v

∂η = ∂w

∂η = 0, x ∈ ∂Ω,

u(x, 0) = u₀(x) ≥ 0, v(x, 0) = v₀(x) ≥ 0, w(x, 0) = w0(x) ≥ 0, x ∈ Ω.

(2)

Notice that in Qiao et al. [5] there is a mistake in the undimensional form of the model as β₁ = β₂, which is true only under the assumption sh = µK, leading to an ungeneric case.

3. The model without diffusion The first step in the analysis of the reaction-diffusion model is to consider the case without diffusion. Assuming d₁ = d₂ = d₃ = 0 in Eqs. (2) we obtain the model without diffusion that reads:

˙

u = u(1 − u) − kuv,

˙v = v s − v

u

− β₁vw,

˙

w = β₂vw − dw − w².

(3)

Basic properties of Eqs. (3) like global existence of non-negative solutions are easy to check. Moreover, following the ideas of Qiao et al. [5], depending on the model parameters we get sets that are positively invariant and globally attractive. More precisely, let w_M = β2s−d for β2s > d or wM be an arbitrary positive number for β₂s ≤ d. Then

• the set D = [0, 1] × [0, s] × [0, w_M] is positively invariant;

• for ks < 1 and every δ_u< 1 − ks, the set E = [δ_u, 1] × [0, s] × [0, w_M] is positively invariant;

• for ks < 1, β₁β₂ < 1, d < β₂s(1 − β₁β₂)(1 − ks), and every δ_u< 1 − ks, δv < s(1 − β1β2)δu, δ_w < wM, the set G = [δ_u, 1] × [δv, s] × [δw, wM] is positively invariant.

This means that positively invariant and globally attractive set exists for any model parameters. Moreover, for ks < 1, β₁β₂ < 1, d < β₂s(1 − β₁β₂)(1 − ks) solutions of Eqs. (3) are bounded above from 0.

Now, we turn to the steady states analysis. There are three such states:

disease-free state (DFS) with w = 0, positive steady state usually called

endemic (ES) and (1, 0, 0) reflecting the absence of predators. Analysing local stability we look for Jacobi matrix. For Eqs. (3) this matrix reads

J =





1 − 2u − kv −ku 0

v²

u² s − 2_u^v − β₁w −β₁v

0 β₂w β₂v − d − 2w



.

Notice that for (1, 0, 0) this matrix is diagonal and eigenvalues are equal to

−1, s and −d, which implies instability, so we will focus on the first two states.

3.1. Stability of DFS Let us first focus on DFS. This state has coor- dinates (u, v, 0), where

u = 1

1 + ks, v = su.

It is easy to check that 1 − u − kv = 0, and hence the Jacobi matrix takes the form

J_DFS=





−u −ku 0

s² −s −β₁v 0 0 β₂v − d



.

We see that local stability depends on the sign of β₂v − d. It is easy to check that β₂v−d < 0 ⇐⇒ β₂ < dk+^d_s and this is the condition for local asymptotic stability of DFS. To investigate global stability we define a Lyapunov function

L1(u, v, w) = u − u − u lnu

u + v − v − v lnv v +β1

β₂w.

This function is obviously nonnegative and L₁(u, v, w) = 0 only at DFS.

Moreover, L₁(u, v, w) → ∞ as ||(u, v, w)|| → ∞. By a lot of tedious calcula- tions, under the assumption β₂ < dk +^d_s (i.e. DFS is locally stable) we also obtain

dL1

dt (u, v, w) ≤ −



1 −k 2 − s

2u



(u − u)²+ 1 u −k

2 − s 2u



(v − v)²



in the set E ; cf. calculations in Qiao et al. [5]. Notice that if 1 − ^k₂ − _2u^s > 0 for any u, then the derivative of L₁ along solutions of Eqs. (3) is negatively defined yielding global stability of DFS. For k < 2 in the set E this inequality is equivalent to δ_u ≥ _2−k^s . On the other hand, it is clear that k < 2 is the necessary condition as for k ≥ 2 we have 1 − ^k₂ ≤ 0. Hence, we are able to formulate the stability result similar to those presented in Qiao et al. [5].

Corollary 3.1 If β2 < dk +^d_s, k < 2 and 1 − ks > δ_u ≥ _2−k^s , then the state DFS is globally asymptotically stable in the set E .

Notice that the inequality 1 − ks > s

2 − k for k < 2 is necessary to find a proper set E . This inequality is equivalent to

sk²− (1 + 2s)k + 2 − s > 0,

for which the discriminant ∆_k= (2s−1)²+4s²is always positive. Hence, there are two zeros of the left-hand side, k₁ = 1 + ¹⁻

√∆_k

2s < 1 + ¹⁺

√∆_k

2s = k2, both positive. Moreover, as √

∆_k > |2s − 1|, it is easy to check that k₁ < 2 < k₂ and k₁ < ¹_s. Now, we are in a position to propose a new stability result, which is more precise comparing to those stated in Corollary3.1.

Theorem 3.2 If k < k1 and β₂ < dk + ^d_s, then the state DFS is globally asymptotically stable in (R⁺)³.

Proof Under the assumptions we have k < 2 and 1 − ks > _2−k^s . Hence, there exists δ_u = _2−k^s for which we define the set E for which Corollary 3.1 implies global stability. However, this set attracts all solutions, meaning global

stability of DFS. _

3.2. Stability of ES Endemic steady state (u^∗, v^∗, w^∗) has all positive coordinates which are expressed as follows

u^∗ = 1 2β₁β₂



β1β2− kdβ₁− ks − 1 +p

(β1β2− kdβ₁− ks − 1)²+ 4β1β2

 , v^∗= 1 − u^∗

k , w^∗ = β2v^∗− d.

Notice that v^∗ > 0 ⇐⇒ u^∗ < 1. However, it is easy to check that u^∗ ∈ (0, 1) independently of the model parameters. Moreover,

w^∗ > 0 ⇐⇒ v^∗ > d β2

⇐⇒ u^∗ < 1 −dk β2

⇐⇒ β₂ > kd +d s,

which means that the state ES exists if and only if β₂ > kd +^d_s and then the state DFS is unstable.

Using the relation s − _u^v∗∗ − β₁w^∗ = 0 we can write the Jacobi matrix as

J_ES=





−u^∗ −ku^∗ 0

v^∗ u^∗

2

−^v_u^∗∗ −β₁v^∗ 0 β₂w^∗ −w^∗



. Moreover,

det (λI − JES) =

= λ³+



u^∗+ v^∗ u^∗ + w^∗

 λ²+



u^∗w^∗+v^∗w^∗

u^∗ + v^∗+k(v^∗)²

u^∗ + β1β2v^∗w^∗

 λ + v^∗w^∗+k(v^∗)²w^∗

u^∗ + β1β2u^∗v^∗w^∗,

and because coordinates and parameters are positive, it is easy to see that the Routh-Hurwitz criterion yields local stability of ES. It is important to notice that we do not need to pose any assumptions so if ES exists, then it is always locally stable. Now, we would like to investigate global stability. As before we define a Lyapunov function

L2(u, v, w) = u − u^∗− u^∗ln u

u^∗+ v − v^∗− v^∗ln v v^∗+β1

β₂



w − w^∗− w^∗ln w w^∗

 . Again, a lot of calculations give us a final inequality

dL2

dt (u, v, w) ≤

−



1 − v^∗ 2uu^∗ −k

2



(u − u^∗)²+ 1 u − v^∗

2uu^∗ −k 2



(v − v^∗)²+β1

β₂(w − w^∗)²

 . Using the same line of reasoning as in the DFS case we obtain the following corollary.

Corollary 3.3 If β2> dk +^d_s, k < 2 and 1 − ks > δ_u ≥ _k(2−k)u^1−u∗ ∗, then the state ES is globally asymptotically stable in the set E .

Now, we would like to check for which parameter values the condition 1 − ks > _k(2−k)u^1−u∗ ∗ necessary for stability of ES is satisfied. This condition is equivalent to

M

β₁β₂− N +p

(β₁β₂− N )²+ 4β₁β₂

> β₁β₂+N −p

(β₁β₂− N )²+ 4β₁β₂, (4) where M = k(1 − ks)(2 − k) > 0 under our assumptions (k < 2 and ks < 1), and N = kdβ₁+ ks + 1 > 1. Moreover, we can check that M < 1. Clearly, if we consider M = M (k, s), then the assumptions imply that this function is defined on O = (0, 2) × 0,_k¹
. Looking for maximal value of M on the closure of O we calculate partial derivatives

∂M

∂k = 2 − 2k − 4ks + 2ks², ∂M

∂s = −k²(2 − k).

Hence, ∇M 6= 0 in O, and therefore M achieves its maximal value on the boundary. Only for s = 0 the function M is non-zero on the boundary. Clearly, for s = 0 there is M = k(2 − k) and we see that the maximum is equal to 1.

This implies that M < 1 in O.

It is easy to see that the right-hand side of Ines. (4) is positive, such that this inequality is not satisfied for all parameter values. Let us rewrite Ineq. (4) in the following way:

p(β1β2− N )²+ 4β1β2(1 + M ) > β1β2(1 − M ) + N (1 + M ) ,

and after squaring both sides we obtain

(β₁β₂−N )²+4β₁β₂(1 + M )² > β₁²β₂²(1 − M )²+N²(1 + M )²+2β₁β₂N 1 − M²
, which is equivalent to

β₁ (β₂− 2kd)(1 + M )²− (β₂+ 2kd)(1 − M )²

> 2 (1 + ks)(1 − M )²− (1 − ks)(1 + M )²
, and eventually to

2β₁ M β₂− kd(1 + M²)
> ks(1 + M²) − M. (5)

Theorem 3.4 If k < 2, ks < 1 and 1. either ^1+ks_1−ks >

1+k(1−ks)(2−k) 1−k(1−ks)(2−k)

2

, β2 > max

kd(^{1+k2(1−ks)2(2−k)2})

k(1−ks)(2−k) , kd + ^d_s

 and β₁ > ^ks(^{1+k2(1−ks)2(2−k)2})−k(1−ks)(2−k)

2(k(1−ks)(2−k)β2−kd(1+k²(1−ks)²(2−k)²)); 2. or ^1+ks_1−ks <1+k(1−ks)(2−k)

1−k(1−ks)(2−k)

2

, kd +^d_s < β₂ < ^kd(^{1+k2(1−ks)2(2−k)2})

k(1−ks)(2−k) and β1 < k(1−ks)(2−k)−ks(^{1+k2(1−ks)2(2−k)2})

2(kd(1+k²(1−ks)²(2−k)²))−k(1−ks)(2−k)β2; then ES is globally asymptotically stable in (R⁺)³.

Proof We want to show that under the assumptions the statement of Corol- lary3.3holds. It could be checked that the assumptions of this theorem cover the assumptions of Corollary3.3. Hence, there exists δ_ε= _k(2−k)u^1−u∗ ∗ for which the state ES is globally stable in E . However, as E is globally attracting, then the state ES is globally stable in the whole phase space. _ 4. Model with diffusion Following the ideas presented in the original article [5] we are able to prove the following inequalities

lim sup

t→∞

max

x∈Ω

u(x, t) ≤ 1, lim sup

t→∞

max

x∈Ω

v(x, t) ≤ s, lim sup

t→∞

max

x∈Ω

w(x, t) ≤ wM, and

lim inf

t→∞ max

x∈Ω

u(x, t) ≥ 1 − ks, lim inf

t→∞ max

x∈Ω

v(x, t) ≥ (s(1 − β₁β₂) − β₁d)(1 − ks), lim inf

t→∞ max

x∈Ω

w(x, t) ≥ β₂(s(1 − β₁β₂) − β₁d)(1 − ks) − d.

Notice, that in [5] it is assumed that our space domain has dimension n = 1.

However, the results could be naturally expanded for arbitrary n.

Inequalities above imply that System (2) is persistent if ks < 1, β₁s > d, β₂s > d and β₂(s(1 − β₁β₂) − β₁d)(1 − ks) > d. The assumption β₁s > d has not occurred before but is necessary to get positivity of the bottom bound of lim inft→∞max_x∈Ωv(x, t).

Similar calculations to those presented in [5] show that appearance of diffusion does not change conditions of local stability of ES, i.e. if ES ex- ists (β₂ > dk + ^d_s), then it is locally asymptotically stable (by linearization and Routh-Hurwitz criterion). To investigate global stability we define the following Lyapunov function

Z(t) = Z

Ω

L₂(u(x, t), v(x, t), w(x, t))dx,

where L₂is the Lyapunov function for the case without diffusion. Using Neu- mann boundary conditions and earlier estimations we obtain

dZ(t) dt =

Z

Ω

∇L₂(u, v, w) × ∂u(x, t)

∂t ,∂v(x, t)

∂t ,∂w(x, t)

∂t

T

dx

= Z

Ω

 u − u^∗

u d₁∆u + v − v^∗

v d₂∆v +β₁ β₂

w − w^∗ w d₃∆w

 dx +

Z

Ω



(u − u^∗)(1 − u − kv) + (v − v^∗) s − v

u − β₁w + +β₁

β2

(w − w^∗)(β2v − d − w)

 dx

≤ − Z

Ω

 d1

u^∗

u²|∇u|²+ d2

v^∗

v²|∇v|²+β1

β₂d3

w^∗ w²|∇w|²

 dx

− Z

Ω

 

1 − v^∗

2(1 − ks)u^∗ −k 2



(u − u^∗)² +



1 − v^∗

2(1 − ks)u^∗ −k 2



(v − v^∗)²+β₁ β2

(w − w^∗)²

 dx.

On that basis we obtain:

Corollary 4.1 If β2 > dk + ^d_s, k < 2, 1 − ks > 0 and _k(2−k)u^1−u∗∗ < 1 − ks then the state ES is globally asymptotically stable.

Analogously, conditions guaranteeing local and global stability of DFE do not change.

Corollary 4.2 If β2 < dk + ^d_s, k < 2, 1 − ks > 0 and _2−k^s < 1 − ks then the state DFS is globally asymptotically stable.

5. Discussion Using the model analysis presented in this paper we would like to discuss the biological meaning of the conditions for the stability of ES versus the stability of DFE. In this discussion we focus on the case without diffusion as main conditions for the diffusive case are the same. Moreover, we only discuss the local stability, as the condition for local stability of ES are easy to interpret, while conditions for global stability are very complicated.

To take a look what our results mean we need to reverse the change of the variables and parameter naming. We start from invariant sets which after reverse reads:

• D = [0, K] ×0,^K_h × [0, ˜w_M],

• E = [˜δ_u, K] × [0,^K_h] × [0, ˜w_M] for kK < rh and every ˜δ_u< 1 − ^kK_rh
K,

• G = [˜δu, K] × [˜δv,^K_h] × [˜δw, ˜wM] for kK < rh, β²K < µsh, d < ^βK_h (1 −

β²K

µsh)(1 − ^kK_rh) and every ˜δ_u < 1 −^kK_rh
K, ˜δ_v < _h¹(1 − ^β_µsh^2K)˜δ_u and δ˜w < ˜w_M,

where ˜wM = ^βK−dh_µh if it is positive or arbitrary small positive number other- wise. This is associated with the fact that for βK ≤ dh (that is the rate of infection and the carrying capacity for preys are sufficiently small) there is

˙

w ≤ −µw² implying that the size of the infected subpopulation of predators will go to extinction. Notice that upper bounds in these sets essentially depend on the carrying capacity for prey K, that is the larger the parameter K, the larger the invariant sets are. The upper bound for healthy predators depends also on the competition coefficient h – with increasing h this bound decreases.

The lower bound also depends on K. Moreover, to have all solutions bounded above from 0, the inequality guaranteeing that w does not tend to 0 must hold.

Hence, in the set G the parameters K and β must be sufficiently small, but on the other hand, the death rate for infected predators d is small comparing to β and K.

The most important inequality that governs the existence and local sta- bility of the positive steady state ES changes in the following way:

β₂> dk + d

s 7−→ β > d k r + h

K



, (6)

which means that the state DFS is stable if the rate of infection is sufficiently small comparing to the death rate of the infected predators, while the state ES exists and is stable in the reverse case. Let us assume that Ineqs. 6holds and try to suggest some actions preventing the transmission of the disease.

This means that we would like to increase the right-hand side. It could be achieved in one of the following ways:

• decreasing the carrying capacity K or the growth rate r for preys;

• increasing the death rate d of the infected preys;

• increasing the hunting coefficient k;

• increasing the competition coefficient h for the healthy preys.

From the biological point of view, the coefficients that would possibly be easiest to control are r and d, as we can try to control birth of preys or hunt the infected preys. These parameters are potential goals for the optimal control problems.

References

[1] C. Holling. The functional response of predators to prey density and its role in mimicry and population regulation. Memoirs of the Entomo-Logical Society of Canada, 45:5–60, 1965. ISSN 0071-075X (Print), 1920-3047 (Online). doi: 10.4039/entm9745fv. URL https://doi.org/10.4039/

entm9745fv.

[2] S. B. Hsu and T. W. Huang. Global stability for a class of predator-prey systems. SIAM J. Appl. Math., 55(3):763–783, 1995.

[3] R. M. May. Stability and Complexity in Model Ecosystems. Princeton University Press, 1973. ISBN 978-0-691-08861-7.

[4] J. D. Murray. Mathematical Biology I. An Introduction. Springer, third edition, 2002. ISBN 978-0-387-22437-4. doi: 10.1137/S0036139993253201.

URL https://doi.org/10.1137/S0036139993253201.

[5] M. Qiao, A. Liu, and U. Foryś. Qualitative analysis for a reaction- diffusion predator-prey model with disease in the prey species. J.

Appl. Math., pages Art. ID 236208, 9, 2014. ISSN 1110-757X. doi:

10.1155/2014/236208. URL https://doi.org/10.1155/2014/236208.

Analiza modelu drapieżnik-ofiara z chorobą w populacji drapieżników.

Piotr Radziński i Urszula Foryś

Streszczenie Przeanalizowaliśmy dyfuzyjny model drapieżnik-ofiara z chorobą w populacji drapieżników zaproponowany przezQiao et al.(2014). W orginalnym arty- kule popełniono błąd podczas ubezwymiarowienia modelu. W niniejszej publikacji poprawiliśmy błędy w ubezwymiarowieniu i pokazaliśmy, że małe zmiany są po- trzebne aby uzyskać wyniki podobne do wyników uzyskanych przezQiao et al.(2014).

Zaproponowaliśmy poprawione warunki na globalną stabilność jednego z istniejących stanów stacjonarnych – stanu wolnego od choroby lub stanu endemicznego, zarówno w przypadku z dyfuzją jak i bez niej. Na podstawie poprawionej analizy zapropono- waliśmy nowe warunki na stabilność.

Klasyfikacja tematyczna AMS (2010): 92D40; 92D30.

Słowa kluczowe: model drapieżnik-ofiara, eko-epidemiologia, lokalna i globalna sta- bilność.

Urszula Foryś is a professor of mathematics at the University of Warsaw. She is an expert in mathematical modelling of biomedi- cal phenomena, mainly related to tumour growth and treatment, and eco-epidemiological modelling. However, she is also inter- ested in other applications, like modelling dyadic interactions, and lastly – neuroscience. She wrote a text-book “Matematyka w biologii” which is now one of the most important basis for courses of biomathematics in Polish universities.

Piotr Radziński holds Bachelor degree in Mathematics, currently studies Mathematics and writes Master’s thesis about proteolytic activity in cancer. His interests lie in biostatistics and biomedical data analysis.

Piotr Radziński University of Warsaw

Faculty of Mathematics, Informatics and Mechanics Institute of Applied Mathematics and Mechanics ul. Banacha 2, PL-02-097 Warszawa, Poland E-mail: pmradzinski@gmail.com

Urszula Foryś University of Warsaw

Faculty of Mathematics, Informatics and Mechanics Institute of Applied Mathematics and Mechanics ul. Banacha 2, PL-02-097 Warszawa, Poland E-mail: urszula@mimuw.edu.pl

Communicated by: Mirosław Lachowicz

(Received: 21th of April 2018; revised: 19th of June 2017)