An Epidemiological Framework for Modelling Fungicide Dynamics and Control

Defining appropriate policies for controlling the spread of fungal disease in agricultural landscapes requires appropriate theoretical models. Most existing models for the fungicidal control of plant diseases do not explicitly include the dynamics of the fungicide itself, nor do they consider the impact of infection occurring during the host growth phase. We introduce a modelling framework for fungicide application that allows us to consider how “explicit” modelling of fungicide dynamics affects the invasion and persistence of plant pathogens. Specifically, we show that “explicit” models exhibit bistability zones for values of the basic reproductive number () less than one within which the invasion and persistence threshold depends on the initial infection levels. This is in contrast to classical models where invasion and persistence thresholds are solely dependent on . In addition if initial infection occurs during the growth phase then an additional “invasion zone” can exist for even smaller values of . Within this region the system will experience an epidemic that is not able to persist. We further show that ideal fungicides with high levels of effectiveness, low rates of application and low rates of decay lead to the existence of these bistability zones. The results are robust to the inclusion of demographic stochasticity.


Introduction
Fungicide use is an essential aspect of disease management in modern agriculture [1]. There are, however, increasing restrictions upon the use of fungicides (and other forms of chemical control) [2,3] and the trend is for reduced usage, to avoid undue release to the environment or to minimise the risk of fungicide resistance. Reduction in the availability and use of fungicides imposes greater demands upon efficient use of the limited resources available. One way to approach this is by the use of mathematical models to investigate the effects of a given fungicide on the host crop and the pathogen to aid the design of appropriate disease management strategies [4][5][6]. The benefits of modelling diseases are well acknowledged: it allows theoretical optimal control strategies to be developed that minimise economic costs and maximise crop returns [7][8][9][10], which can then be tested experimentally [11]. Furthermore the continued widespread use of fungicides is threatened by the emergence of resistant pathogen strains, often as a direct consequence of the application strategy itself [1,2,12] and in order to implement effective resistance management strategies it is often necessary to model those strategies first [9,[13][14][15][16].
Most models that include fungicides often account for fungicide dynamics by simply modifying the parameters of the underlying epidemiological models (i.e. reducing the infection rates and/or increasing the host recovery rates) [9,14,17]. Much work has been done to investigate the consequences of using different underlying models for both the host population and pathogen dynamics (see [6] for a review) but relatively little work has been carried out to investigate the dynamics of the fungicides themselves and how the timing of initial infection relative to host population growth affects invasion and persistence in a chemically controlled system. Generally, previous work has implicitly made one or more of the following assumptions; that there is complete coverage by the fungicide for either the entire host population or a fixed subset of the host population [10,14,16,18,19]; that a generic, multipurpose, fungicide has been applied to the hosts (i.e. they have not readily separated out the effects of different fungicide types such as protectants, curatives or eradicants) [9,20,21]; or that the fungicides are permanent (i.e. that the chemicals do not decay or that their effects do not change over time) [18].
Here we propose a parsimonious model framework that is designed to integrate epidemiological and fungicide dynamics. Specifically we consider purely protectant fungicides, and we distinguish the effects of the rate of application, the decay in activity and the partial effectiveness of the fungicide on the epidemic dynamics. We first use the framework to construct a deterministic compartmental model that incorporates host growth and explicitly allows for the dynamic application of a purely protectant fungicide through the use of multiple susceptible host classes. We compare this model with a conventional model that implicitly incorporates the protectant fungicide dynamics through modification of the underlying infection rates for a single susceptible host class. For brevity, we refer to these two models as the explicit and implicit models respectively. We then extend the explicit model to include demographic stochasticity.
Given that the application of fungicides is usually designed to prevent invasion of a pathogen or to eliminate the pathogen, we use the two models to ask the following questions: N does the explicit inclusion of fungicide dynamics affect the criteria for invasion and persistence? N does the timing of initial infection affect the criteria for invasion and persistence in a chemically controlled system? N how are the differences manifested? N does failure to account for fungicide dynamics lead to erroneous predictions of control effectiveness in epidemiological models?
N are there critical regions of parameter space concerned with coverage, decay and effectiveness of fungicides for which these differences are most exaggerated?
N are the differences maintained when allowance is made for stochastic variability in the transmission of infection?

Deterministic Model Derivation
Consider a pathogen spreading through a population of hosts such as a single field of a crop, in which the epidemiological host unit is an amount of plant tissue i.e. a leaf, stem, or root [4] and there is a simple density dependent growth of the host up to a carrying capacity. In the absence of any form of chemical control the population is divided into two classes, susceptible (Ŝ S) and infected (Î I). We assume that infected host tissue consumes resources but does not contribute to new host growth [22], that infection is adequately described as a mass action process (this assumption holds well for pathogens within areas smaller than their average dispersal distance [23]) and that hosts cannot be reinfected after they cease to be infectious, and are removed. We investigate the effects of applying a purely protectant fungicide to the system. Here we assume that a protectant fungicide affects susceptible hosts, reducing their capacity to become infected.
We first consider a conventional model that only includes implicit fungicide dynamics, we assume that all susceptible hosts are protected. All factors affecting the effectiveness of the protectant are subsumed into a single parameter (r) which represents the average decrease in the infection rate for the host population. The dynamics of this model are represented in Figure 1 and by the equations: Here, density dependent host growth is governed by the ratel l and a carrying capacityk k.b b represents the infection rate in the absence of chemical control. The parameter r (0ƒrƒ1) is a measure of fungicide effectiveness. It reduces the infection rate and represents the combined effects of the partial coverage, decay in activity and incomplete effectiveness of the protectant fungicide (possible expressions for r that reflect different combinations of these effects are summarised in Table 1). Infected hosts are removed at a per capita rate (m m).
We choose an expression for r that takes account of the uninfected population being an aggregation of completely unprotected and partially protected hosts. Given a protectant application rate (p p) and a protectant decay rate (d d) the long term  Table 2 As such we take r~d dzep p d dzp p , which yields our conventional, 'implicit' model: dÎ I dt t~d In order to incorporate fungicide dynamics explicitly we propose an alternative model with an additional protected host class (P P). Here, the per capita protectant application rate (p p) and the per capita protectant decay rate (d d) allow individuals to move between the susceptible (Ŝ S) and protected (P P) states. The protected state experiences reduced infection rates (eb b) whilst the unprotected susceptible state experiences normal infection rates (b b). New host units are assumed to be unprotected which is consistent with a non-systemic protectant fungicide. The dynamics of this model are represented in Figure 1 and by the equations: dÎ I dt t~b Non-dimensionalisation. We introduce the dimensionless variables and parameters The conventional model equations (1) and (2) with implicit fungicide dynamics become Implicit Model and the alternative model equations (5), (6) and (7) with explicit fungicide dynamics become Henceforth we refer to the conventional model as the implicit model and the alternative model as the explicit model. The parameters, their definitions, and typical values used for the numerical simulations are summarised in Table 2. The default values used hereafter are chosen primarily to highlight the characteristic properties of the model system, however the parameter ratios are consistent with those used in disease management practices for both barley and wheat crops [24][25][26][27].

Stochastic Model
In order to demonstrate that existence of the bistability zone is robust to demographic stochasticity and that it is not just a property of the deterministic nature of the explicit model, a continuous-time Markov process version of the explicit model is constructed.
AgainŜ S,Î I andP P represent the actual numbers of hosts in each of the susceptible, infected and protected classes. Given a carrying capacity,k k, the infinitesimal transition probabilities that describe the Markov process are given in Table 3. For computational efficiency a continuous-time Gillespie algorithm [28] is used to generate the sequence of transition event times for each simulation and thus obtain the trajectories for each class.

Equilibrium Analyses
Both the implicit model and the explicit model have the same basic reproductive number for the pathogen: with the disease free equilibria (S ?:0 ,I ?:0 ) and (S ?:0 ,P ?:0 ,I ?:0 ) for the models given by It is trivial to show that both disease free equilibria are stable for R 0 v1 and unstable for R 0 w1. The implicit model only admits a single endemic equilibrium (S ? ,I ? ), given by: whereas the explicit model admits multiple endemic equilibria, given by: where I ? is any root to the cubic equation: that satisfies I ? [½0,1.
It can be shown that the endemic equilibrium for the implicit model is both biologically meaningful and stable only for R 0 w1. This model gives rise to a classical epidemiological bifurcation diagram (a graph of I ? as a function of R 0 showing a transcritical bifurcation at R 0~1 ) given by Figure 2 and we see a single invasion threshold at R 0~1 as expected.
Allowance for fungicide dynamics in the explicit model yields zero, one, or two biologically realistic endemic equilibria depending on the values of the model parameters e, p, d and l.
In the region where R 0 w1 there is only ever one single, stable, endemic equilibrium. However, in the region where R 0 v1 the other parameters affect the number of possible equilibria. Let If cw0 then there is no biologically meaningful (stable or unstable) endemic equilibrium for R 0 v1. This again gives rise to a classical epidemiological bifurcation diagram. If, however, cv0, then a second threshold R c exists and it is possible for the system to have two endemic equilibria for R 0 between R c and 1 (Figure 2). The lower branch of the endemic equilibrium curve is always unstable and the upper branch stable [29]. It follows that failure to account explicitly for fungicide dynamics can lead to an erroneous understanding of the nature of the system. In particular, explicitly accounting for fungicide dynamics means that for values of R 0 less than 1, it is still possible for infection to persist within the system.
The threshold c~0 can be interpreted as a surface in the fungicide parameter (pde) space with l fixed (see Figure 3). Here we note that if the fungicide parameters lie in the region below this surface then the model exhibits bistability. The surface encloses the region of fungicide parameter space near to the origin, i.e. low values of the scaled parameters (p, d and e). It is a key result to note that these values correspond to highly effective fungicides (low e), with long lifespans (low d) that are applied infrequently (low p) and so it follows that a move towards using fungicides with these properties can have extremely undesirable consequences in terms of effective control.

Invasion and Persistence Criteria
We now consider how invasion and persistence criteria depend on the choice of model and on the host growth state of the system at the time of initial infection. In particular we determine how these criteria depend on both the basic reproductive ratio R 0 and the initial levels of infection (I 0 ). Invasion is related to the immediate behaviour of the system, and a pathogen is considered to have invaded if there is an immediate increase in the initial   Both models assume the same logistic growth of uninfected hosts up to a disease free equilibrium carrying capacity. However for the explicit model the relative levels of susceptible and protected hosts will vary during this growth phase (see Figure 4). Analytic solutions exist to both of these sets of equations.
Infection Choice. For each value of I 0 we must choose which uninfected hosts become infected. For the implicit model there is a single uninfected class and so the initial infected individuals must come from this class. If, prior to infection, the system has S i susceptible individuals, where 0vS i v1 then immediately post infection the system must be in the state: For the explicit model, with two uninfected classes, a choice is required. We assume that the proportion of I 0 that comes from each uninfected class is proportional to their relative susceptibilities (i.e. more initial infections come from the susceptible class than the protected class) and their relative densities. i.e. If, prior to infection, the system is in state (S i ,P i ,0) and immediately post infection the system is in state (S 0 ,P 0 ,I 0 ) then proportion of I 0 from susceptible class~S i {S 0 ! 1|S i proportion of I 0 from protected class~P i {P 0 ! e|P i Consequently, for an initial level of infection, I 0 , the initial post infection state (S 0 ,P 0 ,I 0 ) is given by: This gives us a method for moving from a pre-infection system state (S i ,0) or (S i ,P i ,0) to an initial post-infection state (S 0 ,I 0 ) or (S 0 ,P 0 ,I 0 ) for the implicit and explicit models respectively.
Invasion and Persistence Thresholds. For a pathogen to invade we require that the level of infection increases immediately from the post-infection system state. For the implicit model this requires We use the relationship between pre-infection levels and postinfection levels to create a threshold in R 0 {I 0 parameter space for each pre-infection susceptible host level S i . The threshold is given by: This reduces to the classic invasion threshold if S i is taken to be the disease free equilibrium level S ?:0~1 . However for the explicit model we require: This gives the following threshold in (R 0 {I 0 ) parameter space for each pre-infection host state (S i ,P i ): In the absence of analytical tractability persistence thresholds for the two models were determined numerically for a range of pre-infection states. The parameter values used for the numerical calculation are summarised in Table 2.
Using the results from the growth curves in section 0 and Figure 4 we illustrate how both the invasion and persistence thresholds depend upon the inclusion of fungicide dynamics for two distinct pre-infection states: the disease free equilibrium (points E ex and E im in Figure 4) and a pre-infection state during where the uninfected hosts are still growing (points G ex and G im in Figure 4). The thresholds obtained are shown in Figure 5. We conclude that the explicit inclusion of fungicide dynamics leads to a zone of bistability for values of R 0 less than 1, and that this zone is not simply an artefact of the equilibrium analysis but is robust to the system's transient dynamics. Furthermore we show that for preinfection states where the host population is still growing the zone of bistability is extended and lower initial levels of infected hosts (I 0 ) are able to invade and persist within the same range of R 0 values. In addition an extra invasion zone now exists for a range of values of R 0 below 1. Within this invasion zone, a pathogen is able to invade but not persist (see Figure 5).
Deterministic Realisations. Invasion trajectories for both models are shown in Figure 6. For initial conditions at point A (R 0~0 :4, I 0~0 :06), a value of R 0 is chosen so that R i vR 0 vR c for the explicit model lies within the invasion zone for a preinfection state with growing hosts. By considering the trajectories that result from this point it can clearly be seen that for the implicit model the choice of pre-infection state affects the pathogen's ability to invade. For initial conditions at points B (R 0~0 :8, I 0~0 :06) and C (R 0~0 :8, I 0~0 :05), a value of R 0 is chosen so that R c vR 0 v1 for the explicit model parameters. By considering the trajectories that result from these two points it can clearly be seen that for the implicit model the choice of initial infection level does not affect the ability of the pathogen to persist, whereas for the explicit model the initial infection level does affect the ability of the pathogen to persist. For initial conditions at point D (R 0~1 :8, I 0~0 :05) we can see that both models predict invasion for R 0 w1, but that they do not agree on the final endemic level of infection, with the explicit model predicting higher levels.
Stochastic Realisations. Using a nominal carrying capacity of 2000 (k k~2000), 5000 simulations are performed for every initial condition in the I 0 {R 0 parameter space (I 0 [ 0,0:01, . . . ,1 ½ , R 0 [ 0,0:05, . . . ,2:50 ½ , starting from a disease-free equilibrium level of uninfected hosts only) and the proportion of simulations that result in pathogen persistence is recorded and shown in Figure 7. The bistability zone still exists when allowance is made for demographic stochasticity. Figure 8 shows probability density plots obtained from the 5000 simulations for the initial conditions corresponding to points B (I 0~0 :06, R 0~0 :8), C (I 0~0 :05, R 0~0 :8) and D (I 0~0 :05, R 0~1 :8). It can be seen that for a stochastic framework, starting at initial conditions near to the deterministic threshold, individual trajectories may tend towards either equilibrium (endemic or disease free) regardless of whether the initial condition is above or below the deterministic persistence threshold and so we see distinctly bimodal probability density functions as a result. This is due to demographic stochastic effects allowing individual trajectories to enter different basins of attraction and so tend towards either equilibrium. Crucially, this means that by creating a stochastic version of the explicit model we have shown that not only is the bistability zone still present but the range of initial infection levels that can lead to persistence is extended. The robustness of the effects was tested for a range of population sizes (k k[ 200, . . . ,20000 ½ ) and shown to be sustained.

Discussion
In this paper we develop a simple model to investigate the effects of fungicide dynamics on pathogen invasion and persistence. Specifically we incorporate protectant fungicide dynamics into a conventional model in a way that separates out the effects of partial coverage, incomplete effectiveness and decay in activity of a protectant fungicide to create an alternative model that takes  Figure 4): for given R 0 and I 0 values the plot shows the long-term behaviour of both models in terms of the ability of a pathogen to invade and persist. For the implicit model R 0 completely characterises the long-term behaviour of the system (with R 0~1 being the threshold). For the explicit model, the ability of a pathogen to invade and persist is also determined by initial level of infection, I 0 , for a range of R 0 values (R c vR 0 v1). (B) Pathogen invades system during growth phase (points G ex and G im in Figure 4): For the implicit model R 0 still completely characterises the long-term behaviour of the system. For the explicit model the bistability zone is extended and lower initial levels of infected hosts (I 0 ) are able to invade and persist within the same range of R 0 values (R c vR 0 v1). In addition, an extra invasion zone now exists for a range of values of R 0 (R i vR 0 v1) where here R i vR c . Within this zone the pathogen is able to invade but is not able to persist. Initial conditions for selected numerical simulations are shown. A:(R 0~0 :4, I 0~0 :06) is within the invasion zone for the explicit model starting during the growth phase. B:(R 0~0 :8, I 0~0 :06) is just above the invasion and persistence threshold for the explicit model staring from the disease free equilibrium. C:(R 0~0 :8, I 0~0 :05) is just below the same threshold. D:(R 0~1 :8, I 0~0 :05) is significantly above the invasion and persistence thresholds for both models. The default parameter values are given in Table 2. doi:10.1371/journal.pone.0040941.g005 explicit account of the fungicide dynamics. Previous work has utilised simple models of fungicide application to investigate a number of issues e.g. fungicide resistance [9,13,19,30] and optimal control strategies [8][9][10], but the models used have not incorporated explicit fungicide dynamics (as the convention is to take implicit account of fungicide dynamics by mapping the effect onto the transmission rate).
The allowance for explicit fungicide dynamics markedly changes the inferences on both the invasion and persistence of the pathogen. Both frameworks appear superficially to have the same epidemiological properties and thresholds (R 0 is identical for both) but the explicit model exhibits a bistability zone for values of R 0 less than 1, a range that the implicit model considers completely safe from invasion. Previous work on human and animal vaccination models [29,31], re-infection models [32] and models of sexually transmitted diseases (with high and low transmissibility groups) [33] have exhibited similar bistability properties, but this is the first time that this property has been demonstrated for chemical control in agricultural systems. In addition previous epidemiological research has only considered the existence of a bistability zone for pathogens invading a system already at equilibrium. The extension of this result to systems where a pathogen invades a growing host state is a novel result with some striking repercussions. Not only does the inclusion of fungicide dynamics affect the criteria for invasion and persistence, (increasing the risk of a pathogen successfully invading and persisting) when compared with models that only implicitly include fungicide dynamics, but if the initial infection occurs at a time before the host population has reached its equilibrium state then this risk is exacerbated. Now it is even more likely that a pathogen will be able to invade and persist, and also there is an opportunity for a pathogen to simply invade in the short term, causing an early epidemic before eventually dying out.
The existence of both a bistability zone and an invasion zone below the traditionally accepted invasion threshold has several important implications both in the use of models to guide disease control, and in the selection of fungicide traits to promote effective disease control. Firstly, for a non-negligible, initial level of infection, a pathogen may be able to invade a system with a fungicidal control regime that an implicit model would predict to be adequate. For example, protectant application rates, fungicide efficacies and fungicide activity decay rates may be chosen to reduce the R 0 value of the system below 1 yet unknowingly still remain above the critical thresholds R 0~Rc , or R 0~Ri . The system lies either in the bistability zone and is at risk from invasion by a high enough level of initial infection or it lies in the invasion zone and is at risk from an invasion during the host growth phase. The second implication is that for values of R 0 just above 1, the two models predict very different endemic levels of infection; for the implicit model, having an R 0 value just above 1 would always result in an invading pathogen achieving a negligible endemic equilibrium (see Figure 2) whereas for the explicit model an equivalent R 0 value will always lead to a much higher endemic equilibrium (see Figure 2). Consequently it would be possible for a system to experience sudden large invasions from small initial infections as a result of a small change in value of R 0 (from just  It can be seen that persistence is still determined by I 0 for a range of R 0 values (R c vR 0 v1) and furthermore that simulations with initial conditions below the deterministic threshold are still able to persist (see inset). Initial conditions for selected numerical simulations are shown. A:(R 0~0 :8, I 0~0 :06) is just above the deterministic invasion and persistence threshold for the explicit model. B:(R 0~0 :8, I 0~0 :05) is just below the same threshold. C:(R 0~1 :8, I 0~0 :05) is significantly above the invasion and persistence thresholds for both models. The default parameter values are given in Table 2. doi:10.1371/journal.pone.0040941.g007 below 1 to just above 1). Finally, in order to control an existing outbreak, it is necessary to reduce R 0 not just to below 1 but to below the critical threshold value R c . The predicted effort required to eliminate an established pathogen from a system therefore depends on the choice of model. A model without explicit fungicide dynamics will underestimate the effort required and consequently allow the pathogen to persist.
Of particular concern for disease management strategies is the drive in the agrochemical industry to produce longer lasting, more effective fungicides, which necessitate lower application rates. These aspirational fungicide properties are exactly those that correspond to the parameter region for the explicit model in which bistability zones occur (Figure 3).
A continuous-time, Markov process version of the explicit model is constructed and used to demonstrate that the bistability result is not just a property of the deterministic models. The simulations show that the probability distributions for the stochastic model exhibit bimodal behaviour in the region of bistability, with one mode at the endemic equilibrium and one at the disease free equilibrium. It is clear that the effects of fungicide dynamics on invasion and persistence criteria are still maintained when allowance is made for stochastic variability in the various transmission processes.
Taking all of this into account it is reasonable to conclude that failure to account for fungicide dynamics naturally leads to erroneous predictions of control effectiveness.   Figure 7. Here the host carrying capacity,k k is 2000. The left hand figures show the probability density plots for all classes with initial conditions given by Point A (Î I 0~1 20, R 0~0 :8) i.e. just above the deterministic invasion threshold, the middle column shows the probability density plots for Point B (Î I 0~1 00, R 0~0 :8), just below the invasion threshold, and the right hand column shows the probability density plots for Point C (Î I 0~1 00, R 0~1 :8), which is significantly above the invasion threshold. The deterministic trajectory for the corresponding class and initial condition is overlaid on each plot. The numbers within each plot indicate the proportion of realisations that reach those final distinct host levels. Note that the density scale is non-linear. The default parameter values are given in Table 2. doi:10.1371/journal.pone.0040941.g008