A Generic Model Accounting for the Interactions among Pathogens, Host Plants, Biocontrol Agents, and the Environment, with Parametrization for Botrytis cinerea on Grapevines

Although the use of biocontrol agents (BCAs) to manage plant pathogens has emerged as a sustainable means for disease control, global reliance on their use remains relatively insignificant and the factors influencing their efficacy remain unclear. In this work, we further developed an existing generic model for biocontrol of foliar diseases, and we parametrized the model for the Botrytis cinerea–grapevine pathosystem. The model was operated under three climate types to study the combined effects on BCA efficacy of four factors: (i) BCA mechanism of action, (ii) timing of BCA application with respect to timing of pathogen infection (preventative vs. curative), (iii) temperature and moisture requirements for BCA growth, and (iv) BCA survival capability. All four factors affected biocontrol efficacy, but factors iii and iv accounted for > 90% of the variation in model simulations. In other words, the most important factors affecting BCA efficacy were those related to environmental conditions. These findings indicate that the environmental responses of BCAs should be considered during their selection, BCA survival capability should be considered during both selection and formulation, and weather conditions and forecasts should be considered at the time of BCA application in the field.


Introduction
Biocontrol of plant pathogens has emerged as a sustainable method of disease management and as a viable way to reduce the application of chemicals in agriculture [1][2][3][4]. The reasons for increasing restrictions on the use of chemicals and for increasing interest in biocontrol include the negative effects of chemicals on human health and the environment [5,6], pathogen-acquired resistance to commonly applied chemicals, and the lack of replacement products [7]. Biocontrol involves the use of fungi, bacteria, yeasts, or viruses (together referred to as biocontrol agents or BCAs) that may suppress plant pathogens via competition for nutrients or space, antibiosis, parasitism, and induced host plant resistance [1].
Despite the extensive research on biocontrol and the potential of using BCAs as alternatives to chemicals, the global reliance on BCA use remains relatively insignificant [4]. Many BCAs have been reported to suppress plant pathogens under controlled conditions in laboratories and greenhouses, but only a few have performed consistently in the field [8,9]. A possible reason for the lack of  [38]. The core of the model is based on a classic susceptibleinfected-removed (SIR) model, with tissue evolving from healthy-susceptible (HS) to infectious (I) and removed (R). The rate of infection of tissue (RI) depends on primary (STARTP) and secondary infections (I). The rate of resistance induction by a BCA (RRES) depends on BCA application (STARTB) and the total amount of healthy-susceptible tissue (HS). The rates of BCA colonization (RCOLH, RCOLI, RCOLR, and RCOLHr) depend on BCA application (STARTB) and the total amount of colonized tissue (BSUM). The structure incorporates host growth (RG) and physiological senescence (RS). Symbols for state variables, rates, and parameters are explained in Table 1.  , and the environment. The diagram uses the symbols developed by Forrester [38]. The core of the model is based on a classic susceptible-infected-removed (SIR) model, with tissue evolving from healthy-susceptible (HS) to infectious (I) and removed (R). The rate of infection of tissue (RI) depends on primary (STARTP) and secondary infections (I). The rate of resistance induction by a BCA (RRES) depends on BCA application (STARTB) and the total amount of healthy-susceptible tissue (HS). The rates of BCA colonization (RCOLH, RCOLI, RCOLR, and RCOLH r ) depend on BCA application (STARTB) and the total amount of colonized tissue (BSUM). The structure incorporates host growth (RG) and physiological senescence (RS). Symbols for state variables, rates, and parameters are explained in Table 1.

