A spatially explicit model for estimating risks of pesticide exposure to bird populations

Pesticides are used widely in agriculture and have the potential to affect non-target organisms, including birds. We developed an integrated modeling system to allow for spatially-explicit evaluation of potential impacts to bird populations following exposures to pesticides. Our novel methodology builds upon three existing models: the Terrestrial Investigation Model (TIM), the Markov Chain Nest Productivity Model (MCnest), and HexSim to simulate population dynamics. We parameterized the integrated modeling system using information required under the Federal Insecticide, Fungicide, and Rodenticide Act, together with species habitat and life history data available from the scientific literature as well as landcover data representing agricultural areas and species habitat. Our case study of the federally threatened California Gnatcatcher (Polioptila californica) illustrates how the integrated modeling system can estimate the population-scale consequences of pesticide applications. We simulated impacts from two insecticides applied to wheat: one causing mortality (survival stressor), and the other causing reproductive failure (reproductive stressor). We observed declines in simulated gnatcatcher abundance and changes in the species’ distribution following applications of each pesticide; however, the impacts of the two pesticides were different. Our methodology attempts to strike a balance between biological realism and model complexity and should be applicable to a wide array of species, systems, and stressors.


Introduction
Pesticides are used widely in urban, suburban, and rural environments and have the potential to affect organisms other than targeted pests [1][2][3][4]. From 2008-2012, the most recent years for which statistics are available, over 1 billion lbs of pesticide active ingredient (a.i.) were sold a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 annually in the United States [5]. Insecticides accounted for about 64 million lbs a.i. sold, of which approximately 64%, or 38 million lbs a.i./yr were used in the agricultural sector. At least 29 different modes of action (MOA) have been identified for agricultural insecticides, which target different aspects of the biology and physiology of target pests, including the nerve and muscular systems, growth, respiration, and gut function [6].
To identify and characterize risk of pesticides to birds and other non-target organisms, the United States Environmental Protection Agency (USEPA) Office of Pesticide Programs (OPP) uses a tiered risk assessment process [7]. Current USEPA/OPP tiered risk assessment for birds starts with the T-REX model [8], which provides a conservative estimate of exposure through diet. Tier I assessments are based on risk quotients (RQs), which are calculated by dividing a conservative estimate of exposure by a threshold toxicity value representing mortality or sublethal effects (e.g., the median lethal dose or LD50; the median lethal concentration or LC50; the No Observed Adverse Effect Concentration or NOAEC; the No Observed Adverse Effect Level or NOAEL). RQs are compared to levels of concern (LOCs) to determine whether a pesticide use poses risks of concern for mortality, growth or reproductive effects. If RQs are below the acute and chronic LOCs, it is concluded that the insecticide does not pose a risk of concern to the species or taxon. In cases where a RQ exceeds a LOC, there is potential risk of effects, and a higher tier assessment may be warranted. For higher tier avian risk assessments, one option available for OPP is the Terrestrial Investigation Model (TIM, v.3.0) [9], a probabilistic model that focuses on acute exposures to birds. For chronic risks, the Markov Chain Nest Productivity Model (MCnest) [10,11], a mechanistic model of avian breeding, is available to assess potential declines in the annual reproductive success of exposed bird populations. The T-REX, TIM, and MCnest models are focused on field-level exposures to birds and resulting effects. The models assess risks to individuals or small groups of birds (e.g., flocks) that are assumed to be exposed to a pesticide on a use site (e.g., wheat field) or adjacent area receiving spray drift. Under Section 7 of the U.S. Endangered Species Act (ESA) [12][13][14] EPA makes effects determinations to assess whether a pesticide's use may affect a listed (i.e., threatened or endangered) species, and, if so, the US Fish and Wildlife Service (FWS) or the National Marine Fisheries Service (NMFS) may evaluate if a pesticide's use jeopardizes the persistence of that species. Among the tools and methods under development as part of this process are population models, which have long been recognized as a means for understanding the risks and consequences of adverse effects of pesticides on birds and other animals [15,16]. A report from the US National Research Council (NRC) strongly endorsed the use of population models for assessing the risks to threatened and endangered species from exposure to pesticides (NRC 2013). Several population modeling approaches have been applied to avian chemical risk assessment, including matrix projection models [17] and individual-based models [18].
The California Gnatcatcher (Polioptila californica, hereafter CAGN) is an insectivorous songbird that occupies scrub habitats in the Baja peninsula of Mexico and southern California, ranging as far north as Ventura County, CA, USA [19]. CAGN was officially listed as threatened under the ESA in 1993. Principal threats to CAGN include habitat loss and fragmentation, increased frequency of fire, and nest parasitism by brown-headed cowbirds (Molothrus ater) [20]. The extent to which CAGN may be affected by exposure to insecticides has not been studied, though the California Department of Pesticide Regulation lists over 60 crops on which pesticides are used within the range of the species [21]. CAGN typically breed from late March through mid-July and their seasonal reproductive output usually includes multiple nesting attempts by a given breeding pair [19].
In this article, we describe a method for conducting spatially explicit population level risk assessment for birds. This work builds on that of Etterson and colleagues [4] who combined TIM and MCnest into a single unified model (TIM/MCnest) that considers both acute and chronic effects resulting from insecticide exposure during the northern temperate breeding season (approximately March-August). Here we show how breeding season predictions of insecticide effects from TIM/MCnest can be used in the HexSim modeling environment [22] to inform a spatially explicit population model for avian response to insecticide exposure. Hex-Sim is user-friendly development platform within which researchers can construct spatially explicit and individual-based simulation models. Importantly, HexSim allows users to develop models that conform to data limitations, rather than requiring a fixed set of input parameters. We parameterize the models for two insecticides and apply them to study the potential effects of insecticide exposure on the California Gnatcatcher. Our objectives are 1) to create an integrated workflow that allows TIM/MCnest predictions of pesticide effects on avian reproductive success to be implemented within the HexSim environment; and 2) to apply the model to a federally listed species under realistic spatially referenced insecticide use patterns. Our presentation of the California Gnatcatcher simulations is intended as an example of how spatially explicit population level risk assessment for pesticides might be performed for a federally listed bird. Because the simulations we present are not, in their current form, intended to influence policy or management, the pesticides examined in this analysis are unnamed; however, their parameters are based on data available for specific insecticides. Thus, our third objective is to invite review from the risk assessment and modeling communities, in hopes that the best possible tools can be made available for use in avian pesticide risk assessments.

