A Simplified Population-Level Landscape Model Identifying Ecological Risk Drivers of Pesticide Applications, Part One: Case Study for Large Herbivorous Mammals

Environmental risk assessment is a key process for the authorization of pesticides, and is subjected to continuous challenges and updates. Current approaches are based on standard scenarios and independent substance-crop assessments. This arrangement does not address the complexity of agricultural ecosystems with mammals feeding on different crops. This work presents a simplified model for regulatory use addressing landscape variability, co-exposure to several pesticides, and predicting the effect on population abundance. The focus is on terrestrial vertebrates and the aim is the identification of the key risk drivers impacting on mid-term population dynamics. The model is parameterized for EU assessments according to the European Food Safety Authority (EFSA) Guidance Document, but can be adapted to other regulatory schemes. The conceptual approach includes two modules: (a) the species population dynamics, and (b) the population impact of pesticide exposure. Population dynamics is modelled through daily survival and seasonal reproductions rates; which are modified in case of pesticide exposure. All variables, parameters, and functions can be modified. The model has been calibrated with ecological data for wild rabbits and brown hares and tested for two herbicides, glyphosate and bromoxynil, using validated toxicity data extracted from EFSA assessments. Results demonstrate that the information available for a regulatory assessment, according to current EU information requirements, is sufficient for predicting the impact and possible consequences at population dynamic levels. The model confirms that agroecological parameters play a key role when assessing the effect of pesticide exposure on population abundance. The integration of laboratory toxicity studies with this simplified landscape model allows for the identification of conditions leading to population vulnerability or resilience. An Annex includes a detailed assessment of the model characteristics according to the EFSA scheme on Good Modelling Practice.


Introduction
Environmental risk assessment is a key process for the authorization of pesticides, subjected to continuous challenges and updates [1]. In the regulatory context, different agencies have developed models and scenarios for pre-marketing and re-evaluation assessments. The proposed pesticide uses patterns, named "Good Agricultural Practices (GAPs)", which are tested according to these models and scenarios, checking whether or not the estimated level of risk is in line with the protection goals defined by risk managers. Exceeding the accepted level of risk means non-approval or mandatory risk mitigation options, such as maximum application dose, limited number of applications within the season, untreated buffer zones, etc. In the EU system, the protection goals are generic, although the European Food Safety Authority (EFSA) has suggested the development of specific protection goals [2], based on the concept of ecosystem services. For terrestrial mammals, current protection goals include acute lethality and markers for population level pesticide risk assessment [11,16]; additional requirements were considered for ensuring the model capability for supporting environmental management decisions [17,18].

Model Conceptualization
The conceptual model includes two modules: (a) a simplified population dynamics model, and (b) a landscape-based model for integrating the impact of the pesticide exposure on the population model parameters.
The population dynamics is modelled at an individual level for groups of individuals, "nests", located in a defined feeding area. Each "nest" is defined by its location, the initial number of individuals distributed in up to four age groups, and the seasonal reproductive period within the year. Each age group has an associated background mortality rate. Females from the reproductive age groups have an associated background reproduction rate. All parameters can be selected and modified by the user, and adjusted to geoecological conditions. The modifications allow modelling pesticide effects for populations in expansion, recession, or steady-state conditions. At daily intervals, the mortality and reproduction rates are used for setting the number of deaths in each age group and the number of new-borns. Deaths are randomly allocated to the individuals in the group, and the sex of the new-borns is also randomly assigned. Then, the age of all individuals is increased by one day (or other selected time period) and those reaching the age-threshold moved to the next age group. The population dynamics can be represented as the evolution of the total number of individuals in each "nest", the evolution of the distribution of individuals per age group, and the combined distribution of several nests in the area.
The pesticides landscape-based model integrates the impact of the pesticide exposure on the population model parameters for each age-group, as a feeding area around the nest is defined. Each nest is placed in a defined location within the landscape scenario. Each feeding area is then connected to the land use to estimate the percentages for the fraction of the diet obtained for each zone within the feeding area. Following the application of a pesticide, the pesticide concentration in the treated crop and other food items is estimated daily, according to the expected environmental fate. Pesticide concentrations are estimated for each commodity and weighted according to the percentage in the diet. The daily or time-weighted exposures are used as input values for the dose-response curves to estimate the expected impact of the pesticide exposure on the mortality and reproduction rates. These values are selected from the available toxicity data using expert judgement.

