# A compartmental model for Xylella fastidiosa diseases with explicit vector
seasonal dynamics

Alex Giménez-Romero  Instituto de Física Interdisciplinar y Sistemas Complejos
(IFISC, CSIC-UIB), Campus UIB, 07122 Palma de Mallorca, Spain.  Eduardo
Moralejo  Tragsa, Passatge Cala Figuera 6, 07009 Palma de Mallorca, Spain.
Manuel A. Matías  Instituto de Física Interdisciplinar y Sistemas Complejos
(IFISC, CSIC-UIB), Campus UIB, 07122 Palma de Mallorca, Spain.

###### Abstract

The bacterium Xylella fastidiosa (Xf) is mainly transmitted by the spittlebug,
Philaenus spumarius, in Europe, where it has caused significant economic
damage to olive and almond trees. Understanding the factors that determine
disease dynamics in pathosystems that share similarities can help design
control strategies focused on minimizing transmission chains. Here we
introduce a compartmental model for Xf-caused diseases in Europe that accounts
for the main relevant epidemiological processes, including the seasonal
dynamics of P. spumarius. The model was confronted with epidemiological data
from the two major outbreaks of Xf in Europe, the olive quick disease syndrome
(OQDS) in Apulia, Italy, caused by the subspecies pauca, and the almond leaf
scorch disease (ALSD) in Majorca, Spain, caused by subspecies multiplex and
fastidiosa. Using a Bayesian inference framework, we show how the model
successfully reproduces the general field data in both diseases. In a global
sensitivity analysis, the vector-plant and plant-vector transmission rates,
together with the vector removal rate, were the most influential parameters in
determining the time of the infected host population peak, the incidence peak
and the final number of dead hosts. We also used our model to check different
vector-based control strategies, showing that a joint strategy focused on
increasing the rate of vector removal while lowering the number of annual
newborn vectors is optimal for disease control.

Keywords: Epidemiology, Population Biology, Modeling.

Mathematical and computational modeling in Ecology and, in particular,
Epidemiology have been recently recognized as powerful approaches to guide
empirical work and provide a framework for the synthesis, analysis and
development of conservation plans and policy-making ([31, 37, 46, 9]). Plant
epidemics, mainly plant-virus diseases, have been often described by
compartmental models, which deal with the overriding importance of
transmission mechanisms in determining epidemic dynamics ([29, 27, 34]). These
models have contributed to providing answers to some questions related to the
ecology of plant diseases and have led to direct applications in disease
control while guiding research directions ([26]).

The emergence of vector-borne plant pathogens in new areas causing huge
economic impacts, such as Xylella fastidiosa and the Candidatus Liberibacter
spp. (Huanglongbing or citrus greening), has sparked interest in modeling
vector-transmitted plant-disease epidemics ([10, 26]). The vector-borne
bacterium X. fastidiosa (Xf) is a multi-host pathogen endemic to the Americas
that causes economically important diseases, mostly in woody crops ([25]). Xf
is a genetically diverse species with three evolutionary well-defined clades
forming the pauca, fastidiosa, and multiplex subspecies, native from South,
Central, and North America, respectively ([51]). Within each subspecies,
diverse genetic lineages with different host ranges are found. Xf is
transmitted non-specifically by xylem-sap-feeding insects belonging to the
sharpshooter leafhoppers (Hemiptera: Cicadellidae, Cicadellinae) and
spittlebugs (Hemiptera: Cercopoidae) ([41]).

Recently, Xf has gained renewed interest due to the massive mortality of olive
trees in Apulia, Italy ([45]). The first focus of the olive quick decline
syndrome (OQDS) was detected in 2013 around Gallipoli ([44]) and since then
has spread throughout the region by the meadow spittlebug, Philaenus
spumarius. Although this was the first official detection of Xf in Europe, it
has recently been demonstrated that the pathogen arrived much earlier in
Corsica ([48]) and the Balearic islands ([36]). Around 1993, two strains of
the subspecies fastidiosa (ST1) and multiplex (ST81) were introduced from
California to Mallorca (Spain) with infected almond plants ([36]). To date,
over 80% of the almond trees in Mallorca show leaf scorch symptoms and the
outbreak has changed the iconic rural landscape of this Mediterranean island
([38]).

Several epidemic models have been already developed for Xf-diseases, but they
lack a realistic description of some relevant processes ([26]). Furthermore,
some of these models assume a simple general form for infected host dynamics
([53, 1, 16]) or use a simplified S-I compartmental scheme for hosts,
disregarding important features such as the latent period or the host
mortality rate ([48]). Models that do take these features into account,
however, do not explicitly model the population of vectors responsible for
disease transmission ([54]). Other more recent models have taken a step
further in explicitly modeling the vector population ([7]), but the
characterization of its dynamics is still relatively simple, as it overlooks
the known seasonal patterns of vector abundance. Thus, there is a need to
continue advancing in the modeling of Xf diseases by developing more realistic
models that can elucidate the fundamental processes involved in vector-host-
pathogen interactions and help to design effective control strategies.

In this work, we develop a deterministic continuous-time compartmental model
to describe the general epidemiological dynamics of diseases produced by Xf in
Europe, explicitly accounting for key biological aspects of the main vector P.
spumarius, including its seasonal dynamics. Our model is able to describe
field data from the two major European outbreaks: the olive quick disease
syndrome (OQDS) in Apulia, Italy, caused by the subspecies pauca, and the
almond leaf scorch disease (ALSD) in Majorca, Spain, caused by subspecies
multiplex and fastidiosa. We aimed to find the most influential parameters in
the model with respect to incidence and mortality in both diseases by
performing a global sensibility analysis. With this information, the next goal
was to explore control strategies acting especially on the vector population.

## MATERIALS AND METHODS

### Epidemic model: the SEIR-V model

We developed a deterministic continuous-time compartmental model that
incorporates the specific biological features of Xf diseases in Europe,
including the dynamics of the main relevant vector P. spumarius ([8]). To
build the model we took the following considerations: (i) we assume there is
no winter recovery of infected hosts and thus they die sometime after
infection; (ii) hosts show an asymptomatic period in which they are non-
infectious in practice (exposed compartment) because the bacteria are not yet
systemically extended ([50, 49]), while vectors are infectious immediately
after acquiring the bacterium ([18]); (iii) vectors have an annual life cycle
without mother-to-offspring disease transmission ([19, 39]), so we consider
the annual emergence of susceptible newborn vectors and a constant death rate
for both susceptible and infected vectors; (iv) infected vectors carry the
bacterium during their entire lifespan without affecting their fitness ; and
finally, (v) we do not consider host recruitment or natural death given that
the typical development time of Xf-epidemics is faster than the typical host’s
life cycle.

Altogether, our deterministic continuous-time compartmental model consists of
six compartments, four describing the host population (susceptible, $S_{H}$,
exposed, $E_{H}$, infected, $I_{H}$, and removed, $R_{H}$), and two describing
the vector population (susceptible, $S_{V}$, and infected, $I_{V}$). The model
is defined according to the following processes,