Models and integration
Salient features of each of the three simulation models are summarized in Table 1 and described in more detail below. Here we provide enough information to understand the basic workings of each model and supply additional references for those who wish to dig deeper.
TIM v.3.0 (hereafter referred to as TIM) [9] is an individual-based exposure and effects model that predicts avian mortality attributable to acute pesticide exposure resulting from a realistic time-dependent pesticide use scenario. Detailed information on TIM has been

PLOS ONE
provided elsewhere [4,9] and we include only a brief summary of important features here. TIM accounts for exposure through diet, drinking water, inhalation, and dermal contact following insecticide spray application to crops using a 1-hour timestep. Bird foraging on treated fields is stochastic, following a correlated movement on and off treated areas, and cumulative exposure considers avian frequency of usage of field and adjacent habitats that receive spray drift. Exposure in adjacent areas receiving spray drift decreases with distance from the edge of the field, where exposure is estimated using an approach adapted from the AgDRIFT model [23]. TIM runs for a single growing season, which may or may not overlap with avian breeding. Pesticide application methods that may be modeled in TIM include aerial, airblast, ground broadcast, ground banded, and ground in furrow. For all these application methods, exposure can be assessed on the treated field and edge habitat where spray drift is transported. Important assumptions of TIM include: 1. Toxicity to simulated birds is represented by surrogate test species (often the most sensitive of species tested); 2. Dose is a function of diet, inhalation, drinking water, and dermal uptake, as well as elimination; 3. Intake rates are allometrically scaled to bodyweight; 4. Birds move on and off treated fields with varying fidelity based on species characteristics; 5. Birds follow a bimodal feeding pattern where peak feeding occurs soon after sunrise and before sunset; 6. Pesticide concentrations on food items conform to a lognormal distribution; 7. Acute toxicity (LD50) is scaled using body weight and, where available, empirical pesticidespecific scaling factors [24].
The Markov chain nest productivity model (MCnest) is a mechanistic Monte Carlo simulation model that estimates declines in reproductive success for temperate zone birds during the breeding season [10,13]. The model uses a Markov transition matrix [25] to follow female birds as they complete a breeding cycle [26] of nest establishment, survival, nest failure, and re-nesting up to the end of a typical breeding season. The number of successful nests per breeding female in a single breeding season is tracked, along with the cause for each failed nest. Insecticide induced failures are modeled by the inclusion of daily exposure values from TIM. These values are compared to phase-specific NOAELs for different kinds of adverse effects that might be induced by insecticide exposure such as eggshell thinning, reduced egg viability, reduced hatching success, and increased abandonment [10,27]. Populations with identical parameter values are simulated with replication to provide estimates of variability around model predictions. MCnest output consists of estimates of the expected number of successful young per female in a population exposed to a given insecticide use scenario. Bennett and Etterson [28][29][30] provide detailed information on MCnest, including its use, technical background, and species library. Important assumptions underlying MCnest are: 1. Females necessarily attempt to renest if there is sufficient time remaining in the breeding season; 2. Demographic parameters, such as nest survival rates, clutch size, and waiting periods post failure and post fledging are fixed for all individuals; 3. Nest failures occur as a unit (i.e., not on a per-egg basis), whether natural or induced by insecticide exposure; 4. Exceedance of phase-specific toxicity endpoints results in complete nest failure; 5. NOAELs from the avian reproduction test are generally applicable across species.
Assumptions (3) and (4) have been criticized for causing MCnest to over-predict pesticide impacts on birds [31]. Although these assumptions do usually result in conservative predictions, they do not always do so. Their importance in light of the results presented herein is reviewed in the Discussion below.
TIM and MCnest were integrated within the Matlab [32] programming environment by having MCnest call TIM as a subroutine [4]. These integrated simulations first run TIM on adults to generate exposure and mortality results for each simulated breeding female; reproduction is then simulated using MCnest conditional on the mortality and exposure profiles generated from TIM. The integrated model is available online at: https://www.epa.gov/ chemical-research/markov-chain-nest-productivity-model.
HexSim is a spatially explicit, individual-based model (IBM) designed for simulating terrestrial wildlife population dynamics and interactions. For this work, we used HexSimPLE (Hex-Sim Populations Linked by Emigration), which is a customizable hybrid matrix/IBM constructed within HexSim and intended for use in developing parsimonious rapid response models (more detail provided below). HexSim and HexSimPLE are both available from www. hexsim.net. Our HexSimPLE gnatcatcher model consisted of a spatial array of 2-stage projection matrix population models [33] linked by spatially-explicit individual-based densitydependent dispersal (Fig 1). Model inputs included habitat and stress (insecticide) maps, a patch-map used to distribute individual population projection matrices across the landscape, stage-specific vital rates, and a few additional parameters used to specify carrying capacity and movement ability. Insecticide impacts on gnatcatcher vital rates were captured in the stress maps, which in-turn impacted individual matrix elements.
In the simulations described below, our patch-maps contained >24K 79 ha patches (and thus individual projection matrices), except along the range boundaries where patches were clipped and therefore somewhat smaller in size. We used two age-classes, juvenile (j) and adult (a), with the juvenile period lasting from the time of fledging up to the first breeding season at 1 year. Within patches, habitat quality (0 � q � 1) modified vital rates via a coefficient (c) that ranges from 0-1 according to the following function: c = 1−(1−q) α , where α is a user specified slope for the coefficient. For our simulations below, we used values of α = 3 for fecundity and α = 5 for survival, which were determined by tuning the model in the absence of simulated insecticides to produce range-wide population sizes at steady state that were roughly consistent with previous census data [34,40]. Population size estimates obtained this way represent a hypothetical upper limit reflecting assumptions that detection rates were perfect and the pesticides were absent.

