PREDATORS AS A POSSIBLE STRATEGY FOR CONTROLLING A XYLELLA EPIDEMIC?

. In Southern Italy, since 2013, there has been an ongoing Olive Quick Decline Syndrome (OQDS) outbreak, due to the bacterium Xylella fastidiosa , which has caused a dramatic impact from both socio-economic and environmental points of view. Current agronomic practices are mainly based on uprooting the sick olive trees and their surrounding ones, with later installment of olive cultivars more resistant to the bacterium infection. Unfortunately, both of these practices are having an undesirable impact on the environment and on the economy. Here, a spatially structured mathematical model has been proposed to include a predator Zelus renardii as a possible biocontrol agent of the Xylella epidemic. The fact that Z. renardii has been reported to be a generalist predator implies that its introduction is not an eﬃcient control strategy to eradicate a Xylella epidemic. Instead, a specialist predator, whenever identiﬁed, would lead to the eventual eradication of a Xylella epidemic. In either cases it has been conﬁrmed that a signiﬁcant reduction of the weed biomass can lead to the eradication of the vector population, hence of a Xylella epidemic, independently of the presence of predators.


Introduction
The etiological agent of the olive quick decline syndrome (OQDS), is the vector borne bacterial disease Xylella fastidiosa.
The main vector of X. fastidiosa in Southern Italy has been identified in the so called meadow spittlebug, i.e.Philaenus spumarius (Hemiptera, Aphrophoridae), a xylem sap-feeding specialist.Their juvenile forms (nymphs) develop on weeds or ornamental plants, confined in a foam shelter produced by themselves for protection from predators and temperature, while their adult forms move to the tree canopy [1,10,12,19,20].
Once a plant is infected, bacteria multiply within the xylem vessels and the flux of water results inhibited by a bacterial microfilm.Typical symptoms are leaf scorch, dieback of twigs, branches and even of the whole plant [9,21].Possible sanitation of infected olive trees has not been validated yet.Additional details have been reported in [3].
In a series of papers by the same research team [3,5,7], motivated by the outbreak of the Olive Quick Decline Syndrome (OQDS) occurring in Southern Italy, an epidemiological model has been proposed describing the evolution of the main three players, i.e. the insect vector, P. spumarius, the olive trees, and weeds.Since the feeding behaviour and metabolic processes are qualitatively similar for both nymphs and adults [13], we have considered only one stage of active vectors of the infection.
A detailed mathematical analysis and related numerical simulations, reported in the above mentioned papers, have shown the key role played by the removal of weeds (food resources for juvenile insects) from olive orchards and surrounding areas.In addition, as expected, the adoption of more resistant olive tree cultivars has been shown to be a good strategy, though less cost effective, in controlling the pathogen.
In recent papers [14,16] Zelus renardii (Hemiptera, Reduviiidae) has been identified as a predator of P. spumarius for a possible control of a Xylella epidemic.This fact has motivated an extension of the previous spatially structured mathematical model, presented in [3], to include the dynamical behavior of the predator.The fact that Z. renardii has been reported to be a generalist predator implies the choice of an Holling type III functional response of predation in the mathematical model.As a consequence, it has been shown that the introduction of Z. renardii as a predator of P. spumarius is not an efficient control strategy to eradicate a Xylella epidemic.Instead, a specialist predator or of a parasitoid, whenever identified, would lead to the eventual eradication of a Xylella epidemic; as a matter of fact, in this case the appropriate choice for the predation functional response would be an Holling type II.
In either cases it has been confirmed that a significant reduction of the weed biomass can lead to the eradication of the vector population, independently of the presence of predators.A relevant contribution of our approach consists of the restriction of measures of intervention (control) only to a subregion of the whole habitat of interest.
All of the above has been illustrated by a set of computational experiments, within a variety of different possible parameter scenarios.
A warning has emerged concerning weed cut in the relevant habitat: it may negatively affect the efficiency of a predator, due an expected modification of the ecosystem.
On the other hand the infusion of an alien predator may have a catastrophic impact on the relevant ecosystem (see e.g.https://www.cno.it/zelus-renardii-il-suo-impiego-preoccupa-gli-scienziati).
This work may suggest possible model-driven real experiments for providing accurate estimates of all parameters of the model.Once these would be made available, adequate computational experiments may provide predictions on the time behavior of the relevant agroecosystem, under various experimental scenarios.
Our final goal is to analyze optimal control problems which may possibly lead to the identification of an optimal choice of the control parameters and possibly an optimal subregion of intervention.This would require the definition of suitable cost functionals, including all active and passive costs of possible agronomic practices.This is left to future investigations.The paper is organized as follows.In Section 2 the mathematical model is presented.In Section 3 a possible bio-control strategy has been presented, based on the inflow of a predator in the relevant ecosystem.In Section 3.1 the impact of a specialist predator has been analyzed, while the impact of a generalist predator has been analyzed in Section 3.2.The mathematical difference between the two cases is represented by a different predation functional response: Holling type II for a specialist predator; Holling type III for a generalist predator.
In Section 4 numerical simulations are presented which confirm the analytical results.In the numerical simulations the relevant parameters have been taken from [7].Section 5 contains some concluding remarks.A comparison result for parabolic equations is presented in the Appendix.

The mathematical model
We will consider a spatially structured model including the population of vector insects (Philaenus spumarius), the population of predator insects (Zelus renardii) and the population of infected olive trees.
The relevant ecoagronomical habitat will be described by a spatial domain Ω ⊂ R 2 .
As far as the vector insect population is concerned, we will denote by s 1 (x, t) the spatial density of susceptible individuals, and by i 1 (x, t) the spatial density of infected individuals.The spatial density of the total population of P. spumarius will then be In absence of predators, in [3] the evolution dynamics of s 1 (x, t) and i 1 (x, t) has been assumed to be described by the following equations, respectively, subject to suitable initial and boundary conditions (parameters are described in Tab.1): for x ∈ Ω and t > 0.
As a consequence, the evolution equation for the total insect population C I (x, t), in absence of predators, is given by We have assumed that the carrying capacity of all insects depends upon the environmental parameter M (x), x ∈ Ω, that had been already introduced in [3,5,7] to describe the biomass of all other plants, such as bushes, ornamental plants, etc. that we will simply call weeds, either healthy or infected; they are required for the proper development of all insects in the relevant ecosystem.
For the explanation of all other terms in the above equations we refer to [3] (see also [2,8], and references therein).
• PREDATORS As far as the population of predators is concerned, it is denoted by Z(x, t); its evolutionary dynamics, in absence of a population of P. spumarius, may be described as follows (parameters are described in Tab.1): for x ∈ Ω, t > 0, subject to suitable boundary and initial conditions.We have assumed that the carrying capacity of the predators, in addition to a general parameter K(x), depends upon the additional environmental parameter M (x), x ∈ Ω, as mentioned above.
In order to describe mathematically the contribution of P. spumarius to the dynamics of the predator population, we have to take into account that according to what reported in [14], Z. renardii is a generalist predator, and its efficiency in the predation P. spumarius is rather low.In [22, p. 83] (see also [4], [15, p. 38]) it was evidenced that the functional response g(C I ) of a generalist predator to a specific prey follows an S-shaped curve (Holling type III), according to which the predator dedicates its predation only to large quantities of any specific prey C I .In our case this means that in equation (2.4) as additional term describing predation it should with γ 1 , γ 2 > 0.
Remark 2.1.We may remark that in case the relevant predator acts as a specialist one for the P. spumarius, the functional response of predation may be taken of the Holling type II form (see [22]) as follows with γ 1 , γ 2 > 0.
Unfortunately, to our knowledge specialist predators for the P. spumarius have not been identified yet.Anyhow, for a matter of completeness, in this paper both cases have been analyzed.
As a consequence of the above, if we include the impact of predation in equation ( 2.3) we have where g is either g II or g III introduced in (2.5) and (2.6).
Then equations (2.1), (2.2) can be rewritten as follows (2.8) • OLIVE TREES For the olive trees it is better to refer to their canopies, so that we may consider pruning and regrowth.Healthy trees (canopy) are produced at constant net regrowth rate q, get infected by contact with infected insects at rate ζ, or by human activities such as budding and grafting, at rate b.For trees, in view of their long survival, we can neglect natural mortality.Infected trees experience disease-related mortality µ and human-induced mortality due to pruning and logging.
We shall denote by s 2 (x, t) the spatial density of the biomass of healthy trees, and by i 2 (x, t) the spatial density of the biomass of infected trees.As a consequence if C T (x, t) denotes the spatial density of the total biomass of olive trees.

Control, eradicability
Our aim is to analyze the behavior of the insect vector P. spumarius in presence of a predator, which may act as a possible control strategy for the eradication of a Xylella epidemic in an ecosystem.For the habitat of the relevant ecosystem, we take a bounded domain Ω ⊂ R 2 , with a sufficiently smooth boundary ∂Ω.
In the evolution equation for predators we introduce an additional term u(x) ≥ 0, possibly acting only in a suitable subregion ω ⊂ Ω.This control term may be understood as a continuous infusion of predators at position x ∈ ω ⊂ Ω.When we do not specify otherwise we shall take u ∈ L ∞ (ω), u(x) > 0 a.e.x ∈ ω.
Actually the control might also be time dependent, but for the time being we will take it constant in time.
Hence the dynamics of the predator population would be described by Denote Q = Ω × (0, +∞) and Σ = ∂Ω × (0, +∞).Let us consider only the populations of insects, C I for the P. spumarius, and Z for its predator, for a general functional response of predation: Via Banach's fixed point theorem and taking into account the comparison principle for parabolic equations (see the Appendix), we may state the following theorem.
x ∈ ω, we get that (3.2) has a unique nonnegative solution.
Consider now the following eigenvalue problem.
By Rayleigh's principle [6, p. 516] we have that λ1 = inf is the principal eigenvalue for the above mentioned problem to which we may associate a positive eigenfunction φ1 .It is not difficult to show that, if α is a positive real number, the function is the unique solution to the following linear parabolic problem: Now, if we take α a sufficiently large constant, by comparison results for parabolic equations, we may state that As a consequence, we get that C I (•, t) −→ 0 in L ∞ (Ω), as t → +∞.
On the other hand, if λ1 = 0, then there exist two positive constants α, α such that As a consequence, we get that C I (•, t) −→ 0 in L ∞ (Ω), as t → +∞.
We may then state the following theorem.
Theorem 3.2.If λ1 ≥ 0, then for any u ∈ L ∞ (ω), u(x) ≥ 0 a.e.x ∈ ω, we get that This result expresses, in mathematical terms, the fact that for λ1 ≥ 0 the total population of vector insects goes extinct, independently of the presence of predators in the ecosystem.
It is clear that in additional presence of predators, C I (t) will decay to 0 even faster.We may further notice that, according to (3.3) we may make λ1 ≥ 0 by reducing M, the weed biomass, which confirms the results of our previous papers [3,5,7].
Let now λ1 < 0. In this case, in absence of predators we cannot claim that in general the vector population tends to zero, for large times.Additional investigations are required.In particular the qualitative behavior of System (3.2) is different in the two cases of specialist and generalist predators.

Specialist predators
For a specialist predator System (3.2) becomes Consider the following nonlinear parabolic problem Since u(x) > 0 a.e.x ∈ ω, we get that for any t 0 > 0, there exists m t0 > 0 such that Let us prove the following result Proof.Since u(x) > 0 a.e.x ∈ ω, we get that for any t 0 > 0, there exists m t0 > 0 such that Consider a positive t 0 (fixed).Let Z10 , Z 10 ∈ (0, +∞), Z10 sufficiently small and Z 10 sufficiently large.The solutions Z1 and Z 1 to We may infer that as t → +∞, and that 0 < Z * 1 (x) ≤ Z * 2 (x) a.e.x ∈ Ω.Moreover, we may infer that Z * i are positive solutions to (3.5).We immediately get that and so we conclude that Z * 1 (x) = Z * 2 (x) a.e.x ∈ Ω.So, we may infer that there exists a unique positive solution to (3.5).A regularity result for the solutions to the elliptic equations implies that Z * ∈ H 2 (Ω) ⊂ C(Ω).

It follows that the principal eigenvalue λ
is positive.On the other hand, by comparison we get that 0 ≤ Z * (x) − Z1 (x, t) ≤ w(x, t) a.e.x ∈ Ω, ∀t ≥ t 0 , where w is the solution to (at the rate of exp(−λ * 1 t)), and so Analogously, we get that We may then conclude that as t → +∞.It follows that, for any ε > 0, there exists t(ε) ≥ 0 such that On the other hand it can be shown that ∀ε > 0, ∃t 2 (ε) ≥ 0, such that Hence, for any ε > 0 sufficiently small we have that If we consider C ε the solution to then the comparison principle for parabolic equations implies that If the principal eigenvalue λ 1 (ε) to is positive, then as t → +∞, and so On the other hand, by Rayleigh's principle we get that lim To summarize, we get the following result Theorem 3.4.If λ1 < 0 and if λ 1 (0) > 0, then Remark 3.5.According to Theorem 3.2, if λ1 ≥ 0 (which may occur for a small value of M, the weed biomass) the vector population vanishes for large times, independently of the presence of predators.While, according to Theorem 3.4, for λ1 < 0 eradication of the vector population is possible upon the additional information that λ 1 (0) > 0.
The form of λ 1 (0) given by Rayleigh's formula (3.6) shows that there is a monotonicity with respect to Z * ; on the other hand Z * is monotonically increasing with respect to I ω u (and so in particular to ω, by inclusion).
This means that, if the inflow of specialist predators in the relevant habitat is sufficiently large (or we extend the region of intervention ω), eradication of the total vector population is possible.
We may additionally state the following.

Generalist predators
In this case system (3.2) has to be changed into the following one, in which an Holling type III response to predation occurs.In this case (C I , Z) is the solution to We are considering the case in which λ1 < 0, as in Theorem 3.4 concerning specialist predators.In that case we have been able to show that, under the additional assumption of a large inflow of predators in the relevant habitat, the vector population will eventually go extinct.For generalist predators the result is negative, as we will show in the following theorem.
Let us notice that the comparison principle for parabolic equations implies that 0 ≤ Z(x, t) ≤ Z(x, t) a.e.x ∈ Ω, ∀t ≥ 0, where Z is the unique solution to We have that Z ∈ L ∞ (Q) and so there exists M 0 , a positive constant such that defined on [0, +∞)) is bounded.We may conclude (via the comparison principle) that where C I1 is the unique solution to Since C I (•, t) −→ 0 in L ∞ (Ω), as t → 0, we get that λ1 + pεM 0 , the principal eigenvalue to is positive.Since λ1 + pεM 0 > 0 for any ε > 0, we conclude that λ1 ≥ 0, which leads to a contradiction.So, we may conclude that there is no control u such that Let us evaluate, however, how much the population density C I decays if a control u is implemented.
The comparison principle implies that where C I is the unique solution to Hence, C I cannot decay below C I .On the other hand, it can be shown that as t → +∞, where C * is the positive solution to Hence, we may conclude that for any ε > 0, there exists a T (ε) ≥ 0 such that The above results, in particular Theorem 3.2 has confirmed the peculiar role of the parameter M, representing the weed biomass existing in the relevant habitat: independently of the introduction of predators, cutting weeds has been confirmed as a cost effective practice for controlling a Xylella epidemic in olive orchards (see the outcomes of [3,5,7]).
It would be of great interest, from an agronomical point of view, to relate the mathematical threshold parameters, λ1 and λ 1 (0), to some relevant parameters of intervention, such as M (x) (the spatial distribution of the weed biomass) and u(x) (the spatial distribution of the predators inflow).Unfortunately this is a nontrivial mathematical problem, from a purely analytical point of view.A practical response to this issue can be obtained via adequate computational experiments, as it has been reported in Section 4.
On the other hand, the values of λ1 and λ 1 (0) depend, in addition to the possible parameters of intervention M (x) and u(x), upon intrinsic parameters of the relevant agro-ecosystem, which are not yet known numerically.As mentioned in the concluding remarks, it is hopeful that our results may stimulate additional field experiments aimed at the identification of real values of all relevant parameters.