$\displaystyle
S_{H}+I_{V}\stackrel{{\scriptstyle\beta}}{{\rightarrow}}E_{H}+I_{V},\quad
E_{H}\stackrel{{\scriptstyle\kappa}}{{\rightarrow}}I_{H},\quad
I_{H}\stackrel{{\scriptstyle\gamma}}{{\rightarrow}}R_{H}$ (1) $\displaystyle
S_{V}+I_{H}\stackrel{{\scriptstyle\alpha}}{{\rightarrow}}I_{V}+I_{H},\quad
S_{V}\stackrel{{\scriptstyle\mu}}{{\rightarrow}}\varnothing,\quad
I_{V}\stackrel{{\scriptstyle\mu}}{{\rightarrow}}\varnothing\ ,$

which are illustrated in Fig. 1, being the birth of new susceptible vectors
described as a source term. Thus, the host-vector compartmental model is
written as,

$\displaystyle\dot{S}_{H}$ $\displaystyle=-\beta S_{H}I_{v}/N_{H}$ (2)
$\displaystyle\dot{E}_{H}$ $\displaystyle=\beta S_{H}I_{v}/N_{H}-\gamma I_{H}$
$\displaystyle\dot{I}_{H}$ $\displaystyle=\kappa E-\gamma I$
$\displaystyle\dot{R}_{H}$ $\displaystyle=\gamma I_{H}$
$\displaystyle\dot{S}_{v}$
$\displaystyle=N_{v}(0)\sum_{n=1}^{\infty}\delta(t-nT)-\alpha
S_{v}I_{H}/N_{H}-\mu S_{v}$ $\displaystyle\dot{I}_{v}$ $\displaystyle=\alpha
S_{v}I_{H}/N_{H}-\mu I_{v}\ .$

Figure 1: Schematic representation of the model Eq. 2. Boxes are the
compartments in which the population is divided, solid curved arrows represent
changes in state, i.e. transitions between compartments, dashed arrows depict
the crossed interaction between hosts and vectors and solid straight arrows
represent demographic changes in vector population.

The model describes the exposure of susceptible hosts, $S_{H}$, at a rate
$\beta$ through their interaction with infected vectors, $I_{v}$, while
susceptible vectors, $S_{v}$, get infected immediately at a rate $\alpha$
through their interaction with infected hosts $I_{H}$. Exposed hosts get
infected at rate $\kappa$, being the mean latent period $\tau_{E}=1/\kappa$,
while infected hosts die at rate $\gamma$, having a mean infectious period of
$\tau_{I}=1/\gamma$. Infected vectors stay infected for the rest of their
lifetime. Regarding the seasonal dynamics of vectors, we assume that new
adults emerge synchronously each year in fields being all susceptible. This is
represented by the term $N_{v}(0)\sum_{n=1}^{\infty}\delta(t-nT)$ in Eq. 2,
where $T=$1\text{\,}\mathrm{y}\mathrm{r}$$ is the period and $\delta(t-nT)$ is
Dirac delta function, with the property: $\delta(x-x_{0})=0$ for all $x\neq
x_{0}$ and $\int g(t)\delta(t-t_{0})\mathrm{d}t=g(t_{0})$. Vectors die or exit
the field at a given rate $\mu$, which we consider identical for susceptible
and infected vectors. For simplicity, we consider that the quantity of annual
newborn adults, $N_{v}(0)$, is constant. This outburst of new adults followed
by an exponential decay resembles the temporal patterns on the abundance of P.
spumarius observed in crop fields ([2, 4, 12, 33]) (see Fig. S1).

In Eq. 2 the crossed nonlinear terms, $S_{H}I_{v}$ and $S_{v}I_{H}$ are
written divided by the total host population, $N_{H}$. Thus, the plant-to-
vector infection process is modeled using standard incidence, which is
frequency dependent, while the vector-to-plant infection process is modeled
using mass action incidence, which is density dependent ([35]). This implies
that doubling the number of vectors in the crop field would double the number
of resulting infected hosts, as this process is population-dependent, while
doubling the number of hosts would not result in more vectors being infected,
as this process only depends on the contact probability. We think this is the
most reasonable assumption because increasing the number of hosts is expected
to also increase the area of the field, while the number of vectors is
independent.

### Basic reproductive number

The basic reproductive number, $R_{0}$, of the model cannot be trivially
computed using standard methods such as the Next Generation Matrix (NGM)
([17]), as there is no pre-pandemic fixed point in the system of differential
equations Eq. 2. For periodic vector populations, rigorous methods have been
developed ([3]) but not for the case of growing or decaying vector
populations. Here we use the simple method developed in the work of [22] (see
Appendix II), which effectively computes the average number of secondary
infections produced by an initially infected individual in one generation.
Thus, the effective basic reproductive number is given by

$R_{0}=\frac{\beta\alpha}{\mu\gamma}\frac{S_{H}(0)}{{N_{H}}^{2}}\frac{N_{v}(0)}{\mu\tau}\left(1-e^{-\mu\tau}\right)\
,$ (3)

where $\tau$ corresponds to the time length of one generation, in our case one
year. This $R_{0}$ is calculated using the initial susceptible host
population, $S_{H}(0)$. Below we will also use a time-dependent $R_{0}(t)$
using $S_{H}(t)$.

### Epidemiological data

Epidemiological data from an ALSD outbreak in the island of Mallorca Balearic
Islands, Spain were taken from [36]. Dated phylogenetic analysis and estimates
of disease incidence showed that the introduction of both subspecies occurred
around 1993 and $\sim 79$% almond trees were infected in 2017 ([36]). The
annual proportion of infected individuals in the almond tree population
between 1993 and 2017 was estimated by analyzing through qPCR the presence of
Xf-DNA in the growth rings of $34$ sampled trees (cf. Fig. 3 in ([36])). The
disease progression curve was estimated without distinguishing whether
infections were caused by multiplex or fastidiosa subspecies. In addition, a
two-sided bootstrap confidence interval for each data point was set using the
SciPy bootstrap function in Python ([52]). On the other hand, epidemic data
for OQDS were retrieved from ([54]). The data consisted of 2 to 3 yearly
censuses of symptom prevalence in 17 olive groves infected with Xf subsp.
pauca in Apulia, Italy, which were aggregated to fit our model as shown in
Fig. 4 in ([54]). Because the compartments of our model are not in one-to-one
correspondence with those shown in the work of White et al. ([54]), we used
the sum of the symptomatic and desiccated infected trees in the dataset
($I_{S}+I_{D}$) to fit the sum of the infected and dead trees ($I+R$) and the
sum of susceptible and asymptomatic hosts ($S+I_{A}$) to fit the sum of
susceptible and exposed hosts ($S+E$). The processed data used to fit the
model can be found in ([21]), while the raw data can be found in the
supplementary data accessible online of the cited articles ([36, 54]).

### Model fitting through Bayesian Inference