Demography and environmental stochasticity
Demographic processes were modeled separately within each patch by using matrix multiplication to determine annual patch-specific population size, as a vector of age classes. Dispersal was limited to juveniles and modeled as a separate individual-based process (see Dispersal Model below). Therefore, aside from the spatially explicit dispersal step, patch-specific dynamics in ideal habitat (q = 1) can be described by the following equation:

PLOS ONE
In Eq (1) n j (t) is the number of juvenile birds (age-class 1) at time t and n a (t) is the number of adult birds (age-class 2) at time t in the habitat patch. Fecundity (f), the number of female offspring fledged per female is the same for both age classes, s a is the annual survival rate of adult birds and s j is the annual survival rate of juvenile birds. When habitat is not perfectly suitable, the demographic parameters s a , s j , and f are multiplied by patch-specific vital rate coefficients (c) that reflect expected reductions (from their maximum) due to lower habitat quality. Eq (1) assumes a post-breeding census, which proved convenient for the integration of TIM/MCnest output, but HexSimPLE can also employ a pre-breeding census (or multiple censuses) as well. When a census is taken just before the breeding pulse, the youngest individuals will have been alive for almost a full year, and thus have already been subjected to mortality. This is referred to as a pre-breeding census. When the census is taken just after the breeding pulse, the youngest individuals will be at the beginning of their lives, and not yet have experienced any mortality. This is referred to as a post-breeding census.

