Analysis of a Mathematical Model Arising in Plant Disease Epidemiology

In this work we study from the mathematical and numerical point of view a problem arising in vector-borne plant diseases. The model is written as a nonlinear system composed of a parabolic partial differential equation for the vector abundance function and a first-order ordinary differential equation for the plant health function. An existence and uniqueness result is proved using backward finite differences, uniform estimates and passing to the limit. The regularity of the solution is also obtained. Then, using the finite element method and the implicit Euler scheme, fully discrete approximations are introduced. A discrete stability property and a main a priori error estimates result are proved using a discrete version of Gronwall’s lemma and some estimates on the different approaches. Finally, some numerical results, in one and two dimensions, are presented to demonstrate the accuracy of the approximation and the behaviour of the solution.

1. to investigate the relationships existing between the abiotic and biotic disease system components influencing the disease growth in the host plant (see, e.g., [11,18,25,28]). 2. to quantify the temporal and the spatial spread of the disease [13,22]. 3. to assess the efficacy of management options for disease and vector control [10,17]. 4. to develop scenario analysis for the evaluation of different disease management policy (see, for instance, [2,7,12,19,37]).
Furthermore, epidemiological models development and use can identify knowledge or data gap and hence prioritize further research effort. This topic has received an increasing interest and the number of publications studying different mathematical issues is really large. We can highlight the works of Hebert and Allen [14], where vector aggregation and vector dispersal were simulated by using deterministic and stochastic models, Jackson and Chen-Charpentier [16], who considered an ODE system to model the virus propagation, Neofytou et al. [27], who proposed a mathematical model for the interactions between two competing viruses with particular account for the RNA interface, Shi et al [33], where the asymptotic behavior of an epidemic model describing a vector-borne plant disease was studied by using a Lyapunov technique, Zhao et al. [40], who formulated a stochastic plant infectious disease model with impulsive toxicant input by using Markov conversion, and Zhao and Xiao [41], where a method to eradicate plant disease or to maintain the number of infected plants below the economic threshold was studied. We can also refer the recent works [5,8,15,21,23,24,26,29,31,38,39,42] and the references cited therein.
In the present contribution, we propose a general eco-epidemiological model describing the spatio-temporal dynamics of vector borne plant diseases, which constitutes an important extension of the model considered in [16]. It provides a biological meaningful representation of the disease growth in the host plant, the pathogen transmission from the vector to the plant, the vector acquisition of the pathogen from the infected plant, the demography of the vector population, and the susceptibility of the host plants. The model considers the influence of environmental drivers on disease dynamics through the identification of periods of the year characterised by specific rates of disease growth and vector disease transmission. Future model developments may include a time-dependent definition of environmental drivers. The model can consider any possible spatial arrangement of the host plants and any possible plant communities composed by plants of diverse susceptibility to the disease. The model can also describe the metapopulational epidemiology, to account for a patchy distribution of susceptible host plants.
In the paper we show the well-posedness of a nonlinear evolutionary problem (i.e., the eco-epidemiological model) coupling a parabolic PDE for the vector abundance function with an ODE equation for the disease dynamics in the plant. The resulting system is highly nonlinear, so in order to prove existence, uniqueness and regularity of the strong solution, we introduce a backward finite differences scheme and argue on this by proving uniform estimates and passing to the limit on the time step. Parabolic PDEs systems and the related optimal control problems have been widely studied in the literature: we refer, without any sake of completeness, to [1,20,34] and references therein for the analysis of optimal control problems. Then, we introduce a fully discrete approximation of the corresponding variational formulation by using the classical finite element method and the well-known Euler implicit scheme. A discrete stability property and a priori error estimates are proved by using a discrete version of Gronwall's lemma and some rather technical estimates. Finally, we also perform numerical simulations in one and two dimensions to demonstrate the accuracy of the approximation and the behaviour of the solution.

