Using sensitivity analysis to identify key factors for the propagation of a plant epidemic

Identifying the key factors underlying the spread of a disease is an essential but challenging prerequisite to design management strategies. To tackle this issue, we propose an approach based on sensitivity analyses of a spatiotemporal stochastic model simulating the spread of a plant epidemic. This work is motivated by the spread of sharka, caused by plum pox virus, in a real landscape. We first carried out a broad-range sensitivity analysis, ignoring any prior information on six epidemiological parameters, to assess their intrinsic influence on model behaviour. A second analysis benefited from the available knowledge on sharka epidemiology and was thus restricted to more realistic values. The broad-range analysis revealed that the mean duration of the latent period is the most influential parameter of the model, whereas the sharka-specific analysis uncovered the strong impact of the connectivity of the first infected orchard. In addition to demonstrating the interest of sensitivity analyses for a stochastic model, this study highlights the impact of variation ranges of target parameters on the outcome of a sensitivity analysis. With regard to sharka management, our results suggest that sharka surveillance may benefit from paying closer attention to highly connected patches whose infection could trigger serious epidemics.


Supporting Information
Article title: Using sensitivity analysis to identify key factors for the propagation of a plant epidemic Authors: Loup Rimbaud, Claude Bruchou, Sylvie Dallot, David R. J. Pleydell, Emmanuel Jacquot, Samuel Soubeyrand and Gaël Thébaud The following Supporting Information is available for this article: Video S1. scenario B) Single PPV introduction through 1 infectious tree in a moderately connected patch Video S2. scenario C) Single PPV introduction through several infectious trees in a moderately connected patch.
Video S3. scenario D) Multiple PPV introductions through several infectious trees, with 1 st introduction in a weakly connected patch Video S4. scenario E) Multiple PPV introductions through several infectious trees, with 1 st introduction in a moderately connected patch (reference scenario) Video S5. scenario F) Multiple PPV introductions through several infectious trees, with 1 st introduction in a highly connected patch Figure S1. Landscape used in this study, composed of 553 patches in a peach-growing area in south-eastern France. It roughly consists in a 4.5×4.5 km square, with a total area of 524 ha (average patch area of 0.95 ha). Figure S2. Computation of the weights of the exponential dispersal functions. The weights are computed from the integration of the probability density function of the dispersal weighting variable (W) on J=20 sub-intervals of [0; 1]. In this example, Wexp=0.486 and Wvar=0.0434.
Dispersal weighting variable (W) Figure S3. Variation of probability distributions used as input of the model in the broad-range (a-c-e) and the sharka-specific (b-d-f) sensitivity analyses. a-b: Distribution of the probability for a tree to be infected in an orchard contaminated at planting (τ, mixture of beta distributions) for different values of the relative probability of massive introduction (pMI). c-d: Probability distribution of the dispersal weighting variable (W, beta distribution) for different expected values of the dispersal weighting variable (Wexp). e-f: Probability distribution of the latent period duration (gamma distributions) for different expected durations of the latent period (θexp). Figure S4. Variation of the mean dispersal distance in the broad-range sensitivity analysis, for different expected values of the dispersal weighting variable (Wexp). Inset: the corresponding dispersal kernels (mixture of exponential). The colour code matches that of Figure S1c. Figure S5. First-order and total Sobol's sensitivity indices of the 6 target parameters on the standard deviation of the output (σY, standard deviation of the equivalent number of fully productive trees per hectare and per year) of the stochastic replicates, assessed in the broadrange analysis (a, ∑ 1 = 0.32), or informed by knowledge acquired on sharka in peach orchards (b, ∑ 1 = 0.49). θexp: expected duration of the latent period; Φ: probability of introduction at orchard planting; β: transmission coefficient; pMI: relative probability of massive introduction; Wexp: expected value of the dispersal weighting variable; qκ: quantile of the outgoing connectivity of the patch of first introduction. Horizontal bars: CI95 (assessed using 10,000 bootstrap replicates). *: the estimates of the 1 st -order Sobol indices were slightly negative, which may occur when the indices do not significantly differ from 0 [49]. Figure S6. Sensitivity indices of the target parameters and their interactions up to the order 2 on the standard deviation of the output (σY, standard deviation of the number of fully productive trees per hectare and per year) of the stochastic replicates, assessed through a degree-3 polynomial optimised by BIC in the broad-range analysis (a, R 2 =0.48), or informed by knowledge acquired on sharka in peach orchards (b, R 2 =0.62). θexp: expected duration of the latent period; Φ: probability of introduction at orchard planting; β: transmission coefficient; pMI: relative probability of massive introduction; Wexp: expected value of the dispersal weighting variable; qκ: quantile of the outgoing connectivity of the patch of first introduction. Figure S7. Effect of each target parameter on the standard deviation of the output (σY, standard deviation of the equivalent number of fully productive trees per hectare and per year) in the broad-range analysis (a), or informed by knowledge acquired on sharka in peach orchards (b). Each point represents the standard deviation of 50 stochastic replicates of model simulations performed with the same combination of parameters. Blue line: prediction using the degree-3 polynomial optimised by BIC (all other parameters fixed at their mean value); this line is extremely close to the moving average of the number of productive trees, which is not represented for legibility. Green dashed line: standard deviation of the number of productive trees in the disease-free scenario. In (a) the thick line on the x-axis corresponds to the variation range used in the sharka-specific analysis; in (b) the full circle on the x-axis indicates the reference value of each parameter. Wexp: expected value of the dispersal weighting variable (higher values correspond to longer dispersal distances); β: transmission coefficient; θexp: expected duration of the latent period; qκ: quantile of the outgoing connectivity of the patch of first introduction; Φ: probability of introduction at orchard planting; pMI: relative probability of massive introduction. Figure S8. Robustness of the broad-range sensitivity analysis to changes in the upper bound of the expected duration of the latent period (from 5 to 33 years). The curves indicate 1 st -order (solid lines) and total (dashed lines) Sobol's sensitivity indices of the 6 target parameters with regard to the mean output (μY, mean number of fully productive trees per hectare and per year) of the stochastic replicates. The estimated 1 st -order index may be higher than the estimated total index when the two indices are close and the number of simulations used to compute the indices is low. Figure S9. Sharka-specific sensitivity analysis on the mean output (μY, mean of the equivalent number of fully productive trees per hectare and per year) but performed with a wide variation range (as in the broad-range analysis) for Wexp. a: 1 st -order and total Sobol's sensitivity indices of the 6 target parameters (∑ 1 = 0.88). Horizontal bars: CI95 (assessed with 100 bootstrap replicates). b: Effect of each target parameter. Red line: polynomial regression; green dashed line: number of productive trees in the disease-free scenario; full circle on the x-axis: reference value of each parameter. Wexp: expected value of the dispersal weighting variable; θexp: expected duration of the latent period; Φ: probability of introduction at orchard planting; β: transmission coefficient; pMI: relative probability of massive introduction; qκ: quantile of the outgoing connectivity of the patch of first introduction.