PLOS ONE
To introduce insecticide effects in HexSimPLE we modified Eq (1) to incorporate mortality induced by the stressor and informed by the output of TIM/MCnest as an additional survival term, s p = probability of mortality due to the patch-specific insecticide exposure: When there is no insecticide exposure, or no mortality due to insecticide exposure (s p = 1), Eq (2) is equivalent to Eq (1). The insecticide survival term (s p ) and fecundity term, f, are TIM/ MCnest outputs reflecting expected survival and reproductive success during the breeding season under given insecticide exposure conditions [4,11]. The survival terms (s j and s a ) are background rates of survival for juveniles and adults, respectively. By taking the product of the these and s p we assume that reduced survival due to insecticide exposure is an independent competing risk. As in Eq (1), the vital rates in Eq (2) would be modified by habitat quality coefficients (c) to reflect further reductions due to habitat quality. Finally, for all models, an additional 10% variability was introduced to the three basic vital rates (s a , s j , and f) to account for environmental stochasticity. A value of 10% was chosen because it is large enough to represent actual stochasticity, but not so large as to mask the model's response to our treatments. For simplicity of illustration, we assumed perfect correlation between the effects of stochasticity on survival and fecundity, but this assumption is not required by the model.

Life-history data
To obtain life history data for CAGN to parameterize TIM, MCnest, and HexSimPLE we consulted Akçakaya and Atwood [34] and Atwood and Bontrager [19] and followed citations therein to the original sources. We also consulted primary literature and USFWS documents [20]. Table 2 gives the vital rates we used for CAGN in perfect habitat (q = 1).

Habitat suitability model
To model habitat suitability, we took the estimated U.S. range of CAGN from USFWS [35] (since modified slightly) and used the habitat model of Akçakaya and Atwood [34], which requires estimates of four descriptive habitat variables, CSS = amount of coastal sage scrub, ELV = meters elevation above sea level, DGR = distance to grassland, and DTR = distance to nearest pixel (900 m 2 ) containing at least 10% tree cover. We summarized these habitat descriptors using National Gap Analysis Program Land Cover Data [36] from which we aggre-  Table 3) in ArcMap 10.3.1 [37] to produce a map of habitat quality (0 � q � 1) per pixel across the entire U.S. range of

PLOS ONE
CAGN. After exporting the vector habitat map as a raster image, we imported it into HexSima process that involved resampling the data at the scale of an individual hexagon (0.86 ha in this case). We then superimposed a regular grid of habitat patches across the U.S. range of CAGN, with each habitat patch containing up to 91 hexagons (~79 ha), unless clipped by the range boundary. Finally, per-patch habitat quality was set equal to the mean value taken over all patch hexagons. Per patch carrying capacity was set to K = 0.062 CAGN/hexagon � 14 ha/ territory, which is slightly larger than average territory size [19]) but well within the range of reported values (2-18 ha) [38].

Dispersal model
CAGN dispersal was limited to juveniles, who moved according to the following rules. If patch-specific population size (N) > patch-specific carrying capacity (K), then N − K juveniles must disperse (otherwise they remain in their natal patch). For dispersers, maximum dispersal distance is drawn from a uniform distribution ranging from 0.5 to 10 km. Individuals disperse according to a highly autocorrelated random walk, stopping when they reach a patch with N < K. If no such patch is encountered by the time they reach maximum dispersal distance, then they stop anyway. This dispersal mechanism produced simulated dispersal path lengths that were lognormally-distributed. In unusual cases where path length constraints result in patches with N > K after dispersal, then juvenile CAGN are culled until N = K. By design, the number of adults will never exceed K in any patch, so culling of adults is not required. The overall CAGN simulation was calibrated by modifying the dispersal distribution, the autocorrelation parameter, carrying capacity, and habitat exponent (α in Eq 1) to get to a plausible spatial distribution with equilibrium number of CAGN pairs < 5,000, as per Partner's-in-Flight US population estimate for CAGN [39].

