www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Bi-stable dynamics of a host-pathogen model Roumen Anguelov∗†, Rebecca Bekker∗, Yves Dumont∗‡ ∗Department of Mathematics and Applied Mathematics, University of Pretoria, South Africa roumen.anguelov@up.ac.za, rebeccaabekker@gmail.com † Associate Member of the Institute of Mathematics and Informatics Bulgarian Academy of Sciences, Sofia , Bulgaria ‡ CIRAD, Umr AMAP, Pretoria, South Africa AMAP, University of Montpellier, CIRAD, CNRS, INRA, IRD, Montpellier, France yves.dumont@cirad.fr Received: 8 August 2018, accepted: 2 January 2019, published: 23 January 2019 Abstract—Crop host-pathogen interaction have been a main issue for decades, in particular for food security. In this paper, we focus on the modeling and long term behavior of soil-borne pathogens. We first develop a new compartmental temporal model, which exhibits bi-stable asymptotical dynamics. To investigate the long term behavior, we use LaSalle Invariance Principle to derive sufficient conditions for global asymptotic stability of the pathogen free equilibrium and monotone dynamical systems theory to provide sufficient conditions for perma- nence of the system. Finally, we develop a partially degenerate reaction diffusion system, providing a numerical exploration based on the results obtained for the temporal system. We show that a traveling wave solution may exist where the speed of the wave follows a power law. Keywords-Host-pathogen; bi-stability; monotone dynamical system; LaSalle Invariance Principle; partially degenerate reaction-diffusion; traveling wave I. INTRODUCTION The global food supply is currently experiencing pressure from climate change and ever increasing demand. Another major concern is the increasing impact of pathogens. An estimated 16% of the global crop yield is lost to various pathogens an- nually [3], [10], [14]. Consequently there has been an increase in research of botanical pathogens and the resulting diseases, with foliar pathogens being the focus of the majority of published work. One important difference between foliar and soil-borne pathogens is the environment wherein each occurs. Foliar pathogens have to contend with external factors such as wind, radiation and varying tem- peratures. However, the soil environment dampens the effects of such factors, although the inherent opacity of soil poses a number of challenges of its own. Added to these are challenges relating to capturing the direct and indirect influence of the environment on processes such as survival, Copyright: c©2019 Anguelov et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Roumen Anguelov, Rebecca Bekker, Yves Dumont, Bi-stable dynamics of a host-pathogen model, Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 1 of 17 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model dispersal and germination of pathogens, tissue growth, spatial distribution and susceptibility of hosts [12]. To some extent this paper is motivated by an early work of Gilligan [6], [7], [8], where he proposed a SEIR type compartmental model for the propagation of a soil-borne plant disease. This model includes a diffusion term, suggest- ing movement of infectious individuals. Although correct in certain contexts, this is unlikely in plant populations. In the model proposed here we consider compartments for the pathogen, where in- fection/infestation occurs when pathogen attaches to susceptible roots. In the spatio-temporal model we consider diffusion of the unattached pathogen, which we believe is a more biologically realistic assumption. The paper is organized as follows. In the next section we construct the temporal host-pathogen model highlighting the assumptions on which it is based. Section 3 deals with the equilibria of the system. Sufficient conditions for extinction and persistence of the pathogen are presented in Section 4 and 5 respectively. The spatio-temporal model is numerically considered in Section 6. II. THE HOST-PATHOGEN MODEL We consider a population of susceptible host plants with a constant recruitment rate Λ, and a pathogen present in the soil. As usual, the compartments of susceptible and infective/infested hosts are denoted by S and I, respectively, and N = S + I is the total host population. The natural decay rate of the host is d per time unit, and infected hosts have an addition decay rate of α per time unit. We assume that the pathogen is dependent on its host for nutrients or energy, and as such has an expected off–host death rate δ per time units. After coming into contact with a susceptible host it attaches at rate ρ, and grows at intrinsic growth rate of λ, restricted by the carrying capacity γI of the infected/infested roots. The attached pathogen (compartment A) detaches from their hosts at a rate of σ per time unit. The unattached or free pathogen (compartment F ) is responsible for new infections/infestations. It is assumed that if the population of free pathogen is large, the transmission rate from S to I de- pends solely on a constant β and the level of susceptible hosts present. This type of incidence is called saturation incidence and we use the specific form βF M+F . Using a saturating infestation rate is motivated by biological observations that increas- ing the free pathogen beyond a certain level no longer increases infestation proportionally. From a mathematical point of view, if only mass action principle is applied e.g. βFS, then, since F can potentially be very large, S would decrease very rapidly, which is unrealistic. For simplicity, the attachment rate (transfer from F to A is just mass action principle, namely ρFS. However, the growth in the A compartment is limited through the carrying capacity γI = γ(N − S). Since S cannot decrease unrealistically quickly then A cannot increase unrealistically quickly. The flow Fig. 1. Flow chart of the host-pathogen model chart is given in Figure 1. The model is a system of ODEs presented below: dA dt = λA(γI −A) −σA + ρFS dF dt = −δF + σA−ρFS dS dt = Λ −dS − βF M + F S dI dt = βF M + F S − (α + d)I (1) Using the notation x = (A,F,S,I)T the model Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 2 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model (1) is written as ẋ = f(x) (2) with f(x) =   λA(γI −A) −σA + ρFS −δF + σA−ρFS Λ −dS − βF M+F S βF M+F S − (α + d)I   . The local existence and uniqueness of solutions of (1) in R4+ = {x ∈ R4 : x ≥ 0} follows from the fact that f is continuously differentiable on R4+. The vector field defined by f points inwards on the boundary of R4+. Hence, R 4 + is positively invariant. In order to obtain global existence of solutions it remains to show that all solutions initiated in R4+ are bounded. Adding the equations for dS dt and dI dt we have the inequalities Λ−(d+α)(S+I)≤ d(S+I) dt ≤ Λ−d(S+I), (3) which do not depend on the other coordinates of x. Hence, for any solution we have that the interval[ Λ α + d , Λ d ] is a global attractor for S(t) + I(t). More precisely, we have Λ α + d ≤ lim inf t→+∞ (S(t) + I(t)) ≤ lim sup t→+∞ (S(t) + I(t)) ≤ Λ d . (4) Since S and I are also nonnegative, they are bounded. Using (4) it is easy to obtain that the rest of the coordinates of any solution x(t) are also bounded. In fact, since the bounds in (4) do not depend on the initial condition, one obtains that (1) defines a dynamical system on R4+, which is dissipative. III. EQUILIBRIA In the absence of pathogen the population of the host is S and it is governed by the third equation in the model (1). It has an asymptotically stable equi- librium at Λ d with basin of attraction S ∈ [0, +∞). The resulting equilibrium of the model (1) we call Pathogen Free Equilibrium (PFE), that is, PFE = ( 0, 0, Λ d , 0 ) . (5) The basic reproduction number/ratio, R0, is a threshold quantity, which is often used to char- acterize the properties of epidemiological models. It is popularly defined as the number of new infections caused by a single infectious individual in a wholly susceptible population. Its precise definition is that it is the spectral radius of the next generation matrix calculated at an asymptotically stable equilibrium of the population in the absence of disease, [17]. Such equilibrium is usually re- ferred to as Disease Free Equilibrium. Due to the nature of the model in this paper, and as mentioned above, we use the term Pathogen Free Equilibrium (PFE). The model (1) has a unique Pathogen Free Equilibrium given by (5). Following the method in [17] for the compu- tation of R0, the compartment vector x in the model (1) is decomposed into three-dimensional vector of pathogen related compartments y = (A,F,I)T and one pathogen free compartment S and we have ẏ = F(y,S) −V(y,S), where F(y,S) =  λA(γI −A) + ρSF0 βFS M+F   , V(y,S) =   σAδF −σA + ρFS (α + d)I   . The next generation matrix is ∂F ∂y ( 0, Λ d )( ∂V ∂y ( 0, Λ d ))−1 = 1 M(δ + ρΛ d )  ρM Λd ρM Λd 00 0 0 β Λ d β Λ d 0   . Thus R0 = ρΛ dδ + ρΛ . Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 3 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model It is clear that R0 < 1 regardless of the values of the parameters (as long as they are positive). Therefore, it follows from [17, Theorem 2] that the PFE is always asymptotically stable. The asymp- totic behavior of the solutions of (1) depends on the existence of pathogen-endemic equilibria or other invariant sets and their properties. It is interesting to remark that for the model under consideration R0 is not a threshold quantity at all. In order to obtain all equilibria, we set right hand side in model (1) equal to zero, λA(γI −A) −σA + ρFS =0, −δF + σA−ρFS =0, Λ −dS − βF M + F S =0, βF M + F S − (α + d)I =0. (6) From the third equation in (6) we obtain an expression for S in terms of F : S = Λ(M + F) (d + β)F + dM . (7) Then, using (7), from the second equation in (6) we obtain an expression for A: A = σ−1(δ + ρS)F = σ−1 ( δ + ρΛ(M + F) (d + β)F + dM ) F. (8) Similarly, using (7), the fourth equation in (6) gives an expression for I in terms of F : I = βFS (α + d)(M + F) = ΛβF (α + d)((d + β)F + dM) . (9) Substituting these expressions into the first equa- tion and excluding the case F = 0, we obtain a cubic equation about F in the form −a1F 3 + a2F 2 + a3F −a4 = 0, (10) where a1 = λ(α + d)[δ(d + β) + ρΛ] 2, a2 = −2λ(α + d)[δ(d + β) + ρΛ](δd + ρΛ)M −δ(α + d)(d + β)2σ2 +λγΛσβ[δ(d + β) + ρΛ], a3 = λγσΛβ(δd+ρΛ)M−2δ(α+d)(d+β)dMσ2 −λ(α + d)(δd + ρΛ)2M2, a4 = δ(α + d)σ 2d2M2. Clearly, a1 > 0 and a4 > 0, while the signs of a2 and a3 may vary depending on the values of the parameters. However, it is easy to see that for any signs of a2 and a3 there are always either two sign changes in the sequence of the coefficients of (10) or no sign changes at all. Hence the equation has either two positive roots or no positive roots. When these roots exist we denote them by F1 and F2, with F1 ≤ F2. The respective equilibria of the model (10) are denoted by EE1 = (A1,F1,S1,I1) T and EE2 = (A2,F2,S2,I2)T . From the expressions (7), (8) and (9) we see that EE1 > 0 and EE2 > 0. Further, we can also see from the expressions (7)–(9) that S is a decreasing function of F , while I and A are increasing functions of F . Hence we have A1 ≤ A2, I1 ≤ I2 and S1 ≥ S2. (11) The two positive roots of (10) appear simultane- ously as a double root F1 = F2, which then splits into two simple roots. Hence, in the bifurcation state when F1 = F2, the equilibrium EE1 = EE2 appears and then splits into the two distinct equi- libria EE1 and EE2. Since the constant term of (10) is strictly positive, this bifurcation is bounded away from the PFE. The PFE does not undergo any bifurcation and, as mentioned above, it is always asymptotically stable. We perform numerical simulations using non- standard finite difference schemes [2] to solve systems (1) and (33). In all numerical simulations of model (1) we observe two qualitatively different cases and the transition (bifurcation) from one to the other. The Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 4 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model one case is when PFE is the only equilibrium of the system. In this case we observe that all solution converge to PFE, that is PFE is globally asymp- totically stable on R4+. An illustrative example is given Figure 2 with value of the parameters given in Table I. TABLE I PARAMETER VALUES USED IN FIGURE 2 Parameter Value Parameter Value Λ 1.1000 λ 0.5000 γ 0.5000 β 1.0000 σ 0.4000 ρ 0.2000 δ 0.1000 M 10.000 α 0.5100 d 0.5000 The second case is when the model has two positive equilibria. In the simulations presented in Figure 3, EE2 is stable and attracting, while EE1 is unstable (saddle point). Table II contains the pa- rameter values used for the simulations in Figure 3. The solutions that are initiated below EE1 in the (A,F,I)-space converge to the PFE. Solutions φ1 and φ2 are initiated at EE1 with altered value of S, S0 = 2 and S0 = 20 respectively. These values are below and above the PFE value of S, respectively. We observe that φ1 converges to the PFE, while φ2 increases and eventually converges to EE2. The unstable equilibrium is typically very close to the PFE, so that the basin of attraction of the PFE is relatively small. Nevertheless, it contains the whole nonnegative S-axis. TABLE II PARAMETER VALUES USED IN FIGURE 3. Parameter Value Parameter Value Λ 1.0000 λ 0.5000 γ 0.9000 β 5.0000 σ 0.4000 ρ 0.1000 δ 0.1000 M 100.000 α 0.0450 d 0.1000 Due to the complexity of the model, we could not obtain theoretically a general result for the observed properties of the positive equilibria, or alternatively the global asymptotic stability of the PFE. However, we derive in the next two sections sufficient conditions for the two practically impor- tant properties: extinction and persistence of the pathogen. IV. SUFFICIENT CONDITIONS FOR GLOBAL ASYMPTOTIC STABILITY OF PFE We prove sufficient conditions for global asymptotic stability of PFE using LaSalle’s Invari- ance Principle [11, Theorem 2]. Theorem 1. The PFE of model (1) is globally asymptotically stable on R4+ if either condition a) or condition b) below hold: a) λβγ2Λ2 16d2(α + d)δM ≤ 1 (12) b) β < α + d and λβ2γ2Λ2 4d2δM(β + α + d)2 ≤ 1. (13) Proof: Taking into account the inequalities (4), it is sufficient to consider the system (1) on the domain Ω = { x ∈ R4+ : S + I ≤ Λ d } . We consider on Ω the function V (x) = A + F + ξ 2 I2, where ξ is a positive constant with value yet to be determined. We have V̇ (x) = Ȧ + Ḟ + ξIİ = λA(γI −A) − ξ(α + d)I2 + ξ βF M + F SI −δF. Since the first term is a quadratic of A it obtains it largest value when A = 1 2 γI. Using also that F ≥ 0 and SI ≤ ( S + I 2 )2 ≤ Λ2 4d2 , we have V̇ (x) ≤ ( λγ2 4 − ξ(α + d) ) I2 + ( ξβΛ2 4d2M −δ ) F (14) Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 5 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 2. Illustration of the case of global asymptotic stability of the PFE. The coefficient of I2 is nonpositive if and only if ξ ≥ λγ2 4(α + d) (15) Similarly, the coefficient of F is nonpositive if and only if ξ ≤ 4d2Mδ βΛ2 (16) A constant ξ satisfying (15) and (16) exists if and only if λγ2 4(α + d) ≤ 4d2Mδ βΛ2 , or, equivalently, λβγ2Λ2 16d2(α + d)δM ≤ 1, (17) that is if and only if condition a) holds. Hence, if a) holds, we can select ξ such that V̇ (x) ≤ 0 for x ∈ Ω. Let E = {x ∈ Ω : V̇ (x) = 0}. According to LaSalle Invariance Principle, all solutions con- verge to the largest invariant set M of (1) which is contained in E. If the inequality in (17) is strict, then ξ can be selected in such a way that the coefficients in (14) are negative. Hence, V̇ (x) = 0 only if I = F = 0 and then A = 0 as well. Hence, E = { (0, 0,S, 0) : 0 ≤ S ≤ Λ d } . (18) The largest invariant set in E given in (18) is M = {PFE}. If the inequality in (17) holds as equality, then the right hand side of (14) is zero irrespective of the values of I and F . However, if I or F is not zero, then the respective inequalities leading to (14) should be satisfied as equalities. This process leads enlargement of E in (18) by an isolated point ( γΛ 4d , 0, Λ 2d , Λ 2d ) and we again have M = {PFE}. Therefore, by the LaSalle Invariance Principle all solutions converge to PFE. b) The condition β < α+d provides for slightly sharper upper bound of SI and hence slightly weaker condition on the parameters. Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 6 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 3. If solutions of model (1) are initiated ‘close’ to the PFE convergence to this equilibrium occurs. Solutions that are initiated from the ‘smaller’ endemic equilibrium converge to either the PFE or the ‘larger’ equilibrium, dependent on the initial density of the susceptible compartment. From the last equation of the model (1) we have dI dt ≤ βS − (α + d)I ≤ β(S + I) − (β + α + d)I ≤ βΛ d − (β + α + d)I. Therefore lim sup t→+∞ I(t) ≤ βΛ d(β + α + d) . Hence, in the investigation of the asymptotic be- havior of the system we can consider only the subset of Ω where I ≤ βΛ d(β + α + d) . The inequality β < α + d implies that βΛ d(β + α + d) ≤ Λ 2d . Hence, SI ≤ ( Λ d − I ) I ≤ ( Λ d − βΛ d(β + α + d) ) βΛ d(β + α + d) = β(α + d)Λ2 (β + α + d)2d2 . Then, following the same method as for a), in place of (17) we obtain the second inequality in b). The parameter values given in Table I satisfy both condition (a) and condition (b) of Theorem 1. Hence, the global asymptotic stability of PFE observed earlier in Figure 2 can be deduced from either (a) or (b). Another illustration is given on Figure 4. The model parameters, given in Table III, satisfy condition (a), but not condition (b) in Theorem 1. This is sufficient to deduce the global asymptotic stability of PFE seen in Figure 4. We need to remark the conditions (a) and (b) are each only sufficient, but not necessary for the Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 7 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model TABLE III PARAMETER VALUES USED IN FIGURE 4 Parameter Value Parameter Value Λ 0.9000 λ 0.5000 γ 0.5000 β 1.0000 σ 0.2000 ρ 0.3000 δ 0.2900 M 5.000 α 0.0600 d 0.1500 global stability of the PFE. The parameter values given in Table IV satisfy neither of the conditions in Theorem 1. Yet, global asymptotic stability of PFE can be observed in the simulations presented in Figure 5. TABLE IV PARAMETER VALUES USED IN FIGURE 5 Parameter Value Parameter Value Λ 1.1000 λ 0.6900 γ 0.6000 β 1.5000 σ 0.4000 ρ 0.5000 δ 0.2000 M 5.000 α 0.0100 d 0.3000 V. SUFFICIENT CONDITIONS FOR PERSISTENCE We derive sufficient conditions for persistence using the theory of monotone systems. The system (1) is not monotone. We construct an auxiliary system about the vector y = (A,F,I)T which is monotone. From the third equation of (1) we have Λ − (d + β)S ≤ dS dt ≤ Λ −dS (19) Then, it follows that for every solution of (1) it holds Λ β + d ≤ lim inf t→+∞ S(t) ≤ lim sup t→+∞ S(t) ≤ Λ d . Hence, for the asymptotic properties of the solu- tions of (1) it is sufficient to consider the subset of Ω where Λ β + d ≤ S ≤ Λ d . (20) In this subdomain we consider the following sys- tem for y = (A,F,I)T : dy dt = h(y) :=   λA(γI −A) −σA−δF + σA− ρΛ d F βF M+F Λ β+d − (α + d)I   . (21) Let us recall that function h is said to satisfy the Kamke condition if hi is increasing in yj for i 6= j. The Kamke condition implies that the respective system of ODEs is monotone with respect to the initial condition or shortly monotone, [16, Section 3.1]. The Jacobian of h is Jh =  λγI−2λA−σ 0 λγAσ −δ − ρΛd 0 0 βMΛ (β+d)(M+F)2 −(α+d)   . Since the nondiagonal entries of Jh are nonnega- tive, the system (21) is monotone. Moreover, since the system is irreducible in the interior of R3+, it is strongly monotone, [16, Theorem 4.1.1]. Let us recall that a system of the form (21) is called strongly monotone if for any two solution y(1) and y(2) y(1)(0) < y(2)(0) =⇒ y(1)i (t) < y (2) i (t), t > 0, i = 1, 2, 3. Our interest in the system (21) is motivated by the fact that its solutions provide lower bounds for the coordinates A, F , and I of the solutions of (1). This will be shown later by using differential inequalities given in [18] for systems of ODEs with quasi-monotone right hand side. Hence, we carry out first the asymptotic analysis of (21). To find the equilibria of (21), we set the right hand side to zero: λA(γI −A) −σA = 0 (22) −δF + σA− ρΛ d F = 0 (23) βF M + F Λ β + d − (α + d)I = 0 (24) The origin, 0 is an equilibrium. To find the nonzero equilibria, we multiply the first equation by σ λ and add it to the second one to eliminate A. We Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 8 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 4. Illustration of the global stability of the PFE of the host-pathogen model when only condition (a) in Theorem 1 holds. multiply the third equation by γσ α+d and add it to the second one to eliminate I. Hence, we obtain ϕ(y) := σ λ h1(y) + h2 + γσ α + d h3(y) = (25) − ( δ+ ρΛ d ) F− σ2 λ + γσβΛF (α+d)(β+d)(M +F) = 0 (26) Equation (26) is equivalent to the quadratic equa- tion ( δ + ρΛ d ) F 2+( M ( δ + ρΛ d ) + σ2 λ − γσβΛ (β + d)(α + d) ) F + Mσ2 λ = 0. (27) The equation (27) has two positive real roots if and only if the coefficient of F is negative and the discriminant is positive, that is M ( δ+ ρΛ d ) + σ2 λ − γσβΛ (β+d)(α+d) < 0, (28) ∆ = ( M ( δ+ ρΛ d ) + σ2 λ − γσβΛ (β+d)(α+d) )2 − 4Mσ2 λ ( δ + ρΛ d ) > 0 (29) Assuming that conditions (28)–(29) hold, we denote by F̃1 and F̃2, F̃1 < F̃2, the roots of (27) and by Ẽ1 = (Ã1, F̃1, Ĩ1)T and Ẽ2 = (Ã2, F̃2, Ĩ2)T the corresponding equilibria of (21). where Ai = 1 σ ( δ + ρΛ d ) Fi, i = 1, 2, (30) Ii = 1 γ (Ai + σ) , i = 1, 2. (31) Considering the expressions (30) and (31, we have 0 << Ẽ1 << Ẽ2. Theorem 2. Let conditions (28)–(29) hold. Then for every solution y(t) of (21) we have y(0) > Ẽ1 =⇒ lim inf t→+∞ y(t) ≥ Ẽ2. Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 9 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 5. Illustration that PFE of the host-pathogen model may be globally asymptotically stable when neither of the condition in Theorem 1 holds. Proof: It is easy to see that the equilibrium 0 is asymptotically stable, indeed the eigenvalues of Jh(0), ξ1 = −σ, ξ2 = − ( δ + ρΛ d ) and ξ3 = −(α + d) are all negative. We consider the order interval [0, Ẽ1]. It follows from [16, Theorem 2.2.2], that the solutions ini- tiated in this interval, excluding the end points, either all converge to 0 or all converge to Ẽ1. Since 0 is asymptotically stable, this implies that all solutions converge to 0. The Jacobian of h at Ẽ1 after some simplifications is: Jh(Ẽ1) =  −λA ∗ 0 λγA∗ σ −δ − ρΛ d 0 0 βMΛ (β+d)(M+F̃1)2 −(α + d)   , Since the nondiagonal entries of Jh(Ẽ1) are non- negative and the matrix is irreducible, it follows from the theory of nonnegative matrices [4, The- orem 2.1.4] that Jh(Ẽ1) has eigenvector v with positive coordinates and associated eigenvalue ξ, which is a real number. Since Ẽ1 is repelling in [0, Ẽ1], we have that ξ ≥ 0. We will show that ξ > 0. Assume the opposite, namely ξ = 0. Then it is easy to compute that w = ( σ γÃ1 , 1, σγ α + d )T is a left eigenvector. Then, using the expression for ϕ in (25) and that h(Ẽ1) = 0 we have ∇ϕ(Ẽ1) = wTJh(Ẽ1) = 0. The first and the third coordinates of ∇ϕ(y) are always zero since ϕ does not depend on A and I. The fact that ∂ϕ(Ẽ1) ∂F = 0 implies that the equation (26), or, equivalently (27, has a double root. This contradicts (29). Therefore ξ > 0. Next, we consider the order interval [Ẽ1, Ẽ2]. Again following [16, Theorem 2.2.2], that the solutions initiated in this interval, excluding the end points, either all converge to Ẽ1 or all converge to Ẽ2. Considering that Ẽ1 is repelling in the direction of the positive vector v, we conclude that all solutions converge to Ẽ2. Let y(t) be any solution of (21) such that y(0) > Ẽ1. Consider the solution of (21) with initial condition z(0) = min{y(0), Ẽ2}. Since z(0) ∈ [Ẽ1, Ẽ2] and z(0) > Ẽ1, we have lim t→+∞ z(t) = Ẽ2. Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 10 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Then using that y(0) ≥ z(0) and the monotonicity of the system (21) we have y(t) ≥ z(t), t ≥ 0. Therefore lim inf t→+∞ y(t) ≥ lim t→+∞ z(t) = Ẽ2. Using Theorem 2 for the auxiliary system (21) we prove the following theorem for the original model (1). Theorem 3. Let conditions (28) and (29) hold. Then, for any solution of (1) we have that if A(0) > Ã1, F(0) > F̃1, I(0) > Ĩ1, (32) then lim inf t→+∞ A(t) ≥ Ã2, lim inf t→+∞ F(t) ≥ F̃2, lim inf t→+∞ I(t) ≥ Ĩ2. Proof: As discussed earlier, it is sufficient to consider the subdomain, where (20) holds. Let x(t) be a solution of (2) such that A(0), F(0), I(0) satisfy (32) and S(0) ∈ [ Λ β+d , Λ d ] . From (19) it follows that S(t) ∈ [ Λ β+d , Λ d ] for t ≥ 0. Then the coordinates A(t), F(t), and I(t) satisfy d dt  A(t)F(t) I(t)   =  f1(x)f2(x) f4(x)   ≥ h  A(t)F(t) I(t)   Using that h is quasi-monotone and applying [18, Chapter 2, Section 12.X], we obtain that the vector function (A(t),F(t),I(t))T is bounded below by the solution of (21) with the same initial con- dition at t = 0. Then, Theorem 2 implies that Ẽ2 is a lower bound for the limit inferior of (A(t),F(t),I(t))T as t → +∞, which proves the theorem. Theorem 3 shows that if the initial invasion is sufficiently large the pathogen establishes itself at a level above Ẽ2. We should remark that it is not necessary to have initially all compartments A, F and I above Ẽ1. It is sufficient that at some future time they all exceed Ẽ1. For example, and as it can be also expected from biological point of view, the initial infection/infestation is brought in the compartment F . If this initial value of F is sufficiently large that A and I increase above Ã1 and Ĩ1, while F is still above F̃1, the pathogen will persist eventually at least at a level of Ẽ2. Theorem 3 is illustrated numerically by Fig- ure 6. The initial conditions of the solutions in this figure were chosen specifically so that (A0,F0,I0) ≥ Ẽ1. Clearly any solution initiated at or above the level of Ẽ1 ≈ (4.091, 1.7241, 0.2164) persists at a non-zero level above Ẽ2 for all time; and in fact converges to an equilibrium of model (1), at least for the parameter values in Table V. We verify that the data in Table V satisfies conditions (28)-(29) of Theorem 3. Indeed, the left hands sides of (28) and (29) are respectively negative and positive, so that inequalities in (28) and (29) hold. Theorem 3 gives only sufficient conditions for persistence of pathogen. In Figure 7, we show that the pathogen persist, in fact, the solutions converge to EE2, even though the initial conditions do not satisfy the requirements of Theorem 3. TABLE V PARAMETER VALUES USED IN FIGURES 6 AND 7 Parameter Value Parameter Value Λ 0.9000 λ 1.0000 γ 23.52536 β 5.0000 σ 1.0000 ρ 0.9000 δ 0.9000 M 10.000 α 0.0010 d 0.5500 Our numerical investigations of model (1) have revealed that in addition to the PFE which is al- ways locally attracting, there exists an equilibrium close to the PFE which is repelling, and an attract- ing equilibrium which is removed from the PFE. Although we have not proven this mathematically, the important properties of the model, namely extinction and persistence, have been proven under certain conditions. VI. THE SPATIO-TEMPORAL HOST-PATHOGEN MODEL The model in Section II looks only at the tem- poral progression of an infection, and although this Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 11 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 6. Illustration of persistence of the infection as proven in Theorem 3. Observe that solutions of model (1) originating at the level of Ẽ1 remain non-zero for all time. approach is acceptable under certain assumptions (such as a pathogen entering an entire field in a uniform manner), the model can be modified slightly to accommodate the spatial movement of pathogens through the field. Diffusion has been used to model spatial spread in theoretical ecology since the latter half of the twentieth century [9], [15], with its use for modelling fungal growth being justified by the “observation that tip growth occurs to fill space and to capture nutrients” [5]. Davidson (1998) also noted that fungal growth “is, in the main, directed from areas of high hyphal density to areas of low hyphal density”, and included diffusion in his model with the warning ‘that this flux should not be viewed as the movement of existing biomass, but rather the propensity of new biomass to grow away from high density areas’. We reiterate this warning, and include diffusion to model the spatial growth of off-host pathogen in search of new hosts, with µ denoting the diffusion constant. Our Host-Pathogen spatio-temporal model is defined as follows:  ∂A ∂t = λA(γI −A) −σA + ρFS ∂F ∂t = −δF + σA−ρFS + µ∆F ∂S ∂t = Λ −dS − βF M + F S ∂I ∂t = βF M + F S − (α + d)I with A(x, 0) ≥ 0, F(x, 0) ≥ 0, S(x, 0) ≥ 0, I(x, 0) ≥ 0, ∂F ∂x (−L,t) = 0 = ∂F ∂x (L,t). (33) If the initial condition is spatially uniform, and taking the boundary conditions into account, then the solution is also spatially uniform, reducing it to a solution of the corresponding temporal system. The properties and long-term behaviour of the temporal model have been theoretically proven in Section II, and this chapter devotes itself to numerical investigation of the behaviour of Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 12 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 7. Illustration of persistence of the infection when initial conditions do not satisfy conditions (32) in Theorem 3. solutions of system (33). Our interest is mainly in the practically relevant case where the pathogen is introduced in one location, and studying the dynamics of its propagation. In order to solve model (33) numerically, we use non-standard dis- cretization coupled with a second order central- space discretization [2]. A. Numerical investigations 1) Under the conditions for asymptotic stability obtained by application of LaSalle’s Invariance Principle: The parameter values given in Table VI satisfy the conditions (12) which ensures the PFE of the temporal model is globally asymptotically stable. The details are given in Section IV. We investigate whether the solutions of the spatio- temporal model behave in a similar fashion, using the diffusion constant µ = 0.01. Indeed, although convergence occurs over a long time period, the addition of diffusion does not result in observable change in the asymptotic properties of the steady state. In Figure 8, even assuming that the initial population has free pathogen over a quarter of the field, this does not result in the infection spreading. TABLE VI PARAMETER VALUES USED IN FIGURE 8 Parameter Value Parameter Value Λ 1.0000 λ 0.4000 γ 0.2000 β 1.0000 σ 0.0100 ρ 0.4000 δ 0.1000 M 100.000 α 0.0205 d 0.1000 TABLE VII PARAMETER VALUES USED IN FIGURE 9 Parameter Value Parameter Value Λ 1.0000 λ 0.7000 γ 0.4000 β 1.0000 σ 0.1000 ρ 0.2000 δ 0.2000 M 100.00 α 0.0100 d 0.2000 In Figure 9, we also provide an example where the conditions (12) are not satisfied and yet there is convergence to PFE (see also Table VII). 2) Parameter values for persistence of the in- fection: A monotone system, constructed to ap- proximate the temporal host pathogen model from below was proven to admit two interior equilibria Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 13 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 8. The time evolution of pathogen and disease through the field at different times, using the parameter values that satisfy the conditions for stability of the PFE that were obtained by the application of LaSalle’s Invariance Principle. Fig. 9. The parameter values in Table VII do not satisfy condition for asymptotic stability of the PFE, however convergence to PFE is observed. Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 14 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 10. When solutions of model (33) are initiated with pathogen and infectious hosts at the level of EE2, on the left boundary, a field of completely susceptible hosts will experience a travelling infection front. This front connects EE2 and the PFE. in Section V, denoted Ẽ1 and Ẽ1, with Ẽ1 < Ẽ1. Additional conditions were derived, under which the pathogen persists. Indeed, it was found that solutions initiated at or above Ẽ1 and satisfying conditions (28) and (29), page 9, remain non-zero for all t ≥ 0. For the parameter values in Table VIII on page 15, these conditions are satisfied. Indeed, we have M γσ ( δ + ρΛ d ) + σ γλ − βΛ (β + d)(α + d) ≈−0.4204 < 0, ∆ ≈ 0.0053 > 0. The equilibrium, Ẽ1 of the lower approximating system is: Ẽ1 ≈ (4.091, 1.7241, 0.2164). The persistence of pathogen is illustrated in Fig- ure 10. In fact, the solutions converge to EE2, although the stability properties of EE1 and EE2 have not been proven. How does the inclusion of diffusion affect this phenomenon? Solutions initiated at the level of EE2 at the boundary exhibit a travelling infection front, the movement of which is driven by the increase in attached pathogen and infested hosts by the diffusion of the free pathogen (Figure 10). This behaviour suggests a possible control strategy: if the speed of the front can be sufficiently decreased, a percentage of the field would be saved from disease. TABLE VIII PARAMETER VALUES USED IN FIGURE 10. Parameter Value Parameter Value Λ 0.9000 λ 1.0000 γ 23.52536 β 5.0000 σ 1.0000 ρ 0.9000 δ 0.9000 M 10.000 α 0.0010 d 0.5500 To this end, we investigate the relationship be- tween µ and the wave speed c. The parameter val- ues in Table VIII were again used, and a solution with (A0,F0,I0) taking the value of EE2 on the Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 15 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Fig. 11. The speed of the infection front for different values of µ, for solutions of model (33) initiated with levels of inoculum and disease at the level of EE2 on the left boundary. The parameter values in Table VIII were used. Clearly the equation c(µ) = 0.01088µ0.4189 fits the data well. left boundary was considered. The diffusion con- stant µ was taken to be in the interval [10−7, 10−1], which results in c ∈ (0, 4.5 × 10−3]. An equa- tion of the form c(µ) = aµb was fitted to the data in Figure 11, and the fitting process reveals a ∈ (0.010770, 0.011) and b ∈ (0.416, 0.4218) with 95% confidence. In fact, a = 0.01088 and b = 0.4189. Literature indicates that the value of b should be higher, with Gilligan [8] and Metz, Mollison and van den Bosch [13] finding the wave speed to be proportional to the square root of the diffusion constant; that is c ∝ √ µ. Although b < 0.5 the equation fits the data well, and since SSE = 8.59 × 10−6 its use in making predictions would be justified. The coefficient of determination r2 = 0.9933 indicating that 99.33% of the variance of the data is explained by the equation. VII. CONCLUSION In this work, we have derived a Host-Pathogen model where the PFE is always LAS and may co- exist with endemic equilibria. We provided suffi- cient conditions for PFE being globally asymptot- ically stable and for persistence of the pathogen, using two different approaches, LaSalle Invari- ance Principle approach and monotone system ap- proach. We show that these results can be extended to the spatio-temporal system, where we add diffusion in the free pathogen compartment. We also show numerically that a bi-stable travelling wave solution can exist between PFE and a stable endemic equilibrium, here EE2. We show that the speed c of this traveling wave is of the form aµb, where µ is the diffusion parameter. Further theoretical investigations are needed, in order to be able to derive appropriate control strategies to avoid invasion of the pathogen in the whole crop. REFERENCES [1] R. Anguelov, Y. Dumont, J. Lubuma (2012), Mathemat- ical modelling of sterile insect technology for control of anopheles mosquito, Computers & Mathematics with Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 16 of 17 http://dx.doi.org/10.11145/j.biomath.2019.01.029 R. Anguelov, R. Bekker, Y. Dumont, Bi-stable dynamics of a host-pathogen model Applications 64(3), 374-389. https://doi.org/10.1016/j.camwa.2012.02.068 [2] R. Anguelov, Y. Dumont, and J. M.-S. Lubuma (2012). On nonstandard finite difference schemes in biosciences, AIP Conf. Proc. 1487(1), 212-223. https://doi.org/10.1063/1.4758961 [3] D.P. Bebber, T. Holmes, S.J. Gurr (2014), The global spread of crop pests and pathogens, Glob. Ecol. Bio- geogr 23(12), 1398-1407. https://doi.org/10.1111/geb.12214 [4] A. Berman, R.J. Plemmons, Nonnegative matrices in the mathematical sciences, SIAM, 1994. [5] N.J. Cunniffe (2007), Dispersal of soil-borne plant pathogens and efficacy of biological control, PhD The- sis, University of Cambridge. [6] C.A. Gilligan (1985), Construction of temporal models: III. Disease progress of soil-borne pathogens. In: C.A. Gilligan (ed.), Mathematical modelling of crop disease, London, Academic Press. [7] C.A. Gilligan (1990), Mathematical modeling and anal- ysis of soilborne pathogens. In: J. Kranz (ed.) Epidemics of plant diseases: Mathematical analysis and modeling, Heidelberg, Springer-Verlag, 96-142. [8] C.A. Gilligan (1995), Modelling soil-borne plant pathogens: reaction-diffusion models, Canadian Journal of Plant Pathology 17, 96-108. [9] E.E. Holmes, M.A. Lewis, J.E. Banks, R.R. Veit (1994), Partial differential equations in ecology: spatial interac- Partial differential equations in ecology: spatial interac- tions and population dynamics, Ecology 75, 17-29. [10] State of the World’s Plants Report 2016, Royal Botanic Gardens, Kew, https://stateoftheworldsplants.org/2016/ [11] J. LaSalle (1960), Some Extensions of Liapunov’s Sec- ond Method, IRE Transactions on Circuit Theory 7(4), 520-527. [12] J.D. MacDonald, The Soil Environment. In: C.L. Camp- bell, D.M. Benson (eds), Epidemiology and Manage- ment of Root Disease, Springer-Verlag, 1994, 82-116. [13] J.A.J. Metz, D. Mollison, F. van den Bosch (2000), The Dynamics of Invasion Waves. In: U. Dieckmann, R. Law, J.A.J Metz (eds), The Geometry of Ecological In- teractions: Simplifying Spatial Complexity, Cambridge University Press, 482-512, https://doi.org/10.1017/CBO9780511525537.027 [14] E.C. Oerke (2006). Crop losses to pests, J.Agric. Sci. 144(1), 31-43. [15] A. Okubo, S.A. Levin, Diffusion and Ecological Prob- lems: Modern Perspectives, Springer, 2001. [16] H.L. Smith, Monotone Dynamical System, AMS, 1995. [17] P. van den Driessche, J. Watmough (2002), Reproduc- tion numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathe- matical Biosciences 180, 29-48, https://doi.org/10.1016/S0025-5564(02)00108-6 [18] W. Walter, Differential and Integral Inequalities, Springer, 1970. Biomath 8 (2019), 1901029, http://dx.doi.org/10.11145/j.biomath.2019.01.029 Page 17 of 17 https://doi.org/10.1016/j.camwa.2012.02.068 https://doi.org/10.1063/1.4758961 https://doi.org/10.1111/geb.12214 https://stateoftheworldsplants.org/2016/ https://doi.org/10.1017/CBO9780511525537.027 https://doi.org/10.1016/S0025-5564(02)00108-6 http://dx.doi.org/10.11145/j.biomath.2019.01.029 Introduction The host-pathogen model Equilibria Sufficient conditions for global asymptotic stability of PFE Sufficient conditions for persistence The spatio-temporal host-pathogen model Numerical investigations Under the conditions for asymptotic stability obtained by application of LaSalle's Invariance Principle Under the conditions for persistence of the pathogen Conclusion References