The Mathematical Problem: Existence and Uniqueness
In this section, we present a brief description of the model and we obtain its mathematical formulations.
Let Ω ⊂ R d , d = 1, 2, 3, be a bounded domain and denote by [0, T ], T > 0, the time interval of interest. The boundary of the body Γ = ∂Ω is assumed to be C 1,1 at least, with outward unit normal vector ν = (ν i ) d i=1 . Moreover, let x ∈ Ω and t ∈ [0, T ] be the spatial and time variables, respectively. In order to simplify the writing, we do not indicate the dependence of the functions on x = (x j ) d j=1 and t, and a subscript after a comma under a variable represents its spatial derivative with respect to the prescribed variable, i.e.
The time derivative is represented as a point over each variable. Finally, as usual the repeated index notation is used for the summation. We denote by ϕ the health status of the plant, requiring that ϕ varies between 0 and 1. ϕ = 0 refers to a non-infected plant, while for each value greater than zero the plant is infected. The disease severity increases with the increase of ϕ, up to a maximum value of ϕ = 1. The function γ describes at each point the abundance of infected vectors; hence, it is confined between 0 and the local vector population carrying capacity κ.
Therefore, the problem modelling the evolution of an infected plant community is the following.
Problem P. Find the health function ϕ : Ω × [0, T ] → R and the abundance function γ : Ω × [0, T ] → R such thaṫ In the previous problem, a is a positive scalar parameter for the diffusive Laplacian operator acting on γ , and it is related to the dispersal behaviour of the vector, b is the pathogen acquisition rate of the vector, M describes the mortality rate of the vector population due to natural mortality, F describes the progression of the disease in the infected plant and S is a function describing the host plant susceptibility in each point of the domain Ω. Finally, γ 0 and ϕ 0 represent the initial abundance of vectors and the initial diffusion of the disease in Ω, respectively. In order to obtain the variational formulation of Problem P in the next section and to prove the existence and uniqueness result, , and denote by (·, ·) Y , (·, ·) H and (·, ·) V the respective scalar products in these spaces, with corresponding norms · Y , · H and · V . Moreover, let us define the variational space: We state the following assumptions on the problem data: Theorem 1 Assume (5)- (7). Then Problem P stated in (1)-(4) has a unique solution (γ , ϕ) with the following regularity: Moreover, if γ 0,i and ϕ 0,i , i = 1, 2, are given as in (7) and (γ i , ϕ i ), i = 1, 2, are the corresponding solutions to Problem P, then the estimate holds true for a positive constant C > 0 which depends only on meas(Ω), on the final time T , on the shape of the nonlinearities and on the constants and the norms of the functions involved in Assumptions (5)- (7).
Proof First, we will consider an approximation of Problem P. In order to introduce a backward finite differences scheme, we assume that N is a positive integer and Z a normed space. By fixing the time step τ = T /N , we introduce the interpolation maps from Z N +1 into either L ∞ (0, T ; Z ) or W 1,∞ (0, T ; Z ). For (z 0 , z 1 , . . . , z N ) ∈ Z N +1 , we define the piecewise constant functions z τ and z τ , and the piecewise linear functions z τ as follows: if 0 < s < 1 and i = 0, . . . , N − 1. By a direct computation, we have For i = 1, . . . , N , we set and we look for a couple (γ τ , ϕ τ ) satisfying at least the regularity requirements: and solving the following approximation problem that we denote by Problem P τ : We observe that Problem P τ has a unique solution (γ i , ϕ i ). Indeed, from (19) we infer a linear elliptic equation with homogeneous Neumann boundary conditions (21), namely In the following (see (41)), we will prove some uniform estimates for γ i and ϕ i , which ensure that the right hand-side of Eq. (23) is bounded in L ∞ (Ω) and thus in Y . Hence, problem (23)-(24) admits a unique solution γ i ∈ W (see, e.g., [3, Section 9.6-Theorem 9.26, p. 299]), while ϕ i can be explicitly obtained from (20) as Since γ i ∈ W and Assumptions (5)-(7) hold, it immediately follows that the right hand-side of Eq. (25) belongs to L ∞ (Ω), whence ϕ i ∈ L ∞ (Ω). In order to prove that ϕ i ∈ H 1 (Ω) ∩ L ∞ (Ω), we check that the gradient of the right hand-side of (25) belongs to H . Let us focus, e.g., to the not trivial term: Since it follows that the right hand-side of Eq. (26) belongs to H . The gradients of the other terms on the right hand-side of equation (25) can be treated in a similar way. Hence, Now, let us rewrite the discrete Eqs. (19)-(22) by using the piecewise constant and piecewise linear functions defined in (12): Now, we will obtain some a priori error estimates. In the rest of the paper, we will use the Hölder and Young inequalities; that is, for every x, y > 0, α ∈ (0, 1) and δ > 0 there hold x y ≤ δx 2 + 1 4δ Without loss of generality, we assume that According to the bounds of the initial data stated in (7), if rewriting Eq. (20) and using (5)- (7) and (33)- (34), then the following estimates are obtained: Therefore, we conclude that, for every i ≥ 1, Moreover, rewriting discrete Eq. (19) we find that In what follows, we will use the classical notation Due to (5)- (7), (34) and (35), multiplying Eq. (36) by −(γ i ) − and integrating over Ω, we have in Ω, that is, for every i ≥ 1 we have Now, from (19) it follows that Multiplying now Eq. (38) by (γ i −κ) + , integrating over Ω and using (5)-(7), (33)- (34) and (35) we obtain and we conclude that (γ i − κ) + = 0 a.e. in Ω, which leads, for every i ≥ 1, Combining (37) with (39) it follows that and so, we find that Finally, thanks to (35) and (40)- (41), we find that Multiplying Eq. (19) by τ γ i , integrating over Ω and summing up for Due to Assumption (5) the last term on the left-hand side of (44) is nonnegative, while the first term on the left-hand side can be rewritten as follows: Applying Young's inequality (32) to the right-hand side of (44) and using (41), we obtain From now on, without loss of generality, let us assume that τ ≤ 1 4C . Hence, combining (44) with the previous estimates, we have whence, recalling (7) and applying the discrete version of Gronwall's inequality, it leads Combining (41) with (46), and noting that the same estimates hold true both for γ τ and γ τ , we conclude that Moreover, by comparison in (28), we also find that , integrating over Ω and summing up for The second term on the left-hand side of (50) can be rewritten as while, estimating the first term on the right-hand side of (50), and using (5)-(6), (47)-(48) and the above Young's and Hölder inequalities, we have Applying a similar technique to the last term on the right-hand side of (50), we obtain Combining (50) with (51)-(53) and using (47)-(49), it follows that and we conclude that Moreover, noting that from (55) we infer that Finally, by comparison in (27), we also obtain Applying the gradient operator to Eq. (20), multiplying the resulting equation by τ ∇ϕ i , integrating over Ω and summing up for i = 1, . . . , N , we have (59) The left-hand side of (59) can be rewritten as while the first term on the right-hand side of (59) can be bounded as follows Now, we estimate each term of the right-hand side of (61) separately. Applying the Hölder and Young inequalities and using (5)- (7), (41)-(43), (47)-(49), (55)-(58), we find that Finally, estimating the last term on the right-hand side of (59), we have Combining (59) with (60)-(66) we find that From now on, without loss of generality, let us assume that τ ≤ 1 4C . Then, the last term on the right-hand side of (67) can be rewritten as whence we obtain that Due to (7), applying the discrete version of Gronwall's inequality we conclude that Therefore, it follows that Moreover, noting that from (70) we also obtain Finally, from (27), we have thaṫ whence, thanks to (5), (47)-(48), (55)-(58) and (70), we infer that and thus Now, we summarize the above a priori error estimates. According to (41), (47)-(49), (55)-(58) and (70)-(72), we conclude that there exists a positive constant C, independent of τ , such that Now, we will consider the passage to the limit as τ → 0.
Since there exists a unique solution (γ τ , ϕ τ ) to Problem P τ satisfying the regularity requirements (17)-(18), we can pass to the limit as τ → 0 and we can prove that the limit of subsequences of solutions (γ τ , ϕ τ ) to Problem P τ yields a solution (γ , ϕ) to Problem P.
Thanks to (73)-(74) and to the well-known weak or weak* compactness results, we deduce that, at least for a subsequence of τ → 0, there exist four limit functions γ , γ , ϕ and ϕ such that First, we observe that γ = γ . Indeed, thanks to (14) and (75)-(77), we have whence, passing to the limit as τ → 0, we conclude that (γ τ − γ τ ) → 0 strongly in L 2 (0, T ; Y ). Similarly, due to (14) and (78) Now, using (75)-(80) and (81)-(82) and passing to the limit in (27) and (28), we arrive at (1) and (2), respectively. Therefore, we conclude the existence of a solution to Problem P. Finally, we will prove the continuous dependence of the solution to Problem P, from which we will derive the uniqueness of solution.
Thus, the solution to Problem P is unique and it concludes the proof of the theorem.