Methods S1. Relation between the mean and standard deviation of the latent period duration
In order to assess the relation between the mean and standard deviation of the latent period duration, we first studied the analogous relation for the incubation period, for which more data are available. Many studies dealing with the incubation period of infectious diseases suggest a correlation between the mean or median duration of the incubation period and its variance. In order to investigate this hypothesis, we estimated a proportionality coefficient between the mean (μ) and standard deviation (σ) of the incubation period duration for several human, animal, and plant infections.
When the duration of the incubation period was described by a lognormal distribution [1][2][3], we used the following approach: Then, μ and σ were estimated using the equations previously presented.
When the incubation period was described by a gamma distribution [5,6], the mean and standard deviation were calculated from the shape (α) and rate (β) parameters using: = / and = √ / .
We then obtained proportionality coefficients from 0.12 to 0.86 between μ and σ (Table S1). Additionally, in a study focused on a plant fungus [6], age-dependent incubation periods were measured and showed a significant correlation between the mean and standard deviation of this period, with an estimate proportionality coefficient of 0.23 ( Figure S8).
In our model, the latent period is considered synchronised with the incubation period. The proportionality coefficient of 0.35 estimated in [7] is within the range of the coefficients estimated for the incubation period of other infections. Therefore, in each simulation of our model, the standard deviation of the latent period duration is calculated by multiplying its mean by 0.35. Figure S10. Relation between the mean (μ) and standard deviation (σ) of the incubation period duration in Rhizoctonia solani. Data: [6]. Inset: relation for mean incubation periods between 0 and 2, which are better spread along the x-axis. Black lines: linear regressions. Table S1. Proportionality coefficient between the mean (μ) and standard deviation (σ) of the incubation or latent period durations in several human, animal, and plant infections.

