Information retrieved from selected papers was organized based on systems analysis (Leffelaar and Ferrari, 1989) and was used for model development. The model structure is shown in Figure 1. The diagram was drawn following the system representation of Forrester (1997) as used in STELLA^® (abbreviation of Systems Thinking, Experimental Learning Laboratory with Animation), a visual programming language for modeling system dynamics. The diagram combines state variables (rectangles), flows (solid arrows), rates (valves), parameters, and coefficients (circles).

The core of the model’s structure is based on the epidemiological model developed by Zadoks (1971), in which the crop is thought to consist of a large but finite number of infection sites that have equal dimensions and equal probability of becoming infected. A site is defined as a fraction of the host tissues where an infection may occur and where a lesion may develop (Zadoks, 1971; Savary et al., 2015). A site occurs under one of the following mutually exclusive conditions: (i) heathy (H); (ii) latent (L), i.e., without symptoms of stem rust; (iii) infectious (I), i.e., with visible lesions producing spores; and (iv) removed (R), i.e., lesions are older and non-sporulating (Figure 1 and Table 1). The structure used here represents the host in a very simple manner and does not incorporate host growth or senescence; it also does not account for lesion expansion from infected sites. Sites becoming infected at an infection rate (RI), i.e., they change from healthy sites into latent sites. Latent sites become infectious at the end of a latency period (LP) at a transfer rate (Rtrans), and infectious sites become removed at the end of an infectious period (IP) with a rate of removal (Rrem) (Figure 1 and Table 1).

At the beginning of model calculations, H = 1 (i.e., the whole crop is healthy), and the model represents the flow from one state to the following one as a proportion of the H, i.e., on a 0 to 1 scale. The model works with a daily time step, with the following exceptions: the onset of primary infection (see Section “Model Structure”) and the modifiers of RI (see Section “Modifiers of RcOPT”) are calculated hourly and are reported on the daily scale.

The effects of temperature and wetness duration on epidemics are incorporated into the model as mathematical equations. The effects of temperature on germination, penetration, the length of the latency period, and the length of the infectious period are described by the following equations: f(T)_GER, f(T)_PEN, f(T)_LP, and f(T)_IP (see Sections “Modifiers of RcOPT,” “Latency Period,” and “Infectious Period”). The effects of wetness duration on germination and penetration are described by the equations f(WD)_GER and f(WD)_PEN, respectively (see Section “Modifiers of RcOPT”).

Equations were developed and parameterized based on the papers that were collected through the systematic literature review. Useful data were retrieved from tables, figures, and texts of the selected papers. The shape of the data and the Akaike information criterion (Burnham and Anderson, 2002) were used to choose the best-fitting mathematical functions. Equation parameters were estimated using the function nls of the “stats” package of R software (Team, R Core. R: A Language and Environment for Statistical Computing. 2019, which is available at⁴). The parameterized equations were evaluated for goodness-of-fit based on the adjusted R², the concordance correlation coefficient (CCC), the root mean square error (RMSE), and the coefficient of residual mass (CRM) (Nash and Sutcliffe, 1970; Lin, 1989). The adjusted R² was estimated by conducting a linear regression between the observed values and the model predicted values; the linear regression was conducted with the lm function of the R “stats” package (Wickham, 2019). The CCC is a measure of model accuracy (Madden, 2006), and is calculated as the product of the Pearson correlation coefficient and the Cb coefficient; the CCC indicates the difference between the best fitting line and the perfect agreement line (Lin, 1989). The CCC was obtained using the CCC function of the R “DescTools” package (Signorell, 2020). The RMSE, which represents the average distance of real data from the fitted line (Nash and Sutcliffe, 1970), was obtained using the rmse function of the R “modelr” package (Wickham, 2019). The CRM is a measure of the tendency of the equation to overestimate or underestimate the observed values (a negative CRM indicates a tendency of the model toward overestimation) (Nash and Sutcliffe, 1970).