Model Parameterization
The population dynamic parameters for each species are: distribution of age groups, background mortality rates, feeding range, and background reproduction rates, and reproductive period. The information for the selected species was retrieved from targeted searchers in PubMed and WebOfScience. The model allows the user to define all these parameters adapting the results to local habitats and regional conditions. The expected variability of these ecological parameters under different habitats and regions was considered for comparative assessments and sensitivity analysis. The pesticide exposure model implements, at a landscape level, the scenarios of the EFSA guidance for assessing the risk of pesticides to birds and mammals [3]. This guidance is currently under review. The model is based on the last update (guidance published on 26 July 2010) which is mandatory in the EU for regulatory risk assessment of pesticides and will be updated in the future to account for possible revisions. The guidance provides the scenarios, equations, and default values for estimating the level of pesticide residues in the different food items for a large set of crops, and consumption amounts for a set of representative birds and mammals. As the model utterly implements the EFSA exposure assessment for herbivorous mammals, the exposure equations are not duplicated here. The reader is referred to the EFSA guidance [3].
The pesticide effects are modelled through three complementary parameters, that must be selected by expert judgement from the available toxicokinetic and toxicodynamic information: • Change in acute mortality rate • Change in chronic mortality rate • Change in reproduction rate The extrapolation of the experimental mammalian toxicity data observed in the lab to the expected consequences in the field required several assumptions. In repeated-dose studies, animals are dosed every day at the same level, which is not the case in the field. The change in exposure can be modelled using the time-weighted exposure concentrations [3]; some specific recommendations were provided in the opinion of the EFSA Panel [19] that served as a background for the Guidance. The time required for observing relevant symptoms is usually not considered in regulatory assessments but is essential for modelling population effects. This requires access to the daily or weekly observations for each study; using the time between initiation of the dosing and the observation of the symptoms for setting the length of the period to be used in the time-weighting exposure estimation. Regarding the dose-response curve, the ideal situation is the use of the benchmark dose approach [20], but this is rarely available, and in most cases, only the No Observed Adverse Effect Level (NOAEL) and the Lowest Observed Adverse Effect Level (LOAEL) are available. These parameters are selected from statistical analysis, regardless of the magnitude of the effect, the estimation of the actual magnitude at the time window used for the assessment is problematic in most cases, thus a pragmatic assumption was applied when verifiable quantitative information could not be retrieved from the data, or when the assessment was the result of combining several studies. The default level of effect for the LOAEL was assumed as 25%; the default curve was a linear dose-response, provided that the level of effect for the NOAEL was within the range 5-10%. Linearity was considered appropriate as the relationship models the average effect for the average exposure of each population subgroup per nest, which includes individuals with different exposure levels.
The impact of pesticides on the mortality rates was estimated by combining survival rates according to the following equation: where Subgroup SR is 1 minus the background mortality rate for the population group and Px SR is 1 minus the mortality rate associated to exposure to pesticide x (for x = a to i). Equation (1) is used for estimating the accumulated pesticide impact on the acute and on the chronic mortality. Acute mortalities rates are estimated on daily basis according to the actual exposure value and assuming that there is no background acute mortality. Chronic mortality rates are estimated monthly, considering the background monthly mortality for the population group and time window averaged (TWA) exposure levels. The time for estimating TWA is pesticide specific according to their toxicity profile, the maximum TWA exposure value from the previous month is used for the calculations in Equation (1). For reproductive effects, both reproduction and developmental studies are considered; combining the relevant observations for maternal toxicity, effects on fertility, and foetal and neonate mortality. In addition to the dose and magnitude of the effect, the time between application and observation of the effects is critical to establish the time period to which the pesticide effect should be applied. For the selected time period, the impact on reproduction was estimated as follows: Final monthly reproduction rate = (Background rate) ( where the reproduction rate represents the average number of offspring per reproductive female and month during the reproductive window, and Px RE is 1 minus the effect on reproduction associated with pesticide x. The average pregnancy time is considered for selecting the time between exposure and actual observation of the effects on reproduction, i.e., for a pregnancy time of 30 days the effects on the reproductive toxicity rate are observed in the following month. The maximum exposure value from the previous month is used for estimating the effects on mortality and from the two previous months for estimating the effects on reproduction.

Computer Implementation
All calculations were implemented in a computer model developed using the programming language Python version 3. The model was specifically developed for this purpose.
In addition, an executable version with a user-friendly interface was developed. This version offers large flexibility to the user and can be made available under request to potential users for non-commercial purposes. The implemented features include: • replicability assessment through a number of iterations (within nest replications) selected by the user and expression of results as mean with 95th confidence intervals or maximum/minimum, • landscape with field areas defined by the user, • inner and outer adjustable bands for each field, and • selection of crops, crop rotation and pesticide applications at time dates defined by the user, • selection of the lagomorph species and location and all ecological model parameters of each nest defined by the user.
The software implements several graphic and tabular outputs for presenting different types of results according to the user's needs.

Case Study with Lagomorphs
The European rabbit (Oryctolagus cuniculus) and the brown hare (Lepus europaeus) are the representative species of large herbivorous mammals included in the EFSA guidance [3] and have been implemented for this case study.
The population dynamic parameters were selected from a review of the ecology of these two lagomorph species. Default values and ranges were selected based on information retrieved from targeted searchers in PubMed and WebOfScience for the following parameters: distribution of age groups, background mortality rates, feeding range, and background reproduction rates, and reproductive period.
For this case study, we have selected ecological parameters according to previous reviews representing Mediterranean conditions in the EU. The rabbit reproductive season (six months per year in winter and spring) and rate (average of four litter per month for a reproductive female) are based on the reviews by Tablado and co-workers [21,22] and mortality rates were extracted from reviews from the same authors [21,23]. The brown hare reproductive season in central Europe (seven to eight months per year, covering part of the summer) and rate (average 1.75 litter per month for a reproductive adult female, 0.6 litter per month for a reproductive young female) and mortality rates were extracted from several reviews [24][25][26]. Under Mediterranean conditions, the hare reproductive season may be extended to the full year, but accompanied by a reduction in the litter number [27], resulting in an estimated reproduction rate of around 1 litter per month for a reproductive adult female. Table 1 summarizes the selected default values; different values were used in some model estimations for assessing the influence of population patterns, details are provided in each figure caption.
The third element is the toxicodynamic assessment of each pesticide, selecting the relevant exposure time windows and the dose-response relationships for acute lethality, chronic toxicity, and reproduction. Two herbicides with potential risk for mammals according to the EFSA conclusions, glyphosate [28] and bromoxynil [29], have been used as model chemicals in this study. All studies available on mammals and reported in the supporting information for the EFSA Conclusions are available on the EFSA website and were searched for the toxicodynamic assessment.

Model Calibration and Validation for Context of Use Qualification
The model simplifies the ecological complexity to estimate the impact of pesticides. It is based on regulatory scenarios and exposure estimations, and consequently is not suitable for a direct field validation. Instead, the model was calibrated and verified for regulatory use qualification by comparing the outcomes of the specific model subroutines (population dynamics, exposure assessment, time-weight estimations, and impact on survival and reproduction rates) with manual calculations using the equations implemented in the model. Following a set of iterative calibrations, a final verification confirmed the consistency between the model and the manual calculations. This was further confirmed with the case study with lagomorphs. The largest variability is observed during the peak period (end of the reproductive season) for the conditions representing a stable population (C and D). Figure 2 confirms that the proposed default parameters for reproduction and mortality rates represent conditions close to stable populations with large seasonal variability, particularly for initial densities of 140 and 240 rabbits per nest.

Model Replicability and Flexibility
The population peaks at the end of the reproductive seasons, almost triplicating the initial number of individuals and goes back to background levels during the nonreproductive seasons. Figure 3 confirms that under these conditions, the populations are stable for a very long period, while the variability associated with the aleatory allocations increases with time for the first 8-9 years and then stabilises. The model allows the estimation of not stable populations and all kinds of combinations, as presented in Figure 4. The model reproduces rabbits under Mediterranean conditions (reproduction season from December to June) for starting populations of 100 (A and B) and 140 (C and D) individuals per nest. The observed variability is linked to aleatory allocations of deaths within the age group and sex of new-borns and reflects the expected natural variability. The largest variability is observed during the peak period (end of the reproductive season) for the conditions representing a stable population (C and D). Figure 2 confirms that the proposed default parameters for reproduction and mortality rates represent conditions close to stable populations with large seasonal variability, particularly for initial densities of 140 and 240 rabbits per nest.   The model reproduces rabbits under Mediterranean conditions (reproduction season from December to June) for starting populations of 100 (A and B) and 140 (C and D) individuals per nest. The observed variability is linked to aleatory allocations of deaths within the age group and sex of new-borns and reflects the expected natural variability. The largest variability is observed during the peak period (end of the reproductive season) for the conditions representing a stable population (C and D). Figure 2 confirms that the proposed default parameters for reproduction and mortality rates represent conditions close to stable populations with large seasonal variability, particularly for initial densities of 140 and 240 rabbits per nest.  The population peaks at the end of the reproductive seasons, almost triplicating the initial number of individuals and goes back to background levels during the non-reproductive seasons. Figure 3 confirms that under these conditions, the populations are stable for a very long period, while the variability associated with the aleatory allocations increases with time for the first 8-9 years and then stabilises. The model allows the estimation of not stable populations and all kinds of combinations, as presented in Figure 4.   The population peaks at the end of the reproductive seasons, almost triplicating the initial number of individuals and goes back to background levels during the non-reproductive seasons. Figure 3 confirms that under these conditions, the populations are stable for a very long period, while the variability associated with the aleatory allocations increases with time for the first 8-9 years and then stabilises. The model allows the estimation of not stable populations and all kinds of combinations, as presented in Figure 4.

Case Study with Lagomorphs
The toxicological information published by EFSA was used for the data extraction. This information includes the "Conclusion on Pesticides" published in the EFSA Journal and the final Rapporteur Member State assessment with comprehensive summaries of each study published as a background document in the "Registry of Questions" supporting the Conclusion.

Glyphosate
The acute lethality estimation was based on a single dose LD 50 higher than 2000 mg/kg b.w. and the associated LD 0 higher than 200 mg/kg b.w. As these values are expressed as "higher than" values, no additional interspecies correction was applied. A linear relationship (linear regression for pairs (200,0) and (2000,0.5) expressing mortality as rates) provides the following equation: Acute P glyphosate SR = 0.000278ETE t − 0.055556; values below 0 are corrected to 0, (3) where ETE t represents the average estimated theoretical exposure in mg/kg b.w. of day t for each subpopulation group and nest.
The effect on chronic mortality rate considers LOAELs from subacute oral studies (Section B.6.3.1), long-term toxicity and carcinogenicity (Section B.6.5), and relevant effects from reproductive studies (Section B.6.6). A conservative approach was used in line with regulatory assessments. The lowest LOAEL was selected, including not only mortality but also morbidity under the assumption that in the field these effects will affect the survival capacity of the individual. The selected value was a LOAEL of 175 mg/kg b.w. day for maternal toxicity in developmental studies on rabbits. The selected NOAEL is 50 mg/kg b.w. A default effect value of 25% was allocated to the LOAEL. Although the study is on rabbits, the default interspecies variability factor of 5 for chronic exposures was applied for accounting for differences between the domestic rabbit subspecies Oryctolagus cuniculus domesticus, wild rabbits, and hares. Considering the time between initiation of exposure and observation of the effects, a time window of 12 days was fixed for the exposure estimations. A linear relationship (linear regression for pairs (0,0) and (35,0.25) expressing the effects as rates) met the default range effect value (NOAEL ranging between 5 and 10% of effect) and provided the following equation: where ETEtwa (t−12,t) represents the time weight average estimated theoretical exposure in mg/kg b.w. for the 12 previous days to day t for each subpopulation group and nest. Considering that the maternal toxicity was more sensitive than developmental toxicity, a similar equation covers the reproductive effects: where ETEtwa (t−12,t) represents the time weight average estimated theoretical exposure in mg/kg b.w. for the 12 previous days to day t for each subpopulation group and nest. Figure 5 present an example of the evolution of the daily and twa ETEs, the use of the maximum monthly values for the estimations, and the resulting effects on mortality and reproduction rates for two applications of glyphosate at 4.0 kg/ha in cereals during the rabbit reproduction period.

Bromoxynil
The acute lethality estimation was based on a rat single dose LD50 of 130 mg/kg b.w. and the associated LD0 10 mg/kg b.w.; with an interspecies uncertainty factor of 5 applied to the LD50 but not to the LD0 as the lack of mortality is confirmed by other studies in rabbits. A linear relationship (linear regression for pairs (10,0) and (26,0.5) expressing mortality as rates) provides the following equation: Acute Pbromoxynil SR = 0.0313ETEt − 0.3125; values below 0 are corrected to 0, (6) where ETEt represents the average estimated theoretical exposure in mg/kg b.w. of day t for each subpopulation group and nest.
The effect on chronic mortality rate considers LOAELs from subacute oral studies (Section B.6.3.1), long-term toxicity and carcinogenicity (Section B.6.5), and relevant effects from reproductive studies (Section B.6.6). A conservative approach was used in line with regulatory assessments. The lowest LOAEL was selected, including not only mortality but also morbidity under the assumption that in the field these effects will affect the survival capacity of the individual. The selected value was a LOAEL of 17.1 mg/kg b.w. day from a subacute toxicity study on mice, and the default interspecies variability factor of 5 for chronic exposures. Rats and dogs have lower NOAELs for some sublethal effects but also higher NOAELS in 90 days and 1 year studies. The selected value would be also in line with the benchmark approach proposed in the assessment. Considering the time between initiation of exposure and observation of the effects, a time window of 5 days was fixed for the exposure estimations. A linear relationship (linear regression for pairs (0,0) and (3.42,0.25) expressing the effects as rates) met the default range effect value for the NOAEL and provided the following equation: where ETEtwa(t−5,t) represents the time weight average estimated theoretical exposure in mg/kg b.w. for the 5 previous days to day t for each subpopulation group and nest. Regarding reproduction effects, the most sensitive and ecologically relevant observed effect corresponded to a LOAEL of 12.5 mg/kg b.w. day with increased post-implantation loss and malformation, with the default interspecies variability factor of 5 for chronic exposures. Considering the time between initiation of exposure and observation of the effects, a time window of 5 days was fixed for the exposure estimations. A linear

Bromoxynil
The acute lethality estimation was based on a rat single dose LD 50 of 130 mg/kg b.w. and the associated LD 0 10 mg/kg b.w.; with an interspecies uncertainty factor of 5 applied to the LD 50 but not to the LD 0 as the lack of mortality is confirmed by other studies in rabbits. A linear relationship (linear regression for pairs (10,0) and (26,0.5) expressing mortality as rates) provides the following equation: Acute P bromoxynil SR = 0.0313ETE t − 0.3125; values below 0 are corrected to 0, (6) where ETE t represents the average estimated theoretical exposure in mg/kg b.w. of day t for each subpopulation group and nest.
The effect on chronic mortality rate considers LOAELs from subacute oral studies (Section B.6.3.1), long-term toxicity and carcinogenicity (Section B.6.5), and relevant effects from reproductive studies (Section B.6.6). A conservative approach was used in line with regulatory assessments. The lowest LOAEL was selected, including not only mortality but also morbidity under the assumption that in the field these effects will affect the survival capacity of the individual. The selected value was a LOAEL of 17.1 mg/kg b.w. day from a subacute toxicity study on mice, and the default interspecies variability factor of 5 for chronic exposures. Rats and dogs have lower NOAELs for some sublethal effects but also higher NOAELS in 90 days and 1 year studies. The selected value would be also in line with the benchmark approach proposed in the assessment. Considering the time between initiation of exposure and observation of the effects, a time window of 5 days was fixed for the exposure estimations. A linear relationship (linear regression for pairs (0,0) and (3.42,0.25) expressing the effects as rates) met the default range effect value for the NOAEL and provided the following equation: where ETEtwa (t−5,t) represents the time weight average estimated theoretical exposure in mg/kg b.w. for the 5 previous days to day t for each subpopulation group and nest. Regarding reproduction effects, the most sensitive and ecologically relevant observed effect corresponded to a LOAEL of 12.5 mg/kg b.w. day with increased post-implantation loss and malformation, with the default interspecies variability factor of 5 for chronic expo-sures. Considering the time between initiation of exposure and observation of the effects, a time window of 5 days was fixed for the exposure estimations. A linear relationship (linear regression for pairs (0,0) and (2.5,0.25) expressing the effects as rates) met the default range effect value for the NOAEL and provided the following equation: where ETEtwa (t−5,t) represents the time weight average estimated theoretical exposure in mg/kg b.w. for the 5 previous days to day t for each subpopulation group and nest. Figure 6 present an example of the evolution of the daily and twa ETEs, the use of the maximum monthly values for the estimations, and the resulting effects on mortality and reproduction rates for one application of bromoxynil at 1.0 kg/ha in cereals during the rabbit reproduction period. Figure 6 present an example of the evolution of the daily and twa ETEs, the use of the maximum monthly values for the estimations, and the resulting effects on mortality and reproduction rates for one application of bromoxynil at 1.0 kg/ha in cereals during the rabbit reproduction period.
The model allows predictions related to the effect of pesticide applications at a population level. For similar population settings, these effects depend on the toxicity of the pesticide and the application rate, but also on the time of application within the reproductive period. Figure 7 offers an example for the application of bromoxynil at different doses and times for conditions representing a stable population. In line with the information extracted from the EFSA assessment, the exposure to bromoxynil from the treated crop and in-field grass is expected to be a short-term phenomenon, due to the rapid dissipation of the pesticide residues. If the acute lethality threshold is exceeded, some impact may be observed immediately, while most of the effects on the population dynamics will be delayed for several weeks accounting for effects on the mortality rate associated with chronic morbidity and on the reproduction of exposed females. The effects of application rates as low as 0.05 kg/ha of bromoxynil are visualised for applications occurring at the beginning of the reproductive season, but are least evident or even disappear for later applications. In the case of very late applications, only minor effects are observed within the season even at 0.2 kg/ha, but the consequences are very clear in the next season even without additional pesticide applications. Figure 6. Relationship between daily exposure (A1), 5 days, time-weighted average exposure (A2), the maximum monthly twaETA used for assessing mortality (B1) and reproductive (C1) effects, and the associated effects on the mortality (B2) and reproduction (C2) rates, following one application of bromoxynil on cereals at 1.0 kg/ha. Figure 6. Relationship between daily exposure (A1), 5 days, time-weighted average exposure (A2), the maximum monthly twaETA used for assessing mortality (B1) and reproductive (C1) effects, and the associated effects on the mortality (B2) and reproduction (C2) rates, following one application of bromoxynil on cereals at 1.0 kg/ha. The model allows predictions related to the effect of pesticide applications at a population level. For similar population settings, these effects depend on the toxicity of the pesticide and the application rate, but also on the time of application within the reproductive period. Figure 7 offers an example for the application of bromoxynil at different doses and times for conditions representing a stable population. In line with the information extracted from the EFSA assessment, the exposure to bromoxynil from the treated crop and in-field grass is expected to be a short-term phenomenon, due to the rapid dissipation of the pesticide residues. If the acute lethality threshold is exceeded, some impact may be observed immediately, while most of the effects on the population dynamics will be delayed for several weeks accounting for effects on the mortality rate associated with chronic morbidity and on the reproduction of exposed females. The effects of application rates as low as 0.05 kg/ha of bromoxynil are visualised for applications occurring at the beginning of the reproductive season, but are least evident or even disappear for later applications. In the case of very late applications, only minor effects are observed within the season even at 0.2 kg/ha, but the consequences are very clear in the next season even without additional pesticide applications. The decision on the timings selected for connecting the observed experimental results with the expected impacts on rates in the field is also relevant for assessing multiple treatments. The proposal selected for glyphosate and bromoxynil in these estimations is to link max twaETE observed in a month with the rates of the following month. This approach allows for the generic estimations intended for a simplified model, as the actual breeding season will change, but should be considered by the risk assessors in the assessment of multiple applications as exemplified in Figure 8. As fixed month intervals are selected for the monthly rate adjustment if the pesticide dissipation half-life is much shorter than the interval during treatments, the selection of the actual application dates will change the maximum monthly twaETAs and consequently the impact on the rates and population abundance. These differences should be included in the sensitivity analysis. The decision on the timings selected for connecting the observed experimental results with the expected impacts on rates in the field is also relevant for assessing multiple treatments. The proposal selected for glyphosate and bromoxynil in these estimations is to link max twaETE observed in a month with the rates of the following month. This approach allows for the generic estimations intended for a simplified model, as the actual breeding season will change, but should be considered by the risk assessors in the assessment of multiple applications as exemplified in Figure 8. As fixed month intervals are selected for the monthly rate adjustment if the pesticide dissipation half-life is much shorter than the interval during treatments, the selection of the actual application dates will change the maximum monthly twaETAs and consequently the impact on the rates and population abundance. These differences should be included in the sensitivity analysis.
The EFSA guidance considers the brown hare as the focal species for grassland and vineyards. Figure S1 in the Supplementary Material simulates the combined effect of two treatments (end of winter and spring/summer) on a stable brown hare population in central Europe (see Table 1 for details). The results confirm the need for assessing the combined effects of different applications through the season in population-based risk assessments. allows for the generic estimations intended for a simplified model, as the actual breeding season will change, but should be considered by the risk assessors in the assessment of multiple applications as exemplified in Figure 8. As fixed month intervals are selected for the monthly rate adjustment if the pesticide dissipation half-life is much shorter than the interval during treatments, the selection of the actual application dates will change the maximum monthly twaETAs and consequently the impact on the rates and population abundance. These differences should be included in the sensitivity analysis.

Discussion
The examples provided in Figures 5 and 6 indicate the expected exposure estimation (based on EFSA guidance equations [3]) and how this exposure is translated into effects on survival and reproduction rates according to the equations implemented. These effects are implemented as impacts on the population dynamics in Figures 7 and 8. To facilitate the verification of the results, these examples assume maximum exposure from the treated field according to work-cased conditions of the EFSA guidance. The landscape component of the model allows exposure from different fields, crops, and crop rotations according to the user's design.
Topping et al. [30] have addressed the shortcomings of current pesticide environmental risk assessment strategies and suggested an integrated system approach. The proposal implements opinions from the EFSA Panels and includes better use of available ecotoxicity data, the incorporation of landscape tools, and the need for integrating the risks from different pesticides, stressing the limitation of the single product and single crop approach [30].
The ecotoxicity assessment of wild mammals is rarely based on dedicated ecotoxicity studies. Generally, the assessment is based on the re-evaluation of the safety human health studies. This has the benefit of a relatively large number of studies with different designs, as required to cover the human hazard assessment, but with the limitations of studies and endpoints targeted to humans. As confirmed by the glyphosate and bromoxynil examples, instead of focusing on selecting the lowest relevant NOEL as the Reference Point, the information can be evaluated from a different perspective and used to establish dosetime-response relationships. The model addresses the impact of pesticides on mammalian populations [30] by combining all morbidity effects into the combined effect on the monthly mortality rate. The interpretation of the time patterns from the toxicological studies is a critical element for assessing the effects of pesticides. The information available is based on standard toxicity testing for assessing repeated dose effects with continuous dosing. This situation requires the reanalysis of the daily data from each test to identify the time required for observing the effects. The comparison of the NOAELs and LOAELs from studies with different duration may offer relevant information [31]. For glyphosate and bromoxynil, the lowest LOAELs were observed already at short exposure durations, suggesting the use of exposure time windows of 12 and 5 days for glyphosate and bromoxynil, respectively. Achieving the maximum toxicity level within a relatively short time period is a frequent feature in chronic toxicity testing, observed for about 50% of pesticides [32] as well as chemicals in general [33], being more frequent for chemicals with high toxicity [31]. Based on this information, substance and endpoint specific values are required regarding the number of days to be used for estimating the time-weighted average exposure.
While standard risk assessment approaches focus on setting thresholds, our proposal for population modelling is to use the full dose-response curve and to integrate the time component, selecting for each pesticide the relevant duration of effects associated with chronic exposure according to the information extracted from the different studies. This is in line with the use in ecotoxicology of novel dose-response tools such as the benchmark dose approach [31]; it also considers that at the population level, the pesticide exposure will vary both over time and among individuals, and consequently the effect on mortality and reproduction rates is better modelled by a continuous rather than a binary response. Extracting information from the toxicity studies and constructing the dose-time-response models requires expert knowledge and has an associated level of uncertainty. In addition to general tools for reducing expert bias such as expert knowledge elicitation [32], the advantage of modelling versus testing is the possibility for conducting sets of assessments covering the variability of the different parameters and presenting sensitivity analysis estimations.
In line with the EFSA Scientific Committee and PPR Panel [2,33] specific protection goals linked to the provision of ecosystem services can be established through five complementary dimensions: ecological entity, attribute, magnitude, temporal, and spatial scales. For wild mammals, the ecological entity is the population, and the most relevant attributes of survival and reproduction. Our model provides a tool for assessing the magnitude of the expected effects on these attributes and their temporal and spatial scale. In addition, these models offer the possibility for implementing in practice the "ecological recovery option" [34,35], which is proposed in the EFSA guidance for aquatic organisms [36], but not implemented for terrestrial vertebrates yet. Seasonal recovery is directly estimated by the model and results from the combination of the duration of the reproductive season and the background reproductive rate. The estimations presented in Figures 1-8 assume that the background reproductive and mortality rates are independent of the number of animals per nest; this is a simplification as for most animal species population's growth rate decreases with animal density [37]. The use of the recovery option requires the inclusion of this dependency in the rate equations; from the modelling perspective, this inclusion can be easily incorporated once the dependency is clarified. However, from the ecological perspective, quantifying this dependency requires specific local information elements on factors, triggering the link between population density and population dynamics, such as predation pressure, competition on habitat, and availability of food. In fact, recovery may be relevant under some conditions but not under others.
In regulatory risk assessments, assessing the effects of single or multiple applications of a pesticide on a specific crop and field is standard practice. However, the combination of effects from treatments in different fields and crops is currently lacking. The inclusion of this possibility is essential for ensuring that the risk characterisation results are informative and can be considered in environmental impact assessments [14]. Most approaches for assessing the effect of pesticide mixtures focus on mixture toxicity principles [38], assuming co-exposure. Significant research in this area is on-going [39]; however, in landscape modelling at the population level, in addition to the same individual being exposed to several pesticides, the need for combining effects of different individuals nesting on different crops is also very relevant. Offering this capacity is a key feature of this simplified model.
Predictive models reproduce sets of patterns observed at different scales and levels; the more patterns a model can reproduce simultaneously, the more reliably it may capture the essential features of a real system's organization; however if a model is too complex, its analysis will be difficult to interpret; thus there is a need for finding the optimal zone of model complexity, the so-called "Medawar zone" [40]. For regulatory use, the level of complexity should be adapted to the "regulatory question". In our model this need is covered by offering the user the capacity to decide the level of complexity adequate for each problem formulation; e.g., from a simplified landscape (single crop and edge of crop scenario), single pesticide and default parameters; to a complex landscape simulating realistic local conditions, an unlimited number of pesticides and applications, and ecological parameters adapted to actual conditions and environmental pressures. The landscape is selected by the user, based on fields specified by polygons. For local and regional assessments, the real landscape structure can be reproduced by the model, e.g., imported from GIS databases, such as those developed in the EU for implementing the Common Agricultural Practice. A simplified model with four rectangular fields is suggested for generic assessments, but the landscape can be defined in complex polygons resembling an actual or hypothetical location. The crop and crop rotation are defined for each field. The user also defines, for each crop, the pesticides and their application patterns; defining the day of application, the application rate, and the edge-of-field drift according to the intended risk management options. The nests can be located within or outside the fields. The feeding area for each group is directly estimated by the model, as well as the average twaETE according to the EFSA guidance calculations. A key element of this model is that it can refine current risk assessments using only existing standard regulatory toxicity test results. This is a key feature for population models to be used in the regulatory context [41].
There are clear challenges for incorporating higher tier data in the risk assessment of pesticides [42]. The use of the model offers a set of benefits. The exposure to and the impact of pesticides is highly dependent on the landscape agroecological conditions, and models allow the assessment of this variability. Models can also identify the most sensitive factors and prioritize further data collection and refinement [43]; this is also implemented in this model, which supports tiered approaches, moving from simple to complex assessments, which is a key principle for regulatory risk assessments.

Conclusions
Although the use of population models has been identified as a clear opportunity, the incorporation of ecological models in the regulatory context is still minimal, and most examples focus on the aquatic compartment. Based on the EFSA, good modelling practices, and other international recommendations we have developed a simplified, flexible, and versatile population model for assessing the risk of pesticides to herbivorous mammals. This landscape model combines (a) regulatory exposure scenarios, (b) standard regulatory toxicity results, and (c) ecological principles, to model the effect of pesticides on the population dynamics of rabbits and hares. The model reproduces simplified and realistic agricultural landscapes, considers crop rotation, and can run estimations for single and multiple applications of one or several pesticides per field. Background reproduction and mortality conditions (seasonality and rates) can be modified by the user and adjusted for each sub-population group (nest). A case study with estimations for two pesticides glyphosate and bromoxynil, confirms the capacity of the model for identifying the risk drivers for pesticide applications under different ecological and landscape conditions. Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ijerph18157720/s1, Figure S1: Effect of two herbicide treatments in vineyards (end of winter and spring/summer) on the evolution of brown hare abundance (total population number).  Data Availability Statement: All data used for the model calculations are publicly available at the cited references or at the EFSA webpage (www.efsa.europa.eu, (accessed on 13 December 2020)). An executable version with a user-friendly interface was developed. This version offers large flexibility to the user. The code and the executable version may be available under request to potential users interested in checking the usability of the model, with no cost for non-commercial purposes.

Conflicts of Interest:
The authors declare that they have no actual or potential competing financial interests. The author J.V.T. is employed with the European Food Safety Authority (EFSA) in the Scientific Committee and Emerging Risks Unit, which provides scientific support to EFSA's Scientific activities. However, the present article is published under the sole responsibility of the authors and may not be considered as an EFSA scientific output. The positions and opinions presented in this article are those of the authors alone and do not represent the views of EFSA. The author G.T. is employed with PharmaMar, however this work has not been conducted as part of this current affiliation. The author D.T. has conducted this work on personal capacity as independent researcher, not linked to current or past affiliations.
available. This information offers a sufficient level of detail for extracting the information to be used for adjusting the mortality and reproduction rates to different levels of pesticide exposure. Information, although with different levels of detail, is also published by other jurisdictions world-wide, including the US EPA Office of Pesticides Programs, the Health Canada Pest Management Regulatory Agency, the Australian Pesticides and Veterinary Medicines Authority, or the Food and Agricultural Materials Inspection Center in Japan.  Animals are grouped in nests, and the population dynamics is estimated for each nest, which can also be grouped for several nests. The number, sex, and age of each individual in the nest is estimated with daily iterations, adjusted by the associated survival rate for each group, the number of new-borns during the reproductive period, and pass to the next age-class when reaching the selected threshold. There are random allocations for mortality within the group and the sex of new-borns, providing environmental variability to the model estimation.
For each age-group, a feeding area around the nest is defined. Each nest is placed in a defined location within the landscape scenario. Each feeding area is then connected to the land use to estimate the percentages for the fraction of the diet obtained for each zone within the feeding area. Following the application of a pesticide, the pesticide concentration in the treated crop and other food items is estimated daily, according to the expected environmental fate. Pesticide concentrations are estimated for each commodity and weighted according to the percentage in the diet. The daily or time-weighted exposures are used as input values for the dose-response curves to estimate the expected impact of the pesticide exposure on the mortality and reproduction rates.

. Formal Model
The model variables for the different modules are summarised below. Landscape variables: Figure A1. Representation of the conceptual model. Animals are grouped in nests, and the population dynamics is estimated for each nest, which can also be grouped for several nests. The number, sex, and age of each individual in the nest is estimated with daily iterations, adjusted by the associated survival rate for each group, the number of new-borns during the reproductive period, and pass to the next age-class when reaching the selected threshold. There are random allocations for mortality within the group and the sex of new-borns, providing environmental variability to the model estimation.
For each age-group, a feeding area around the nest is defined. Each nest is placed in a defined location within the landscape scenario. Each feeding area is then connected to the land use to estimate the percentages for the fraction of the diet obtained for each zone within the feeding area. Following the application of a pesticide, the pesticide concentration in the treated crop and other food items is estimated daily, according to the expected environmental fate. Pesticide concentrations are estimated for each commodity and weighted according to the percentage in the diet. The daily or time-weighted exposures are used as input values for the dose-response curves to estimate the expected impact of the pesticide exposure on the mortality and reproduction rates.

Appendix A.4. Formal Model
The model variables for the different modules are summarised below. Landscape variables: Total area. Defined by the parcels defined by the user, accommodating any actual size. Land sectors. Number, form, and location defined by the user. Each sector (field) includes:

•
Location and form (polygons) are defined by setting coordinates for the vertices. • Associated crop, or a different land use (e.g., grassland, bare soil, forest), defined by the user, which can be changed to represent crop rotation. • Associated inner and outer edge. The inner edge allows setting in-field risk mitigation options, such as uncropped or unsprayed buffer zones. The outer edge is used for estimating pesticide contamination due to spry drift in adjacent areas. The edges can be selected for each land sector. • Residual concentration of pesticides at time zero, allowing the consideration of residues from previous exposures • Pesticide applications, the number is defined by the user and each application is defined by:

Pesticide
Date of application Dose Interception, to estimate the dose that reaches the feeding items for the species Dose ratio for the inner edge, to account for unsprayed Dose ratio for the outer edge A default four-land scenario, distributing the land in four equal squares, is available and allows most estimations for regulatory purposes. Specific local scenarios can be easily included allowing GIS modelling.
Population dynamics is modelled through groups of animals or "nests" defined by the following variables: • Species • Location, defined by coordinates • Initial population, defined by a number of individuals, average age, and age standard deviation • Breeding months • Four age groups, with defined age thresholds and food commodities, and the following associated variables that can be modified by the user: Initial male/female ratio Age range Average mobility range Maximum mobility range Background mortality rate for males Background mortality rate for females Background reproduction rate Attractiveness factor for the associated food commodities.
The toxicodynamic parameters to be established for each pesticide are: • Time variables: Use of actual or time-weighted averages for the pesticide exposure, and the averaging period when relevant.

•
The exposure-response curves, for acute mortality, chronic mortality, and reproduction • The time-scheme for updating the chronic mortality and reproduction rates The pesticide exposure estimations are based on the parameters and equations of the EFSA Guidance.