State Variables and Connecting Flows
The site of the system consists of K units of plant tissue that can be potentially occupied (i.e., affected) by the pathogen during the epidemic. The K units represent the state variables of the model and belong to one of the following non-overlapping categories of tissue: (i) healthy and susceptible to infection (HS); (ii) affected by the pathogen and infectious, i.e., can generate new, secondary infections (I); (iii) affected by the pathogen and removed, i.e., no longer infectious (R); (iv) healthy and colonized by the BCA, i.e., resistant to infection by the pathogen (H r ); (v) healthy and colonized by the BCA, i.e., which is protected from the pathogen (H b ); (vi) infectious and colonized by the BCA, i.e., unable to generate new infections (I b ); and (vii) removed and colonized by the BCA (R b ). The seven state variables are mutually exclusive so that: (1) The model considers that, during the epidemic and as a consequence of BCA application, the K units move from one state variable to another by means of rates.
At the beginning of a simulation, all of the plant tissue is in the state variable HS. The size of HS is dynamic and increases over time as a consequence of plant growth (in such a way that HS = 1 at the time of maximum plant size) or decreases as a consequence of senescence (which is relevant for those diseases in which the senescent plant tissue is no longer susceptible to infection). Inflow (rate of growth, RG) and outflow (rate of senescence, RS) of host tissue with respect to HS is calculated as follows: in which t is the current day, t−1 is the day before, and RRG t and RRS t are relative rates of host growth and senescence on day t, respectively. The host tissue in the state variable HS moves to state variable I as a consequence of infection by the pathogen; this flow is regulated by RI, the rate of infection, which is calculated as follows: (4) in which STARTP is the initial inflow of the pathogen into the system, b is the relative rate of infection, I is as previously defined, and COFR is the correction factor for occupied tissue.
In Equation (4), STARTP is calculated by assuming that the pathogen enters the system starting on day PIN (the day of the first seasonal infection) and continues to enter at a constant rate RPIN for a period of PDUR days; PIN, RPIN, and PDUR are all model parameters that are defined for each situation.
In Equation (4), COFR is calculated as follows: in which K and HS are as previously defined. The host tissue in the state variable I moves to state variable R when the infectious period (i.e., the period during which the pathogen continues producing inoculum on affected tissue) is over; this outflow is regulated by RR, the rate of removal, which is calculated as follows: (6) in which h is the relative rate of removal and I is as previously defined. The model considers that, at any time during the simulation period, a BCA enters the system because of human intervention (i.e., a treatment with the BCA); this can be before, at the same time as, or after the pathogen. The BCA inflow is regulated by STARTB, which is calculated for a period of BDUR days (i.e., the period during which the BCA is applied), starting from day BIN (i.e., the day on which the BCA is applied) at a constant rate equal to RBIN; BIN, RBIN, and BDUR are all model parameters that are defined for each situation.
The introduction of the BCA generates outflows from HS, so that the healthy tissue cannot be infected by the pathogen and, therefore, cannot move to I. The model considers that this outflow can be caused by BCAs that induce resistance in the host tissue and/or that prevent infection due to competition and/or antibiosis.
For BCAs that induce resistance, the outflow from HS (named RRES) is calculated as follows: RRES t = c 0 t × H r t−1 × COFRH r t−1 (7) in which c 0 is the relative rate of change from HS to H r and COFRH r is the correction factor for plant resistant tissue and is calculated as follows: in which K and H r are as previously described.
For BCAs that prevent infection by the pathogen, the outflow from H r (named RCOLH r ) is calculated as follows: (9) in which c 1 is the relative rate of change from H r to H b ; BSUM is the total of the tissue colonized by the BCA (i.e., BSUM = H b + I b + R b ); and STARTB and COFRH r are as previously defined.
Since induction of resistance in the host tissue is transitory, the model considers that the H r tissue can go back to HS and become susceptible to infection. The flow from H r to HS is calculated as follows: RSUS t = e t × H r t−1 (10) in which e is the relative rate of change from H r to HS and H r is as previously described. The introduction of a BCA that prevents infection by the pathogen also generates an outflow from HS (named RCOLH), which is calculated as follows: (11) in which c 1 is the relative rate of change from HS to H b and BSUM, STARTB, and COFR are as previously defined. The introduction of a BCA also generates an outflow from I. This occurs for those BCAs able to inhibit or reduce the sporulation on affected and infectious plant tissue (i.e., I) because of mycoparasitism and/or antibiosis, so that the infectious tissue reduces its ability to generate new infections. The outflow from I (named RCOLI) is calculated as follows: (12) in which c 2 is the relative rate of change from I to I b , BSUM and STARTB are as previously defined, and COFRI is the correction factor for infectious tissue and is calculated as follows: in which K and I are as previously described.
The introduction of a BCA also generates an outflow from R, even though this does not directly affect the epidemic. This outflow (termed RCOLR) is calculated as follows: (14) in which c 2 is the relative rate of change from R to R b , BSUM and STARTB are as previously defined, and COFRR is the correction factor for removed tissue and is calculated as follows: in which K and R are as previously described.
The model considers that as the plant tissue becomes colonized by the BCA (which is accounted for by Equations (9), (11), (12), and (14)), the plant tissue can revert to BCA-free tissue because of BCA mortality. The flows from H b , I b , and R b to HS, I, and R, respectively, are calculated through a rate of BCA mortality, BMOR (BMORH, BMORI, and BMORR, respectively), as follows: (16) in which f is the relative rate of mortality (i.e., the relative rate of change from H b , I b , or R b to HS, I, or R, respectively).

