Evaluation of predicted Medfly (Ceratitis capitata) quarantine length in the United States utilizing degree-day and agent-based models

Sites and temperature data

Hourly air temperature data for 11 sites was downloaded from NOAA’s publicly available Integrated Surface Database (ISD) dataset^28,29.

The airport sites shown in Figure 1 were chosen for their biological relevance and availability of high quality hourly data over a long time frame.

Sites are referred to here by the last three letters of the callsign shown in Table 1. For 8 sites (SFO, FAT, LAX, RIV, SAN, JAX, TPA, and MIA), temperature data starting on 1950-01-01 was used. The 3 other sites contained large (>14 days) gaps or other problems in the early years of their data, so data starting on 1970-01-01 for IAH and 1973-01-01 for BUR and MCO was used. For all sites, temperature data from the start date through 2017-05-15 was used to generate PQLs for dates ranging from the start date to 2016-01-01.

Data was fetched and parsed using the Fetching and parsing ISH.ipynb program. Records for the same station callsign were merged, since identification, format, and precise location of stations has changed over time. The data was then cleaned using the Cleaning temperatures.ipynb by removing outliers, identifying large gaps (> 3 hours), resampling to every hour on the hour using linear interpolation, and filling the large gaps using day-over-day linear interpolation (interpolating using values for the same hour of day from previous and following days). The resulting temperature datasets are available.

Degree-day calculation

Degree-days were computed by the single-sine method²⁶, using a base development temperature of 12.39°C (53.3°F) and 345.56 degree-days Celsius (DDc; 622 DDf) per generation following the standard set by California Department of Food and Agriculture regulation 3406(b)^23,30. Since hourly temperature data are available, we also calculated degree-days by simple summation for comparison²⁴. For each date, the number of days required for 3 generations of degree-day based life cycles was computed. These calculations are implemented in Temperature functions.ipynb.

Agent-based simulations: MED-FOES

MED-FOES^25,27 is an agent-based simulation explicitly modeling the eradication of a population of Medflies under inundative sterile male releases (sterile insect technique or SIT) and other interventions, such as increased trapping and foliar sprays. A MED-FOES simulation models a single non-spatial population, starting from a given population size and age distribution, tracking the number of individuals through time until the last fly (Agent) dies and the population is eliminated. In addition to hourly temperatures, simulation parameters include: the initial population, additional mortality induced by control efforts, the effectiveness of SIT, and a large number of biological parameters for which ranges are known from the literature including temperature-dependent development and mortality. The simulations were performed using the same hourly time series of temperature values used for degree-day calculations.

Due to the fact that many of the parameters are only known to within a range, 2500 individual MED-FOES simulations were run for each start date at each site, evenly sampling different regions of parameter-space via the Latin Hypercube Sampling³¹ procedure. This set of simulations, encompassing a range of possible elimination outcomes, is referred to as a ‘run’. The number of days from the start date required for 95% of the simulations in a run to be eliminated is taken as a conservative prediction of needed quarantine length and referred to as ABS PQL²⁵.

Varying the start date for different simulations was achieved by simply starting at different points in the input temperature file; for this study a run was started every 7 days over the range of dates available for each site. Each set of runs for a single site over a range of starting dates is referred to as a ‘runset’. All runsets were conducted with the same input parameters aside from temperature. Initial population numbers were chosen as a “standard outbreak” based on seven real outbreaks modeled previously²⁵. The 7 day interval ABS PQL values were upsampled to daily values using linear interpolation to allow day-of-year aggregations across years and comparisons with daily degree day based PQLs.

MED-FOES version 0.6.2 was run under Open Grid Scheduler/Grid Engine 2011.11 on a CentOS 6.6 HPC cluster. The MED-FOES code, configuration files, helper scripts, and raw results are available. Overall, we created 11 runsets (one for each site). Each runset contained runs starting every 7 days over the input temperature data range for that site, and each run contained 2500 individual simulations sampling different regions of biologically plausible parameter space. This sums to a total of approximately 86×10⁶ simulations.

Statistical analysis