We employed an informative normal $\mathcal{N}(\hat{\mu},\hat{\sigma}^{2})$
prior distribution, with $\hat{\mu}$ and $\sigma$, the mean and standard
deviation, respectively, for previously measured parameters in the literature,
such as the infected and latent periods for ALSD,
$\tau_{I}\sim\mathcal{N}(14,4)$, $\tau_{E}\sim\mathcal{N}(4,1)$ ([50, 36]) and
OQDS, $\tau_{I}\sim\mathcal{N}(3.5,1)$, $\tau_{E}\sim\mathcal{N}(1.75,0.5)$
([18]). The corresponding rates are given by $\gamma=1/\tau_{I}$ and
$\kappa=1/\tau_{E}$, respectively. Similarly, a prior normal distribution was
used for the removal rate of vectors, $\mu\sim\mathcal{N}(0.02,0.0075)$, as
the mean value $\mu=0.02$ already captures the vector dynamics observed in
field-data (Fig. S1). Regarding the prior distribution for the transmission
rates a very wide and uninformative uniform prior distribution,
$\beta\sim\mathcal{U}(0.001,1)$ and $\alpha\sim\mathcal{U}(0.001,1)$, was used
for each parameter. The number of hosts, $N_{H}$, was already provided in the
datasets, while, given the lack of information about the vector population, we
assumed $N_{v}(0)=N_{H}/2$ for the initial vector population of each year.
However, we tested the robustness of our results by changing $N_{v}(0)$.

The posterior distributions of the parameters were approximated using the
Markov Chain Monte Carlo algorithm No U-Turn Sampler (NUTS) with the
recommended target acceptance rate of 65% ([24]). To ensure a proper
convergence, we constructed three independent Markov Chains with $10^{5}$
iterations each after a burn-in of $10^{4}$ iterations and checked that the
results were statistically equivalent. For each chain, we started at the
maximum-likelihood parameters yielded by the Nelder-Mead algorithm with 1000
iterations.

The parameters of our compartmental model were determined by fitting the model
to data by means of a Bayesian Inference framework using the Turing.jl package
([20]) in Julia ([5]). The scripts used to fit the model can be found in
([21]).

### Sensitivity Analysis

We performed a Global Sensitivity Analysis (GSA) ([43]) of the model to assess
the relative contribution of its parameters and their interactions with
different features of the epidemic. In contrast to the Local Sensitivity
Analysis (LSA), the GSA assesses the influence of a large domain of the
parameter space in the desired outputs of the model. We performed GSA by means
of a variance-based analysis, the Sobol method ([47]). This particular method
provides information not only on how a particular parameter alone influences
the model outputs (as happens with LSA), but also due to the nonlinear
interactions among two or more parameters. Briefly, the method considers the
model output, $Y$, as a general function of the inputs, $f(x_{1},...,x_{n})$,
so that the variance of the output, $Var(Y)$ is decomposed as the sum of the
variances given by the variations of the parameters alone and its
interactions,
$Var(Y)=\sum_{i=1}^{n}Var(f(x_{i}))+\sum_{i<j}^{n}Var(f(x_{i},x_{j}))+\cdots$.
This information is organized in what are known as Sobol indices. The total
order indices are a measure of the total variance of the output quantity
caused by variations of the input parameter and its interactions,
$S_{T}=Var(f(x_{1},...,x_{n}))/Var(Y)$. First order (or “main effect”) indices
are a measure of the contribution to the output variance given by the
variation of the parameter alone, but averaged over the variations in other
input parameters, $S_{i}=Var(f(x_{i}))/Var(Y)$. Second-order indices take into
account first-order interactions between parameters,
$S_{ij}=Var(f(x_{i},x_{j}))/Var(Y)$. Further indices can be obtained,
describing the influence of higher-order interactions between parameters, but
these are not going to be considered.

Following the Sobol method, we analyzed the variation of the time at which the
infected population peaks, $t_{peak}$, the magnitude of this peak, $I_{peak}$
and the final number of dead hosts, $R_{\infty}$, relative to variations of
the model parameters. The method was implemented within the Julia high-level
programming language ([5]) using the sub-package DiffEqSensitivity.jl in
DifferentialEquations.jl package ([40]).

## RESULTS

### Model fit and parameter estimates

The posterior distributions of the fitted parameters including their estimated
mean and median for ALSD and OQDS are shown in Figs. 2 and 3, respectively,
together with the assumed prior distributions. We observe that the literature-
driven priors for the latent and infected period, $\tau_{E}$ and $\tau_{I}$,
were already very good guesses and changed slightly converging to the
appropriate distribution that better fitted the epidemic data for both ALSD
and OQDS (Fig. 2(A-B) and Fig. 3(A-B)). Similarly, the prior for the vector
removal rate, $\mu$, obtained from field data, was good enough so that little
changes were needed for convergence (Fig. 2(C) and Fig. 3(C)). On the other
hand, we also observe that the completely uninformative priors for the
transmission rates successfully converged to the posterior distributions (Fig.
2(D-E) and Fig. 3(D-E)).

Figure 2: Posterior (blue histograms) and prior (red line) distributions of
the model parameters for ALSD. Solid and dashed black lines correspond to the
mean and median of the posterior distributions. (A) Host infected period
$\tau_{I}=1/\gamma$. (B) Host latent period $\tau_{E}=1/\kappa$. (C) Vector
removal rate $\mu$. (D) Vector infection rate $\alpha$. (E) Host infection
rate $\beta$. (F) The variance of the field data $\sigma^{2}$.

The latter distributions are far from a Gaussian-like shape (note that the
x-axis is log-scaled), being heavy-tailed. This kind of distribution highly
distorts the statistical measures of mean, median and standard error,
indicating that the estimates for transmission rates are not as robust as the
estimates for the other parameters. These rather uninformative distributions
are probably arising because of the lack of data about the vector, i.e.
$S_{v}(t)$ and $I_{v}(t)$, to constrain the fits. In essence, many
combinations of $\alpha$ and $\beta$ can similarly fit the host data while
yielding quite different time series for $S_{v}(t)$ and $I_{v}(t)$, which
cannot be contrasted due to the lack of field data. Nevertheless, the obtained
best-fit mean and median parameters, although quite different, are able to
perfectly fit the data (Fig. 4). Finally, we also observe that the variance
for the field data also converged to a bell-shaped distribution.

Mean and median parameter estimates, i.e. the best-fit parameter values for
ALSD and OQDS, are summarized in Tables 1 and 2, respectively. As already seen
from the posterior distributions, the best-fit values for $\tau_{E}$,
$\tau_{I}$ and $\mu$ are close to the ones given by literature and field data
for both diseases. Conversely, $\alpha$ and $\beta$ are rather uninformative,
as their 95% confidence intervals cover almost two orders of magnitude. This
again indicates that without some data about the evolution of the vector
states in time, $S_{v}(t)$ and $I_{v}(t)$, it is nearly impossible to find the
proper values for these parameters.

Figure 3: Posterior (blue histograms) and prior (red line) distributions of
the model parameters for OQDS. Solid and dashed black lines correspond to the
mean and median of the posterior distributions. (A) Host infected period
$\tau_{I}=1/\gamma$. (B) Host latent period $\tau_{E}=1/\kappa$. (C) Vector
removal rate $\mu$. (D) Vector infection rate $\alpha$. (E) Host infection
rate $\beta$. (F) Variance of the field data $\sigma^{2}$. Table 1: Estimated
epidemiological parameters from Bayesian model fitting to the disease
progression curve of ALSD in Mallorca.