Methods S2. Bounds for the expected value of the latent period duration and the transmission coefficient in the broad-range sensitivity analysis
The broad-range sensitivity analysis (i.e. the analysis designed to assess the intrinsic influence of the target parameter on model behaviour), required to define appropriate minimal and maximal bounds of the variation ranges of parameters whose definition domain is not finite. Thus, we tried to find values below and above which the epidemic process does not fundamentally change any more. For example, the mean duration of the latent period (θexp) can be left-bounded by 1 day. Indeed, lower values would result in similar processes, since the time step of the model is one week. The estimate of the maximal bound is less easy. We assumed that the epidemic process is unchanged for θexp values above , for which only 5% of the infected hosts at orchard planting become infectious before orchard removal. This value was calculated as follows: Let Ψ describe the duration of orchard cultivation; Ψ~Poisson( = 15).
Let Θ describe the duration of the latent period; Θ~( , 0.35 × ) (parameterising the Gamma distribution with its expectation and variance; see also Methods S1).
We seek such that: (Θ ≤ Ψ) equals 5.43% with θexp=32; 4.77% with θexp=33; and 4.19% with θexp=34. The minimal bound (β=0.08) corresponds to the value for which = 2 after n=30 years, and the maximal bound (β=12.80) corresponds to the value for which = 0.95 × 0 after n=3 years (i.e. 95% of the orchard is infectious when the orchard just becomes productive; in such extreme scenario, prunus cultivation would be economically unfeasible). The broad-range sensitivity analysis confirmed that the epidemic process does not change anymore beyond this maximal value, as shown by the stabilisation of the output variable for β values beyond 5.00 ( Figure 4a).

Methods S3. Computation of Sobol's sensitivity indices
Let = ( 1 , … , ) be a function of a set of input parameters Xi (i=1,…,p) defined on a pdimensional cube. Assuming that these input parameters are independent, the function can be decomposed into uncorrelated terms of increasing dimensionality [1][2][3], Videos. Simulations of different scenarios of sharka epidemics in a French peach production area in the absence of disease control (one year every 0.5 second). Each polygon represents a patch on which orchards are cultivated, with a turnover of 1/15 year -1 . Orchard colour represents the proportion of healthy trees, from 100% (white, all trees in the S state) to 0% (red, all trees in the E or I states). The green area is not cultivated with prunus. Except the introduction parameters, all model parameters were fixed at their reference value. 'Single introduction' means only one pathogen introduction at year 1, whereas 'multiple introductions' means a first introduction at year 1 and a risk of subsequent introductions at every orchard planting.
Video S1. scenario B) Single PPV introduction through 1 infectious tree in a moderately connected patch.
Video S2. scenario C) Single PPV introduction through several infectious trees in a moderately connected patch.
Video S3. scenario D) Multiple PPV introductions through several infectious trees, with 1 st introduction in a weakly connected patch.
Video S4. scenario E) Multiple PPV introductions through several infectious trees, with 1 st introduction in a moderately connected patch (reference scenario).
Video S5. scenario F) Multiple PPV introductions through several infectious trees, with 1 st introduction in a highly connected patch.