Driving Variables for the Pathogen
Driving variables are those functions that determine the relative rate of change of the system as influenced by external variables [39].
For pathogen infections that are influenced by temperature and relative humidity, the relative rate of infection (b) is calculated by Equation (17a): in which γ, ζ, and ν are the equation parameters accounting for the effect of temperature; ρ and ψ are the equation parameters accounting for the effect of humidity; Teq are temperature equivalents calculated as (T t − Tmin)/(Tmax − Tmin), in which T t is the average temperature (in • C) of day t; Agronomy 2020, 10, 222 8 of 21 Tmin and Tmax are minimal and maximal temperatures at which the pathogen can cause infection, respectively; and RH is the average relative humidity (%) of day t. For pathogen infections that are influenced by temperature and the duration of a moist period, the relative rate of infection (b) is calculated by Equation in which α, β, and θ are the equation parameters accounting for the effect of temperature; ϑ and ς are the equation parameters accounting for the effect of moisture; Teq is as previously described; MD is moisture duration (number of wet hours per day or number of hours with high RH, depending on the pathogen). Equation (17a) is a logistic equation, and Equation (17b) is a Gompertz equation, and both describe the S-shaped increase in infection as "moisture" (RH or MD, respectively) increases [40] up to an asymptote that is defined by temperature by means of a bell-shaped beta equation of Analytics [41]. In the beta equation, parameters γ and α define the top of the curve, ζ and β its symmetry, and ν and θ its size.
The relative rate of change from I to R (h) is calculated as follows: in which φ is the duration of the infectious period (in days) at the optimum temperature (Topt, • C), and T t , Tmin, and Tmax are as previously described.
In Equation (18), the temperature response curve is derived from Reed et al. [42] and Wadia and Butler [43].

Driving Variables for the BCA
The model considers four main biocontrol mechanisms: mycoparasitism, competition, antibiosis, and induced resistance. As in Jeger et al. [28], a single BCA can have one or more biocontrol mechanisms, and these may operate additively. The biocontrol mechanisms characterizing an individual BCA are included in the model as the BCA profile (PROF): (19) in which P, C, A, and IR are the relative contribution of mycoparasitism, competition, antibiosis, and induced resistance, respectively, to the overall BCA activity, considering that P + C + A + IR = 1. The relative rates of change from HS to H r (c 0 ), HS and H r to H b (c 1 ), I to I b , and R to R b (c 2 ) are calculated as follows: (22) in which GRO is the BCA growth rate under fluctuating temperature and moisture; and EFF 0 , EFF 1 , and EFF 2 are overall BCA efficacies in preventing the infection of the HS and H r tissue by induced resistance (EFF 0 ), antibiosis, and mycoparasitism (EFF 1 ), and in reducing the sporulation of the I tissue (EFF 2 ). In Equations (20), (21), and (22), GRO is calculated by using Equation (17b), in which α, β, and θ are replaced by χ, δ, and ε (the equation parameters accounting for the effect of temperature); ϑ and ς are replaced by ω and η (the equation parameters accounting for the effect of moisture); and Teq and MD are as previously defined. The relative rate of change from H b to HS, I b to I, and R b to R (f ) is calculated as follows: in which T t and RH t are as previously defined, and Tmin, Topt, Tmax, RHmin, RHopt, and RHmax are the minimal, optimal and maximal temperatures and RHs for BCA survival, respectively. The temperature and RH response curve in polynomial Equation (23), in which T and RH are the independent variables, is derived from Equation (18). The relative rate of change from H r to HS (e) is constant and depends on the duration of the induced resistance in the plant tissue, which depends on the combination of BCA and pathogen.

Model Output
The model output is represented by changes over time of the state variables in the system. An example of model output is shown in Figure 2A for three categories of host tissue: (i) healthy and susceptible (HS, green line); (ii) healthy and occupied by the BCA (H b , purple line); and (iii) occupied by the pathogen and infectious (I, red line). The simulation describes the changes in the proportion of the three categories of host tissue following the application of a preventative BCA on day 1 for a 40-d period during which the host tissue does not change because of plant growth and/or senescence. In Figure 2A, the proportion of HS tissue declines on day 1 because of the introduction of the BCA, which colonizes 60% of the tissue, and declines again at day 4 because of infection by the pathogen. Following infection, the tissue colonized by the BCA remains relatively constant until day 25; during this period, the BCA is effective in controlling the pathogen, which does not colonize additional tissue. After day 25, the tissue colonized by the BCA rapidly decreases, and the tissue occupied by the pathogen increases. Weather conditions ( Figure 2B) are important drivers for these dynamics, with a decrease in air temperature and wetness duration favoring the pathogen more than the BCA.
Agronomy 2019, 9, x FOR PEER REVIEW 9 of 20  An additional state variable, the area under the disease progress curve (AUDPC) [40], was calculated to evaluate the overall effects of BCA characteristics and usage and of environmental conditions on the disease development. The AUDPC was calculated at a daily rate (RAUDPC) as the sum of the total K units of plant tissue occupied by the pathogen (POCC), as follows:

Model Parametrization
The model was parametrized for the biocontrol of BBR in grapevine clusters during ripening, which is between "veraison" (growth stage GS83, [44]) and "berries ripe for harvest" (GS89). The simulation period was set at 40 days, with northern Italy as the reference environment [45]. K units are single berries, which can dynamically belong to one of the seven categories (state variables) of the model. Since the number of berries is already defined at the beginning of the simulation period (K = 1) and does not change because of plant growth from GS83 to GS89, relative rate of growth (RRG) was set to 0 for the entire simulation. Since we assumed that no berries become resistant to B. cinerea infection because of senescence during that growth stage, relative rate of senescence (RRS) was set to 0 for the entire simulation.
During ripening, BBR can develop under favorable weather conditions through three main pathways: (i) latent infections become visible as rotted berries, (ii) air-borne conidia germinate on and An additional state variable, the area under the disease progress curve (AUDPC) [40], was calculated to evaluate the overall effects of BCA characteristics and usage and of environmental conditions on the disease development. The AUDPC was calculated at a daily rate (RAUDPC) as the sum of the total K units of plant tissue occupied by the pathogen (POCC), as follows:

Model Parametrization
The model was parametrized for the biocontrol of BBR in grapevine clusters during ripening, which is between "veraison" (growth stage GS83, [44]) and "berries ripe for harvest" (GS89). The simulation period was set at 40 days, with northern Italy as the reference environment [45]. K units are single berries, which can dynamically belong to one of the seven categories (state variables) of the model. Since the number of berries is already defined at the beginning of the simulation period (K = 1) and does not change because of plant growth from GS83 to GS89, relative rate of growth (RRG) was set to 0 for the entire simulation. Since we assumed that no berries become resistant to B. cinerea infection because of senescence during that growth stage, relative rate of senescence (RRS) was set to 0 for the entire simulation.
During ripening, BBR can develop under favorable weather conditions through three main pathways: (i) latent infections become visible as rotted berries, (ii) air-borne conidia germinate on and infect berries, and (iii) aerial mycelium produced on rotted berries infects adjacent healthy berries (berry-to-berry infection) [46,47]. In the current study, we considered that latent infections and resulting berry-to-berry infections are more common than conidial infections [47][48][49][50][51]. We then assumed that the BBR epidemic starts with the onset of rotted berries that have been latently infected in early growth stages, which constitutes the initial inflow of the pathogen into the system (STARTP). This inflow is assumed to occur on the 4th day of the simulation (i.e., PIN = 4) and to continue for a period of PDUR = 1 day, at a rate RPIN = 0.2 (meaning that 20% of the berries are affected by latent infections).
During the simulation, new berries become affected (i.e., rotted) through the berry-to-berry pathway at the relative rate b, which is calculated by using Equation (17a), following Ciliberti et al. [52].
We assumed that as berries become affected, they begin producing conidia and enter in the I category. Afterwards, the affected berries continue producing conidia until harvest [46]; therefore, there is no outflow from I to R and h = 0.
Two BCAs with different multiple MOA (i.e., having different PROFs) were entered in the system (i.e., BCAs are applied to clusters) in different simulation runs. Specifically, the MOA profile of the first BCA is PROF = P (0.0) + C (0.8) + A (0.2) + IR (0.0). This profile can represent, for example, Aureobasidium pullulans, which is effective against B. cinerea by competing for nutrients at the infection site, which is its main MOA, and also by releasing hydrolytic enzymes that inhibit the pathogen [53,54]. The MOA profile of the second BCA is PROF = P (0.8) + C (0.2) + A (0.0) + IR (0.0). This MOA profile can represent, for example, Pythium oligandrum, which is mainly a mycoparasite but which also competes for nutrients with pathogens [55].
In the model, the overall BCA efficacies (EFF 0 , EFF 1 , and EFF 2 ) in preventing the infection are considered at their maximal (i.e., EFF 0 , EFF 1 , EFF 2 = 1), meaning that tissue colonized by the BCA totally prevents or reduces B. cinerea development.
Both BCAs are applied to clusters as a preventative treatment on the 1st day of simulation (BIN = 1) or as a curative treatment on the 7th day (BIN = 7). These applications constitute the initial inflow of the BCA into the system (STARTB), which has BDUR = 1 day (i.e., the day of BCA application) at a rate RBIN = 0.6 (meaning that the BCA covers 60% of the K units at the time of application).
Parameters of driving functions for calculating b, GRO, and f were derived from the literature and are indicated in Table 2. Rate b is calculated by using Equation (17a), as in Ciliberti et al. [52]. GRO is calculated by using Equation (17b) and by using different parameter values that describe the different responses to temperature and moisture of nine BCA strains (named S1 to S9, see Table 2). Table 2. Parameter estimates of the equations fitting the following relationships: the effects of temperature and relative humidity on b (the relative rate of Botrytis cinerea infection), the effects of temperature and moisture duration on GRO (the relative rate of growth of the BCA), and the effects of temperature and relative humidity on f (the relative rate of BCA mortality).  Finally, rate f is calculated using Equation (23), with different settings of parameter values referring to three temperature and humidity conditions under which the BCA survives (Table 2); these settings simulate different BCA manufacturing processes and/or formulations that result in different survival capabilities under stressful vineyard conditions [56,57].