Insecticides and model applications
We developed model applications by simulating the use of two acetylcholinesterase-inhibiting insecticides on wheat crops within the range of CAGN. Simulations for the two insecticides were independent of each other. The insecticides were assumed not to affect carrying capacity, nor did we model indirect effects (e.g., prey reduction). We obtained spatial coverages of crop distributions from USEPA [12] and assumed that a given insecticide was used on 100% of the wheat crop, though potentially applied on different dates. Wheat was chosen as it is found in areas adjacent to the habitat of CAGN; however, other crops (e.g., lettuce, tomatoes) that potentially receive insecticide exposures also occur near the habitat of CAGN. Data on the effects of the two insecticides on birds were taken from information submitted in support of

PLOS ONE
FIFRA pesticide registrations (Table 3). We excluded the names to emphasize that the focus of this work is on model development and procedure, not the actual effects of any specific products. In this case study, one pesticide was chosen because of individual-level risk concerns for mortality after acute exposures (hereafter survival stressor), and the other has risk concerns for reproductive effects following chronic exposures (hereafter reproductive stressor). Complete parameter sets used for TIM/MCnest are provided as S1 File. We used pesticide labeling information to determine the maximum application rates and assumed that growers could divide the maximum a.i. (lbs/acre) evenly into 2, 3, or 4 applications separated by the minimum interval between applications. For a given combination of insecticide and number of applications we created a hypothetical moving window of application dates by shifting the date of each application by 1 week and running TIM/MCnest to generate a curve of TIM/MCnest-predicted effects for application scenarios indexed by date of first application. We thus generated three such curves for each insecticide, corresponding to 2, 3, or 4 total applications. These curves were then used to generate a spatially explicit distribution of insecticide effects by randomly assigning an application scenario (combination of date of first application and number of applications) to each habitat patch and using linear interpolation from the effects curves (produced by TIM/MCnest, and described above) to estimate fecundity or survival associated with that habitat patch and application scenario. In the above procedure, the proportion of each habitat patch planted in wheat was used to modify the expected effects of the insecticide on patch-specific vital rates using a binomial mixing model: In Eq (3) v e = vital rate (fecundity or survival) predicted by TIM/MCnest under exposure conditions for the supplied patch-specific application date and rate, p = proportion of patch cultivated in wheat, and v u = value of a vital rate that in the absence of exposure in that particular habitat patch (i.e., taking into account habitat quality, q). For CAGN, unexposed survival s u = 1 and fecundity f u = 2.26 female fledglings/female in the best habitat (q = 1). We generated 100 realizations of the random distribution of insecticide usage calculated in this way and then resampled those 100 realizations to produce 10 trials for each insecticide. For each trial we randomly selected 50 realizations of the insecticide use pattern and used these as a fixed 50-year time-series. Each trial consisted of ten replicate simulations of 50 years of CAGN population growth, using the fixed time-series of insecticide use maps. To create a new trial, we randomly selected another 50 realizations (with replacement) from the 100 generated use maps, and repeated the ten replicate simulations.