Numerical simulations
The numerical strategy adopted to approximate the controlled system (Eqs.(2.7)-(3.1))consists of the finite element method for space discretization and the finite difference method for time discretization.This procedure is the state-of-the-art for the solution of parabolic Partial Differential Equations (PDEs); see e.g.[18].
Space discretization.We first apply a standard Galerkin procedure to the weak formulations of the controlled system (Eqs.(2.7)-(3.1)).As computational domain Ω we have taken a two-dimensional slab of size 400 × 80 km 2 , which mimics e.g. the whole region of Apulia in Southern Italy, from South (right hand side of the domain) to North (left hand side of the domain).The rectangular domain has been discretized by a uniform grid of 200 × 40 bi-linear finite elements (Q1), yielding a total amount of 8241 discretization nodes.The stiffness matrix is computed exactly, whereas the mass matrix is obtained by applying the mass lumping technique.
Time discretization.After space discretization by finite elements, we obtain the semi-discrete problem that consists of a system of ordinary differential equations (ODEs).We solve these ODEs by employing a first order semi-implicit finite difference scheme, where the linear diffusion terms are approximated by Backward Euler, whereas the non-linear reaction terms are approximated by Forward Euler.As a result, at each time step it is required the solution of a linear system of algebraic equations of dimension 5 × 8241 = 41, 205 degrees of freedom.The linear system is solved by Gaussian elimination with the built-in function of Matlab.The time step size is ∆t = 0.0002 years.
For further details on the numerical discretization of parabolic problems we refer e.g. to [18].