Model Running
The model was used to study the effect of the following sources of variation on BBR development: (i) MOA of the BCA (2 levels: mainly competition and mainly mycoparasitism); (ii) BCA application time (2 levels: preventative and curative); (iii) BCA strain (9 levels: 3 ranges of temperatures combined with 3 moisture requirements for BCA growth); and (iv) BCA survival capability (3 levels: low, medium, and high). These sources of variation generate 108 combinations (2 MOAs × 2 application times × 9 strains × 3 survival capabilities). In addition, a situation with no BCA application was considered as the untreated control (NT). To study the effects of environmental conditions, each combination was run under nine scenarios that reflect three climate types with three scenarios per type: (i) warm and dry, (ii) mild and semi-arid, and (iii) cold and wet. Scenarios are represented by fluctuating conditions of temperature, relative humidity, and wetness duration (see Table 3). Therefore, 981 model runs were generated: (1 NT + 108 BCA combinations) × 9 climate scenarios. An example of the effect of the previously mentioned sources of variation on the disease dynamics in the three climate types is provided in Figure 3. Each graph shows the simulated proportion of the host tissue occupied by the pathogen (POCC) over the entire simulation period for each climate type (blue lines, cold and wet; green lines, mild and semi-arid; and yellow lines, warm and dry; Figure 3). Simulations of Figure 3 refer to a competitive BCA ( Figure 3A,B) or to a mycoparasitic BCA ( Figure 3C,D) applied as a preventative treatment ( Figure 3A,C) or as a curative treatment ( Figure 3B,D), with the temperature and moisture requirements of S8 and with low survival capability ( Table 2). In the cold and wet climate type, BCA application reduced the final (i.e., at day 40) value of POCC by 11% to 16%. In the mild and semi-arid climate type, BCA application reduced the final value of POCC by 53% to 68%, irrespective of application time or MOA. In the warm and dry climate type, BCA application reduced the final value of POCC by 40% to 45%.
The final values of AUDPC simulated for each of the 981 runs of the model were used to calculate the efficacy (E) of each BCA combination (T) in relation to the untreated control (NT), as follows: E = (NT -T)/NT. A factorial analysis of variance (ANOVA) was carried out for each climate type to determine whether the efficacy of each BCA combination was significantly affected by the main sources of variation (MOA, application time, strain, and survival capability) or their interactions. The three scenarios per climate type were used as replicates. The ANOVA was conducted by using the function anova of R software (v 3.6.0; R core team, [58]). Figure 4 summarizes the efficacy of the BCA in controlling BBR for the different simulation runs. Under the cold and wet climate ( Figure 4A), in which BBR developed rapidly and occupied all of the host tissue in < 20 days (Figure 3), BCA efficacy ranged from 0% to 99% and was significantly influenced by all of the main sources of variation (p < 0.001) and by the following interactions: MOA × application time, application time × strain, application × survival capability, and strain × survival capability (Table 4). MOA and application time (preventative or curative) accounted for 2.6% and 1.7% of total variance, respectively, and all of the interactions accounted for <1.7% of the total variance (Table 4). Those factors that are greatly affected by environmental conditions (the BCA strain and its survival capability) together accounted for 91% of the total variance ( Table 4). The average efficacy was higher for BCAs with medium or high survival capability ( Figure 5A) than for BCAs with low survival capability. The average efficacy in the cold and wet climate was higher for S2 and S8 than for S6 ( Figure 6). proportion of the host tissue occupied by the pathogen (POCC) over the entire simulation period for each climate type (blue lines, cold and wet; green lines, mild and semi-arid; and yellow lines, warm and dry; Figure 3). Simulations of Figure 3 refer to a competitive BCA ( Figure 3A,B) or to a mycoparasitic BCA ( Figure 3C,D) applied as a preventative treatment ( Figure 3A,C) or as a curative treatment ( Figure 3B and 3D), with the temperature and moisture requirements of S8 and with low survival capability (Table 2). In the cold and wet climate type, BCA application reduced the final (i.e., at day 40) value of POCC by 11% to 16%. In the mild and semi-arid climate type, BCA application reduced the final value of POCC by 53% to 68%, irrespective of application time or MOA. In the warm and dry climate type, BCA application reduced the final value of POCC by 40% to 45%. The final values of AUDPC simulated for each of the 981 runs of the model were used to calculate the efficacy (E) of each BCA combination (T) in relation to the untreated control (NT), as follows: E = (NT -T) / NT. A factorial analysis of variance (ANOVA) was carried out for each climate type to determine whether the efficacy of each BCA combination was significantly affected by the main sources of variation (MOA, application time, strain, and survival capability) or their interactions. The three scenarios per climate type were used as replicates. The ANOVA was conducted by using the function anova of R software (v 3.6.0; R core team, [58]). Figure 4 summarizes the efficacy of the BCA in controlling BBR for the different simulation runs. Under the cold and wet climate ( Figure 4A), in which BBR developed rapidly and occupied all of the Dashed lines indicate POCC dynamics when no BCA application is applied (NT), and solid lines indicate POCC dynamics when a BCA is applied. Blue, green, and yellow lines indicate the simulation in a cold and wet, mild and semi-arid, and warm and dry climate type, respectively. Each line corresponds to the POCC dynamics averaged across the three scenarios used as replicates for each climate type.
( Table 4). Those factors that are greatly affected by environmental conditions (the BCA strain and its survival capability) together accounted for 91% of the total variance ( Table 4). The average efficacy was higher for BCAs with medium or high survival capability ( Figure 5A) than for BCAs with low survival capability. The average efficacy in the cold and wet climate was higher for S2 and S8 than for S6 ( Figure 6).  figure); responses of 9 BCA strains to temperature (S1 to S9; X axis, see Table 2 for details); BCA survival capability (low, medium, and high, as indicated on the right side of each plot); and climate (A: cold and wet, B: mild and semiarid, and C: warm and dry). Each point represents the average, and the bars represent the standard errors of three scenarios per each climate type. Red and blue colors indicate that the BCA is applied as a preventative or a curative, respectively, i.e., before or after Botrytis cinerea infection.  figure); responses of 9 BCA strains to temperature (S1 to S9; X axis, see Table 2 for details); BCA survival capability (low, medium, and high, as indicated on the right side of each plot); and climate (A: cold and wet, B: mild and semi-arid, and C: warm and dry). Each point represents the average, and the bars represent the standard errors of three scenarios per each climate type. Red and blue colors indicate that the BCA is applied as a preventative or a curative, respectively, i.e., before or after Botrytis cinerea infection.  Figure 5. BCA efficacy against Botrytis bunch rot in ripening grapevine clusters averaged across nine BCA strains and as affected by BCA survival capability (low, medium, and high; see Table 2) and climate type: (A) cold and wet, (B) mild and semi-arid, and (C) warm and dry. The thick line in the boxes is the median, the lowest value in each box represents the 1st quartile (25th percentile), the highest value of each box represents the 3rd quartile (75th percentile), red points in the boxes are the means, and black points in the graph are outliers.
. Figure 6. BCA efficacy against Botrytis bunch rot in ripening grapevine clusters of nine BCA strains (see Table 2) in the cold and wet climate. The thick line in the boxes is the median, the lowest value in each box represents the 1st quartile (25th percentile), the highest value of each box represents the 3rd quartile (75th percentile), red points in the boxes are the means, and black points in the graph are outliers.
Under the mild and semi-arid climate ( Figure 4B), in which BBR developed gradually and did not completely occupy the host tissue until the end of the simulation period (Figure 3), BCA efficacy ranged from 1% to 71% and was significantly influenced by the MOA of the BCA (P = 0.017), which accounted for 0.8% of total variance (Table 4), and by the growth requirement of the BCA (as indicated by the strain; P = 0.002), which accounted for 0.4% of total variance (Table 4). BCA efficacy was significantly affected (P < 0.001) by the survival capability of the BCA, which accounted for 97.9% of the total variance (Table 4); BCA efficacy increased with the survival capability of the BCA ( Figure  5B).
Under the warm and dry climate ( Figure 4C), in which BBR developed slowly and occupied only 50% of the host tissue at the end of the simulation period (Figure 3), BCA efficacy ranged from 0% to 40% and was significantly influenced by the survival capability of the BCA (P < 0.001) and by the interaction between application time and survival capability (P = 0.007), which accounted for 97.3% and 2.2% of total variance, respectively ( Table 4). The average efficacy was higher for BCAs with a high survival capability than for BCAs with a low or medium survival capability ( Figure 5C). MOA, strain, and application time did not significantly affect BCA efficacy (Table 4).  Table 2) and climate type: (A) cold and wet, (B) mild and semi-arid, and (C) warm and dry. The thick line in the boxes is the median, the lowest value in each box represents the 1st quartile (25th percentile), the highest value of each box represents the 3rd quartile (75th percentile), red points in the boxes are the means, and black points in the graph are outliers.
Agronomy 2019, 9, x FOR PEER REVIEW 14 of 20 Figure 5. BCA efficacy against Botrytis bunch rot in ripening grapevine clusters averaged across nine BCA strains and as affected by BCA survival capability (low, medium, and high; see Table 2) and climate type: (A) cold and wet, (B) mild and semi-arid, and (C) warm and dry. The thick line in the boxes is the median, the lowest value in each box represents the 1st quartile (25th percentile), the highest value of each box represents the 3rd quartile (75th percentile), red points in the boxes are the means, and black points in the graph are outliers.
. Figure 6. BCA efficacy against Botrytis bunch rot in ripening grapevine clusters of nine BCA strains (see Table 2) in the cold and wet climate. The thick line in the boxes is the median, the lowest value in each box represents the 1st quartile (25th percentile), the highest value of each box represents the 3rd quartile (75th percentile), red points in the boxes are the means, and black points in the graph are outliers.
Under the mild and semi-arid climate ( Figure 4B), in which BBR developed gradually and did not completely occupy the host tissue until the end of the simulation period (Figure 3), BCA efficacy ranged from 1% to 71% and was significantly influenced by the MOA of the BCA (P = 0.017), which accounted for 0.8% of total variance (Table 4), and by the growth requirement of the BCA (as indicated by the strain; P = 0.002), which accounted for 0.4% of total variance (Table 4). BCA efficacy was significantly affected (P < 0.001) by the survival capability of the BCA, which accounted for 97.9% of the total variance (Table 4); BCA efficacy increased with the survival capability of the BCA ( Figure  5B).
Under the warm and dry climate ( Figure 4C), in which BBR developed slowly and occupied only 50% of the host tissue at the end of the simulation period (Figure 3), BCA efficacy ranged from 0% to 40% and was significantly influenced by the survival capability of the BCA (P < 0.001) and by the interaction between application time and survival capability (P = 0.007), which accounted for 97.3% and 2.2% of total variance, respectively ( Table 4). The average efficacy was higher for BCAs with a high survival capability than for BCAs with a low or medium survival capability ( Figure 5C). MOA, strain, and application time did not significantly affect BCA efficacy (Table 4).  Table 2) in the cold and wet climate. The thick line in the boxes is the median, the lowest value in each box represents the 1st quartile (25th percentile), the highest value of each box represents the 3rd quartile (75th percentile), red points in the boxes are the means, and black points in the graph are outliers.
Under the mild and semi-arid climate ( Figure 4B), in which BBR developed gradually and did not completely occupy the host tissue until the end of the simulation period (Figure 3), BCA efficacy ranged from 1% to 71% and was significantly influenced by the MOA of the BCA (P = 0.017), which accounted for 0.8% of total variance (Table 4), and by the growth requirement of the BCA (as indicated by the strain; p = 0.002), which accounted for 0.4% of total variance (Table 4). BCA efficacy was significantly affected (p < 0.001) by the survival capability of the BCA, which accounted for 97.9% of the total variance (Table 4); BCA efficacy increased with the survival capability of the BCA ( Figure 5B).
Under the warm and dry climate ( Figure 4C), in which BBR developed slowly and occupied only 50% of the host tissue at the end of the simulation period (Figure 3), BCA efficacy ranged from 0% to 40% and was significantly influenced by the survival capability of the BCA (p < 0.001) and by the interaction between application time and survival capability (P = 0.007), which accounted for 97.3% and 2.2% of total variance, respectively ( Table 4). The average efficacy was higher for BCAs with a high survival capability than for BCAs with a low or medium survival capability ( Figure 5C). MOA, strain, and application time did not significantly affect BCA efficacy (Table 4).