Results
Across the range of CAGN, modeled habitat quality varied considerably (Fig 2). Except for some small inland patches, the largest and highest quality habitat patches were predicted to occur in coastal areas. Approximately 13% of the land within the CAGN range [35] was classified as potentially suitable (q > 0.5) and mean habitat quality among suitable hexagons (79 ha ea.) was 0.63 (+/-0.09 SD, n = 275,110). Habitat quality was negatively correlated with frequency of occurrence (Fig 2, lower inset), and high-quality habitat (q > 0.75) accounted for only about 1.5% of suitable habitat (q > 0.5).
When the vital rates of Table 2 are used in Eq (1), they give an estimate for density-independent finite rate of increase (λ) of 1.495, clearly demonstrating that, in ideal habitat, the species would be expected to exhibit strong positive population growth. With habitat exponents (α) of 3 and 5 for fecundity and survival, respectively, the transition between source and sink (λ = 1) would occur when habitat quality is approximately q � 0.32.
Predicted survival (s p ) from TIM/MCnest contrasted sharply between the two modeled insecticides. All survival stressor scenarios resulted in predicted pesticide mortality probability (1-s p ) of 1 for all exposed individual CAGN. Conversely, all reproductive stressor scenarios resulted in a predicted pesticide mortality probability of 0 for exposed individuals. In contrast, predicted seasonal productivity for exposed CAGN was eliminated by exposure to both insecticides (Fig 3) because all birds died under the survival stressor and all birds experienced zero fecundity under the reproductive stressor. Therefore, all TIM/MCnest predictions of effects of the reproductive stressor on fecundity are attributable to sublethal reproductive effects observed in the avian reproduction test results (Table 3) whereas all TIM/MCnest predictions  Table 2, S1 File), but any effects on fecundity are made irrelevant by the prediction of 100% mortality, given exposure, in the modeled scenarios.
Date of first application had a strong effect on predicted fecundity, whereas the number of applications affected predicted fecundity little (Fig 3A, reproductive stressor) or not at all ( Fig  3B, survival stressor). Seasonal patterns in model predictions for reductions in fecundity differed strongly between the two insecticides (Fig 3). For the reproductive stressor, TIM/MCnest predicted reductions in fecundity associated with applications as early as 1 January (Fig 3A), but not complete loss of reproductive output for the full season. This result is due to the presumed persistence of this chemical in the environment (lacking an estimate of the degradation half-life, the default value of 35 days was used). In contrast, TIM/MCnest predicted a complete loss of reproductive success from early applications of the survival stressor due to the lethal nature of exposure. For both insecticides, as the date of first application moved beyond the window of nest initiation for CAGN, which ends in roughly mid-July, predicted fecundity increased because modeled birds were able to successfully reproduce before experiencing any exposure to the insecticide.
HexSimPLE simulations of CAGN over 50 years showed differences between the two insecticides and control (Fig 4). Under control simulations (no insecticide) CAGN persisted at approximately 5,000 breeding pairs. Model runs predicted declines under both insecticides, with a new equilibrium being reached after approximately 15-20 years (Fig 4). Predicted

PLOS ONE
population size at year 50 under the survival stressor was lower than under the reproductive stressor due to the high mortality associated with the former insecticide. Individuals exposed to the survival stressor were predicted to be permanently removed from the population, whereas individuals exposed to the reproductive stressor, though they may have lost some or all expected fecundity for a given year, could nevertheless survive and attempt to breed the following season (or later the same season, if time permits). Fig 5 shows the predicted change in distribution of CAGN with exposure to both insecticides compared to control. The spatial distribution of the bird remains relatively constant, with some losses of habitat patches, particularly in the inland sites. However, most patches show some predicted contraction under each insecticide and even large contiguous patches show some reduction in cumulative abundance.

Discussion
Our results suggest that if both chemicals are used as described in this manuscript, then the survival stressor is predicted to have a greater impact on CAGN than the reproductive stressor. This result makes sense from a strictly mathematical point of view. The finite rate of increase (λ = dominant eigenvalue) for the matrix presented in Eq (1) is λ = s a + s j f (other symbols as previously defined). Parameterized with the data from Table 2 (i.e., for ideal habitat), the

