Bi-stable dynamics of a host-pathogen model

—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 ﬁrst 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 sufﬁcient conditions for global asymptotic stability of the pathogen free equilibrium and monotone dynamical systems theory to provide sufﬁcient 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. pathogen compartment. We also show numerically that a bi-stable travelling wave solution can exist between PFE and a stable endemic equilibrium, here EE 2 . We show that the speed c of this traveling wave is of the form aµ b , where µ is the diffusion parameter. Further investigations in order to be able to derive appropriate control strategies to avoid invasion of in the crop.


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 annually [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 temperatures. 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, 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, suggesting 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 infection/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 depends 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 increasing 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. βF S, 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 ρF S. 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  Figure 1. The model is a system of ODEs presented below: (1) Using the notation x = (A, F, S, I) T the model (1) is written asẋ The local existence and uniqueness of solutions of (1) in R 4 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 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 R 4 + , 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 equilibrium at Λ d with basin of attraction S ∈ [0, +∞).
The resulting equilibrium of the model (1) we call Pathogen Free Equilibrium (PFE), that is, The basic reproduction number/ratio, R 0 , is a threshold quantity, which is often used to characterize 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 referred 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 computation of R 0 , 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ẏ The next generation matrix is Thus It is clear that R 0 < 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 asymptotic 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 R 0 is not a threshold quantity at all.
In order to obtain all equilibria, we set right hand side in model (1) equal to zero, From the third equation in (6) we obtain an expression for S in terms of F : Then, using (7), from the second equation in (6) we obtain an expression for A: Similarly, using (7), the fourth equation in (6) gives an expression for I in terms of F : Substituting these expressions into the first equation and excluding the case F = 0, we obtain a cubic equation about F in the form where Clearly, a 1 > 0 and a 4 > 0, while the signs of a 2 and a 3 may vary depending on the values of the parameters. However, it is easy to see that for any signs of a 2 and a 3 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 F 1 and F 2 , with F 1 ≤ F 2 . The respective equilibria of the model (10) are denoted by From the expressions (7), (8) and (9) we see that EE 1 > 0 and EE 2 > 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 The two positive roots of (10) appear simultaneously as a double root F 1 = F 2 , which then splits into two simple roots. Hence, in the bifurcation state when F 1 = F 2 , the equilibrium EE 1 = EE 2 appears and then splits into the two distinct equilibria EE 1 and EE 2 . 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 nonstandard 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 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 asymptotically stable on R 4 + . 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 The second case is when the model has two positive equilibria. In the simulations presented in Figure 3, EE 2 is stable and attracting, while EE 1 is unstable (saddle point). Table II contains the parameter values used for the simulations in Figure 3. The solutions that are initiated below EE 1 in the (A, F, I)-space converge to the PFE. Solutions φ 1 and φ 2 are initiated at EE 1 with altered value of S, S 0 = 2 and S 0 = 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 EE 2 . 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. 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 important 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 Invariance Principle [11,Theorem 2]. (1) is globally asymptotically stable on R 4 + if either condition a) or condition b) below hold: Proof: Taking into account the inequalities (4), it is sufficient to consider the system (1) on the domain We consider on Ω the function where ξ is a positive constant with value yet to be determined. We havė 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  The coefficient of Similarly, the coefficient of F is nonpositive if and only if A constant ξ satisfying (15) and (16) exists if and only if or, equivalently, that is if and only if condition a) holds. Hence, if a) holds, we can select ξ such thatV (x) ≤ 0 for x ∈ Ω. Let E = {x ∈ Ω :V (x) = 0}. According to LaSalle Invariance Principle, all solutions converge 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, The largest invariant set in E given in (18) 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)  b) The condition β < α+d provides for slightly sharper upper bound of SI and hence slightly weaker condition on the parameters.
Hence, in the investigation of the asymptotic behavior of the system we can consider only the subset of Ω where .
The inequality β < α + d implies that Hence, 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  Figure 5.

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 Then, it follows that for every solution of (1) it holds Hence, for the asymptotic properties of the solutions of (1) it is sufficient to consider the subset of Ω where In this subdomain we consider the following system for y = (A, F, I) T : Let us recall that function h is said to satisfy the Kamke condition if h i is increasing in y j for i = 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 Since the nondiagonal entries of J h are nonnegative, the system (21) is monotone. Moreover, since the system is irreducible in the interior of R 3 + , 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  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: 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 multiply the third equation by γσ α+d and add it to the second one to eliminate I. Hence, we obtain Equation (26) is equivalent to the quadratic equation 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 Assuming that conditions (28)-(29) hold, we denote byF 1 andF 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 Proof: It is easy to see that the equilibrium 0 is asymptotically stable, indeed the eigenvalues of J h (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 initiated 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 E 1 after some simplifications is: Since the nondiagonal entries of J h (Ẽ 1 ) are nonnegative and the matrix is irreducible, it follows from the theory of nonnegative matrices [4, Theorem 2.1.4] that J h (Ẽ 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 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 .
Theorem 3. Let conditions (28) and (29) hold. Then, for any solution of (1) we have that if 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 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 condition at t = 0. Then, Theorem 2 implies that E 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 abovẽ A 1 andĨ 1 , while F is still aboveF 1 , the pathogen will persist eventually at least at a level ofẼ 2 .
Theorem 3 is illustrated numerically by Figure 6. The initial conditions of the solutions in this figure were chosen specifically so that (A 0 , F 0 , I 0 ) ≥Ẽ 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 EE 2 , even though the initial conditions do not satisfy the requirements of Theorem 3. Our numerical investigations of model (1) have revealed that in addition to the PFE which is always locally attracting, there exists an equilibrium close to the PFE which is repelling, and an attracting equilibrium which is removed from the PFE. 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: 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  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 discretization coupled with a second order centralspace 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 spatiotemporal 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. 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 infection: A monotone system, constructed to approximate the temporal host pathogen model from below was proven to admit two interior equilibria    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. Table VIII on page 15,  these conditions are satisfied. Indeed, we

For the parameter values in
The equilibrium,Ẽ 1 of the lower approximating system is:Ẽ 1 ≈ (4.091, 1.7241, 0.2164). The persistence of pathogen is illustrated in Figure 10. In fact, the solutions converge to EE 2 , although the stability properties of EE 1 and EE 2 have not been proven. How does the inclusion of diffusion affect this phenomenon? Solutions initiated at the level of EE 2 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. To this end, we investigate the relationship between µ and the wave speed c. The parameter values in Table VIII were again used, and a solution with (A 0 , F 0 , I 0 ) taking the value of EE 2 on the left boundary was considered. The diffusion constant µ was taken to be in the interval [10 −7 , 10 −1 ], which results in c ∈ (0, 4.5 × 10 −3 ]. An equation 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 r 2 = 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 sufficient conditions for PFE being globally asymptotically stable and for persistence of the pathogen, using two different approaches, LaSalle Invariance Principle approach and monotone system approach. 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 EE 2 . 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.