Parameter calibration
In the computational experiments, we have considered the following five scenarios: The predation inflow u has been applied only in the region ω = {x 1 > 250}.
The initial conditions are set as: Coherently with our mathematical analysis, in both cases A1 and A2, in which the level of weed biomass has been kept sufficiently small (M=0.4), the epidemic is stopped, independently of the specificity of the predators, either specialist or generalist, i.e. either Holling II or Holling III response functional, respectively.See Figure 1 for the case A1, while the case A2 is not reported because it is analogous to case A1.
In cases B11 and B12, Holling II has been adopted as predation response functional.They differ only by the level of predation inflow u : smaller in case B11 with respect to case B12.In case B11 the epidemic propagates, see Figure 2; in case B12 the epidemic is stopped (not shown, analogous to case A1, see Fig. 1).
This result highlights the fact that even in the case of predators with a narrow range of possible preys, or in the case of parasitoids, the amount of released control agents represents a crucial aspect that has to be taken into account when planning a control strategy based on their release in the environment.Here it is confirmed that a release of control agents below a certain threshold frustrates the possibility of success of the control strategy.
Given accurate experimental values of all parameters of the mathematical model of the relevant ecosystem, a suitable iteration of the computational experiments B11 and B12 may lead to the identification of the above mentioned threshold level for u; it is clear that such threshold value will depend upon the selected region of intervention ω, too.
In last case B22, Holling III has been adopted as predation response functional: even for a "large" predator inflow u the epidemic persists, see Figure 3.

Theorem 3 . 3 .
The parabolic equation (3.4) admits a unique positive steady state Z * , i.e. a unique positive solution of the following semilinear elliptic problem

Figure 1 .
Figure 1.Case A1.Spatial distribution of the total insect population C I , and of the infected trees i 2 , at three time frames.The epidemic is stopped.

Figure 2 .
Figure 2. Case B11.Spatial distribution of the total insect population C I , and of the infected trees i 2 , at three time frames.The epidemic propagates along the directions of the red arrows.

Figure 3 .
Figure 3. Case B22.Spatial distribution of the total insect population C I , and of the infected trees i 2 , at three time frames.The epidemic propagates along the directions of the red arrows.