Parameter | Definition | Units | Posterior Mean | Posterior Median | 95% C.I.
---|---|---|---|---|---
$\tau_{I}$ | Host infected period | yr | $13.84$ | $13.82$ | $[7.12,20.47]$
$\tau_{E}$ | Host latent period | yr | $4.46$ | $4.47$ | $[2.88,5.99]$
$\beta$ | Host infection rate | $\text{\,}\frac{\\#host}{\\#vector\cdot day}$ | $0.062$ | $0.02$ | $[0.0061,0.3013]$
$\alpha$ | Vector infection rate | $\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}^{-1}$ | $0.15$ | $0.086$ | $[0.0047,0.54]$
$\mu$ | Vector removal rate | $\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}^{-1}$ | $0.0222$ | $0.0221$ | $[0.015,0.030]$
$R_{0}$ | Basic reproductive number | - | 133 | 25 | -
|  |  |  |  |

Table 2: Estimated epidemiological parameters from Bayesian model fitting to
the disease progression curve of OQDS in Apulia.

Parameter | Definition | Units | Posterior Mean | Posterior Median | 95% C.I.
---|---|---|---|---|---
$\tau_{I}$ | Host infected period | yr | $3.61$ | $3.60$ | $[2.06,5.20]$
$\tau_{E}$ | Host latent period | yr | $1.24$ | $1.25$ | $[0.70,1.75]$
$\beta$ | Host infection rate | $\text{\,}\frac{\\#host}{\\#vector\cdot day}$ | $3.44$ | $2.60$ | $[0.55,8.79]$
$\alpha$ | Vector infection rate | $\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}^{-1}$ | $0.0031$ | $0.0022$ | $[0.0005,0.0084]$
$\mu$ | Vector removal rate | $\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}^{-1}$ | $0.0240$ | $0.0240$ | $[0.014,0.035]$
$R_{0}$ | Basic reproductive number | - | 33 | 21 | -
|  |  |  |  |

Overall, the estimated mean and median parameters consistently fit the
provided data within the 99% confidence interval for the ALSD and OQDS
outbreaks (Fig. 4(B,D)). We also computed the instantaneous reproductive
number, $R_{0}(t)$, by using Eq. 3 with $S_{H}(t)$ instead of only $S_{H}(0)$
along the simulation. Noteworthy, $R_{0}(t)=1$ coincides with the stopping of
new infections being produced, i.e. the number of exposed hosts does not
increase (Fig. 4(A,C)). This supports our approximate method for computing the
reproductive number for Xf diseases (Appendix II, Eq. 3). Due to the different
time scales of both epidemics
($\tau_{I}^{ALSD}+\tau_{E}^{ALSD}>\tau_{I}^{OQDS}+\tau_{E}^{OQDS}$), the OQDS
outbreak dies out earlier than the one for ALSD.

We notice that for ALSD a large proportion of the vector population gets
infected every year (Fig. 4(A)), while a very small proportion is needed in
OQDS to produce a lethal outbreak (Fig. 4(C)). However, this last statement is
rather unrealistic, as around 50% of the vectors that are captured in Apulia
are indeed infected by Xf ([8, 13]). Thus, the evolution of the infected
vector population should be qualitatively similar to that obtained for ALSD
(Fig. 4(C)). As previously explained, different suitable values of these
parameters should give rise to similar progression curves for the hosts while
different ones for the vectors cannot be contrasted due to the lack of data.

Figure 4: (A) Simulation of the model with the best-fit parameters for ALSD.
(B) Model fit to field data by means of the mean and median values of the
posterior distributions of the parameters for ALSD. (C) Simulation of the
model with the best-fit parameters for OQDS. (D) Model fit to field data by
means of the mean and median values of the posterior distributions of the
parameters for OQDS. The gray-shaded area corresponds to the 99% confidence
interval. The error bars for the field data correspond to their 95% confidence
interval obtained with a bootstrapping technique.

Nevertheless, by manually exploring other values for $\alpha$ and $\beta$
parameters, we can obtain a more biologically plausible scenario for the OQDS
that is still able to fit the available data for the hosts. Fig. 5(A) shows a
simulation of the model with previously inferred best-fit median parameters
for OQDS. By changing the values of $\alpha$ and $\beta$, we obtain a more
realistic scenario, i.e. around a 50% of the vector population getting
infected during the outbreak (Fig. 5(B)) ([8, 13]). This change in the
transmission parameters only affects the progression curve of the infected
vector population, being the progression of the host compartments practically
unchanged (Fig. 5(C)). Anyway, both sets of parameter values for $\alpha$ and
$\beta$can properly fit the field data (Fig. 5(D)).

Figure 5: (A) Simulation of the model with the original best-fit parameters
for OQDS. (B) Simulation of the model with the original best-fit parameters
for OQDS but with different $\alpha$, $\beta$ values. (C) Comparison of the
progression curves. Note that the curves for the hosts are very similar while
the curve for the infected vector population is very different. (D) Comparison
of the model fit to the data with both simulations. Solid lines correspond to
results with the original best-fit parameters while dash-dot lines correspond
to the results of the more realistic scenario with different $\alpha$ and
$\beta$.

The model adjusted to the progression curves of both diseases indicates that
the transmission rate $\alpha$ must be greater than $\beta$ when the
proportion of infected vectors is relatively high ($>30\%$). We checked if the
relation between $\alpha$ and $\beta$ held when changing the assumed
$N_{v}(0)=N_{H}/2$, obtaining that it kept approximately the same for very
different values of the initial vector population.

### Global Sensitivity Analysis

We computed the sensitivity indices for the model parameters with respect to
the more relevant quantities of interest, namely, the time at which the number
of infected hosts is maximal, $t_{\textrm{peak}}$, the maximum number of
infected hosts, $I_{\textrm{peak}}$ and the final number of dead hosts,
$R_{\infty}$. The results were obtained exploring the parameter space
constrained to the intervals $\\{\beta\in(0.001,0.1),\ \tau_{E}\in(3,7),\
\tau_{I}\in(5,25),\ \alpha\in(0.001,1),\ \mu\in(0.01,0.04)\\}$ using $10^{4}$
Quasi-Monte Carlo samples and are summarized in Fig. 6.

Figure 6: Global Sensitivity Analysis of the model parameters performed with
the Sobol method with respect to the time at which the infected population
peaks, $t_{peak}$ (A-B), the magnitude of this peak, $I_{peak}$ (C-D) and the
final number of dead hosts, $R_{\infty}$ (E-F). The left column (A,C,E) shows
the total and first-order indices and the right column (B,D,F) shows the
second-order indices.