The main results reported here are ‘normals’ in a meteorological sense of the term, but without the typical running mean smoothing which would complicate interpretation. For a variable of interest (eg. temperature or PQL), all values for the same calendar day irrespective of year (eg. 20-July) are aggregated, and summary statistics such as mean, minimum, maximum, and standard deviation are computed for each aggregation. Temperature functions.ipynb contains the code used to perform normal calculations, and the code generating figures as well as all statistical analysis is Summary Figures.ipynb (Jupyter Notebook³², module and version information documented in the file).

The results reported here are the normals of PQL, computed using the full temperature time series as opposed to computing PQL from the normal of the temperature time series. While the latter is fairly common practice, it is not mathematically proper since, as with means, the normal of a function of X is not generally equal to the function applied to the normal of X. Additionally, by computing the normals of the predicted quarantine durations, we can investigate properties of the distribution of values as shown in Figure 3 and Figure 4 and the “supernorm” supplementary figure S1, supplementary figure S2, and supplementary figure S3.

Results

Figure 2 shows the mean of the normal PQL based on 3 generation degree day accumulation and MED-FOES 95% elimination along with the minimum and maximum of the normals for temperatures. Figure 3 and Figure 4 show the standard deviations (σ) of the normals for the degree day and ABS based PQL.

Figure 2. Summary of normal predicted quarantine length for each site and start date (last fly detection).

Year range of input temperature data used is inclusive. All panels have identical limits except SFO quarantine.

Figure 3. Variation in predicted quarantine length for each site and start date based on 3 generations of single-sine degree day accumulation.

Figure 4. Variation in predicted quarantine length for each site and start date based on 95% of MED-FOES agent-based simulations showing elimination.

There is significant variation in PQL across both time and location. The temporal variation in PQL is dominated by a yearly cycle, characterized by the normal values shown in Figure 2. Table 2 shows the percentage of variance in quarantine length predictions captured by the mean of the normal yearly cycle (R²) for each site. At all but one site, greater than 75% of the variance in both degree day and ABS based PQLs is accounted for by the mean normal, and the majority exceed 90%. SFO is an exception to this common trend, with the mean normal accounting for only 9.1% of the variation in degree day based PQL and 28.0% of the ABS based PQL. This is also reflected in supplementary figure S2 and supplementary figure S3.

Seasonal dependence

Seasonal variation, evidenced by the general shape of the curves shown in Figure 2, is doubtless familiar to anyone engaged in Medfly pest management. Outbreaks starting in the late summer, autumn, or early winter will extend through relatively cold periods, when thermal dependent development will be slow and therefore extend the duration of quarantine required for 3 generations of degree days to accumulate (referred to as DD PQL hereafter). Similarly, outbreaks starting in the spring or early summer often lead to short quarantines due to the relatively high temperatures.

This familiar pattern is also seen in the ABS PQLs despite it being quite different in nature from simple degree day accumulation. However, the ABS predictions show a smaller seasonal swing. The ABS generally produces a smaller overall range of PQLs, with longer quarantines than DD PQL for spring and early summer outbreaks, and shorter quarantines for late summer through early winter in almost all cases.

A particular feature of interest, shown most dramatically at FAT in Figure 2, is that ABS PQL often flattens out or even dips for quarantines starting in the late autumn or early winter. This can be due to relatively rare and brief cold-snaps, normally lasting only a few hours, which increase mortality. Since DD PQL does not account for mortality, it misses the effect of cold-snaps entirely. This effect is most clearly seen at more northern and inland sites where cold-snaps are more likely: particularly FAT and RIV, but also BUR, LAX, JAX, and IAH.

Geographic dependence

PQL generally shows a positive correlation with latitude, and sites are ordered by latitude in the figures and tables here. As seen in Figure 2, higher latitude sites tend to have longer PQLs as well as larger seasonal swings for both degree day and ABS based predictions.

Figure 5 shows the relationship between PQL and latitude. An ordinary least squares fit to the median PQL at each site shows a significant slope for both DD PQL (F =14.08, p=0.005) and ABS PQL (F =10.55, p=0.010), but the degree day based predictions are more sensitive to latitude than the ABS (coefficients of 17.39 and 4.78 respectively). Additionally, the ABS predictions are more stable for SFO, and to a lesser extent FAT, where the degree day model for Medfly produced PQLs that appear either unrealistically long (SFO) or are subject to rapid and extreme seasonal variation in the mid year (FAT).