Fully Discrete Approximations: An a Priori Error Analysis
By using Green's formula and boundary condition (3), we write the variational formulation of Problem P. γ 0 , and, for a.e. t ∈ (0, T ) and for all w, r ∈ V ,

Problem VP. Find the health function ϕ : [0, T ] → V and the abundance function
In the rest of this section, we will consider a fully discrete approximation of Problem V P. This is done in two steps. First, we assume that the domain Ω is polyhedral and we denote by T h a regular triangulation in the sense of [6]. Thus, we construct the finite dimensional space V h ⊂ V given by where P 1 (T r) represents the space of polynomials of degree less or equal to one in the element T r, i.e. the finite element space V h is composed of continuous and piecewise affine functions. Here, h > 0 denotes the spatial discretization parameter. Moreover, we assume that the discrete initial conditions, denoted by ϕ h 0 and γ h 0 are given by where P h is the classical finite element interpolation operator over V h , respectively (see, e.g., [6]). Secondly, we consider a partition of the time interval [0, T ], denoted by 0 = t 0 < t 1 < · · · < t N = T . In this case, we use a uniform partition with step size k = T /N and nodes t n = n k for n = 0, 1, . . . , N . For a continuous function z(t), we use the notation z n = z(t n ) and, for the sequence {z n } N n=0 , we denote by δz n = (z n − z n−1 )/k its corresponding divided differences.
Therefore, using the implicit Euler scheme, the fully discrete approximations are considered as follows.
Problem VP hk . Find the discrete health function ϕ hk = {ϕ hk n } N n=0 ⊂ V h and the discrete abundance function and, for n = 1, . . . , N and for all w h , r where R 1 : R → [0, κ] and R 2 : R → [0, 1] are truncation operators defined as We point out that these truncation operators are introduced for mathematical reasons because we are not able to guarantee that the discrete solutions ϕ hk n and γ hk n are bounded. The arguments used in the continuous case are not valid here, so it remains an open problem even if, in the numerical simulations, we always find that ϕ hk n ∈ [0, κ] and γ hk n ∈ [0, 1]. Moreover, the existence of a unique discrete solution to Problem V P hk is obtained proceeding as in the continuous case.