Parameters $\alpha,\ \beta$ and $\mu$ are the most influential with regard to
the time at which the infected host population peaks, $t_{peak}$, the
magnitude of the peak, $I_{peak}$, and the final number of dead hosts,
$R_{\infty}$. The total output variance (total order indices) cannot be
explained by the variances of the parameters alone (first order indices) (Fig.
6). Therefore, higher-order interactions among the parameters importantly
affect the sensitivity of the quantities under study. Indeed, the contribution
to the total output variance of $\gamma$ and $\kappa$ for $t_{peak}$ and
$R_{\infty}$ come notably from higher-order interactions. This can be checked
in panels (B,C,F) of Fig. 6, in which the contribution to the output variance
from interactions between pairs of parameters (second order indices) is
represented. Interactions among the parameters contribute to increasing the
output variance with respect to $t_{peak}$ and $I_{peak}$, while the effect is
more heterogeneous in the case of $R_{\infty}$. In particular, the
interactions between $\alpha-\beta$ and $\alpha-\mu$ produce the main
contributions to the increase of output variance in all cases, while
$\kappa-\beta$, $\kappa-\alpha$ and $\kappa-\mu$ interactions decrease the
output variance.

### Epidemic control through vector management

The sensitivity analysis clearly indicates that acting on $\alpha,\beta$ and
$\mu$ is the best strategy to lower disease incidence and mortality. However,
controlling transmission rates is cumbersome so a different control strategy
based only on vector control is considered in this section. In our model,
there are two ways of implementing vector-population control: (i) decreasing
the typical time, $1/\mu$, that vectors spend between crops each year by some
mechanism (thus increasing $\mu$) and (ii) reducing the initial number of
vectors that invade crops each year (e.g. lowering $N_{v}(0)$ via egg or nymph
control ([30]).

We analyzed the effect of vector management by simulating epidemic outbreaks
using different values of $\mu$ and $N_{v}(0)$, and keeping the rest of
parameters as fitted for both ALSD and OQDS outbreaks (Fig. 7). In both
epidemics, decreasing the presence time as well as the number of vectors
contribute to controlling the epidemic by lowering $R_{0}$ and, consequently,
the final size of the epidemic, $R_{\infty}$. Furthermore, we observe that
decreasing vector presence is more efficient than decreasing its annual
initial population, i.e. we further reduce $R_{\infty}$, the final size of the
epidemic, by applying a similar reduction in the residence time $1/\mu$. This
could also be anticipated as $R_{0}$ depends quadratically on $1/\mu$ while
only linearly on $N_{v}(0)$ (Eq. 3). However, the minimal intervention
strategy, starting from the current situation in the $(1/\tau,N_{v}(0))$
parameter space that yields an absolute control of the epidemic, $R_{0}<1$,
involves a mixed strategy of lowering both $1/\mu$ and $N_{v}(0)$.

Figure 7: Epidemic control through vector management for ALSD in Mallorca (A)
and OQDS in Apulia (B). The white shaded line denotes $R_{0}=1$, the white
diamond corresponds to the parameter values of the fitted model. The white
star is the closest disease-free state to the current situation in this
representation.

## DISCUSSION

In this work, we have developed a deterministic continuous-time compartmental
model for Xylella fastidiosa vector-borne diseases in Europe. The model
attempts to characterize the main biotic processes that lead to the
development of epidemics, including the seasonal dynamics of the main vector,
P. spumarius. We show how the model is sufficiently general to represent with
some accuracy the parameters that determine the ALSD in Mallorca (Spain) and
the OQDS in Apulia (Italy), both transmitted by P. spumarius. To our best
knowledge, this is the first mathematical model describing Xf epidemics that
considers the temporal pattern of vector abundance observed in field data,
faithfully representing the known biological information about the
pathosystem. It includes a dynamic approximation of the non-stationary
populations of P. spumarius, mathematically represented by a sporadic source
term through which vectors are born every year, and an exponential decay term.
Due to the non-stationarity of the vector dynamics, $R_{0}$ in the model
cannot be computed with standard methods such as the Next Generation Matrix
([17]). To circumvent this problem, we applied an approximate method to
compute it as previously proposed by ([22]). We show that this approximate
$R_{0}$ correctly characterizes the epidemic, further validating the method
proposed by [22].

Nonlinear mathematical models of disease transmission enhance our
understanding of the different mechanisms operating in an epidemic, especially
compared with correlative or machine learning methods, often very useful in
practice but offering very little understanding. A key aspect to render these
models useful is the determination of the parameters from available data. If
this step can be properly performed, these models become very predictive and
especially helpful to design disease control strategies. However, an
appropriate calibration of the model relies on access to good-quality field
data, which is often the bottleneck for the application of this kind of model.
In the present study, the parameters have been obtained using a Bayesian
inference framework, which relies on probability distributions rather than
point-like measures. This way, mean or median values can be considered
together with their confidence intervals able to characterize the robustness
of the obtained parameters.

One of the conclusions of the study is that the available data for both
diseases is not enough to obtain robust estimates for all of the model
parameters. The lack of data about the vector population compartments yields
many possible values for the parameters that regulate transmission, $\alpha$
and $\beta$, provided that the progression of the host compartments correctly
fits the field data. In other words, very infectious vectors (high $\beta$)
that hardly ever get infected (low $\alpha$) can produce a similar outbreak
within the host population to that produced by very low infectious vectors
(low $\beta$) that get infected very often (high $\alpha$). The great
difference in these situations would be that, in the former, the infected
vector population would be very low, while in the latter, it would be quite
high. This is a manifestation of parameter unidentifiability from the fit
([11, 42]), which stresses the importance of transmission and calls for
detailed measurements of the vector population, and not just of the hosts.

In any case, our model shows that the vector-plant transmission process,
mediated by $\beta$, is somehow different than the plant-vector one, mediated
by $\alpha$. This heterogeneity can be caused by several factors: differences
in the efficiency of plant-vector transmission with respect to vector-plant
transmission, differences in contact rates, i.e. susceptible vectors contact
trees at a different rate than infected vectors; vector feeding preferences,
i.e. differences in the probability of contacting a susceptible host compared
to an infected host, etc. Indeed, our mathematical model assumes constant
contact rates with no preferences over any host state, so that under these
assumptions, it indicates that the probability of effectively transmitting the
pathogen from plant to vector is greater than from vector to plant.
Furthermore, we found that the timing and magnitude of the infected host peak
and the final number of dead hosts are mostly controlled by the vector-plant
transmission rate, $\beta$, the plant-vector transmission rate, $\alpha$ and
the vector removal rate $\mu$. Because these parameters are strongly related
to the vector, the analysis makes clear that enhancing the knowledge about the
vector, as well as obtaining precise data, is crucial to improve the modeling
of Xf diseases and pose important questions to be solved in specifically
designed experiments.

The fact that the most influential parameters of the model are those related
to the vector can be used to design appropriate disease control strategies.
Because acting on transmission rates is rather cumbersome, we argue that
control strategies should focus on reducing the vector population in crop
fields. In our model, this depends on two parameters, $\mu$, the rate at which
vectors exit the field and $N_{v}(0)$, the number of newborn susceptible
vectors every year (assumed constant in this study). Our results show that a
mixed strategy acting on both parameters is optimal to lower disease
prevalence and, eventually, eradicate the disease. Interestingly, we also show
that acting on the vector removal $\mu$ is more effective than controlling the
newborn vector population $N_{v}(0)$. In fact, most control strategies carried
out in practice for Xf diseases focus on the latter factor, reducing
$N_{v}(0)$ via egg or nymph control ([14, 32]). However, our results indicate
that alternative strategies based on increasing the removal (or dispersal)
rate of vectors should be explored. Furthermore, the evolution of the
population compartments of the hosts and vectors provides relevant information
on the epidemiology of both diseases. In both cases, the newly defined basic
reproductive number that accounts for a decaying vector population is very
predictive of the moment in which new infections are not produced anymore,
coinciding approximately with the peak of infected hosts. Therefore, any
intervention with control measures after this peak would have marginal effects
on future disease progression.

Our mathematical model is still rather simple, implementing only a few
relevant epidemic processes in contrast to the high complexity of the
pathogen-vector-host interactions occurring in plant epidemics. Indeed, the
model itself raises some questions about these interactions, for example,
whether or not contact rates are homogeneous. Another simplification of the
model is the fact that the spatial constraints and the intrinsic stochasticity
of the transmission processes are neglected. A straightforward extension of
the model would be to include a specific spatial setting and implement the
explicit motion of the vector within a stochastic framework, such as
Individual Based Models ([23]). With this, the effectiveness of current and
further control strategies could be tested and improved controlling for the
motion of the vector. For instance, the control strategy based on the removal
of symptomatic trees together with their surrounding trees at a given distance
could be implemented in the model, evaluate the current effectiveness
according to the present protocols and even provide improved parameters to be
implemented in the field. Of course, implementing a model in which the spatial
degrees of freedom are explicitly represented would require access to further
information about vector mobility and spatially resolved data to confront the
model, which is not currently available.

Mathematical models tested against experimental data increase our
understanding of the system under study. They also help to identify critical
parameters that require better prior information to adjust functions relating
to different variables and make the model predictions more accurate to suggest
and test control strategies ([15, 28]). Our mathematical model suggests a
certain lack of knowledge of the transmission processes and reveals that the
currently available data is not enough to fit complex models dealing with the
explicit dynamics of the vector population.

Funding. AGR and MAM acknowledge financial support from Grants No.
RTI2018-095441-B-C22 (SuMaEco) and No. PID2021-123723OB-C22 (CYCLE) funded by
MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe” and from
Grant MDM-2017-0711 (María de Maeztu Program for Units of Excellence in R&D)
funded by MCIN/AEI/10.13039/501100011033.

The authors declare no conflict of interest.

## References

  * [1] C. Abboud, O. Bonnefon, E. Parent and S. Soubeyrand  “Dating and localizing an invasion from post-introduction data and a coupled reaction–diffusion–absorption model”  In _Journal of Mathematical Biology_ 79.2, 2019, pp. 765–789  DOI: 10.1007/s00285-019-01376-x
  * [2] S. Antonatos et al.  “Seasonal Appearance, Abundance, and Host Preference of Philaenus spumarius and Neophilaenus campestris (Hemiptera: Aphrophoridae) in Olive Groves in Greece”  In _Environmental Entomology_ 50.6, 2021, pp. 1474–1482  DOI: 10.1093/ee/nvab093
  * [3] N. Bacaër  “Approximation of the Basic Reproduction Number $R_{0}$ for Vector-Borne Diseases with a Periodic Vector Population”  In _Bulletin of Mathematical Biology_ 69.3, 2007, pp. 1067–1091  DOI: 10.1007/s11538-006-9166-9
  * [4] D.. Beal, M. Cooper, M.. Daugherty, A.. Purcell and R… Almeida  “Seasonal Abundance and Infectivity of Philaenus spumarius (Hemiptera: Aphrophoridae), a Vector of Xylella fastidiosa in California Vineyards”  In _Environmental Entomology_ 50.2, 2021, pp. 467–476  DOI: 10.1093/ee/nvaa178
  * [5] J. Bezanson, A. Edelman, S. Karpinski and V.. Shah  “Julia: A fresh approach to numerical computing”  In _SIAM Review_ 59.1 SIAM, 2017, pp. 65–98  DOI: 10.1137/141000671
  * [6] F. Brauer, C. Castillo-Chavez, A. Mubayi and S. Towers  “Some models for epidemics of vector-transmitted diseases”  In _Infectious Disease Modelling_ 1.1, 2016, pp. 79–87  DOI: 10.1016/j.idm.2016.08.001
  * [7] M. Brunetti, V. Capasso, M. Montagna and E. Venturino  “A mathematical model for Xylella fastidiosa epidemics in the Mediterranean regions. Promoting good agronomic practices for their effective control.”  In _Ecological Modelling_ 432, 2020, pp. 109204  DOI: 10.1016/j.ecolmodel.2020.109204
  * [8] V. Cavalieri et al.  “Transmission of Xylella fastidiosa Subspecies Pauca Sequence Type 53 by Different Insect Species”  In _Insects_ 10.10, 2019, pp. 324  DOI: 10.3390/insects10100324
  * [9] Y.. Chew et al.  “Mathematical Models Light Up Plant Signaling”  In _The Plant Cell_ 26.1, 2014, pp. 5–20  DOI: 10.1105/tpc.113.120006
  * [10] C. Chiyaka, B.. Singer, S.. Halbert, J.. Morris Jr and A… Bruggen  “Modeling huanglongbing transmission within a citrus tree”  In _Proceedings of the National Academy of Sciences_ 109.30 National Acad Sciences, 2012, pp. 12213–12218  DOI: 10.1073/pnas.1208326109
  * [11] Gerardo Chowell  “Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts”  In _Infectious Disease Modelling_ 2.3, 2017, pp. 379–398  DOI: 10.1016/j.idm.2017.08.001
  * [12] D. Cornara et al.  “Spittlebugs as vectors of Xylella fastidiosa in olive orchards in Italy”  In _Journal of Pest Science_ 90.2, 2017, pp. 521–530  DOI: 10.1007/s10340-016-0793-0
  * [13] D. Cornara et al.  “Transmission of Xylella fastidiosa by naturally infected Philaenus spumarius (Hemiptera, Aphrophoridae) to different host plants”  In _Journal of Applied Entomology_ 141.1-2 Wiley Online Library, 2017, pp. 80–87  DOI: 10.1111/jen.12365
  * [14] D. Cornara et al.  “EPG combined with micro-CT and video recording reveals new insights on the feeding behavior of Philaenus spumarius”  In _PLoS One_ 13.7 Public Library of Science, 2018, pp. 1–20  DOI: 10.1371/journal.pone.0199154
  * [15] N.. Cunniffe et al.  “Thirteen challenges in modelling plant diseases”  In _Epidemics_ 10 Elsevier, 2015, pp. 6–10  DOI: 10.1016/j.epidem.2014.06.002
  * [16] M.. Daugherty and R… Almeida  “Understanding How an Invasive Vector Drives Pierce’s Disease Epidemics: Seasonality and Vine-to-Vine Spread”  In _Phytopathology_ 109.2, 2019, pp. 277–285  DOI: 10.1094/PHYTO-07-18-0217-FI
  * [17] O. Diekmann, J… Heesterbeek and M.. Roberts  “The construction of next-generation matrices for compartmental epidemic models”  In _Journal of The Royal Society Interface_ 7.47, 2010, pp. 873–885  DOI: 10.1098/rsif.2009.0386
  * [18] A. Fierro, A. Liccardo and F. Porcelli  “A lattice model to manage the vector and the infection of the Xylella fastidiosa on olive trees”  In _Scientific Reports_ 9.1, 2019, pp. 8723  DOI: 10.1038/s41598-019-44997-4
  * [19] J.. Freitag  “Host range of the Pierce’s disease virus of grapes as determined by insect transmission”  In _Phytopathology_ 41.10, 1951, pp. 920–934
  * [20] H. Ge, K. Xu and Z. Ghahramani  “Turing: a language for flexible probabilistic inference”  In _International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain_ , 2018, pp. 1682–1690  URL: http://proceedings.mlr.press/v84/ge18b.html
  * [21] À. Giménez-Romero  “A compartmental model for Xylella fastidiosa related diseases”  In _GitHub Repo_ , 2022  URL: https://github.com/agimenezromero/A-compartmental-model-for-Xylella-fastidiosa-related-diseases
  * [22] À. Giménez-Romero, R. Flaquer-Galmés and M.. Matías  “Vector-borne diseases with nonstationary vector populations: The case of growing and decaying populations”  In _Physical Review E_ 106 American Physical Society, 2022, pp. 054402  DOI: 10.1103/PhysRevE.106.054402
  * [23] V. Grimm and S.. Railsback  “Individual-based Modeling and Ecology”  Princeton (NJ): Princeton University Press, 2005
  * [24] M.. Homan and A. Gelman  “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo”  In _J. Mach. Learn. Res._ 15.1 JMLR.org, 2014, pp. 1593–1623  URL: https://dl.acm.org/doi/10.5555/2627435.2638586
  * [25] D.. Hopkins and A.. Purcell  “Xylella fastidiosa : Cause of Pierce’s Disease of Grapevine and Other Emergent Diseases”  In _Plant Disease_ 86, 2002, pp. 1056–1066  DOI: 10.1094/PDIS.2002.86.10.1056
  * [26] M.. Jeger and C. Bragard  “The Epidemiology of Xylella fastidiosa; A Perspective on Current Knowledge and Framework to Investigate Plant Host–Vector–Pathogen Interactions”  In _Phytopathology_ 109.2, 2019, pp. 200–209  DOI: 10.1094/PHYTO-07-18-0239-FI
  * [27] M.. Jeger, J. Holt, F. Van Den Bosch and L.. Madden  “Epidemiology of insect-transmitted plant viruses: modelling disease dynamics and control interventions”  In _Physiological Entomology_ 29.3, 2004, pp. 291–304  DOI: https://doi.org/10.1111/j.0307-6962.2004.00394.x
  * [28] M.. Jeger, L.. Madden and F Van Den Bosch  “Plant virus epidemiology: Applications and prospects for mathematical modeling and analysis to improve understanding and disease control”  In _Plant disease_ 102.5 Am Phytopath Society, 2018, pp. 837–854  DOI: 10.1094/PDIS-04-17-0612-FE
  * [29] M.. Jeger, F. van den Bosch, L.. Madden and J. Holt  “A model for analysing plant-virus transmission characteristics and epidemic development”  In _Mathematical Medicine and Biology: A Journal of the IMA_ 15.1, 1998, pp. 1–18  DOI: 10.1093/imammb/15.1.1
  * [30] C. Lago et al.  “Degree-day-based model to predict egg hatching of Philaenus spumarius (Hemiptera: Aphrophoridae), the main vector of Xylella fastidiosa in Europe”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2022  DOI: 10.1101/2022.11.10.515963
  * [31] “Mathematics and Biology: The interface, challenges and opportunities”  Lawrence Berkeley Lab., CA (United States), 1992  eprint: https://escholarship.org/uc/item/870966k3
  * [32] J. López-Mercadal, P. Mercadal-Frontera and M.. Miranda  “Mechanical management of weeds drops nymphal density of Xylella fastidiosa vectors”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2022  DOI: 10.1101/2022.10.18.512680
  * [33] J. López-Mercadal et al.  “Collection of data and information in Balearic Islands on biology of vectors and potential vectors of Xylella fastidiosa (GP/EFSA/ALPHA/017/01)”  In _EFSA Supporting Publications_ 18.10 Wiley Online Library, 2021, pp. 6925E  DOI: 10.2903/sp.efsa.2021.EN-6925
  * [34] L.. Madden, M.. Jeger and F. van den Bosch  “A Theoretical Assessment of the Effects of Vector-Virus Transmission Mechanism on Plant Virus Disease Epidemics”  In _Phytopathology_ 90.6, 2000, pp. 576–594  DOI: 10.1094/PHYTO.2000.90.6.576
  * [35] M. Martcheva  “An Introduction to Mathematical Epidemiology”  New York (USA): Springer, 2015  DOI: 10.1007/978-1-4899-7612-3
  * [36] E. Moralejo et al.  “Phylogenetic inference enables reconstruction of a long-overlooked outbreak of almond leaf scorch disease (Xylella fastidiosa) in Europe”  In _Communications Biology_ 3.1, 2020, pp. 560  DOI: 10.1038/s42003-020-01284-7
  * [37] J.. Murray  “Mathematical Biology”  Berlin (Germany): Springer, 1989
  * [38] D. Olmo et al.  “Landscape Epidemiology of Xylella fastidiosa in the Balearic Islands”  In _Agronomy_ 11.3, 2021, pp. 473  DOI: 10.3390/agronomy11030473
  * [39] A.. Purcell and A. Finlay  “Evidence for noncirculative transmission of Pierce’s disease bacterium by sharpshooter leafhoppers”  In _Phytopathology_ 69.4, 1979, pp. 393–395  DOI: 10.1094/Phyto-69-393
  * [40] C. Rackauckas and Q. Nie  “DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia”  In _The Journal of Open Research Software_ 5.1, 2017, pp. 15  DOI: 10.5334/jors.151
  * [41] R.. Redak et al.  “The biology of xylem fluid-feeding insect vectors of Xylella fastidiosa and their relation to disease epidemiology”  In _Annual Review of Entomology_ 49.1, 2004, pp. 243–270  DOI: 10.1146/annurev.ento.49.061802.123403
  * [42] Kimberlyn Roosa and Gerardo Chowell  “Assessing parameter identifiability in compartmental dynamic models using a computational approach: application to infectious disease transmission models”  In _Theoretical Biology and Medical Modelling_ 16.1, 2019, pp. 1  DOI: 10.1186/s12976-018-0097-6
  * [43] A. Saltelli, S. Tarantola, F. Campolongo and M. Ratto  “Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models”  Hoboken (NJ), USA: Wiley, 2004
  * [44] M. Saponari, D. Boscia, F. Nigro and G.. Martelli  “Identification of DNA sequences related to Xylella fastidiosa in oleander, almond and olive trees exhibiting leaf scorch symptoms in Apulia (Southern Italy)”  In _Journal of Plant Pathology_ 95.3 Edizioni ETS, 2013  DOI: 10.4454/JPP.V95I3.035
  * [45] M. Saponari, A. Giampetruzzi, G. Loconsole, D. Boscia and P. Saldarelli  “Xylella fastidiosa in olive in Apulia: Where we stand”  In _Phytopathology_ 109.2 Am Phytopath Society, 2019, pp. 175–186  DOI: 10.1094/PHYTO-08-18-0319-FI
  * [46] S. Sarkar et al.  “Biodiversity conservation planning tools”  In _Annual Review of Environment and Resources_ 31, 2006, pp. 123–59  DOI: 10.1146/annurev.energy.31.042606.085844
  * [47] I.. Sobol  “Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates”  In _Mathematics and Computers in Simulation_ 55.1, 2001, pp. 271–280  DOI: 10.1016/S0378-4754(00)00270-6
  * [48] S. Soubeyrand et al.  “Inferring pathogen dynamics from temporal count data: the emergence of Xylella fastidiosa in France is probably not recent”  In _New Phytologist_ 219.2, 2018, pp. 824–836  DOI: https://doi.org/10.1111/nph.15177
  * [49] J.. Stevenson, M.. Matthews, L.. Greve, J.. Labavitch and T.. Rost  “Grapevine Susceptibility to Pierce’s Disease II: Progression of Anatomical Symptoms”  In _American Journal of Enology and Viticulture_ 55.3 American Journal of EnologyViticulture, 2004, pp. 238–245  URL: https://www.ajevonline.org/content/55/3/238
  * [50] B.. Teviotdale and J.. Connell  “Almond Leaf Scorch”  In _ANR University of California_ , 2003, pp. 8106  DOI: 10.3733/ucanr.8106
  * [51] M. Vanhove et al.  “Genomic diversity and recombination among Xylella fastidiosa subspecies”  In _Applied and Environmental Microbiology_ 85.13 Am Soc Microbiol, 2019, pp. e02972–18  DOI: 10.1128/AEM.02972-18
  * [52] P. Virtanen et al.  “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”  In _Nature Methods_ 17, 2020, pp. 261–272  DOI: 10.1038/s41592-019-0686-2
  * [53] S.. White, J.. Bullock, D… Hooftman and D.. Chapman  “Modelling the spread and control of Xylella fastidiosa in the early stages of invasion in Apulia, Italy”  In _Biological Invasions_ 19.6, 2017, pp. 1825–1837  DOI: 10.1007/s10530-017-1393-5
  * [54] S.. White, J.. Navas-Cortés, J.. Bullock, D. Boscia and D.. Chapman  “Estimating the epidemiology of emerging Xylella fastidiosa outbreaks in olives”  In _Plant Pathology_ 69.8, 2020, pp. 1403–1413  DOI: https://doi.org/10.1111/ppa.13238

Appendices

##  Appendix I Vector population dynamics

Fig. S1 shows a time series for the population of Philaenus spumarius in
Mallorca, taken from ([33]) (in blue). Superimposed (in orange) is the
assumption used in our model Eq. 2, the $\delta(t-nT)$, i.e. every year
susceptible vectors appear in the system.

Figure S1: Vector dynamics produced by the model compared to field-data from
([33]).

##  Appendix II Determination of $R_{0}$

The handicap of determining the basic reproductive number of the model Eq. 2
is that the pre-pandemic fixed point given by $I_{H}=I_{v}=0$ and
$S_{H}=S_{H}(0)$ is not a fixed point of the system of differential equations,
because vector population decays, so that the standard methods to compute
$R_{0}$ such as the Next Generation Matrix ([17, 22]) do not apply. In ([22])
a method was suggested to determine the basic reproductive number in the case
of compartmental models of vector-borne transmitted diseases in which the
vector population grows or decays. It consists in averaging the instantaneous
basic reproductive number over the time of a generation.

To proceed we consider that $I_{H}=I_{v}=0$, $S_{H}=S_{H}(0)$ is indeed a
fixed point of the system. Then, the basic reproductive number could be
determined, e.g. as shown in ([6]). First, an infectious host infects vectors
at a rate $\beta S_{H}(0)/N_{H}$ for a time $1/\gamma$. This produces $\beta
S_{H}(0)/\gamma N_{H}$ infected vectors. The second stage is that these
infectious vectors infect hosts at a rate $\alpha N_{v}(0)/N_{H}$ for a time
$1/\mu$, producing $\alpha N_{v}/\mu N_{H}$ infected hosts per vector. The net
result of these two stages is

$\tilde{R}_{0}=\frac{\alpha\beta}{\mu\gamma}\frac{S_{H}(0)}{N_{H}^{2}}N_{v}(0)=R_{0}^{*}\cdot
N_{v}(0)\ .$ (1)

This result coincides with the value of $R_{0}$ obtained using the standard
NGM method, that can be applied in this case because we are assuming that we
use a nongeneric initial condition that sits at the fixed point of the model.

In practice, our initial condition will never be a fixed point of the model,
and, as mentioned above, we will obtain an approximate basic reproductive
number, to which we will refer as $R_{0}$ using the method suggested in
([22]), that consists in ([22]). calculating the average number of secondary
infections produced by an infected host in one generation. First defines an
instantaneous basic reproductive number,

$R_{0}^{(i)}(t)=\frac{\beta\alpha}{\mu\gamma}\frac{S_{H}(0)}{{N_{H}}^{2}}N_{v}(t)=R_{0}^{*}N_{v}(t)\
,$ (2)

from which the average is simply computed as

$R_{0}=\left<R_{0}^{(i)}(t)\right>\Big{\rvert}_{0}^{\tau}=R_{0}^{*}\left<N_{v}(t)\right>\Big{\rvert}_{0}^{\tau}=R_{0}^{*}\frac{1}{\tau}\int_{0}^{\tau}N_{v}(t)\,\mathrm{d}t\
,$ (3)

In our model, the time-dependent vector population can be obtained from Eq. 2,

$\dot{N}_{v}=\dot{S}_{v}+\dot{I}_{v}=-\mu N_{v}\Longrightarrow
N_{v}(t)=N_{v}(0)e^{-\mu t}\ ,$ (4)

and introducing this expression for $N_{v}(t)$ in Eq. 3 the integral can be
solved

$R_{0}={{\beta\alpha S_{H}(0)}\over{\mu\gamma
N_{H}^{2}}}\,\frac{N_{v}(0)}{\mu\tau}\left(1-e^{-\mu\tau}\right)=R_{0}^{*}\
\frac{N_{v}(0)}{\mu\tau}\left(1-e^{-\mu\tau}\right),$ (5)

that is an approximated expression to the basic reproductive number for our
model, in which the vector population is nonstationary, where, in Eq. 2 and
Eq. 5 it has been defined, $R_{0}^{*}=(\beta\alpha S_{H}(0))/(\mu\gamma
N_{H}^{2})$.

Note that in our model one generation correspond to one year and that
$N_{v}(0)$ is reset every year.