Controlling the Spatial Spread of a Xylella Epidemic

In a recent paper by one of the authors and collaborators, motivated by the Olive Quick Decline Syndrome (OQDS) outbreak, which has been ongoing in Southern Italy since 2013, a simple epidemiological model describing this epidemic was presented. Beside the bacterium Xylella fastidiosa, the main players considered in the model are its insect vectors, Philaenus spumarius, and the host plants (olive trees and weeds) of the insects and of the bacterium. The model was based on a system of ordinary differential equations, the analysis of which provided interesting results about possible equilibria of the epidemic system and guidelines for its numerical simulations. Although the model presented there was mathematically rather simplified, its analysis has highlighted threshold parameters that could be the target of control strategies within an integrated pest management framework, not requiring the removal of the productive resource represented by the olive trees. Indeed, numerical simulations support the outcomes of the mathematical analysis, according to which the removal of a suitable amount of weed biomass (reservoir of Xylella fastidiosa) from olive orchards and surrounding areas resulted in the most efficient strategy to control the spread of the OQDS. 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 this paper for a more realistic description and a clearer interpretation of the proposed control measures, a spatial structure of the epidemic system has been included, but, in order to keep mathematical technicalities to a minimum, only two players have been described in a dynamical way, trees and insects, while the weed biomass is taken to be a given quantity. The control measures have been introduced only on a subregion of the whole habitat, in order to contain costs of intervention. We show that such a practice can lead to the eradication of an epidemic outbreak. Numerical simulations confirm both the results of the previous paper and the theoretical results of the model with a spatial structure, though subject to regional control only.


Introduction
The etiological agent of olive tree disease known as olive quick decline syndrome (OQDS) is the plant pathogenic bacterium Xylella fastidiosa, which is a vector-borne bacterium.
The main vector of Xylella fastidiosa in Southern Italy has been identified in the so-called meadow spittlebug, i.e., the Philaenus spumarius, a xylem sap-feeding specialist. Their juvenile form (nymphs) develops on weeds or ornamental plants, confined in a foam produced by themselves for protection from predators and temperature, while their adult form moves to olive tree canopies. Experiments have shown a larger infection prevalence of adults on olive trees than on weeds; this fact may lead to the assumption of infection of adults from infected olive trees more than from weeds.
X. fastidiosa transmission is the result of four events [see e. g. Almeida et al. (2005), Redak et al. (2004) and references cited therein]: (a) acquisition from a source plant; (b) attachment and retention to the vector's foregut; (c) inoculation into a new host; (d) development of the infection in a plant after inoculation. A successful transmission also includes bacterial multiplication. Once a plant is infected, bacteria multiply within the xylem vessels inducing the production of a gel in the plant xylem, which occludes the xylem vessels, thus inhibiting the flux of water through the lymph vessels eventually blocking the nutrition of the plant. Typical symptoms are leaf scorch, dieback of twigs, branches and even of the whole plant (Carlucci et al. 2013).
Sanitation of infected olive trees is unfeasible; the scope of our research is the mathematical modeling of the population dynamics of the ecosystem in the presence of infection. The availability of a sound mathematical model may lead to predictive analysis of the relevant populations, so as to suggest possible eradication strategies, or at least possible optimal control strategies.
In a previous paper (Brunetti et al. 2020), based on the outbreak of OQDS in Southern Italy, a model describing the epidemic was presented. It consists of a system of ordinary differential equations (ODEs) describing the evolution of the main three players, i.e., the insect vector, Philaenus spumarius, the olive trees, and weeds. A preliminary mathematical analysis and related numerical simulations have shown that "the removal of a significant amount of weeds (acting as a reservoir for juvenile insects) from olive orchards and surrounding areas has resulted in the most efficient strategy to control the spread of the OQDS. 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." The scope of the present paper is to extend the above results to a spatially structured model and to show in a more rigorous way that the spatial expansion of an OQDS outbreak can indeed be stopped by acting either on the weed biomass or on the choice of the olive cultivar. We have adopted a deterministic reaction-diffusion model. We recall that reaction-diffusion models can be interpreted as mean-field approximations of individual based models, which are more appropriate at the micro scale. Of course, in our case individual behaviors and possible randomness are lost. On the other hand, our approach allows a mathematical "qualitative" analysis of the system, including the derivation of eradicability theorems. Such a "qualitative" analysis has driven on the one side, our numerical experiments and, on the other side, anticipates future investigations on optimal control problems in a variety of scenarios.
In order to keep mathematical technicalities to a minimum, here the weed biomass is taken as a given field parameter and as in Brunetti et al. (2020), the insect life cycle has not been taken into account [more information about this last aspect can be found in Rossini et al. (2020)]. As possible control actions, insect traps, weed cut, choice of the cultivar, and reduction of contact rates have been taken into account. In the numerical simulations, particular attention has been paid to weed cut and choice of the cultivar, the other control measures being more difficult to implement. It is worth mentioning that, in recent investigations presented in Schneider et al. (2020), the authors, by means of a cellular automaton simulator, have confirmed the relevance of the olive cultivar as a possible effective control strategy.
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. Due to diffusion, any point of the domain Ω is strongly connected to any other point of Ω, so that any action taken in a subregion ω ⊂ Ω will eventually propagate to the whole habitat. This is the "leit motiv" of our proposal concerning regional control: "think globally, act locally" [see Aniţa and Capasso (2009)]. This practice may contribute in a significant way to improve the ratio cost/effectiveness of real control strategies. Mathematical analysis and numerical simulations have been carried out showing that it is indeed possible to eradicate the disease by such local action. Our aim is to analyze optimal control problems in future investigations, which may possibly lead to the identification of an optimal choice of the subregion of intervention.
The paper is organized as follows. In Sect. 2, the mathematical model is presented. In Sect. 3, possible control strategies are proposed, based on which an eradicability result is shown. Finally, in Sect. 4, numerical simulations are presented which confirm the analytical results. In the numerical simulations, the relevant parameters have been taken from Brunetti et al. (2020).