Remark 1
We note that a semi-implicit scheme has been also checked. In that case, Problem VP hk has the following form: and, for n = 1, . . . , N and for all w This algorithm has been used for obtaining the numerical results of the real cases because the CPU time increases drastically when the implicit scheme is employed. We note that the discrete solutions of both algorithms are rather similar. Moreover, the numerical analysis performed in this section could be extended to the semi-implicit scheme without difficulty. Now, we have the following stability property.
where > 0 is assumed small enough and we used several times Cauchy-Schwarz inequality and Young's inequality (32), we easily find that Now, taking ϕ hk n as a test function in variational Eq. (103) we find that Taking into account that it follows that Combining the previous estimates and using a discrete version of Gronwall's inequality (see, for instance, [4]) we conclude the desired stability property. Now, we will obtain a priori error estimates on the numerical errors ϕ n − ϕ hk n and γ n − γ hk n . Thus, we have the following a priori error estimates result.
where C is again a positive constant which is independent of the discretization parameters h and τ .
Proof First, we obtain the error estimates for the abundance function. Then, we subtract variational Eq. (98) at time t = t n for a test function w = w h ∈ V h ⊂ V and discrete variational Eq. (102) to obtain, for all w h ∈ V h , Therefore, we find that, for all w h ∈ V h , Taking into account that where δγ n = (γ n − γ n−1 )/k and we used again Cauchy-Schwarz inequality and Young's inequality (32), it follows that Now, we obtain the error estimates on the health function. Therefore, subtracting variational Eq. (99) at time t = t n for a test function r = r h ∈ V h ⊂ V and discrete variational Eq. (103), it leads, for all r h ∈ V h , and therefore, we have, for all r h ∈ V h , Keeping in mind that where δϕ n = (ϕ n − ϕ n−1 )/k and we used again Cauchy-Schwarz inequality and Young's inequality (32), we have, for all r h ∈ V h , Thus, combining the previous estimates on the health and abundance functions, we find that, By induction, it follows that, for all Finally, taking into account that applying again a discrete version of Gronwall's inequality (see [4]) we conclude the a priori error estimates.
We note that a priori error estimates (104) can be used to obtain the convergence order of the algorithm. Thus, if we assume that the continuous solution has the additional regularity: then, using classical results on the approximation by finite elements (see [6]), after straightforward estimates we can show that there exists a positive constant, again independent of the discretization parameters h and k, but depending on the continuous solution, such that

Numerical Scheme
As a first step, given the solution ϕ hk n−1 and γ hk n−1 at time t n−1 , the discrete health and the abundance functions at time t n are obtained by solving the following discrete nonlinear system, for all w h ∈ V h and r h ∈ V h , This numerical scheme was implemented on a 3.2 Ghz PC using Matlab and a typical 1D run (h = k = 0.01) took about 0.52 seconds of CPU time meanwhile a typical 2D run (k = 0.01 and 665 nodes) took about 2.64 seconds of CPU time.

