Abstract
The COVID-19 pandemic exposed challenges of balancing public health and economic goals of infection control in essential industries like food production. To enhance decision-making during future outbreaks, we developed a customizable agent-based model (FInd CoV Control) that predicts and counterfactually compares COVID-19 transmission in a food production operation under various interventions. The model tracks the number of infections as well as economic outcomes (e.g., number of unavailable workers, direct expenses, production losses). The results revealed strong trade-offs between public health and economic impacts of interventions. Temperature screening and virus testing protect public health but have substantial economic downsides. Vaccination, while inexpensive, is too slow as a reactive strategy. Intensive physical distancing and biosafety interventions prove cost-effective. The variability and bimodality in predicted impacts of interventions caution against relying on single-operation real-world data for decision-making. These findings underscore the need for a proactive infrastructure capable of rapidly developing integrated infection-economic mechanistic models to guide infection control, policy-making, and socially acceptable decisions.
Teaser COVID-19 model helps navigate trade-offs between public health and economic impacts of infection control interventions in essential industries.
Introduction
The United States (US) food industry, known for its labor-intensive nature (1), was significantly affected by the Coronavirus disease 2019 (COVID-19) pandemic, alongside other essential industry sectors (2). During the early phases of the pandemic, food facilities/operations abroad and in the US were forced to close or reduce production due to labor shortages (2, 3). US livestock processing, including poultry, pig, and cattle slaughter, was reduced by up to 45%, resulting in job losses, financial impacts, retail shortages and loss of animals (4, 5). During the COVID-19 pandemic in 2020, the combined value of production for beef, pork, broilers, turkeys, eggs and milk was reduced down by $12.8 billion in the US, to 9% below the pre-pandemic forecast for 2020, based on price and production quantity projections (6). Dairy supply chain disruptions caused increased milk dumping in the US, where 2.5% of all federally regulated milk was dumped compared to 0.2-0.5% dumping recorded in the normal course of production (7). By September 2021, nearly 100,000 workers in US meatpacking facilities, food processing facilities, and farms were reported positive for COVID-19, most of which were from meatpacking facilities (65%) and other food processing facilities (21%) (8); importantly, these statistics are likely underestimates (9). These outbreaks drove infection rates in rural communities, as individuals infected at work transmitted their infection to others outside of work (10). By the end of 2020, COVID-19 cases attributable to meatpacking facilities were reported to be the source of an estimated 334,000 infections in the US, with associated mortality and morbidity costs totaling more than US$11.2 billion (11). These impacts affected the functioning of the national food supply chain.
Several mitigation strategies have been considered, encouraged, or enforced to control the spread of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) within the US food industry (12, 13). These control strategies include vaccination, practicing physical distancing, use of face coverings, screening for infection, practicing personal hygiene (e.g., hand washing), cleaning and disinfection of the working environments, ventilation improvements, and minimizing community spread. Produce (i.e., fruit and vegetable) farms and food processing facilities, while having certain features in common, vary greatly in size, physical characteristics, and organizational structure (14). Thereafter, for brevity, we use ‘operation’ when referring to any individual produce farm operation or food processing facility. Food operation’s varying locations, policies, and workforce demographics have resulted in significant differences in worker histories with respect to vaccination, boosting, and past infection (15, 16). The diversity of operations and mitigations have led to strong interest from industry stakeholders in modeling tools tailored to the particular characteristics of their individual operations (17) to aid them in making predictions, such as regarding the expected outbreak dynamics and impacts of possible interventions, and decisions, such as what level of investments to make in biosafety measures or when to start or stop an intervention (personal communication with the Industry Advisory Council for the study). Importantly, the evaluation of COVID-19 mitigation strategies should be based not only on public health metrics but also on economic metrics that account for the production losses in the operation due to worker shortages or strict infection control strategies, as well as considerations of negative societal impacts of food supply chain disruptions and possible food shortages.
Several mathematical models (18–30) have been developed to evaluate and compare COVID-19 mitigation strategies and assess their effectiveness across different levels of compliance. These models are primarily designed for national scale assessments (21, 25–27) but also include more localized communities, encompassing cities (24, 28, 30), closed societies with shared environments (19, 29), and even smaller communities in universities (22), companies (20, 23), and office spaces (18). Some of the evaluated interventions include physical distancing, mask use, vaccination, asymptomatic/symptomatic testing, contact tracing, quarantine, restrictions on travel, isolation, and school closures (18, 19, 26–28). While most of the modeling studies concentrate on health outcomes in the general population, only a handful have considered the health of individuals within workplace settings. These settings include a generalized company building (20), an oil and gas facility (23), a meatpacking plant (29), and a university building (22). These studies have accounted for the complex process of disease transmission between individuals by using agent-based models (ABMs), which can simulate employees’ decisions based on their social and physical profiles. A few studies have also utilized ABMs to simulate the economic impacts of COVID-19 (24, 25, 30, 31). Nevertheless, there remains a need for models that explore COVID-19 spread as well as the health and economic impacts of mitigation strategies in individual food operations to help prevent similar impacts in future infection outbreaks.
Here, we provide Food Industry CoVid-19 Control Tool (FInd CoV Control), a customizable tool based on an ABM developed to simulate COVID-19 transmission dynamics in the food industry work environment. We integrated the ABM with an economic model to predict the direct and (certain) indirect costs of interventions. Using COVID-19 in the food industry workforce as a model system, our objective was to develop a tool that helps policymakers and individual operations navigate tradeoffs between public health and economic impacts of infection control interventions in an essential industry.
Results
Model setup
FInd CoV Control consists of three modules: Employee population, Work environment, and Disease transmission (Figure 1A; definitions in Text S1, further details in Texts S2-S4 and Tables S1-S13). The Employee population includes all employees (agents) in a modeled operation, each of which is characterized by a set of attributes (Table 1). The Work environment module defines the characteristics of the work environment in terms of a produce farm or processing facility setting (thereafter referred to as ‘farm’ and ‘facility’ for brevity), shift schedule (Figure 1B.i), and agent hierarchy and contact network (Figures 1B.ii, 1B.iii, and S1). The Disease transmission module tracks COVID-19 infection spread in a population representing employees of a operation, using an elaborated variant of a “Susceptible-Exposed-Infectious- Recovered-Susceptible” (SEIRS) model (Figure 1C). FInd CoV Control is customized to the population and work environment of a particular operation based on the user-set parameters (Table 2), which are also used to calculate the number of agents with various immunity trajectories and histories (Table 3). We used FInd CoV Control to make predictions about COVID-19 transmission dynamics following the arrival at work of an index case infected outside the workplace, both under a no-intervention “baseline” and under various interventions. The evaluated interventions included: a temperature screening intervention, three virus testing interventions, a total of five primary vaccination and/or boosting-promoting interventions, and three direct basic reproduction number (R0)-reduction (physical distancing and/or biosafety) interventions (Figure 1A). The baseline and interventions were simulated in a way that allows counterfactual comparisons. Interventions were evaluated using two groups of metrics: (i) Public health: the number of employees with symptomatic and asymptomatic infection (and total infected); and the initial effective reproductive number (Reff.) and (ii) Economic: the number of employees unavailable to work; the fraction of shifts with employee shortage; and total direct expenses, production losses, and total costs associated with an intervention (expressed in US$) (Figure 1C). While economic effects are often interrelated and ripple over multiple dimensions, in FInd Cov Control, the economic analysis is limited to the costs directly borne by operations (Text S5) and is meant to serve as a reference, together with the infection model, for employers’ decision-making.
Model validation
FInd CoV Control was validated with publicly available data on outbreaks from early in the pandemic when few, if any, interventions would have been implemented. Specifically, for produce farm operations, FInd CoV Control was validated using two outbreaks, one on a farm with shared (i.e., employer-provided dormitory style) housing and one with a mix of shared and individual housing. For processing facilities, FInd CoV Control was validated using three outbreaks in facilities with individual housing, one in each of the dairy, pork, and produce processing facilities. These outbreaks and validation results are described in Text S6. The results of the validation analysis indicated a reasonable fit between the reported data and model predictions.
Main results
We begin by presenting a representative set of results over a 90-day simulation, for a facility with 103 employees, shared housing, and otherwise default parameters (see Tables 2 and 4), in Figures 2, 3, and 4. Results for the farm model with similar parameters are qualitatively comparable, and outcome estimates for a larger facility with 1,003 employees (Figures S2, S3 and S4) are generally similar with the main differences being later-peaking outbreaks (due to greater incidence growth required to saturate the larger population), and an associated modest increase in the effectiveness of some interventions. There is also a modest reduction in noisiness, which leads to a reduction in the probability of experiencing labor shortages in individual runs, given that the average outcome is not to experience a shortage in either case. Differences associated with shared versus individual housing are covered in the “Scenario analysis” section (and Text S7).
Results pertaining to the number of symptomatic infections, our primary public-health outcome, are presented in Figure 2. There is little qualitative difference between the curves depicting the mean incidence (Figure 2A) and mean prevalence (Figure 2B) of symptomatic infection over time, although there is a difference in scale and a slight difference in location (with peak incidence occurring slightly prior to peak prevalence) and noisiness (with more visual noise in the incidence vs prevalence curves). This result is expected, so to avoid redundancy, we focus on (cumulative) incidence in the remaining panels.
Bimodality and variability in symptomatic infections
For all interventions (including the no-intervention baseline), over 40% of all runs result in no symptomatic infections at all (Figure 2C), and as a result, the number of symptomatic infections for a given intervention is strongly bimodal at baseline and for all interventions except for moderate viral testing (p = 0.3, i.e., 30% of scheduled workers tested each shift, amounting to testing every worker 1.5 times per week), high-intensity viral testing (p = 1, i.e., every worker scheduled for testing each shift), and “high-effectiveness (-80% R0) Physical distancing/Biosafety” (Figure 2D). This reflects that, in the absence of repeated reintroduction from the broader community, a major source of variance is the possibility of early stochastic die- off, even at an initial Reff. well above 1 (e.g., Reff. = 2.52 at baseline for this scenario), which effectively partitions the vast majority of outcomes (for interventions for which Reff. remains appreciably above 1) into two modal regions (groups) with respect to the outbreak size: (i) small or non-existent outbreaks and (ii) large outbreaks. At baseline, over the 90 days since the introduction of an index case, group (i) has 0-12 (median 0, mean 0.60) total infections, not including the index case, and 0-4 symptomatic infections (median 0, mean 0.16), with no infected individuals left by the end of simulation, while group (ii) has 51-97 total infections (median 82, mean 81.9), not including the index case, and 17-49 (median 34, mean 33.1) symptomatic infections. There are also 5 runs (out of 1000) that fell between these two groups, with 21-42 total infections (median 39, mean 35.8) and 5-14 symptomatic infections (median 11, mean 10.6).
Bimodality and variability in counterfactual effects of interventions
We can refine and expand observations detailed in the previous section by taking advantage of the steps we have taken to make counterfactual comparisons as precise as possible (see “Interventions” section); specifically, the i-th run of the model with any intervention corresponds in a meaningful way to the i-th run of the model at baseline, with the difference between the two being attributable solely the intervention. Consequently, it is meaningful to examine the pairwise differences, with respect to a particular outcome, between corresponding runs with and without a given intervention. The distributions of these pairwise differences, for runs that have at least one symptomatic infection at baseline, are presented in Figure 2E. There, we can see that there is not only a great deal of variance in outcomes within an intervention (or within the baseline), but also a substantial variance in the counterfactual effects of an intervention. This is meaningful, considering that the individual ABM runs reflect the real-world variation in epidemiologic outcomes, and the model provides a view into the counterfactual comparisons within individual runs that cannot be observed in the real world. To illustrate this phenomenon, a single intervention can be examined in detail (Figure 2E). The “moderate- effectiveness (-40% R0) Physical distancing/Biosafety” intervention is modestly beneficial on average (mean reduction in number of symptomatic infections = 10.4) and in its typical performance (median reduction = 6). Nevertheless, it can be extremely effective in individual runs (maximum reduction of 45, close to the maximum across all interventions of 48). On the other hand, it can prove entirely ineffective or even counterfactually counterproductive in other individual runs (61 runs with no reduction in number of symptomatic infections, and 68 runs with an increase of 1-20. This is due to the timing and chance effects as explained in Text S8, where results for additional interventions (temperature screening and viral testing) are also illustrated.
We can also note that most of these distributions of pairwise differences are themselves bimodal, reflecting two different ways that an intervention can counterfactually affect a run that produces large outbreak in the no-intervention scenario. In these counterfactual comparisons, on the one hand, an intervention may prevent a large outbreak altogether, producing a data point in the high- effectiveness modal region in Figure 2E (and contributing to the difference in the number of large outbreaks between the intervention and no-intervention scenarios in Figure 2D).
Alternatively, it may produce a smaller difference in the outbreak size (or none at all), producing a data point in the low-effectiveness modal region (and still contributing a large outbreak to both the intervention and no-intervention scenarios in Figure 2D). Depending on the intervention, one of these modal regions may be extremely small, or they may both be substantial.
We can further refine these observations by considering the change in number of symptomatic infections as a fraction of the baseline number of symptomatic infections (Figure 2F; Figure S2F for a large facility). In particular, for those interventions with a significant fraction of runs in the high-effectiveness modal region, the primary way that they shift runs from having large outbreaks in the no-intervention scenario, to not having large outbreaks in the presence of the intervention is, by causing them to have no symptomatic infections at all; the apparently symmetrical lower modes seen in Figure 2E are primarily produced by variation in the number of infections (in a large outbreak) at baseline, not by variation in the number of infections (in a small or (effectively) non-existent outbreak) under the intervention.
“All good things in moderation” may backfire in viral testing
Based on the different patterns of effects in Figure 2F, the most effective reductions in symptomatic infections are seen for moderate- and high-intensity viral testing. Figure 4 (Figure S4 for a large facility) shows that viral testing at a fairly high rate is also generally more costly. In the case of high-intensity viral testing, this cost is overwhelmingly due to direct intervention expenses (primarily the cost of test kits), and this result is robust across a variety of scenarios (not shown). This is also true in most runs for moderate-intensity viral testing, but in some runs, moderate-intensity viral testing can result in significant costs due to both direct intervention expenses and production losses. The latter reflects the ability of testing at a “moderate” rate (p=0.3/working day) to generate a “worst of both worlds” scenario. This scenario generates large numbers of employees who are isolated at the same time, resulting in large numbers of worker- shifts missed due to isolation, yet, infected employees are not identified and isolated fast enough to prevent a large outbreak from occurring. A more frequent version of this (rarely producing > 15% absences on a production shift, which is considered to cause shortages and, thus, production losses) can be seen for low-intensity viral testing (p = 0.05/work day) where mean and median increases in unavailability are both greatest (Figures 3E and 3F; Figures S3E and S3F for a large facility). This observed pattern reinforces and extends the result from prior research that existent but inadequate larger-scale (city-level) non-pharmaceutical interventions can result in what the authors describe as a “dual blow of increased deaths and unemployment,” which fall disproportionately on low-income workers (30).
Health benefits of physical distancing/biosafety interventions at low cost
Finally, the next-most effective intervention in preventing symptomatic infections, after the moderate- and high-intensity viral testing interventions is the “high-effectiveness (-80% R0) Physical distancing/Biosafety” intervention (Figure 2F), for illustration purposes represented by a combination of masking, face shield use, and ventilation improvements. We found this to be much more effective than the “moderate-effectiveness (-40% R0) Physical distancing/Biosafety” intervention (represented by masking and face shield use, without ventilation improvements), but only modestly more costly (and substantially less costly than the more effective viral testing interventions) (Figure 4).
Scenario analysis
For our scenario analysis, we first defined several elements whose effects and interactions with intervention effects we wished to examine (Figure 5). These scenario elements were “setting” (“farm” vs. “facility”), “housing” (“individual” vs. “shared”), “vaccinated” (“high” (based on US national levels in early 2022 (32,33)) vs. “none”), and “recovered” (“high” US national levels (32, 34, 35) vs. “none”) (details in section Text S7 and Table S12). We then conducted a full factorial analysis for all 16 combinations of these four factors (and for all 13 intervention scenarios) in the default facility size of 103 workers over 90-day-long simulation runs. Many combinations of these scenarios are intended to represent limiting cases, rather than realistic scenarios, e.g., a scenario with both “vaccinated” = “high” and “recovered” = “none” represents a limiting case of relatively high vaccination and no history of infection. Results of this analysis, evaluated using regression trees for each of the three primary outcomes (symptomatic infections, worker-shifts unavailable, and total cost) and the two separate contributors to total cost (production losses and intervention expenses), indicated that, for outcomes other than production losses, the effects of “setting” were relatively limited, and mostly pertained to which other effects were strong enough to be included in the pruned partition trees. Because transmission in the two settings is defined by the user-settable value of R0 (Table 2), which was kept the same between the two settings, the observed differences between settings can be attributed to differences in the work environment (Figures 1B and S1). In general, the evaluated outcomes were slightly higher for the “facility” setting than for the “farm” setting. Because of this, and to omit explanations of the different equations used to set default production-per-week in the two different settings, we chose to describe results from the facility model here. Because total cost is simply a sum of intervention expenses and production losses, and because various factors affect each of those components differently, we will focus our discussion on each cost individually (results presented in Figure 5).
Health and economic outcomes are driven by the interaction between the worker infection history and intervention intensity
For almost all evaluated outcomes, the two biggest factors driving outcomes are intensity of virus testing intervention and “recovered” (i.e., whether the employee population has a significant history of natural infection) (Figure 5); the only exception is intervention expenses (Figure 5C), for which intensity of virus testing was the strongest factor, but “recovered” did not produce a sufficient impact to appear in the regression tree. While “recovered” being “high” (rather than “none”) had a desirable impact (i.e., produced lower symptomatic infections, unavailability, and intervention expenses) on all outcomes for which it was relevant, the effects of virus testing were more variable. Symptomatic infections are minimized by viral testing at a rate high enough to reliably control an outbreak before it can get large (testing every worker every shift or roughly every 3 shifts (p = 1 or 0.3/ work day, respectively)) (Figure 5A). Unavailability (Figure 5B), on the other hand, is lowest when either testing is non-existent (and so no workers are isolated as a result of testing, but only as a result of hospitalization) or testing is extremely intensive (and so the outbreak(s) is/are rapidly contained; p = 1/work day); however, even such intense testing may be insufficient to achieve a reasonable level of control in the face of a population with insufficient natural and hybrid immunity (“recovered” = “none”) and constant reintroduction (housing = “individual”). Conversely, unavailability is highest when the testing rate is intermediate (p = 0.05 or, even more so, p = 0.3/work day), resulting in enough asymptomatic and mildly symptomatic cases being detected to increase unavailability, but not enough to rapidly contain the outbreak(s). Similarly, production losses (which are driven by unavailability of ≥ 15% on a production shift) are highest when p = 0.3/work day, to the point that recursive partitioning with default parameters results in further division only of the node with p = 0.3; all scenario-intervention combinations with either a lower (p = 0, p = 0.05) or a higher (p = 1) rate of virus testing are combined in a single node (Figure 5D). Intervention expenses, on the other hand, move in an opposite pattern to symptomatic infections, and are highest when the rate of viral testing per work day is highest (p = 1), and lowest when it is low (p = 0.05) or non- existent (p = 0) (Figure 5C).
Strongest outcome drivers are viral testing intensity and history of natural infection
While “recovered” is consistently more important than “vaccinated” (i.e., splits defined by it occur closer to the root of each tree, where either occurs at all), and both share a uniformly desirable effect (where they show any effect at all), the interaction of the two is more complex: For production losses, there is only a split defined by “vaccinated” in a branch in which “recovered” is set to “none,” but for symptomatic infections, the reverse is true – the only split defined by “vaccinated” occurs in a branch in which “recovered” is set to “high.” In more conceptual terms, this amounts to saying that, at the population level, immunity resulting from recovery from natural infection plays a stronger role in determining a wide range of simulation outcomes than immunity resulting from vaccination - perhaps an unsurprising result late in a pandemic; whether the interaction of these two is sub-additive or super-additive depends on which outcome one is considering. In particular, with respect to symptomatic infections, there may be a synergistic effect of vaccination and recovery from natural infection, likely reflecting the strong protective effect of modeled hybrid immunity. For production losses, on the other hand, vaccination is less influential in the presence of moderate-to-high levels of natural recovery within the past year. This likely reflects the threshold effect in our model of production losses (of 15%) – in the presence of sufficient protection from natural recovery, the probability of suffering production losses at all may be low enough, even in the absence of vaccination, to reduce the importance of vaccination in predicting or determining that outcome.
The only panel from which the “recovered” factor is absent, or even not one of the two strongest factors, is Figure 5C Intervention Expenses; the “vaccinated” factor is absent from this panel as well (although the cost of giving workers time off for vaccination is accounted for). This is unsurprising, given that intervention expenses are driven far more - at least, for the relatively simple interventions that we consider in this analysis - by what interventions one decides to implement than by transmission dynamics; as a result, no scenario parameters appear in it. The only split that does, other than the virus testing splits, is a split by whether there is a Physical distancing/Biosafety (“R0 reduction”) intervention, which raises costs (mean cost = US$7,171 vs. US$1,108; Text S5) over the alternative (a weighted average of no-intervention baseline, temperature screening, and vaccination and/or boosting interventions), in line with what we would expect. Temperature screening, being similar in certain respects to virus testing, but substantially cheaper and generally substantially less effective, appears only in the tree for unavailability (Figure 5B), where its use increases unavailability, more than a maximal rate (p = 1) of viral testing (117 vs. 85=(26+164)/2).
Intervention effectiveness is highly sensitive to the degree of community transmission
Housing affects both unavailability (Figure 5B) and number of symptomatic infections (Figure 5A); both are higher when housing is “individual”. This is not to say that “shared” dormitory housing is a poorer environment for transmission than individual housing; rather, it reflects the role of community transmission in creating opportunities for reintroduction of infection from outside the employee population. This result is confirmed and elaborated by tests in Text S7, where we treat presence or absence of community transmission and presence or absence of dormitory transmission as separate factors. In line with this, additional analyses (Text S7) further indicate that our predictions about intervention effectiveness can be highly sensitive to the degree of community transmission.
“R0 reduction” strategies are cost-effective
Physical distancing/Biosafety interventions (“R0 reduction”) can reduce the number of symptomatic infections (when sufficiently effective) with minor increase in intervention expenses (Figures 5A and 5C, respectively). This suggests that highly effective R0 reduction strategies are cost-effective and, hence, should be prioritized for implementation.
Sensitivity analysis
To test the sensitivity of our model to a variety of parameters that are not user-settable, we visually and numerically examined the results when each of these parameters was varied, using a One Factor at a Time (OFAT) approach, with all other parameters (both user-settable and otherwise) at their default values, except that we examined a facility with “individual” housing. In this analysis, we focused on our three primary outcomes -- total symptomatic infections, total number of worker-shifts missed, and total cost over the simulation length -- and, additionally, on the employee-to-employee (contribution to the) effective reproduction number (Reff.) at the start of simulation. In Figure 6, we present results for 5 representative intervention scenarios, for the 5 parameters that showed the greatest sensitivity across all evaluated parameter-outcome combinations. For each of these parameters, we present the mean value for each outcome when the parameter is halved, when it is left (along with all other parameters) at its default value (Table S14), and when it is doubled. Symptomatic infections are most strongly affected by parameters defining the mean duration (μIM) of mild symptomatic infection and relative per-contact probability of transmission (βIM) during this stage), followed by a parameter (φ) governing the protection from developing symptomatic disease provided by Recovered and Hybrid immunity (Text S8B). Worker-shifts missed are most strongly affected by the relative frequency of severe infection, given any symptomatic infection (ψ), followed by the mean duration of severe infection (μIS), while the effects of these 5 parameters on Reff. and total cost were generally smaller (Text S8C).
Discussion
This study presents an ABM for tracking COVID-19 transmission and control in the food industry workforce to serve as a decision-support tool that can be used to mitigate the impacts of infectious disease outbreaks on essential worker populations and the food supply chain. The model can be customized to produce farm or processing facility settings, the type of employee housing predominantly used, as well as the vaccination and infection history and age characteristics of the workforce. Additionally, the model allows testing of a number of interventions and evaluating them counterfactually with regard to several public health and economic outcomes, and interpretation of predictions at the population and individual operation levels. The two strongest themes in our results are bimodality and trade-offs. Finally, the model also provides insights about effectiveness of different possible interventions and areas requiring further research. The developed model is expected to facilitate the food industry’s resilience and responsiveness to COVID-19 and similar future outbreaks, as well as to help navigate tradeoffs between public health and economic impacts of infection control interventions in essential industries.
Bimodality and variability in outcomes and intervention effectiveness
The simulation runs of FInd CoV Control can be interpreted to represent a population of operations with similar workforces and work environments. The simulation predicts how an outbreak would unfold in each operation (i.e., run) following infection introduction, and in counterfactual versions of the same operation that implemented different interventions. This allows us to interpret the predictions at the population level, answering questions such as: “What fraction of operations would experience certain health and economic outcomes?” Not only are most outcome distributions bimodal, but the counterfactual effects of most potential interventions are bimodal as well. Relatively ineffective (on average) interventions (e.g., temperature screening) not only sometimes appear to produce good outcomes, but also genuinely produce strong positive effects in a counterfactual sense, albeit with low probability. This is a particular consequence of a broader phenomenon: Much of the positive effect (when there is one, and especially when there is a strong one) of effective and ineffective interventions alike comes from their potential to control an outbreak at a very early stage, often before there is a single symptomatic infection. This presents a further challenge for “reactive” interventions (i.e., those implemented after the detection of a first infected case), above and beyond the issue of how quickly they can be deployed; by the time that operation managers (or policy-makers) are aware of an outbreak, the best opportunity to control it has already passed. This further supports the value of tools like FInd CoV Control that can be used as planning/forecasting tools, perhaps quarterly or on a rolling 90-day basis, to proactively prepare for the potential disease introduction.
Conversely, even fairly effective (on average) interventions can result in little or no effect in individual outbreaks. Some may even be capable, albeit with low probability, of producing counterfactually worse outcomes. This can largely be attributed to matters of timing of infection, for individuals who are infected at some point in either case; this can have an impact both through chance occurrence of opportunities for secondary transmission at particular points in time, and through the increase in probability of symptomatic infection that comes with increased time since last immunity event (i.e., last vaccination or last recovery from natural infection). Together, these possibilities, reflected in the bimodality in outcomes, further increase the (already substantial) real-world probability of misleading conclusions from anecdotal observations in individual operations, and thereby reinforce the importance of mechanistic, predictive models such as FInd Cov Control. As an alternative, data-driven evidence-based recommendations would require large-scale data collection, including time-series information on the infection spread and control (e.g., incidence of symptomatic and asymptomatic infections, results of diagnostics tests, vaccination history, infection, and isolation-related absence from work). Importantly, the data collected would need to include detailed metadata that explain the context of infection spread and control in individual work environments, because the characteristics of the work environment may serve as effect modifiers on the infection dynamics and intervention effectiveness. There are limited examples of studies of COVID-19 epidemiology in the food industry work environments (36, 37). These studies highlight the potential for rapid transmission of SARS-CoV-2 in the food industry’s work settings, which they attributed to the close arrangement of workstations and extended contact among employees. The data made available by the food operations where the investigations took place played a crucial role in enabling these studies. Our findings emphasize the need for much larger-scale data collection. However, given the cost, technical, and confidentiality-related obstacles to collecting such data, the prospect of purely data-driven decision-support models is dim. Collecting confidential data on infection spread in a food operation for the operation’s private decision-making is of course encouraged. However, it should be noted that in isolation from data on infection transmission in other comparable operations, such data will have a limited value even for the operations’ private use since the data will represent just one of many possible ways an outbreak (with or without an intervention) has unfolded as demonstrated by the model simulations. This underscores the need to proactively develop infrastructure capable of rapidly building and analysis of mechanistic or hybrid (e.g., combining ABM and machine learning (38)) models to guide infection control and policymaking under urgency and sparse data conditions.
Tradeoffs between health and economic impacts of interventions
One of the biggest trade-offs we see in our results is that of cost vs. effectiveness, particularly with respect to viral testing. Sufficiently intensive testing is highly effective at controlling transmission, but viral testing can be quite expensive, whether from the cost of test kits themselves, the cost of increased unavailability due to isolation of individuals who test positive, or both. Importantly, the intuitive solution of trying to find a moderate level of testing that optimizes this tradeoff is not necessarily a productive approach—testing, but at an insufficient frequency to achieve reliable control can actually be more expensive than either more frequent testing or not testing at all. The level of testing at which this economic “worst of both worlds” occurs is one of multiple aspects of intervention effects that is heavily dependent on the level of community transmission (modeled as individual housing that provided opportunity for acquiring infection within the community), further complicating the effort to select an optimal approach. It is helpful to realize the multiscale nature of infection control in work environments deemed essential to society, where the interest is to control the infection and its effects at the worker individual level (to protect the individual employee’s health) and at the worker population level (to reduce infection spread and cost of control, and increase labor availability). These scale- related trade-offs spill over into the trade-off between costs of control (borne primarily by the company) and effectiveness (borne by both the company and individuals), leading to inefficiencies commonly faced in the private provision of public goods (39). Designing public health policies that align operation incentives with desired public health outcomes is therefore critical to ensure the optimal provision of infection control interventions by operations. This trade-off spills over into a broader challenge around both protecting essential workers and supplying the country with food. Thus, there is a need for more discussion around essential categories of industry and appropriate metrics for evaluating “costs of control”.
Cost-effectiveness of counterfactual interventions
At an individual operation level, FInd CoV Control can be used preemptively to ask questions such as: “Given the characteristics of the workforce and work environment in my operation, if an infected worker enters my facility in the near future, how likely it is that we will experience an outbreak?” “If we have an outbreak, how likely (in terms of the measures of central tendency and variation) are health and economic impacts under different intervention scenarios?” FInd CoV Control evaluates cost-effectiveness of 12 intervention scenarios and a no-intervention scenario. As seen in Figure 2 and Figure 3, even intensification of vaccination in an already moderately vaccinated population can yield modest but meaningful benefits. For proactive control of COVID-19 in the food industry, maintaining a vaccinated and boosted workforce to be prepared for a new outbreak remains a cost-effective intervention (albeit not sufficient to make other interventions unnecessary). Vaccination uptake can be increased by removing convenience and confidence barriers and leveraging workers’ motivation to protect self, family, and community (40). Our model findings seem to resonate with perceptions of the food industry’s leadership. Certain companies, particularly those in the labor-intensive meatpacking sector, took proactive measures by mandating vaccinations for their workforce (41–43) and prevented risk of infection among employees and potential plant closures that were prevalent at the beginning of the pandemic (44). In situations where vaccine mandates were not in place, other strategies such as physical distancing requirements and quarantines were implemented; however, these interventions were reported to lead to high worker absenteeism and hindered the efficient operation of processing plants (1). Strategies like screenings for disease were valuable in controlling workplace transmission, but also had serious limitations regarding reliability of the results and posed challenges due to being labor-intensive and costly compared to simpler strategies like use of face coverings and practicing personal hygiene (45). These observations in the food industry match the findings in our model about the effectiveness of screening, testing and physical distancing/biosafety strategies in preventing symptomatic cases, albeit with high expenses and production losses associated with testing strategies. Effectiveness can also trade off against “costs” that are not strictly monetary. For example, even very intensive physical distancing and biosafety measures may be cost-effective, but some aspects of such measures (e.g., masking and face shields) can be highly unpopular in the long run (and even limit productivity in harsh work environments, such as extreme cold, hot and wet/damp), especially given the previous observation of the limitations of reactive interventions (46). A comparison between the effects of highly and moderately effective interventions in this category supports the idea that, if one is going to implement masking and face shields, the addition of ventilation improvements is likely to be cost-effective. Our model revealed that predictions about intervention effectiveness are highly sensitive to the degree of community transmission. This emphasizes the importance of interpreting effectiveness of work-based interventions in the light of the epidemiology of the disease in the community. This also emphasizes the importance of mitigating disease spread outside of work; however, this is particularly challenging for agricultural worker populations that are not stationary and typically share housing and transportation, allowing for easy employment- related transmission of the virus (47, 48).
Limitations and future directions
The empirical support for model parameter values varies, and this is of particular concern for parameters identified as influential in the Sensitivity analysis (Figure 6). One area of particular concern, given a combination of high sensitivity, moderate support (i.e., a good amount of data, but somewhat coarse-grained, and with significant potential for confounding), and a history of changing strains, is the magnitude of long-lasting immunity provided by boosting (e.g., parameter φ), revealing that this is a critical knowledge gap requiring further research. Some more structural limitations to our model include our relatively simple model of change in infectiousness over the course of infection, our assumption that voluntary self-isolation is rare enough in essential workers to be omitted from the model, our relatively simple model of vaccination and boosting interventions (exponential decay of the eligible but unvaccinated/unboosted), and a simplified binary notion of housing for a given operation as being either shared or individual.
There are numerous possible refinements of this model, many of which are facilitated to a greater or lesser extent by its modular structure. Three areas seem particularly likely to be fruitful: (i) Replacement of the current discrete-staged model of infectiousness over time with continuous infectiousness curves, analogous to the continuous curves that make up aspects of our continuous immunity trajectories; (ii) Continued improvement of our model of immune effects, immune boosting, and immune waning, as well as accounting for the changes in vaccination guidelines; and (iii) Incorporating multiple simultaneous interventions (starting at different times) and incorporating mixed individual and shared employee housing. More broadly, we hope to further increase the modularity, flexibility, and ease-of-use of the model, to facilitate easy modification to address other respiratory pathogens and/or other critical infrastructure sectors to enhance resilience and responsiveness to similar future outbreak events. Towards these goals, we created a user-friendly web interface for an early version of FInd CoV Control described in this article, which allows the user to customize it to the characteristics of their workforce and generate a clear and easily interpretable confidential result (49).
Materials and Methods
Employee population module
We model a heterogeneous population of agents (employees) with a variety of attributes reflecting both their current state and certain aspects of their personal history. Agent attributes set at the simulation start are summarized in Table 1. Attributes that represent past events and current state include age, (directly) immunity-related attributes, vaccination history, and the current state of infection, if any.
Age
Agents are randomly assigned an age category, with probabilities that are derived from industry-wide data about the age of agricultural workers (50). This age category is then used to determine their probabilities, in the absence of immunity, of experiencing symptoms or dying (Table S1).
Immunity-related attributes
In this model, all of an agent’s attributes that are directly relevant to that agent’s immunity, and that are not a consequence of their age, represent acquired immunity (whether complete or partial) to SARS-CoV-2 infection and COVID-19 disease. These attributes are associated with specific past events that created or boosted that agent’s immunity; for ease of reference, we refer to these simply as immunity events (Table S2). Thus, all agents had, conceptually speaking, an immune status of fully Susceptible (S) and a vaccination status of Not Vaccinated (NV, also referred to as unvaccinated) at the start of the COVID-19 pandemic, but some may have different values for one or both of these attributes at the start of the simulation.
Immunity-related attributes include the time (tlast, i) of the most recent immunity event, the time (tR, i) of the most recent recovery from natural infection (if any), an “immunity trajectory” (Ci, representing a trajectory defined by a combination of vaccination status and whether the agent has previously recovered from natural infection), and what the agent’s level of immunity was immediately prior to their most recent immunity event (Plast, E, i, Plast, IP, i, and/or Plast, IS, i) (Table 1). Together, these determine an agent’s Protection against Any Infection (PE, i), Protection against Symptomatic Infection given Any Infection (PIP | E, i), and Protection against Severe Infection given Symptomatic Infection (PIS | IP, i) at the present time (Table 1 and Text S3B). Each immunity trajectory includes curves for PE, i, PIP | E, i, and PIS | IP, i. These may include an initial increase in immunity, referred to as “ramp-up”; an initial period of total immunity; and/or immune waning. During both ramp-up and immune waning, tlast, i is relevant; during complete immunity (for immunity trajectories R, HV1, HV2, and HB only), tR, i is relevant (specifically, in determining that there is complete immunity); and during ramp-up, Plast, E, i, Plast, IP, i, and/or Plast, IS, i are also relevant. These attributes are discussed further in the section “Immune dynamics” and in Text S3B and Tables S4-S7.
Vaccination History
We assume for the sake of simplicity that all vaccination and boosting uses the monovalent Pfizer vaccine. Therefore, vaccination history includes the number of doses of vaccine the agent has received (i.e., whether they are unvaccinated (NV), partially vaccinated (V1), fully (primarily) vaccinated (V2), or boosted (B)), and when they received each of their previous doses, if any (Figure 1B.iv). Generally, only the time of their most recent dose is relevant to the model dynamics (Table S4). The current version does not account for repeated boosting in the vaccination history.
State of Infection
An agent’s current infection status can be: Not Infected (NI), infected, but not yet infectious (E, for “Exposed”), Asymptomatically Infectious (IA), Presymptomatically Infectious (IP), Mildly symptomatic (IM), Severely symptomatic (IS), Critically symptomatic (IC), or Dead (D) (Figure 1B.iv). Additionally, if that state is anything other than Not Infected, we record how long they have been in that state; together with their precalculated duration (see below) for that state, this determines how much longer they will remain in it before progressing, recovering, or dying.
Initial setting of agent infection and vaccination history
Agents’ history of immunity events prior to the start of simulation is important for immunity (and immune ramp-up and waning), and their history of vaccination events specifically is important for their eligibility for future vaccination (both primary and boosting). Therefore, user input (Table 2), is used to determine what fraction of the population has and has not experienced various immunity events and when. We then randomly generate exact times in a simple fashion. This aspect of run initialization involves (i) current infection status initialization, (ii) infection history initialization, and (iii) vaccination history initialization. These aspects do not directly interact with each other, although infection history and vaccination history interact in their effects on immunity. The initial settings are described in Text S2.
Disease transmission module
Course of infection
The main aspects of the course(s) of infection are summarized in Figure 1B.iv. Agents are categorized as Not Infected (NI); Exposed (E); Infectious (I); or Dead (D). Infectious agents are further divided by whether they currently have clinical disease (IA and IP vs. IM, IS, and IC), whether they will subsequently develop clinical disease (IA vs. IP), and/or how severe their symptoms are (IM vs. IS vs. IC). They are further categorized by their immunity trajectory – fully Susceptible (S); one of the Vaccinated trajectories (V1, V2, or B); Recovered (R); or one of the Hybrid (H) immunity trajectories (HV1, HV2, or HB) – and (for all immunity trajectories other than S) the times since their last immunity event (their entry or reentry into that immunity trajectory) and their last recovery, if any (see “Immune dynamics” section).
Infectible agents (i.e., agents whose infection status is NI, and whose susceptibility to infection (see “Transmission model” section) is greater than 0) can acquire infection with SARS- CoV-2 from a potentially infectious contact with an Infectious agent (i.e., IA, IP, IM, IS, or IC). The probability of such a contact resulting in an infection depends on the susceptibility to infection of the infectible agent (1 – PE, i) and the infectiousness of the infectious agent (βIA, βIP, or βIM). How contacts are made is described in the “Work environment module” section.
Upon infection, the formerly infectible agent enters the Exposed state (E), and after a period of time (DE, i), the Exposed agent becomes infectious (I). At this time, they enter one of two states, both of which are continuously infectious, but without clinical disease: either IA or IP, the latter of which is the first stage of the symptomatic path. The distinction between IP and IA reflects a substantial difference in per-contact transmission rates between the two (51).
All IA agents are assumed to recover following a period of time (DIA, i). Agents taking the symptomatic path progress through up to four stages: IP, IM, IS, and/or IC. All IP agents progress to the IM state after a period of time (DIP, i), but agents who are in IM, IS, or IC may, after the corresponding period of time (DIM, i, DIS, i, or DIC, i, respectively), either continue their disease progression to the next state on the symptomatic path or recover.
When an agent recovers, regardless of which infectious state they recovered from, their last immunity event time (tlast, i) and Recovered immunity event time (tR, i) are both set equal to the current time, their infection status is set to NI, and their immune status is set to Recovered (R) if they have never been vaccinated (NV), or to the appropriate Hybrid state (HV1, HV2, or HB) otherwise. For agents in IC that do not recover, the next step is Death (D).
We do not explicitly model individual symptoms such as fever, cough, etc. However, for the temperature testing intervention, we do tacitly assume that fever is only present if the individual is symptomatically infected (IM, IS, or IC). We define “Severe” symptoms as those requiring hospitalization; consequently, we assume that only agents with an infection status of IA, IP or IM can transmit to their fellow employees, because agents with Severe or Critical symptoms are so sick that they require hospitalization. Relative transmissibility (per contact) is set based on the infection stage (Table S3). Absolute transmissibility has no effect in the model, as we set the average expected contact rate in order to achieve a specified basic reproduction number (R0), and our assumptions about contacts per day (Text S3A) make the distinction between twice as high transmissibility per contact and twice as many expected contacts mathematically irrelevant.
For agents in any of the infected states, disease progression is based on the duration of each state (Table 4), age-dependent baseline probabilities of entering each disease state during disease progression (Table S1), and immunity-dependent modification of those baseline probabilities.
Transmission model
Agents who are infectible (i.e., agents whose infection status is NI, and whose susceptibility to infection is greater than 0) can be infected by contacts with either infectious coworkers or infectious people outside of work, in the broader community (if housing is “individual”). Contact structures of the agents while at work were determined by the place of agents in the hierarchical structure of the farm or facility and their work schedule. To be more precise, for each shift type (e.g., weekday Production Shift 1, weekend Cleaning Shift, etc.), we have a matrix of expected contact rates between pairs of agents. For some combinations of a pair of agents and a shift type, a non-zero contact rate may represent contacts made at work; for others, it may represent contacts made in shared housing. Further details are in Text S3A.
Immune dynamics
We distinguish between 8 basic states (immunity trajectories) with respect to immunity: fully Susceptible (S), partially vaccinated (i.e., with one dose of a two-dose primary series) (V1), Fully Vaccinated (for the purpose of this study defined as a 2-dose primary series; V2), Boosted (B), Recovered (R), and Hybrid immunity (H) with partial vaccination (HV1), with full vaccination (HV2), or with full vaccination and a Booster (HB). Non-hybrid vaccinated trajectories (V1, V2, and B) feature a smooth ramp-up from their individual’s previous level of immunity, that lasts for TV1→V2 = 21 days, Tramp, V2 = 14 days, or Tramp, B, 1 + Tramp, B, 2 = 14 days, respectively, counting from the time since the individual’s last immunity event (i.e., first vaccine dose, second vaccine dose, or booster shot, respectively). The non-hybrid Recovered trajectory features an initial Ttotal, R = two months (61 days) period of total immunity, counting from the time of their (most recent) recovery. The Hybrid immunity trajectories (HV1, HV2, and HB) have characteristics of both vaccinated and recovered trajectories, and can be entered either by recovery following vaccination or by vaccination following recovery. Consequently, a particular individual’s experience of one of these trajectories may include either or both of total immunity and ramp-up. The transitions between these immunity trajectories are summarized in Table S2, and the equations for the protection that they offer are given in Table S7. As evidence has mounted that waning immunity, both from natural infection and from vaccination, plays an important role in the dynamics of transmission during the pandemic, we included this in our model. To take full advantage of the agent-based nature of our model, we assigned each agent three variables indicating the major factors influencing their level of susceptibility to both infection and progression: (i) the immune state that they entered at the time of their last immunity event, (ii) the time at which that event occurred (and hence, at any given subsequent point in time, how long it has been since that event), and (iii) the time of their last recovery from natural infection. The exceptions are agents whose immune state is fully Susceptible (S), whose time of last event is not defined, and agents whose immune state is either fully Susceptible (S) or one of the non-hybrid Vaccinated states (V1, V2, and B), whose time of last recovery is not defined. We then created functions (Table S7) giving, for any valid combination of state, time since entry, and time since recovery (and previous immunity, if they are currently in a ramp-up period), their level of relative protection from each of infection, symptomatic infection conditional on any infection, and severe infection conditional on symptomatic infection. We used exponential or exponential- mixture waning for long-term behavior of V2 and B (fitted from data in (52)), and logistic waning for long-term behavior of R, HV1, HV2, and HB (with parameters inferred from the tables in (53)), with some special case behavior at the start of states other than S (ramp-up and/or a period of complete immunity), in order to account for the delay in reaching full protection following vaccination, and to prevent unrealistic cycles of extremely rapid reinfection. This is further explained in Text S3B.
General model of vaccination
Agents’ vaccination status can be unvaccinated (NV), partially vaccinated (V1), fully vaccinated (V2), or boosted (B). This vaccination status directly corresponds to their immune status (with NV corresponding to fully Susceptible (S)) if they have never recovered from a natural infection; if they have, then their immune status is either Recovered (R), if they have never been vaccinated, or one of the Hybrid immunity trajectories (HV1, HV2, and HB).
Vaccination trajectories
Partially vaccinated agents become eligible to receive a second shot TV1 -> V2 = 21 days after receiving their first one, and fully vaccinated agents become eligible to receive a booster shot TV2→B = 5 months (treated as a deterministic 152 days) after their second shot (54). We assume that all V1 agents enter V2 states as soon as they are eligible, but that only a fraction of V2 agents enters the B state as soon as they are eligible. We further assume that agents who are eligible to become V1/B at simulation start, but have not yet done so, will not become V1/B (respectively) during the simulation, in the absence of an intervention to promote primary vaccination/boosting, respectively.
Work environment module
In both farm and facility models, each week consists of 5 work days and 2 non-working days (i.e., a typical “work week” and weekend in the US). Each calendar day is modeled as consisting of 3 eight-hour periods we call “shifts.” Each agent spends two shifts awake, and one asleep. For simplicity, we tie the agent’s sleep cycle to their work cycle, so that each agent spends their first shift awake at work if it’s a work day. This structure is illustrated for a sample agent (one scheduled to work on the first shift of the calendar day) in Figure 1B.i and shown in greater detail in Table S8.
In the facility model, each working day consists of Production Shift 1, Production Shift 2 (which may actually be a non-working shift, if the facility in question only has one production shift per day), and a Cleaning Shift. In the farm model, all agents are scheduled to work on the same shift, which, by analogy with the facility model, we refer to as Production Shift 1. (Hence, by the same analogy, all agents are awake, but not working, during Production Shift 2, and asleep during the shift that would be the “Cleaning Shift” in the facility model.)
We distinguish between available and unavailable agents for the purposes of SARS-CoV- 2 transmission between agents. Available agents are available to work their scheduled shifts – meaning that they can do work (relevant in the economic analyses), and can also potentially infect, or be infected by, other agents. Agents are available by default, but may become unavailable to work (and, subsequently, may become available again).
In the baseline (no intervention) model, unavailable agents are limited to those who are either hospitalized (i.e., those who have an infection status of IS or IC) or dead (those who have an infection state of D). Hospitalized agents become available again upon recovery. Under certain interventions, agents who are not hospitalized may become diagnosed and isolated; these remain unavailable until they are deisolated (Text S3C).
We model a facility or farm with a hierarchical organization, although the details of this hierarchy differ somewhat between the facility model and the farm model. This hierarchy is assumed to be fixed over the time horizon of the model, as is the associated work schedule. The expected number of contacts between pairs of workers who are both “available” (i.e., not isolated, hospitalized, or dead) is likewise held constant (i.e., we do not reassign workers between work crews or production lines based on other workers’ absences). To make the interface more manageable, we assume that the structure of the operation is “regular” in the sense that if there are two or more of the same sub-structure (e.g., two production shifts, two or more teams of work crews, two or more work crews within a team, or two or more production lines), then each of those sub-structures has the same structural characteristics (e.g., if one work crew consists of 10 workers and a foreman, than all work crews consist of 10 workers and a foreman).
Features specific to the farm and facility models are described in Text S4A and S4B. Briefly, we assume the following:
All workers live in the same type of housing (i.e., either individual or employer-provided shared housing).
All workers work a regular 40-hour, 5-day work week (8-hour shifts) and 2-day weekend, with the model simulation starting on a random day of the week, except for a small number of floating workers (e.g., quality assurance technician, mechanic) in the facility model, for whom the work shift(s) on a given day is/are randomly selected.
There are many more contacts within the hierarchical structure than outside it (e.g., more contacts between workers on the same crew/production line than between workers on different crews/production lines, more contacts between foremen and supervisors than between other workers and supervisors, etc.), but that contacts are possible between any two agents who are both present on the same shift.
All contact between workers occurs either (a) while traveling to, at, or traveling home from work, or (b) in shared, employer-provided housing, i.e., that workers in individual housing do not socialize with each other outside of work.
There is homogeneous mixing within employer-provided housing.
Worker contacts on the way to and from work follow the same basic patterns as worker contacts at work (e.g., we tacitly assume that shared transportation is substantially more likely to group together workers who are on the same crew than workers who are on different crews).
Interventions
Apart from the baseline (non-) intervention, we modeled 12 interventions falling into four groups (described in detail in Text S3C):
“Temperature screening” (1 scenario) denotes temperature screening and isolation where all employees arriving for their shift are tested prior to admitting employees to work for that shift. The temperature threshold is set at 38°C.
“Virus test” (3 scenarios) denotes low-, moderate-, and high-intensity of viral testing and isolation, respectively modeled with the probability p = 5, 30, or 100% of each employee being tested each shift upon arrival at the workplace.
“Vaccine” (5 scenarios) denotes primary “vaccination” of unvaccinated workers with 2 doses at a daily probability of p = 2 or 4%; “boosting” of boosting-eligible but unboosted workers at p =2 or 4% per day; or a combination of both primary and boosting vaccination interventions (“vax+boosting”) at p = 2% per day.
“Physical distancing/Biosafety” (3 scenarios) where physical distancing and/or biosafety interventions are modeled as generating a 20, 40, or 80% reduction in R0 at work (i.e., -20, -40 or -80% R0), such as through the application of masks, face shields, physical barriers, and/or ventilation. The exact approach necessary to achieve these desired effects on R0 will vary across work environments. However, to illustrate how higher-effectiveness interventions may be built by “stacking” multiple lower-effectiveness strategies and to be able to quantify their net cost in the subsequent economic analysis, in absence of required data we assumed: (i) For a low-effectiveness (-20% R0) intervention, the use of KN95 masks, one per employee per shift; (ii) For a moderate-effectiveness (-40% R0) intervention, a low-effectiveness intervention combined with the use of face shields, one per employee per 30 days; (iii) For a high-effectiveness (-80% R0) intervention, a moderate-effectiveness intervention combined with ventilation improvement, such as with the use of (a) portable air cleaner(s).
We performed counterfactual comparisons of the no-intervention baseline and 12 interventions (Figure 1A). To make these comparisons more precise, on a run-by-run basis, we reseed the pseudorandom number generator (pRNG) with the same value before processing each intervention and (b) make use of the pRNG in such a way as to ensure that two runs that start with the same pRNG state and that differ only in the factors that we allow to vary between possible interventions under a given scenario, will have the same pRNG state at all analogous points thereafter.
Model running
For each run, the model is first initialized. Then, at each time step, the following processes occur:
Agents eligible for deisolation (Text S3C) are deisolated.
If any testing is being performed (Text S3C), agents who are scheduled to work and (potentially) available are tested.
o If the testing probability per shift is 1, then all (potentially) available agents are tested.
o If the testing probability per shift is < 1, then the number of tests to be performed is determined, and these are performed on the (potentially) available agents in order from least to most recently tested, randomizing ties.
o If any agents test positive, they are isolated, and their isolation time is set to the present time.
Agents to be vaccinated are randomly selected (meaning that they are not infected and either they have just become eligible for a booster and are boosting on time, or they are eligible to receive some form of vaccine, and there is a vaccination-promoting intervention), and their immunity event times, immunity trajectories, and vaccination statuses are updated accordingly.
Transmission (potentially) occurs, with probabilities as described in Text S3A.
Infected agents who are eligible to leave their current state of infection (i.e., the sum of their time of entering that state and their precalculated duration for that state is less than the time at which the current time step ends) do so, randomly selecting which new infection state to enter, if necessary, as described above.
o This step is repeated as necessary, i.e., an agent may in principle progress twice in a single time step if the duration of one of their infection states is sufficiently small.
o If an agent recovers, their last recovery and last immunity event times (tR, i and tlast, i) are set to the present time, and their immunity trajectory is updated.
Outcomes are recorded for use in subsequent analysis.
The ABM model and all analyses were implemented in R software (version 4.0.4) (58).
Model outcomes
Our outcomes of interest can be broadly grouped into those which pertain primarily to public health, and those which pertain primarily to business disruptions and economic impact. We compare these outcomes both across (potential) interventions within a setting scenario and across different scenarios.
Public health outcomes
Our primary public health outcome of interest is the total number of symptomatic SARS- CoV-2 infections that occur across the course of a run. We define this as the number of occasions of an agent transitioning from the IP to IM state between the start of simulation (t = 0) and the end of simulation (t = T). This includes transitions into IM of agents who were in the E state at simulation start, if they transition from E to IP (and, subsequently, to IM) rather than to IA, but does not include the transitions into IM before t = 0 for agents that were in IM at simulation start. It also does not include the transitions into IM after t = T of any agents who are in E or IP at simulation end. If the simulation length and transmission dynamics are such as to result in some agents transitioning into IM more than once during the simulation, then we count each such transition separately; hence, it may be possible, depending on user-set parameters, to have a number of symptomatic infections that exceeds the number of agents. An additional outcome of interest is the number of total SARS-CoV-2 infections (i.e., including asymptomatic infections) that occur across the course of a run. We present public health outcomes, for each intervention, as a curve of the mean over time, as well as violin plots summarizing the cumulative distribution over a defined planning horizon, to allow both comparison of expectations and understanding of the variance.
Economics outcomes
The main dimensions considered for economic analysis are: (i) direct costs of performing an intervention, and (ii) the productivity costs or benefits of the intervention. To assess (i), we add certain details to our hypothetical interventions, based on how interventions are implemented in real life. Relative to a “no-intervention” baseline of doing nothing, the direct costs of interventions can be expected to always be greater than or equal to zero. Each intervention is compared to the baseline to estimate direct costs of performing each intervention. Calculations of the costs of interventions are described in Text S5. To answer (ii), we estimate the productivity loss based on the worker absences from the infection model with a Cobb-Douglass production model. The
Cobb-Douglass production function takes the general form where operations use labor L and capital K to produce output Q, and A is some constant. The output elasticities of labor and physical capital are β and α, respectively. In choosing output level Q, operations face the following cost function (abstracting from fixed costs): where w is the wage rate and r is the rental rate for capital. Profit-maximizing operations operate at the efficient production frontier in equilibrium.
Based on a previous survey of producers (17), one important assumption made in the model is that the facility can maintain full production (Qf) with up to 15% of absenteeism without incurring appreciably higher production costs. If a operation can reduce its “full output” labor input allocation Lf to a short-staffed model where Ls = 0.85Lf without reducing full output (Qs = Qf) in a costless way, the operation must be able to substitute enough additional capital above the current “full output” allocation in the short run, such that the cost of increased capital is completely offset by reduced labor costs. To capture this short-run flexibility, we therefore assume that the current equilibrium is not unique, but instead one in a set of cost-minimizing equilibria where 0.85Lf ≤ L ≤ Lf.
If labor absenteeism exceeds 15%, however, we assume that the capital allocation now remains fixed at the (higher) 15% absenteeism point, and the production paradigm takes the previously described Cobb-Douglass form. Under that framework, since operations cannot substitute capital for labor, they will have to reduce output based on these labor shortages. In an empirical work application of productivity analysis in the food and agricultural sector in the US, Ahmad (59) estimates a Bayesian stochastic Cobb-Douglas production function using US state- level agricultural data from 1960-2004 (n=2,160). He estimates α = 0.316 [0.271, 0.362] and β = 0.437 [0.392, 0.483]. Thus, production is Full production quantity is provided by the users, and estimated production loss is simply The results of the above analyses are summarized by the following economic outcomes: the mean number of employees unavailable to work production shifts over time; violin plots summarizing the distribution of the cumulative number of production worker-shifts missed; fraction of shifts short; and total direct expenses, production losses, and total costs associated with an intervention (US$).
Scenario analysis
We tested scenarios corresponding to factors: “setting” (“farm” vs. “facility”), “housing” (“individual” vs. “shared”), “vaccinated” (“high” vs. “none”), and “recovered” (“high” vs. “nothing”). In a full factorial analysis approach, we tested all 16 combinations of these four factors and for each scenario we ran the no-intervention baseline and all 12 interventions. We then constructed regression trees, using the R package rpart (version 4.1.19; using default values for all control parameters) for each of our three primary outcomes (symptomatic infections, worker-shifts unavailable, and total cost) and the two separate contributors to total cost (production losses and intervention expenses) vs. the four scenario parameters and the 5 intervention parameters defining our 12 interventions (Table S13). In the main text of this paper, we only present a subset of results for which setting is “facility,” for reasons discussed in the “Scenario analysis” subsection of the Results section.
Sensitivity analysis
We evaluated the parameter sensitivity for four outcomes using the OFAT approach. Three outcomes were the total symptomatic infections, total number of worker-shifts missed and total cost over the simulation length. The fourth outcome was the effective reproduction number (Reff.) at the start of simulation. Because of the large number of parameters to be examined, under (initially) a variety of scenarios, we averaged results over batches of 100 runs each, rather than the 1000 runs that we use in most other contexts. In general, for a given parameter whose value in our model was x, we examined results when that parameter was set to 1/2 x, x, and 2x (Table S14). For 2 parameters (TV2→B and T1st V2 -> 0), one of these values was impossible, given a requirement that it be possible for some individuals to be eligible for boosting at simulation start. For these parameters, we set the relevant value to the most extreme possible value in the same direction (i.e., 212.5 days for (increased) TV2→B and 305 days for (reduced) T1st V2 -> 0).
In order to estimate Reff. at the start of simulation (other outcomes were measured without this modification), we modified the model so that, when an employee was infected, instead of their infection status being changed from Not Infected to Exposed, their immunity was set to Recovered, and their tlast, i was reset to the current time. Conceptually, this amounts to their “skipping over” any infected states, and going directly from infection to recovery. This ensured that all employee-to-employee transmissions would necessarily be from employees who were infected at the start of simulation. We then ran the model starting with a single Exposed employee (as is the default), and observed the number of employee-to-employee transmissions that occurred. The average of that number across multiple runs was then used as an empirical estimate of Reff. In scenarios with “individual” housing, the community-acquired infection was allowed to occur (likewise transferring infected employees directly into a Recovered state, meaning that they were not immediately infectible by the index case), but was not included in the count of employee-to-employee transmissions measured by Reff.
We initially examined the results of sensitivity analysis, using default user-settable parameters, except as noted, across the same range of settings as we initially considered for scenario analysis. We found that, for most parameters, sensitivity was higher in “high” settings (i.e., with a history of both vaccination and recovery from natural infection) than in historical or pseudo-historical settings, higher in facility settings than in farm settings, and higher in settings with shared housing than in settings with individual housing. For this reason, we focused our further examination on a setting of a facility, with all parameters default except for housing, which was “shared” instead of “individual” (with physical distancing in shared housing at the default level). Looking back at scenario analysis, this is equivalent to the scenario in which both “vaccination” and “recovered” were “high.”
For each parameter-intervention-outcome combination, we obtained three mean outcome values, as described above. We divided the largest of those three values by the smallest to obtain a simple measure ω of the sensitivity that could address both monotonic and non-monotonic effects. For each parameter, we then chose the largest such measure (ωmax) across all intervention- outcome combinations. We then selected for depiction and discussion in the Results section those parameters for which ωmax ≥ 2 (all results are shown in Figures S5-S8). As the only exception to this process, we omitted the baseline no-intervention scenario in the Total Cost outcome, as the sensitivity measure for this outcome-intervention was generally undefined or (rarely, and seemingly at random) infinite, due to the extreme rarity of non-zero Total Cost at baseline.
Data Availability
All data are available in the main text or the supplementary materials. R code is available under GPL-2.0 license at https://github.com/IvanekLab/covid-journal-code
Funding
Agriculture and Food Research Initiative from the USDA National Institute of Food and Agriculture, Competitive Grant no. 2020-68006-32875 (RI, AA, SA, MW).
Cornell Institute for Digital Agriculture (CIDA) (RI, MW)
Artificial Intelligence Institute for Next Generation Food Systems from the USDA National Institute of Food and Agriculture, grant 2020-67021-32855 (RI, MW)
The funders had no role in the study design, data collection, and analysis, decision to publish, or preparation of the manuscript.
Author contributions
Conceptualization: CH, RI, AA
Methodology: CH, RI, AA, CZ, DW
Investigation: CH, RI, EB, SIM
Visualization: CH, RI
Supervision: RI, AA
Writing—original draft: CH, RI
Writing—review & editing: CH, EB, SIM, CZ, AA, DW, MW, SA, RI
Competing interests
CZ is employed at iFoodDecisionSciences, Inc., which hosted the web-interface developed as an extension of a previous version of FInd CoV Control described here (49) until December 2023; this interface is currently hosted by iDecisionSciences, LLC. FInd CoV Control was licensed to iDecisionSciences, Inc. under the GNU General Public License. DW is the founder of iFoodDecisionSciences, Inc and iDecisionSciences, LLC. All other authors declare they have no competing interests.
Data and materials availability
All data are available in the main text or the supplementary materials. R code is available under GPL-2.0 license at https://github.com/IvanekLab/covid-journal-code
Supplementary Materials
Supplementary Materials enclosed.
Acknowledgments
The authors thank the Industry Advisory Council, which was assembled to guide modeling work and comprised of executive-level leaders in the major US produce farms and food processing companies and an official in the Centers for Disease Control and Prevention, for valuable feedback on the structure of FInd CoV Control. The authors also thank Ms. Yu Tang (Dyson School of Applied Economics and Management, Cornell SC Johnson College of Business, Cornell University, Ithaca, NY, USA) for her contribution to the methodology, investigation and writing of the original draft for the evaluation of economic impacts of interventions.