Our Model
Since the feeding behavior and metabolic processes are qualitatively similar for both nymphs and adults (Janse and Obradovic 2010), we will consider only one stage of active vectors of the infection. We will consider a spatially structured model which includes the population of vector insects and the population of infected olive trees.
We will denote by C I (x, t) the spatial density of the total population of insects, by s 1 (x, t) the spatial density of susceptible insects, and by i 1 (x, t) the spatial density of infected insects: C T (x, t) will denote the spatial density of the total biomass of olive trees, s 2 (x, t) the spatial density of the biomass of healthy trees, i 2 (x, t) the spatial density of the biomass of infected trees: All the parameters in the model are nonnegative quantities.

• INSECTS
Due to the fact that bacteria reside only in the foregut of an adult insect, the latter generate only healthy offspring [see, e.g., Almeida et al. (2005), Redak et al. (2004, and references cited therein]. The parameter r denotes the birth rate of new insects. Reproduction, however, can occur only if bushes and other plants, whether healthy or infected, are present, explaining the presence of a "carrying capacity" M(x) in the logistic term, describing the biomass of all such plants, which we simply call weeds.
Both healthy and infected insects experience natural mortality at a rate n; χ is a scale parameter.
Finally, we assume that insects may diffuse in the habitat (with constant diffusion coefficient to avoid purely technical complications).
As far as the local incidence of infection for insects, at point x ∈ Ω, and time t ≥ 0, as in previous models (Capasso 1984;Aniţa and Capasso 2009), we assume that it arises due to biting of infected olive trees at any point x ∈ Ω of the habitat, within a spatial neighborhood of x represented by a suitable probability kernel k(x, x ), depending on the specific structure of the local ecosystem [see also Shcherbacheva et al. (2018)]; as a trivial simplification, one may assume k(x, ·) as a Gaussian density centered at x; hence, the "local incidence rate (i.r.)" for insects, at point x ∈ Ω, and time t ≥ 0, is taken as Hence, the spatial dynamics of susceptible insects is expressed by the following equation while the spatial dynamics of infective insects is expressed by the following equation  (1) and (2) act in a spatial domain Ω ⊂ R 2 , at times t ∈ (0, +∞). They are subject to suitable initial conditions, and we assume homogeneous Neumann boundary conditions.

• 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.
Hence, the spatial dynamics of trees is expressed by the following equations Both Eqs.

Control and Eradicability
By summing Eqs. (1) and (2), we obtain an equation for the total insect population C I (x, t) For M constant, pure Neumann boundary conditions allow constant in space and time solutions, satisfying the following equilibrium equation

This equation may admit two solutions
This shows that for M ≤ n r , no other nonnegative solutions are feasible, but the trivial one. Otherwise, another equilibrium is feasible provided M > n r .
From the above simple reasoning, we may conjecture that a way to eradicate the disease is to eliminate the insect population by a significant reduction of the carrying capacity M, which may be obtained by eliminating weeds in the relevant olive orchards. Equation (5) shows the quantitative role of the scale parameter χ ; a smaller value of χ allows a larger value of C * * I , so that we may say that M χ plays the role of an "effective carrying capacity" of insects. Numerical simulations, reported in Sect. 4, illustrate these facts. A more accurate analysis follows. Assume that certain constant controls are considered in a non-empty open subset ω ⊂ Ω.
Possible regional controls (acting in ω) are the following: C1: Traps : n → n + γ 1 ; increase insect death rate by insecticides and/or treated nets.
For various reasons, including cost reduction, we may consider the possibility of implementing the proposed control measures only on a suitable subregion ω ⊂ Ω of the whole habitat, including areas already affected by the epidemic, augmented by a preventive confinement area.
Let I ω denote the characteristic function of ω. Then, the controlled system is We assume that there is no flux of insects through the boundary of Ω (the domain is isolated): and that the following initial conditions are satisfied Assumptions: . We postpone, for the time being, the proof of the existence and uniqueness of such a solution.
As far as the trees are concerned, we may easily obtain an upper bound for the total canopy.
If we now add Eqs. (8) and (9) and take into account the initial conditions, we find that C T = s 2 + i 2 satisfies If we denote byC T the solution to we may state that (here, and throughout this paper, · p denotes the usual norm in L p (Ω)); it is obvious thatC T is space independent.
To prove this inequality, we set y( Since y 0 and g are nonnegative, we have that y(x, t) ≥ 0, a.e. x ∈ Ω, ∀t ≥ 0, and the inequality (14) follows.
As a consequence of (13) We note that the quantity C(q − ) + is the carrying capacity of the tree population.
We may now also note that if we denote by C T the solution to by similar arguments as for the previous inequality, we may also state that It is evident that C T is space independent too. Altogether we may state For the insect population, we may add Eqs. (6) and (7) and take into account the boundary and initial conditions, to obtain that (15) By Banach's fixed point theorem and using the comparison result for the solutions to the linear parabolic equations [see Friedman (1964) and Protter and Weinberger (1994)], the existence and uniqueness of a nonnegative solution to (15) follow. Moreover, the solution satisfies where C I is the solution to the linear parabolic equation Turning back to (6)-(11), it follows via Banach's fixed point theorem (and using the comparison of the solutions to the linear parabolic equations and the comparison of the solutions to the linear first order ODEs) that there exists a unique solution (s 1 , i 1 , s 2 , i 2 ) such that j ∈ {1, 2} [for such an argument see Aniţa andCapasso (2009, 2012) and ]. Let C * I be the maximal nonnegative solution to and letλ 1 be the principal eigenvalue for the problem It is obvious thatλ 1 depends on Ω, ω, and the controls γ 1 and γ 21 , and is an increasing function of γ 1 , γ 21 , and ω-by inclusion-via Rayleigh's principle.
By the same methods as in , the following proposition can be shown to hold.
We note that Proposition 1 confirms for the full reaction-diffusion system, the outcomes of the preliminary analysis presented at the beginning of Sect. 3. In particular, it is of interest to observe that, for small values of the weed distribution, M(x), or a large value of the control parameter, γ 21 , or a large domain of intervention, ω,λ 1 may become greater than or equal to zero, so that the only nonnegative solution of (16) is the trivial one, and the whole insect population eventually goes extinct.
The caseλ 1 < 0 requires further investigation, since an additional nontrivial value of C * I is possible; we may then require suitable threshold conditions for the eventual extinction of the epidemic. Now, let λ 1 be the principal eigenvalue for the problem Notice that λ 1 ≥ n and that if γ 1 +∞, then λ 1 λ * 1 , where λ * 1 is the principal eigenvalue for Notice that λ * 1 may be as large as we wish if γ 1 is sufficiently large and/or if ω is sufficiently large.
The following result concerns the eradicability of the disease in the most interesting situation whenλ 1 < 0, q > and C I 0 = 0 L ∞ (Ω) .
Notice that condition (19) holds if, for example, γ 1 is sufficiently large and/or the subset ω is sufficiently large.
Let (ĩ 1 ,ĩ 2 ) be the solution to 20) (the existence, uniqueness, and nonnegativity of the solution follows via a fixed point argument). We have that 2}. This follows from considering the system satisfied by (ĩ 1 − i 1 ,ĩ 2 − i 2 ) and showing in a standard manner that If we multiply the first equation in (20) byĩ 1 and integrate over Ω, we obtain that If we multiply the second PDE in (20) byĩ 2 and integrate over Ω, we obtain that We may infer that This condition holds for any sufficiently small ε > 0 [because (19) holds]. On the other hand, since condition (19) holds, it follows that for any ε > 0 sufficiently small, we have that and we conclude that ĩ 1 (t) 2 2 + ĩ 2 (t) 2 2 converges to 0 as t → +∞, at the rate of exp{−2at}. Since Ω is bounded, it follows that ı 1 (·, t) → 0,ĩ 2 (·, t) → 0 in L 1 (Ω), and that as t → +∞, at least as fast as exp{−at} (i.e., the total number of infected insects and the total number of infected trees tend to 0, exponentially). It means that the disease is eradicable.