Discussion and Conclusions
The model that was developed by Jeger et al. [28] and that was improved by Xu et al. [29], Xu et al. [30], and Xu and Jeger [31], accounts for the biocontrol mechanisms involved and is able to predict the dynamics of pathogen and biocontrol agent (BCA) populations. In Xu and Jeger [31], the significant effects of varying BCA-temperature relationships and application times on BCA efficacy suggested the importance of considering environmental conditions under which the BCA and target pathogen interact. In the present research, the model of Jeger et al. [28] was enlarged to include crop growth and senescence and the environmental effects on the pathogen and on BCA growth and survival. Like the model of Jeger et al. [28], the enlarged model has a generic structure and can be applied to any pathosystem involving fungal pathogens of aerial plant organs, as well as different pathogen-BCA interactions involving different BCA mechanisms of action.
We parametrized the enlarged model for B. cinerea causing Botrytis bunch rot (BBR) on grapevines. The model parametrization was derived from the epidemiological studies performed by Ciliberti and colleagues [52,59,60]. Those epidemiological relationships were incorporated in a mechanistic model for B. cinerea-grapevine developed by González-Domínguez et al. [47], but the latter model did not include a BCA component. The use of BCAs for BBR control has been extensively studied [61][62][63][64][65][66], with emphasis on biocontrol mechanisms and field efficacy; less research has been conducted to understand how environmental conditions affect BCA fitness and efficacy [11]. In the current study, the model parametrization for BCAs used different parameter values represented by nine BCA strains, which differed in their growth and survival in response to temperature and moisture conditions.
The model was run under three climate types to study the combined effects of the following factors: (i) mechanism of action of the BCA, (ii) timing of BCA application with respect to the pathogen (preventative vs. curative), (iii) temperature and moisture requirements for BCA growth, and (iv) BCA survival capability. All of these factors affected, although to different degrees, biocontrol efficacy. Environmental conditions were the most important factors, accounting for > 90% of the variance in simulated biocontrol efficacy; other factors, even though significant under some climate types, accounted for only a minor percentage of the variance. This finding may help explain why the application of BCAs often results in inconsistent control of the target pathogen in the field [66]. In other words, our results suggest that the inconsistent BCA efficacy in repeated experiments [8,16,67] and in the practical biocontrol of diseases [63,[68][69][70] can be caused, at least to some extent, by differences in environmental conditions between experiments or by fluctuations in environmental conditions in the same experiment [1,12,13]. This finding also stresses the importance of considering the environmental response of the BCA during its selection, BCA survival capability during both selection and formulation, and weather conditions and forecasts at the time of BCA application in the field.
Concerning the environmental response of the BCA during its selection, BCAs that are able to grow under a wide range of environmental conditions (i.e., strains S2 and S8 in this study) and that share the temperature and moisture requirements of the target pathogen may be more effective than BCAs with a more limited ability to grow under a range of environmental conditions. BCA responses to temperature and moisture can be evaluated by means of environmentally controlled experiments [71], and the effects of temperature and moisture on the pathogen-BCA relationship can be evaluated by using environmental niches [11]. It is essential that the effects of environments be included when screening BCAs for market development [72,73].
Concerning the BCA survival capability during both selection and formulations, our model simulations indicate that BCAs may be more effective in controlling the target pathogen for long periods and under a range of weather conditions if they have a high rather than a low survival capability. This result confirms previous findings [74][75][76] and also the importance of protective effects provided by additives or adjuvants used in the formulation of the commercial product [56,77,78]. This result also confirms that survival capability should be a key property used to screen microorganisms for biocontrol [73].
Finally, weather conditions and forecasts at the time of BCA application in the field should be considered so as to maximize the probability that the BCA will grow and control the pathogen. Although the current model could be useful in this respect, its utility should be verified with field experiments [79]. On the other hand, developers of BCAs could use the current model to predict the efficacy of candidate organisms under different scenarios of weather conditions and application timings.