In the southeastern U.S., H. zea completes several generations that cycle through non-crop weeds and then crops that include corn (Zea mays L.), upland cotton (Gossypium hirsutum L.), soybean (Glycine max (L.) Merr.), and several other minor crops (12–15). Larvae are the economically important life stage because they consume vegetative and reproductive plant structures that result in crop yield loss (e.g., cotton bolls, soybean pods, and corn kernels). As the season progresses, populations grow when high-quality crop hosts become abundant in the landscape, with some of the largest seasonal populations originating from corn (15). The adult flight originating from corn often coincides with blooming cotton and soybean in early August.

To provide an early warning system for H. zea flights, university extension researchers maintain trap networks to monitor H. zea adults across the region (16). Resulting trap data is used as an early warning system to incentivize intensive in-field H. zea egg and larval monitoring. Importantly, growers often infer H. zea activity for their locations based on a small number of traps distributed across large geographic regions. For example, the North Carolina trap network has included 22 unique black light trap locations in 2024 that span approximately 36,000 km² of land area from central to eastern North Carolina. Traps in this network are concentrated in major agricultural regions with the goal of providing growers with timely pest activity information (data available at: https://www.ces.ncsu.edu/trap-data/). These networks provide useful H. zea activity information each year but may not account for population variation across spatially and temporally heterogenous agricultural landscapes (17, 18).

This study documented variation in H. zea adult activity in an experimental trap network distributed across 21 locations in a major row crop production landscape from 2020-2023. The trap network spans approximately 700 km² of agricultural land in five North Carolina counties (i.e. Northampton, Halifax, Nash, Edgecombe, and Wilson). The fine-scale trap network monitored adult male H. zea using commercially available pheromone lures (PHEROCON^® CEW, Trécé Inc., Adair, OK) and fabricated metal Texas Hartstack traps (11). We collaborated with row crop growers to select 21 trap locations that were distributed across the study extent and located in open areas adjacent to agricultural fields. Each year, traps were installed at the same location late in July and monitored weekly for 6-11 weeks during the peak flight period ( Figure 1 ). Captured moths from individual traps were returned to the laboratory and counted each week.

Our dataset includes n = 705 observations, L = 21 locations, T = 4 years (2020-2023). We have 83 site-year combinations (site 3 did not participate in 2020) and 6 - 11 weeks of observations per year. In 2020 there were 6 weeks of trapping, 2021 had 7 weeks, 2022 had 10 weeks, and 2023 had 11 weeks. In 2023, Trap 13 was missing an observation for week 2, so for this observation we imputed the average value of all other sites for this same week. Let S = { s 1 , … , s L } be the set of all sampling locations. Our response variable of interest is the cumulative pest count for each week. In particular, if Z t s x is the pest count (in hundreds of pests) for location s in year t at week x , then our response variable, Y t s x , is the cumulative pest count for location s in year t at week x , i.e.,

Note that by definition, Y t , s , x ≤ Y t , s , x + 1 .

The logistic growth curve represents biological population growth by including both exponential increase and a density-limiting carrying capacity. It is defined as

for all x , where x represents the week of the observation. The model is defined by three parameters. β is supremum of the function, and here represents the maximum number of cumulative pests. A larger value of β means a greater number of pests. The growth parameter, k , represents how quickly the pest population is increasing, and is inversely correlated with growth rate. If k is large, then the growth is slower, while a small k indicates fast growth. Finally, γ is the midpoint and represents where the function attains half of its supremum, i.e., f ( γ ) = β / 2 .   A larger γ means that the midpoint is later in the growing season, and vice-versa for a smaller γ .

As formulated in Equation 1, the logistic growth curve is the same at each trapping location. Because spatial dependencies among trapping locations are expected, β and γ depend on location. Mathematically, the model becomes

where

The asymptote parameter, β , has been replaced with a spatially-varying parameter, β ( s ) . Similarly, the midpoint γ , is now modeled spatially with γ ( s ) . Spatial terms are modeled as a Gaussian Process with exponential correlation, i.e.,

and

γ ( s 1 ) , … , γ ( s L ) ~ N o r m a l ( 0 ,   σ γ 2   Σ γ ) ,

where Σ β and Σ γ represent exponential correlation structures, a special case of Matern correlation. If Σ is an exponential correlation matrix, then

where ρ is the spatial range parameter, and | | ⋅ | | is the Euclidean distance between trapping locations. In this correlation function, points are more correlated if they are closer together, and the correlation decreases exponentially fast with distance. This model was then fit using all locations and all years (before 2023) via Bayesian inference. Please see the Supplementary Materials for full prior specification.

There are several noteworthy points about this model. First, the model has spatial variation, but is the same for each year, i.e., there is no temporal component. There were not enough years of data to accurately learn a parameter differentiating years, and the model thus overfit on year when such a parameter was included. Next, the model has spatial variation in the supremum and midpoint, but not in the growth parameter. Additional models were tested with spatial variation in k , but there was very little variation among sites. Thus, the growth parameter was fixed across space. Additionally, both β ( s ) and γ ( s ) have exponential correlation structures, but are allowed to have different spatial range parameters, ρ β and ρ γ . Finally, covariates were not included in the model (e.g., landscape composition), because we found that models including candidate covariates tended to overfit to the data and did not provide good out-of-sample prediction. While we might reasonably assume that covariates such as landscape composition weather, etc. do affect the pest population, the current model implicitly accounts for such factors in the spatial processes. Moreover, the focus of this paper is selecting the optimal trapping locations, not necessarily the reason that these sites were chosen. We consider this latter question an important avenue for future work.

After empirically establishing the model form, we proceeded to determine how to select the optimal trapping location sets. We seek a small number of trapping locations, K , with K < L , such that if we fit a model on these K locations, then the model can most accurately predict at the L − K unobserved sites. For now, we consider K as known and fixed but we will generalize this in the experimental results.

We seek to minimize the mean log error (MLE) between the real and predicted cumulative counts at the unobserved locations in 2023. Let S i = { s i 1 , … , s i K } be a set of K locations where S i ⊆ S , and let Y ^ τ s x ( S i ) be the predicted value at site s ∉ S i in year τ = 2023 for week x ∈ { 1 , … , 11 } , using the locations in S i to fit the logistic growth model. Thus, Y ^ τ s x ( S i ) are out-of-sample predictions. Then our loss function is the MLE on the sites we did not train on for 2023, defined as

where log() is the natural logarithm. Optimization thus seeks

In this case, the optimum is the set of sampling locations, S ˜ , such that, if the model is trained on this set of locations from 2020-2023, and then used to predict trap data from other locations in 2023, MLE will be minimized for predictions of these out-of-sample points. 2023 was chosen as the testing data subset to simulate forecasting of future data using past data. MLE was chosen for our loss function due to the following observation. The observed cumulative counts vary in orders of magnitude, i.e., the first week may have tens of counts while the cumulative sum at the end of the season can be in the thousands. By taking the logarithm of the difference between the observed and fitted values, it ensures that all errors are on the same scale, such that each component of the sum in (3) contributes approximately equally to the loss function. Please see the Supplementary Materials for further discussion on the loss-function as well as potential other options.

The loss function in Equation 3 is complex enough to complicate optimization in practice. Indeed, it is not easy to minimize this loss function directly, primarily because it is a combinatorial optimization problem with L choose K many solutions. Because fitting the logistic growth curve model takes a nontrivial amount of time, a brute force approach is not feasible. Thus, we employ Bayesian Optimization (BO) to approximate the optimal trapping locations S ˜ . There are three main steps to a BO schema: the objective function, statistical model, and acquisition function, after Yanchenko (19). Please see Figure 2 for an overview of the modeling workflow.

We now apply the BO method from Section 2.4 to our H. zea trap data. For a fixed K , we follow the method outline in Section 2.4 to obtain our approximate optimum x * . We then use the trap locations S * corresponding to x * to fit a logistic growth model based on the data from 2023 and compute the MLE on the hold-out sites for 2023, i.e. L ( S * ) . We do this for K = 5 , 6 , … , 14 . For comparison, we also randomly select K trap locations and compute the out-of-sample MSE for 2023, which serves as a baseline method. Both the BO algorithm and random sampling are repeated for 50 Monte Carlo (MC) samples, and the average MLE is reported. We use the standard Kriging estimate to predict cumulative counts at new locations, and all Bayesian models are fit using Stan in R (20).

The main results show that the proposed BO approach greatly outperforms randomly sampling locations with a lower MLE for all K ( Figure 3 ). Additionally, the MLE monotonically decreases as the number of training sites increases for the BO algorithm. We plot the proportion of MC samples for which a trap location was chosen in the optimal set for K = 5 , 10 ( Figures 4A, B ), we also plot the proportion of samples a trap was chosen in the optimal set averaged over all K ( Figure 4 ). The central traps are chosen for the optimal set more often when K = 5 , whereas the spread of optimal sites for K = 10 is biased toward the south. When the results are averaged over all K , the southern and central locations have a higher proportion of MC samples for which they are in the optimal set.

While Figures 4 , 5 show the proportion of MC samples for which a site was chosen in the optimal set, Figure 6 considers the binary inclusion/exclusion of sites in the optimal set. For each value of K , the proportion of samples for which each site was chosen in the optimal set was calculated. Then the K sites with the largest proportion are chosen as the optimal sites. This is represented in Figure 6 as a black tile if the trap has one of the K largest proportions, and a white tile otherwise. In case of ties, both sites have a black tile. Traps 1, 5 and 14 are included in the optimal site for every value of K . Conversely, traps 13, 15, 16, 18, and 21 are never selected. The plot also shows a substantial amount of nesting in the optimal sites, i.e., trap i is selected in the optimal set for K , then it is also selected in the optimal set for K + 1 .

As an extension of this work, we could explicitly include money, time, distance, etc. into our loss function if we wanted these to factor into our optimal choice (e.g., 9). For example,

where C ( S i ) is the cost to sample at locations S i and λ is a tuning parameter which controls how much we weigh fit (small λ ) versus cost (large λ ) in the loss function.

The general implication of these empirical findings is that optimization enriches the information provided per unit trap in a multi-site network. First, the magnitude of the MLE associated with an optimized set is appreciably smaller than that for random sampling. For example, a MLE of 1.2 means that the predicted count was within 230 units of the true cumulative count value, on average. In the Supplementary Materials , we also compare the results of the BO algorithm with a brute-force computation for K = 5 and find that the BO results are close to the global optimum with similar trap combinations detected by each method ( Supplementary Table S1 ). This establishes the utility of the method for efficiently selecting trapping locations.

The BO algorithm was run 50 times and the average MLE was recorded. Ideally, the algorithm would be run once and the global optimum would be found, such that running multiple MC iterations would be unnecessary. The necessity of running the BO algorithm multiple times could be because the algorithm is finding a local optimum and/or because the objective function surface is relatively “flat.” This would mean that several different sets of K trap locations yield comparable MLEs. While these considerations and the computational expense of running the BO algorithm multiple times are relevant, running the algorithm multiple times also has advantages. The proportion of MC samples that a trap is included in the optimal set yields a continuous measure of importance as opposed to a binary result. This allows for some uncertainty quantification in the optimal sites and a richer understanding of which sites to choose for maximum precision in pest density estimation, and perhaps which sites to choose for research to understand why they are inconsistent or otherwise different from the optimal set.

In general, the optimal trapping locations appear in the southern and central locations. This could be interpreted one of two ways. It could mean that knowing the number of H. zea at these locations is very informative for knowing the number of moths at other locations. On the other hand, it could mean that these locations are very difficult to predict so we should trap at these sites to ensure an accurate value. For example, trap 1 could be included in every optimal set because knowing its value helps predict at the other sites, or because predicting the number of H. zea at trap 1 is extremely difficult. In the Supplementary Materials , we looked for systematic reasons why some trapping locations were chosen more than others. We found a moderate relationship showing that traps with a high probability of being included in the optimal set tended to have larger amounts of corn and soy planted within a 1 km radius. These are initial findings, however, and require further study.

Together, these steps toward optimization are important, as the costs of monitoring H. zea are immediately clear. Establishing which traps in the network provide the most return on investment would be a significant step toward sustainable monitoring approaches that continue to provide useful information to agricultural stakeholders. In the future, this approach would also apply to invasive species monitoring, which is a common approach when outbreaks occur. The persistent threat of pests with similar ecology as H. zea (i.e., Helicoverpa armigera) is a major focus of regulatory agencies tasked with responding to incipient threats of invasive species in the United States (e.g., USDA-APHIS). Ecological niche modeling suggests that H. armigera could establish in the southern U.S. and would pose a significant threat to agriculture across broad geographies (21). Effective monitoring strategies would be a first step toward effective management of a widespread invasive outbreak. The BO approach is one tool which could be adapted to model effective monitoring networks that are informed by the current distribution and activity of a closely related species.

These results suggest that understanding these contrasting motivating factors will help to assess which locations are informative and/or what about these locations make them informative. Future studies focused on assessing location-wise characteristics from the perspective of H. zea could generate new context for the optimized dataset. For example, overwintering temperatures limit H. zea survival during the pupal stage (22) which contributes to regional trends in population abundance over time (16). Moreover, soil conditions which vary at much smaller spatial scales are also linked to overwintering survival which affects H. zea population size in the spring (23, 24). During the growing season, the abundance of highly suitable host plants for H. zea development across multiple generations may also explain variation among traps across the relatively small geographic region where we established our trap network (12–14, 18). Although these studies independently assessed specific aspects of H. zea ecology, they are unable to account for multiple variables influencing these insects in agricultural systems. Absent an optimized network, it is difficult to learn what factors at key trap sites are responsible for those sites’ providing information that can positively affect grower management decisions. Next steps connecting optimized networks to relevant ecological and environmental drivers would further refine our understanding about H. zea population dynamics in the context of realistic agroecosystems. In short, the improvement of observational efficiency results in the improvement of any study or application that relies on precision estimates of the observed process. Optimizing the provision of evidence used to describe biological systems, to make system management decisions, or to discover novel or incipient events in those system is a valuable step to take early in a study or management program.

This study builds on a deep body of literature focused on improving sampling optimization for insect pests in agriculture. Entomologists apply numerous sampling techniques for the purpose of estimating arthropod population densities, sometimes adapting from methods used in engineering (25) or analysis of data from clinical trials (26, 27). A commonly applied sampling technique in entomological pest management is sequential sampling, first developed in the 1940’s (28), for the purpose of increasing the efficiency of estimation for values important to industrial processes. The primary motivation for the development of sequential sampling techniques was efficiency: the potential to accelerate estimation or increase precision per sample, while minimizing wasted sampling effort. In pest management, this translates to iterative sampling and interim analysis based on observations. Samplers use pre-established stopping conditions to determine when sufficient pest density or crop damage information has been collected to make a management decision. This level of tolerated pest activity is informed by the economic thresholds, which is a pest density or other damage metric above which not managing the pest is expected to result in net economic loss. For example, the threshold at which the value of the expected yield loss to H. zea larval herbivory exceeds the cost of management. Effectively, this effort becomes a statistical hypothesis, and sequential sampling proceeds until the hypothesis is rejected or fails to be rejected under pre-established tolerances for Type I and Type II error rates. While some important developments have extended the rationale of increasing sampling efficiency for decision support (e.g. 25, 29–31, 32), efficiency gains through innovative statistical approaches continue to represent a major opportunity for improving the environmental compatibility and economy of agricultural pest management. The BO algorithm developed in this study enhances sampling in a way parallel to existing practices of sequential sampling or interim analysis, but adds a substantial amount of value to data in the process by showing the aforementioned nesting of high-value sites within larger optimal sets, by providing a measure of uncertainty around the importance of individual sites, and by including correlation structure that allows nonrandom associations between sites to be exploited for optimization. Similar approaches could be developed for other pest species and geographic regions but will require initial investment to collect activity data for model development and validation.