Remark 1
Recall that if the diffusion coefficient is strictly positive, d > 0, then any point of the domain Ω is strongly connected to any other point of Ω, so that actions taken in a subregion ω ⊂ Ω will eventually propagate to the whole habitat.
On the other hand, due to the lack of diffusion for i 2 , it is reasonable to expect the control to be more effective (leading to a faster decay of i 1 and i 2 ) if γ 22 , γ 23 , and γ 3 act on the whole domain Ω. In this case, the controlled system becomes The boundary conditions and the initial conditions are as before.
In this special case, the following result holds: as t → +∞.

Notice that (25) is a weaker assumption than (19).
Proof Using a comparison result for (i 1 , i 2 ), we obtain that (26) If we multiply the first equation in (26) byĩ 1 and integrate over Ω, we obtain that (by Rayleigh's principle).
If we multiply the second PDE in (26) byĩ 2 and integrate over Ω, we obtain that We may infer that This condition holds for ε > 0 sufficiently small (because (25) holds). On the other hand, since condition (25) is satisfied, then for ε > 0 sufficiently small, we have that and we conclude that ĩ 1 (t) 2 2 + ĩ 2 (t) 2 2 converges to 0 as t → +∞, at the rate of exp{−2ãt}. Since Ω is bounded, it follows that and that as t → +∞, at least as fast as exp{−ãt} (i.e., the total number of infected insects and the total number of infected trees tend to 0, exponentially). This means that the disease is eradicable. Notice thatã ≥ a.

Numerical Simulations
The numerical strategy adopted to approximate the controlled system (Equations (6)-(9)) 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., Quarteroni and Valli (1994).
Space discretization We first apply a standard Galerkin procedure to the weak formulations of the controlled system [Eqs. (6)-(9)]. For the computational domain Ω, we have taken a rectangle of size 400 × 80 km 2 , which mimics, for example, 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 bilinear finite elements (Q1), yielding a total number 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 the spatial discretization, we obtain a 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 nonlinear reaction terms are approximated by forward Euler. As a result, at each time step we must find the solution of a linear system of algebraic equations of dimension 4 × 8241 = 32964 degrees of freedom. The system is solved by Gaussian elimination with the built-in function in MATLAB. The time step size is Δt = 0.0002 years. For further details on the numerical discretization of parabolic problems, we refer the reader to Quarteroni and Valli (1994).

Parameter Calibration
The values of the simulation parameters are listed in Table 2. The diffusion coefficient in the PDEs is set to d = 1e − 4 km 2 /year. The total simulation time is T = 10 years. The initial distributions, used in all next numerical simulations, are the following: That is, we assume constant initial distributions for healthy insects and trees, and we localize the initial presence of infected insects on the right-hand side of the domain.

Numerical Results
A couple of numerical experiments are reported here to show the efficacy of the controls on the weed biomass and the olive cultivar.
In Figs. 1 and 2, we have taken M = 1 and χ = 0.03, while in Figs. 3 and 4, M = 1 and χ = 0.01. The colorbar indicates the scaled values of the distribution corresponding to the colors adopted in the map: It goes from blue, associated with low values, to yellow, associated with high values.
In both computational experiments, control measures have been applied only in the right-hand section of the habitat, simulating the subregion of the Apulian region already affected by the xylella epidemic, augmented by a "containment band," i.e., the subregion Experiment 1 χ = 0.03, ζ = 0.8 year −1 : we first run the model without control, and thus, the control parameters (γ 21 , γ 1 , γ 23 , γ 22 , γ 3 ) are set to zero everywhere. The xylella epidemic starts spreading as a travelling wave from the right portion of the domain, where the initial condition for the infected insects is positive; see Fig. 1.  Silva et al. (2015) and Yurtsever (2000) χ   When we apply regional control, setting the control parameters (γ 21 , γ 1 , γ 23 , γ 22 , γ 3 ) = (0.999, 0.2, 0.2, 0.2, 0.2), the xylella epidemic is stopped at the origin and no travelling front occurs; see Fig. 2. This result confirms our conjecture that a significant weed cut in the olive orchards strongly decreases the carrying capacity of insects, yielding a successful blocking of the xylella epidemic spread.
In order to investigate how the diffusion coefficient d influences the results, we have re-run both tests with d = 2e−4, 5e−4, 1e−3 (all in km 2 /year). We do not observe any significant qualitative difference.
Experiment 2. χ = 0.01, ζ = 0.2 year −1 : as above, we first run the model without control, and thus, the control parameters (γ 21 , γ 1 , γ 23 , γ 22 , γ 3 ) are set to zero everywhere. The xylella epidemic starts spreading as a travelling wave from the right portion of the domain, where the initial condition of the infected insects is positive; see Fig. 3.
When we apply a regional control, setting the control parameters (γ 21 , γ 1 , γ 23 , γ 22 , γ 3 ) = (0.6, 0.2, 0.2, 0.6, 0.2), the xylella epidemic is stopped at the origin and no travelling front occurs; see Fig. 4. This result clearly indicates that for more resistant olive cultivar, we may still block an epidemic with a lower level of weed cut in the olive orchards, even though, by using χ = 0.001, we have a larger effective insect carrying capacity. As in the previous experiment, we have also re-run both tests with d = 2e−4, 5e−4, 1e−3 (all in km 2 /year). Even in this case, we do not observe any significant qualitative difference.

Concluding Remarks
The results reported in this paper confirm that the most promising target for an effective and cost-efficient control of the X. fastidiosa epidemic is represented by agricultural management practices consisting of the removal of the weeds in the whole relevant habitat of the olive orchards. A further interesting strategy, as expected, is represented by the use of more resistant cultivar (Brunetti et al. 2020) [see also Schneider et al. (2020) and references cited therein].
Here, we have extended the ordinary differential equation (ODE) model introduced in Brunetti et al. (2020) to an integro-partial differential system which takes into account the spatial structure of the relevant epidemic system, subject to possible control strategies, including weed cut, insect traps, treated nets, more resistant cultivar, etc. As emphasized in Remark 1, theoretical results show that eradication of an epidemic outbreak is possible by the implementation of the above mentioned control measures only on a suitable subregion of the whole habitat, including the area already affected by the outbreak, possibly augmented by a suitable "containment band."   Unfortunately in practice, it might be unfeasible to act on the diffusion parameters; they are related to the ecological structure of the relevant habitat, including the behavior of the insect population. Intervention strategies should then include possible modification of insect behavior, such as by the use of treated nets.
Future investigations will be devoted to the search of an optimal set of control parameters, γ 's, together with an optimal subregion of intervention ω. It is clear that, for a realistic optimal control problem, relevant participating costs need to be included: e.g., production, losses, and management costs.
Once again, validation of the model proposed here represents a key issue: although we have tried to make explicit the assumptions underlying our model, they have not yet been validated by comparison with experimental data. Therefore, we caution that our results are far from being conclusive for X. fastidiosa subsp. pauca -Philaenus spumarius olive tree epidemics. However, it is desirable that, with additional features that make it more realistic, our model might provide the foundations for designing optimal control strategies by public authorities.
It is worth reporting here a statement (abridged) taken from the very recent paper (Matricardi et al. 2020) on COVID-19 modeling, which can be applied to the role of any model: "a model is only an approximate interpretation of reality and it is always wrong in some small or relevant elements. The destiny of the model presented here is to be rapidly improved thanks to novel knowledge coming from new observations and better assumptions. The Authors hope that many and more brilliant minds will read the present pages, will identify and highlight putative mistakes, will get inspiration for their research, and will produce better, more complete, and useful models........ If the speculations presented here on implications for surveillance, control, and therapy of [Xylella] will contribute, even only minimally, to save ..... [olive trees]........., then the Authors have accomplished their small mission. " by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.