The date of the establishment of the epidemic (Start) depends on Onset, which represents the occurrence of weather conditions favorable for the first infection of the season. Onset is calculated hourly and reported at a daily time step as either Onset = 0, i.e., conditions are not favorable for infection, or Onset = 1, i.e., conditions are favorable for primary infection. Information retrieved from Singh et al. (2006, 2002), Roelfs (1985); Burrage (1970), Rowell and Romig (1966), and Sharp et al. (1958) were used to set Onset = 1 as follows: (i) rain event (P) ≥ 1 mm h^–1; (ii) the wetness period (WD) following rain is ≥3 h; and (iii) mean temperature during WD (Tw) is between 15 and 32°C. At Onset, the proportion of sites entering L depends on the parameter INOCP, which is an estimate of RI for the first seasonal infection.

Once the first seasonal infection has established, secondary infections are responsible for the development of epidemics through the transfer of sites from healthy to latent through a rate of infection (RI; i.e., the proportion of newly infected sites per unit of time). RI is modeled as a function of the number of infectious sites (I), a basic infection rate corrected for the removal (Rc; Vanderplank, 1963), and a correction factor for diseased sites (COFR) as follows:

COFR is calculated as follows:

where S is the sum of disease sites (S = L + I + R).

In equation (1), Rc represents the proportion of daughter lesions generated per mother lesion. In the model, Rc depends on the optimum corrected basic infection rate (RcOPT), which is the basic infection rate under optimum environmental conditions on a susceptible wheat variety, and on modifiers (RcMod). RcMod is the daily reported value of modifiers for the effect of temperature (RcT), wetness duration (RcWD), and light regime (RcLR) calculated hourly. Rc is calculated as follows:

RcOPT was estimated following Sun and Zeng (1994) from disease progress curves as follows:

where p is the latency period, which was set as 9 days under favorable conditions for epidemics (Mortensen and Green, 1978); i is the infectious period, which was set as 32 days under favorable conditions for epidemics (Mortensen and Green, 1978); and r is the apparent infection rate (Vanderplank, 1963, 1975), calculated as follows:

where x₁ and x₂ are disease fractions on two successive dates (t₁ and t₂) at the early stage of the epidemic under conditions conducive to the disease. Equation (5) was used to calculate r from published disease progress curves in susceptible and unprotected crops (Romig and Calpouzos, 1970; Nazareno and Roelfs, 1981; McGrath and Pennypacker, 1991). The first two non-zero severity values (expressed on a 0 to 1 scale) were retrieved from these curves and used for the calculation, resulting in an RcOPT value of 4.5.

Modifiers of RcOPT are used in the model to account for major environmental conditions influencing the infection process (Loomis and Adams, 1983). As mentioned before, modifier contributions are calculated by the model on an hourly basis and are subsequently reported as RcMod on a daily basis; this enabled the model to accurately account for the effects of environmental conditions on infection. RcMod is calculated as follows:

Development and parameterization of modifiers of RcOPT were based on the available literature on urediospore germination and germ tube penetration via stomata.

The modifier RcT accounts for the effect of air temperature (T, in°C) on both germination (f(T)_GER) and penetration (f(T)_PEN), as follows:

The effect of T on germination is calculated by a Bete equation (Analytis, 1977) in the following form:

where Teq is the equivalent of temperature, calculated as Teq = (Tw – Tmin)/(Tmax – Tmin); Tw is the mean temperature recorded during the wet period; Tmin is the minimum temperature for uredospore germination, i.e., 4°C (Burrage, 1970); and Tmax is the maximum temperature for uredospore germination, i.e., 35°C (Katsuya and Green, 1967). Germination does not occur when Tw < Tmin or Tw > Tmax. Equation (8) was developed and parameterized using the data of Burrage (1970) and Katsuya and Green (1967); estimates and standard errors of equation parameters were α = 4.375 ± 1.016, β = 1.391 ± 0.419, and γ = 0.396 ± 0.201, with adjusted R² = 0.916, CCC = 0.962, MRSE = 0.1, and CRM = 0.004.

The effect of T on penetration is described by the Duthie equation (Duthie, 1997) in the following form:

with ζ = ((ε + 1)/ε) × ε^{(1/(ε + 1))}; equation (9) was developed and parameterized using the data of Burrage (1970) and Sharp et al. (1958). Estimates and standard errors of equation parameters were ε = 2.699 ± 1.83, δ = 27.02 ± 1.37, and η = 0.683 ± 0.155, with adjusted R² = 0.824, CCC = 0.914, MRSE = 0.139, and CRM = 0.028.

The modifier RcWD accounts for the effect of wetness duration (WD, in h; i.e., cumulative number of h with leaf wetness) on both germination (f ‘(WD)_GER) and penetration (f ’(WD)_PEN), as follows:

The effect of WD on germination is calculated as the first derivative (f ’(WD)_GER) of the following equation:

which was developed and parameterized by fitting the data from Burrage (1969, 1970), and Eversmeyer and Kramer (1989). Estimates and standard errors of equation parameters were ι = 35.224 ± 3.663 and κ = 2.27 ± 0.618, with adjusted R² = 0.971, CCC = 0.982, MRSE = 0.073, and CRM = –0.066.

The effect of WD on penetration is calculated as the first derivative (f ’(WD)_PEN) of the following equation:

which was developed and parameterized by fitting the data from Burrage (1970). Because germinated uredospores form germ tubes on stomata about 3 h after the beginning of a wetness period (Roelfs, 1985; Singh et al., 2002), there is no penetration and f(WD)pen = 0 when WD < 3 h. Estimates and standard errors of equation parameter were λ = 1.207 ± 0.124 and μ = 0.068 ± 0.017, with adjusted R² = 0.951, CCC = 0.978, MRSE = 0.084, and CMR = 0.032.

The RcLR modifier accounts for the light/dark regime during infection. Germination of uredospores is slightly higher in the dark than in daylight (Givan and Bromfield, 1963), but penetration of germ tubes is substantially inhibited in the dark. Therefore, RcLR = 1 in daylight hours and =0.452 in dark hours (Sharp et al., 1958).

Latent sites gradually flow from L to I, with a rate of transfer (Rtrans) governed by a temperature-dependent equation accounting for the latency period (LP). The rate of transfer is calculated daily as the first-order derivative of the following Gompertz equation:

where DD is the degree-days accumulated from the day in which the infection occurred, with 5°C as the minimum temperature (Singh et al., 2002); therefore, when T ≤ 5°C, T = 0; when T > 5°C, T = (T – 5). Equation (13) was developed and parameterized by fitting the data of Katsuya and Green (1967); estimates and standard errors of equation parameters were ν = 2632 ± 255.8 and ρ = 0.07 ± 0.008, with adjusted R² = 0.823, CCC = 0.904, MRSE = 0.084, and CMR = 0.032.

Sites I remain infectious for a period (IP), the length of which is determined by temperature. The flow from I to R is therefore governed by a temperature-dependent rate of removal (Rrem). The length of the infectious period is determined by the following logistic equation:

where DD is the degree-days accumulated from the day of transfer of sites from L to I, with 5°C as the minimum temperature. As a consequence, when T ≤ 5°C, T = 0; when T > 5°C, T = (T – 5). The f(T)_IP ranges from 0 to 1, and sites are removed from I when f(T)_IP = 1. Equation (14) was developed and parameterized by fitting the data of Mortensen and Green (1978); estimates and standard errors of equation parameters were ς = 23.97 ± 4.15 and σ = 0.012 ± 0.0006, with adjusted R² = 0.965, CCC = 0.982, MRSE = 0.063, and CMR = –0.005.

The model calculates disease severity (DS) during the progress of the epidemic as the sum of the proportion of sites that are carrying visible lesions, i.e., infectious and removed sites, as follows:

To evaluate the ability of the model to predict the first seasonal infection of stem rust, model outputs were validated against real data obtained from nine experimental wheat fields in Italy, which were established for a range of both durum (Triticum durum L., cv. Fuego, Farah, Telemaco, Antalis, Brancaleone, Salgado, Acropolis, RGT Aventadur, Ramirez, Gibraltar, Puro, Platone, Tirex, Teodorico, Don Matteo, and Pigreco) and bread (T. aestivum L., cv. Drusilla, Lucilla, Tocayo, Caronte, Nabucco, Bologna, and Nogal) varieties of wheat between 2016 and 2021 (Table 2). All the tested varieties were susceptible to stem rust. All plots (10 m²) in each experimental field were inspected at 5- to 10-day intervals for the onset of stem rust symptoms in any variety. Weather data were recorded by meteorological stations (PESSL iMetos 3.3) located <2 km from the experimental fields, except for USS16 and USS17, where the weather station was about 15 km away.

The model was operated by using weather data in each year and location starting from 9 days before the beginning of the disease assessment period (Table 2). Disease assessments (cases) were categorized as either 0 (disease symptoms were not observed on wheat plants, O–) or 1 (disease symptoms were observed on wheat plants, O+). Similarly, cases were categorized as either 0 or 1 based on whether the onset of the disease was predicted (P+) or not (P–) by the model. Because first seasonal infections are seldom severe (Prescott et al., 1986), predicted lesions were considered to be visible for a field assessor when 50% of the uredia generated by the first seasonal infection were predicted to erupt.

According to Madden (2006), a Bayesian analysis was used to evaluate the correspondence between model predictions and stem rust symptoms observed in the field. We tested the hypothesis that new stem rust symptoms were observed in the field (O+) on those days when the model predicted lesion appearance (P+), and that no disease symptoms were observed (O–) on those days when the model did not predict lesion appearance (P–). A contingency table (2 × 2) was prepared containing (i) the true positive proportion (TPP or sensitivity), (ii) the true negative proportion (TNP or specificity), (iii) the false positive proportion (FPP), and (iv) the false negative proportion (FNP).

Prior and posterior probabilities of predicting the first seasonal stem rust occurrence based on model output were calculated (Madden, 2006). The prior probabilities of disease to occur P(O+) or not to occur P(O–) were compared with the following posterior probabilities: (i) the probability of disease occurrence when predicted by the model, P(P + O+); (ii) the probability of disease occurrence when not predicted by the model (i.e., missed real infections), P(P–O+); (iii) the probability of no disease occurrence when not predicted by the model, P(P–O–); and (iv) the probability of no disease occurrence when predicted by the model (i.e., unjustified alarms), P(P + O–). Positive and negative likelihood ratios were calculated as LR(+) = TPP/FPP and LR(–) = FNP/TNP, respectively. The diagnostic ability of the model was evaluated by means of the Youden index, calculated as J = TPP – FPP. Overall model accuracy was estimated as the ratio between correct and total predictions.

To validate the ability of the model to predict disease development throughout the season, six disease progress curves were retrieved from the literature. Details of experimental sites used for disease progress validation are summarized in Table 3. Initial infections resulted from naturally occurring airborne inoculum in two epidemics and from inoculum originating from artificially and heavily infected “spreader wheat rows” located near the experimental plots in four epidemics (Table 3).

The model was operated beginning at 9 days before the day of the last assessment in which no disease was observed. For epidemics caused by artificially inoculated spreader rows, rain was not considered in the calculation of the variable Onset. Values of INOCP used to initialize model runs are reported in Table 3; because of the use of heavily infected spreader rows, the value of INOCP used for epidemics caused by artificial inoculum was about 10 times greater than the value used for epidemics caused by natural inoculum (Table 3). Because there was no information in the literature for calculating INOCP, values were estimated empirically. For each epidemic, the value of INOCP was determined based on a comparison of predicted and observed final disease severity. The value of INOP that resulted in the closest agreement between the final predicted disease severity and the final observed disease severity was used for model runs.

Predicted disease severities (i.e., sum of infectious and removed sites) and observed disease severities were compared. For the evaluation of model performance, RMSE, CRM, and CCC were calculated (Nash and Sutcliffe, 1970; Lin, 1989).