Figure 5. Predicted quarantine length dependence on latitude.

For each site, the mean, median, and inter-quartile range are shown, similar to a boxplot. An ordinary least-squares linear fit to the median values is shown by the green lines. The left panel is for single sine degree day predictions, and MED-FOES ABS based predictions are in the right panel.

In addition to the variation associated with latitude, large differences in PQLs computed for the same start date can exist between even relatively nearby sites. For example, the differences in both degree day and ABS PQLs for the three sites in the Los Angeles region (LAX, BUR, RIV) (shown in the supplementary figure S4) display a strong seasonal component with a spike in July and/or August. The difference in DD PQL between LAX and BUR is normally about a month (overall median=35 days; overall 25% & 75% quantiles are 28 & 45 days), but the median difference of the normal exceeds 75 days in August with some PQL differences up to 142 days. Differences in ABS PQLs are more seasonally stable, with the LAX minus BUR difference not exceeding 42 days for any start date in the 43 years analyzed here.

Variance and uncertainty

Figure 3 and Figure 4 report the standard deviation (σ) of the normal for DD PQL and the MED-FOES ABS PQL respectively. These indicate the year to year variability of the PQL for outbreaks starting at a given time of the year and can be used to gauge the uncertainty of predictions based on past PQLs relative to the actual quarantine length which will be required. Similar information is represented by the inter-quartile ranges shown in Figure 5 and supplementary figure S2 and supplementary figure S3. The distributions of PQL values for a site and day-of-year (aggregating across years) are generally not highly skewed, making σ a relatively easy to interpret measure of uncertainty.

Excluding SFO, the mean normal is a good predictor of DD PQL with σ values below 20 days except for the late summer and early autumn, where variance increases due to quarantines extending through the cold season. FAT and, to a lesser extent, RIV show this increase more dramatically, presumably due to their more arid/inland climates where both daily and seasonal temperature ranges are larger (also see Figure 2). The standard deviation generally decreases with decreasing latitude, together with reduced means. The standard deviation in DD PQL for SFO shows an inversion of the seasonal trend other sites exhibit. This is due to the colder temperatures leading to extremely long DD PQLs, frequently extending across two winter seasons.

The standard deviations of the ABS PQL normals shown in Figure 4 are generally about ½ as large as for DD PQL. This indicates that the ABS PQL not only shows less dramatic seasonal swings, but is also produces more consistent predictions across years. Values again generally decrease with latitude, but less consistently than DD PQL σ of normals. Also, unlike with the DD PQL, the results for SFO appear consistent with other sites.

A notable feature is that BUR, LAX, and SAN all show an increase in the year to year variation in ABS PQLs starting in July and extending through November, while that increase for all other sites starts in July or August but extends to January or February. Additionally, results for FAT show a sharp increase in uncertainty starting in September, fitting with the more arid/inland climate. RIV shows a significant but more gradual increase.

Historical quarantines

Thirty-four Medfly quarantines in CA dating from 1975 to early 2017 were analyzed (supplementary table S1). The start of all but two of these quarantines was in the latter half of the year (July through December), when DD PQLs are typically relatively long, with 68% (23/34) occurring in September through October, when DD PQLs are longest. August, the month where uncertainty in DD PQL often spikes (see Figure 3), accounts for 30% (7/34) of historic quarantines.

For each historic quarantine start date, the DD PQL and ABS PQL for the closest of the 11 sites analyzed above (see Figure 1 and Table 1) to the actual outbreak location was determined (see supplementary table S1). For this set of hypothetical quarantines, the ABS produced significantly shorter quarantines (mean=169.7 days, σ=21.8 days) than simple 3 generation degree day accumulation (mean=234.2 days, σ=79.2 days) (df =33, t=6.01, p<10⁻⁵). Additionally, the variance in the difference between quarantine lengths using a specific date and the mean of the normal PQL for that day of year was smaller for the ABS (σ=8.2 days) than with degree day (σ=25.9 days) (df =33, F =9.92, p<10⁻⁸).