First Example: Numerical Convergence in a One-Dimensional Case
We will consider the following academic problem: Problem P ex . We note that Problem P ex corresponds to Problem P with the following data: Since the exact solution to Problem P ex is unknown, we take the solution obtained with parameters h = 1/32768 and k = 10 −5 as the reference solution. Thus, the approximation errors estimated by max 0≤n≤N ϕ n − ϕ hk n Y + γ n − γ hk n Y are presented (multiplied by 10 4 ) in Table 1 for several values of the discretization parameters h and k. Moreover, the evolution of the error depending on the parameter h + k is plotted in Fig. 1. We notice that the convergence of the algorithm is clearly observed although the linear convergence, stated in (105), is not achieved.

Second Example: Dynamics of ' and in a 2D Spatial Domain
In this example, the model (Eqs. (1)-(4)) is used to explore the dynamics of a plant disease transmitted by a vector. In this case, the disease epidemiology is highly dependent on the rate of increase of the pathogen population in the host plant and on the vector efficacy in transmitting the disease. This scenario represents a very common case for most of the plant diseases caused by bacteria or viruses. More in detail, we consider the case of an infection caused by a bacterium and vectored by an insect moving in a two-dimensional spatial domain Ω = (0, 1)×(0, 1). The insect vector both acquires and transmits the pathogen when feeding on the leaves or the green parts of the host plant.
We aim to investigate the role of the parameter F in Eq. (2), the intrinsic rate of increase of the bacterial population on the disease growth in the plant, and to analyse the influence of the plant health status on the processes of vector acquisition and transmission of the infection. Since the importance of parameter F in determining the pattern of local disease dynamics, the exploration of the parameter values is highly relevant for the definition of control strategies of the epidemiology of new foci.
We set the initial condition for the plant disease as a gradient of plant infection, ϕ 0 (x, y) = 0.01(x + 1), and we evaluate how the infection of the vector is triggered by the infected plants over the whole spatial domain. To better represent the starting phase of a new disease outbreak, we consider a low disease growth rate (F = 0.01) and we explore how spatial heterogeneity in the disease influences the process of vector acquisition of the infection. The following data have been employed in the simulations: Using the time discretization parameter k = 0.01, the proportions of infected vector (γ ) and the plant health status (ϕ) functions at final time are plotted on the 2-D domain (see Fig. 2). Both functions γ and ϕ seem to have a similar linear increase according to the initial conditions. Furthermore, the spatial heterogeneity in the plant health status and in proportion of infected vector are strictly related.
The simulated scenario represented in Fig. 2 is a quite common case in plant health. In new outbreaks, the disease is imported through infected plants that become the source of infection for vectors responsible for spreading the disease in wider spatial domain, not previously infected.
The model was implemented over the whole domain but, for clarity of exposition, the results are presented for the point (x i , y i ) = (0.27, 0.52), where we assumed there is a plant and a sub-population (i.e., local population) of vectors feeding on it. We compared the growth dynamics of the disease in the selected plant and the local dynamics of infected vectors, using five different values of the bacteria growth parameters (F = 0.1, 1, 2, 5, 10). We note that, in this case, we have used values greater than 1 for this function but the analysis of the previous sections can be extended straightforwardly. The initial conditions in (x i , y i ) are an infected plant with health status ϕ 0 = 0.01 * 0.52 * (0.52 + 1) and a non-infected local vector population (proportion of infected vectors γ 0 equal to 0). Using the time discretization parameter k = 0.01, the estimated temporal dynamics of plant health ϕ(x i , y i , t) and the proportion of infected vector γ (x i , y i , t) are shown in Fig. 3. Consistent with expectations, as F increases, the level of disease infection in the plant grows faster. The same pattern is followed by the local population of infected vectors. For the greatest value of the parameter F, the quadratic shape of the functions changes into a logistic pattern because ϕ and γ have an asymptotic behaviour at their limit value.
In conclusion, the results of this second example show how the proposed model can simulate complex spatial and temporal epidemiological scenarios, considering also how variability in key biological parameters can influence the onset and the spread of an infection in a previously disease-free area. The exploration of these epidemiological scenarios allows to investigate the spread of the disease in a new outbreak and to support the comparative evaluation of management and eradication strategies.