PLOS ONE
sensitivity of λ to changes in s a is 1, whereas the sensitivity of λ to changes in f is 0.4314. Thus, the underlying demographic model suggests that the population growth rate should be more than twice as sensitive to changes in adult survival than to changes in fecundity, and this is reflected in the greater predicted population declines resulting from survival stressor exposure. However, this result should be interpreted with caution. The elasticities [33] of λ to changes in adult survival versus fecundity are reversed (0.35 vs. 0.65) suggesting that proportional changes in fecundity would have a greater effect on population growth than proportional changes in adult survival.
Our simulation results in the absence of insecticides differ strongly from those of Akçakaya & Atwood [34], who predicted decline and a high risk of extinction of CAGN. They further predicted that the studied Orange County population would fall below 60 individuals with 19% probability within 20 years and with 76% probability within 50 years. A more recent study by Winchell and Doherty [40] estimated the combined Orange County and San Diego County population on public lands to be 1,324 (976-1,673, 95% CI) [20,40]. We calibrated our model to equilibrate at around 5,000 pairs across all habitat in the absence of insecticide exposure.
Our results differ from those of Akçakaya & Atwood [34] in part because we used vital rates that are at the high end of their ranges. Thus, for example our survival rate estimates (s j = 0.4314 and s a = 0.52, Table 2) are the highest of those rates reported by Akçakaya & Atwood [34]. Similarly, our fecundity rates derived from MCnest without insecticide applications were 2.26 female offspring/female, which closely matches the highest fecundity value (2.3) reported by Akçakaya & Atwood [34]. Thus, our vital rates are optimistic, but they apply only in habitat receiving a habitat quality score of q = 1, which represents a miniscule proportion of the total PLOS ONE available habitat. The lowest rates reported by Akçakaya & Atwood [34] still occur in our model, but are associated with poor habitat quality scores (q < 0.5). Further, our vital rates ( Table 2) are for parameterizing the density-independent matrix that describes within-patch dynamics. Once patch-specific carrying capacity is reached, juvenile survival declines as juveniles disperse from their natal patch, some of whom will fail to locate a suitable breeding patch, and are culled. Thus, our reported juvenile survival rate (s j ) doesn't include the contribution of density dependence to juvenile mortality, which is an emergent property of habitat quality, carrying capacity, and dispersal in our HexSimPLE model.
Another reason we chose the higher rates reported by Akçakaya & Atwood [34] is that in so doing we minimize the potential for double-counting the effects of insecticides (both of which were already registered and in use when the reported studies were conducted). In any case, our intent here is not to improve upon existing vital rate estimates, previous predictions of population trajectory, or population size for CAGN. Rather, we hope to build upon previous work and demonstrate how the TIM/MCnest/HexSim model, together with available information on vital rates, habitat quality, and population size for CAGN can be used to better understand and interpret the risk of pesticide exposure and effects to this federally threatened species.
Many sources of variability contribute to our simulation results (Table 1), including stochastic survival (TIM) and reproduction (MCnest), variable pesticide exposure and tolerance (TIM) among individuals, habitat heterogeneity (HexSim), probabilistic dispersal (HexSim), and environmental stochasticity (HexSim). With this demonstration we have not tried to quantify the contributions of these sources of uncertainty, though that remains an interesting avenue for further research.
Recently USEPA conducted a biological evaluation of several organophosphates under Section 7 of the ESA [12][13][14]. Thus, it is interesting to contrast those results with ours. Using TIM/MCnest, USEPA concluded that there were risks of mortality and/or reductions in reproductive success across all organophosphates examined in the study. Our results also show reductions in survival with the survival stressor modeled herein (and therefore reproductive success), but continued persistence of California Gnatcatchers. This helps to resolve an apparent paradox pointed out by Moore [31], that few populations could withstand such mortality, yet populations persist, even though these products have been used for many years. An important distinction is that the USEPA [12][13][14] risk characterizations consider exposed individuals only. Thus, when similar scenarios are run in a spatially explicit environment that includes habitat quality, a mix of exposed and unexposed individuals, and refuges from exposure, the populations are predicted to persist for at least 50 years in an apparently stable trajectory (Fig 4).
An interesting result of our simulations is that the number of applications had little or no effect on survival and/or fecundity of CAGN. This result is not general and should be interpreted with great caution. For the survival stressor, repeat applications had no effect because the model predicted 100% mortality upon the first application, due in large part to high predicted inhalation toxicity. For the reproductive stressor, repeat applications had little effect because of the 35d (default) half-life used for residues on dietary items. This had the effect of maintaining residue concentrations at toxic levels throughout the breeding season.
Perhaps the most useful aspect of the modelling framework we evaluate here will be the ability to evaluate the influence of factors surrounding how, when, and where an insecticide is applied on impacts to species of birds. For example, simulated application dates can be modeled in TIM/MCnest to determine when the risk of pesticide use is minimal. This may be most useful for insecticides that do not have high mortality risk but may cause sublethal reproductive effects. Similarly, the spatially explicit nature of HexSim can be used to explore the overlap of pesticide use, animal demography, species distribution, and critical habitat in a way that would protect essential resources, while perhaps allowing insecticide usage in marginal habitats. In addition, this framework allows for consideration of all relevant use sites in assessing exposures and risk to the population. In this example, wheat was simulated; however, exposures from other use sites that are adjacent to the resources used by CAGN could also be integrated into HexSimPLE, along with mortality and fecundity declines associated with those other uses. Although we did not do so in this case study, available crop specific usage data on the percent of treated area and actual application rates (e.g., California's Pesticide Use Reporting data, https://www.cdpr.ca.gov/docs/pur/purmain.htm) could also be incorporated in order to evaluate the impact of actual application practices on risk.
While we believe our model has promise and offers substantial advances to current methods for avian pesticide risk assessment, we must also acknowledge that it has some important limitations. In particular, the assumptions listed in Methods should be carefully noted (we briefly touch on them here-they have been reviewed in more depth elsewhere [4,28]). With some effort, some of those assumptions could be relaxed, though perhaps not all. One MCnest assumption was noted recently as an important limitation [31], which is that complete nest failure occurs with exceedance of phase-specific toxicity endpoints. As noted by Etterson et al. [4], relaxation of this assumption would be difficult, requiring quantitative description of the relationship between specific endpoints from toxicity tests and avian fledging success with exposure to insecticides. However, the toxicity tests, especially the avian reproduction test, simply weren't designed to provide quantitative dose-response curves and would have to be significantly modified [28]. See references [4 and 10] for more in-depth discussion of the MCnest assumptions. The assumption of complete nest failure with endpoint exceedance is not universally conservative. A female who loses a complete nest attempt early to pesticide exposure, but subsequently renests successfully may experience greater reproductive success than a female who suffers low fledging success (e.g., raises one fledgling) from a successful nest attempt if doing so leaves insufficient time to renest.
Other important limitations to the model application described here include lack of connectivity between the Mexican and US CAGN populations. Our model includes no dispersal into or out of Mexico. We do not believe this to be a significant source of error in our simulations presented, or any plausible modifications thereof, primarily due to the relatively limited dispersal of CAGN [20,34]. Our model also requires a habitat description. In the presented example, we had an existing high-quality habitat model available from the scientific literature upon which we could draw. For some listed species this will not be the case and a more qualitative model would have to be substituted. We do not know how such a model would perform in our linked modelling system, but we believe it is an important area for further research. Finally, an important limitation is that we have included only insecticide exposure as a single stressor (or arguably two stressors if poor-habitat quality is also considered), whereas many listed species are facing multiple threats, the combination of which result in the realized population trajectory of the taxon. For example, brown headed cowbird nest parasitism is an important stressor for CAGN that was not included in the example simulations provided herein, except as a component of background rates of nest failure. Similarly, insecticides may also affect birds indirectly by reducing the availability of invertebrate prey [41][42][43] [44]. Our modeling system has the capacity for modeling nest parasitism and indirect effects as well as an arbitrary number of additional stressors, provided their spatial distribution can be described and their effects are known. In some cases, it may be argued (as we did immediately above concerning cowbird parasitism) that the effects of other stressors may already be incorporated into the background vital rates used for modeling. This will not be universally true and should be evaluated on a case by case basis. Another important and difficult example of multiple stressors that would require considerable further research to include would be the effects of exposure to multiple insecticides used simultaneously, or nearly so, at the landscape scale. Thus, while any number of complexities can be considered for inclusion, we believe these are best handled on a caseby-case basis, as most of these processes will require substantial effort to evaluate and parameterize.

Conclusions
We have presented a spatially explicit population-level model for assessing risks of pesticides to passerine and near-passerine birds that use agricultural areas and adjacent habitats. We provided an example using the threatened California Gnatcatcher and showed how the model output could be used to better understand how the use of an insecticide across time and space may influence risk at the population level. However, this demonstration should be viewed as a beginning. A more thorough risk assessment using these models would include more stakeholders and a larger exploration of parameter space. Of course, a prerequisite to such an assessment involves identifying a modeling approach that is up to the challenge, but also scientifically sound and defensible. Here, we have described a new research tool that we believe has these very attributes. Our integrated models represent a major step towards the recent vision laid out by the National Academies [3] for a probabilistic approach to pesticide risk assessment for federally listed species.