Toxicokinetic-Toxicodynamic Modeling of the Effects of Pesticides on Growth of Rattus norvegicus.

Ecological risk assessment is carried out for chemicals such as pesticides before they are released into the environment. Such risk assessment currently relies on summary statistics gathered in standardized laboratory studies. However, these statistics extract only limited information and depend on duration of exposure. Their extrapolation to realistic ecological scenarios is inherently limited. Mechanistic effect models simulate the processes underlying toxicity and so have the potential to overcome these issues. Toxicokinetic-toxicodynamic (TK-TD) models operate at the individual level, predicting the internal concentration of a chemical over time and the stress it places on an organism. TK-TD models are particularly suited to addressing the difference in exposure patterns between laboratory (constant) and field (variable) scenarios. So far, few studies have sought to predict sublethal effects of pesticide exposure to wild mammals in the field, even though such effects are of particular interest with respect to longer term exposure. We developed a TK-TD model based on the dynamic energy budget (DEB) theory, which can be parametrized and tested solely using standard regulatory studies. We demonstrate that this approach can be used effectively to predict toxic effects on the body weight of rats over time. Model predictions separate the impacts of feeding avoidance and toxic action, highlighting which was the primary driver of effects on growth. Such information is relevant to the ecological risk posed by a compound because in the environment alternative food sources may or may not be available to focal species. While this study focused on a single end point, growth, this approach could be expanded to include reproductive output. The framework developed is simple to use and could be of great utility for ecological and toxicological research as well as to risk assessors in industry and regulatory agencies.

Introduction

Before chemicals can be registered for use they undergo ecological risk assessment (ERA), this process is particularly rigorous for agricultural pesticides, which are designed to be toxic to pest species.[1,2] It is not practically possible to determine a chemical’s ecological impact experimentally. Instead, risk assessment relies on extrapolation from summary statistics such as the “no observed adverse effect level” (NOAEL) generated for a few species in standardized laboratory studies.[3] However, such statistics should only be extrapolated with caution as they do not account for the processes that lead to toxic effects and are dependent on duration of exposure.[4,5] Moreover, they are generated in controlled (supposedly optimal) conditions—regulated temperature and freely available food and water[6]—that are unrealistic in the field. The resulting lack of ecological realism in current standard risk assessment methods is problematic.[7]

Mechanistic effects models (MEMs) aim to simulate the mechanisms by which chemicals affect individuals, populations, and communities and therefore enable us to predict how they will respond in untested and more ecologically relevant conditions.[3] This is an appealing prospect with great potential for use in ERA of pesticides.[8,9] By focusing on the underlying processes, modeling techniques can add ecological realism to extrapolations and even reduce animal testing requirements.[10]

Accounting for the mismatch in exposure between laboratory and field[11] was identified as one of five key obstacles to long-term risk assessment of pesticides for mammals (along with selection of suitable toxicity end points, extrapolation of toxicity between species, exposure assessment, and evaluation of population level effects[12]). In chronic toxicity tests rats or mice are exposed to a constant concentration of a pesticide in their diet for periods as long as 2 years.[6,13−15] Such constant exposure is unrealistic in the field where pesticides are not applied at a constant rate all year round. This disparity can be addressed through the use of toxicokinetic-toxicodynamic (TK–TD) modeling.[10] TK–TD models work at the individual level, predicting an internal measure of chemical concentration over time (toxicokinetics) and the stress this places on an organism (toxicodynamics). As such, the effects on a given end point resulting from time varied exposure can be predicted.[16]

The use of TK–TD modeling has now been recommended for certain regulatory purposes, such as predicting survival of aquatic organisms.[17] However, European protection goals for birds and mammals state that there should be no visible mortality associated with pesticide[1] use so sublethal effects are more relevant with respect to realistic exposure. Sublethal effects can be predicted using the “DEBtox” modeling framework,[18,19] combining TK–TD modeling with the Dynamic Energy Budget (DEB) theory.[20] DEB is an established metabolic theory which has been applied to a range of taxa,[21] mathematically describing the processes of energy acquisition and allocation that determine the life history of an organism. Using TK–TD modeling to place stress on these processes can produce predictions of effects on sublethal end points such as growth and reproduction. Very little research has concerned mammals, however, as DEBtox studies have thus far mainly focused on invertebrates[22] and more recently fish.[23,24]

At present DEBtox is limited to research applications as it is not regarded as user-friendly enough for use by regulators.[25] To this end, a simplified version of the theory, “DEBkiss”, was developed in which only structural body mass (bones, muscle, organs, etc.) is considered with no reserve storage.[26] The model retains many DEB principles but with fewer parameters and model equations. It was developed for applications where simplicity and ease of use are important, such as the analysis of toxicity data or for use within individual based population models. The only published study to date in which TK–TD modeling has been used to predict sublethal effects in mammals utilized DEBkiss.[27] Although limited data were available, the model accurately simulated observed effects of environmental toxicants on growth and reproduction in the American mink (Mustela vison). These results suggest that this simplified framework may be sufficient for practical and regulatory applications.

Here we tested the utility of DEBkiss by working with raw data from repeated dietary dose toxicity tests and modeled the effects of several pesticides on rats. As the first study to use regulatory data for this purpose, we adopted the practice of beginning with the simplest possible methods and identifying areas where more complex techniques may be required. Internal pesticide concentration was modeled with a one compartment model and a single end point, body weight, was modeled over time using the DEBkiss growth model. Our study aimed to establish a practical procedure for the parametrization, calibration, and validation of DEBtox models using regulatory data and to assess the quality and utility of predictions. Furthermore, we aimed to improve the interpretation of standard toxicity studies by extracting novel, meaningful information using modeling.

Methods

Data

All data used here were made available from existing regulatory studies (Syngenta, unpublished) carried out according to 94/79/EC,[28] investigating toxicokinetics and chronic toxicity of acibenzolar-S-methyl (benzothiadiazole; fungicide, insecticide and plant activator[29]), azoxystrobin (strobilurin; fungicide[30]), fenpropidin (unclassified; fungicide[31]), fludioxonil (phenylpyrrole; fungicide[32]), mandipropamid (mandelamide; fungicide[33]) and prosulfuron (sulfonylurea; herbicide[34]), and thiamethoxam (neonicotinoid; insecticide[35]) in laboratory rats (Rattus norvegicus).

The toxicokinetics studies followed OECD 417[36] guidelines. Animals were treated with a single oral (gavage) dose of a 14C radiolabeled pesticide with total radioactivity found in various tissues and excreta monitored over a period of days. The animals were allowed free access to a certified standard diet.

Data sets differed between pesticides but followed a common framework. At least two dose levels were studied with typically three male and three female animals in each treatment group. The reports include data detailing the proportion of the initial dose excreted in faeces, urine, and bile over ∼48 h after a single oral dose, providing an average percentage of the dose which was absorbed into the body. Pesticide concentration in the blood of animals were measured over ∼48 h following a single high or low dose. Pesticide concentration were also measured in different body tissues from animals terminated at ∼4 time points following a single high or low dose. Details of dosing, including exact dose (mg(AI)), body weight (g) at the start (and in some cases the end) of testing and achieved dose (mg(AI) × kg(BW)–1) were provided for each individual animal.

Chronic toxicity studies lasting 28 days (OECD 407), 90 days (OECD 408) or 2 years (OECD 414 or 416) were carried out according to OECD guidelines.[6,13−15] Animals of around 5–7 weeks in age were provided with diet containing pesticide and multiple toxicological end points monitored over the study period. Animals were kept in standard conditions with food and water available ad libitum. Each study provided individual weekly observations of body weight (g) and food consumption (g(diet) × day–1). Sample size was typically 5 animals per sex per treatment in 28 day studies, 10 animals per sex in 90 day studies and 50–80 animals per sex in 2 year studies. The measured concentration of pesticide in the diet of each treatment was also reported. Each study comprised a control group and at least 3 treatments fed diets containing different concentrations of pesticide. As study duration was increased dietary doses were generally decreased.

Models

Toxicokinetic Model

For each toxicant (denoted “AI” for active ingredient), the internal concentration was modeled using a one compartment TK model with first order kinetics. The internal concentration here refers only to toxicant present in body tissues at a given time and excludes any in the gut which has not yet been absorbed. As terrestrial mammals are primarily exposed to pesticides via the diet, the toxicant concentration in the gut was also modeled as an intermediate “depot” compartment. This was a very similar approach to Bednarska et al.[37] but we also account for change in body size by including dilution by growth and changes to surface area to volume ratio as per Gergs et al.[38] The model equations are shown below.where ΔC indicates change in the body weight (denoted “BW”) normalized dose, C, of toxicant over time (Mass(AI)× Mass(BW)–1 × t–1) and subscripts Gut and Int denote gut and internal respectively; I is toxicant ingestion rate (Mass(AI) × Mass(BW)–1 × t–1); F represents bioavailability (dimensionless); ka and ke represent the rate constants of toxicant absorption from the gut and toxicant elimination from the system respectively (t–1); W denotes body weight, and ΔW is change in weight over time; L is volumetric length (the cube root of body volume), and L is the ultimate volumetric length of the test species.

Growth Model

Body weight over time was modeled with the DEBkiss[26] growth equation:where ΔW represents the change in total body weight (W) over time, yVA represents the efficiency with which assimilates are converted to structural mass (Mass(BW) × Mass(Assimilates)–1), k is the proportion of assimilates allocated to the soma, f is the scaled feeding rate (unitless), JaAM is the maximum surface area specific assimilation rate (Mass(Assimilates) × Mass(BW)–2/3 × t–1) and JvM the mass specific maintenance rate (Mass(Assimilates) × Mass(BW)–1 × t–1). The DEBkiss model is represented graphically in Figure .

Figure 1

A graphical representation of the DEBkiss model. The value of k determines the proportion of resources assimilated from food allocated to maintenance and growth or maturity and reproduction. Processes outlined in red are those that can be subjected to stress.

Endotherms are also subject to surface area specific maintenance costs, accounting for heat loss to the environment. However, as long as the ambient temperature is within the thermoneutral zone of a species[39] these are assumed to be zero.[40] Laboratory guidelines require rodents to be kept at 22 ± 3 °C, as this was considered to be within the thermoneutral zone of the rat.[41] More recent research has suggested that this temperature range is too low[42] but for simplicity it was assumed that heat loss could be omitted.

The parameter k represents the proportion of assimilates allocated to maintenance and growth with the remainder (1-k) going toward maturation and reproduction. DEBkiss theory states that up to puberty, these resources are used up as animals develop and are stored as “reproduction buffer” during puberty and adulthood.[26] However, as reproduction was not modeled in this study, the reproduction buffer would serve no purpose other than to make the model more complex. Instead, it was assumed that the 1-k branch continues to be used up as the animals develop into sexually mature adults (at 70–90 days of age[43]) and then to maintain maturity, a process which can also be included in DEBkiss.[26] Only unmated animals were included in this study, so body mass was not impacted upon by pregnancy.

If food intake is reduced such that kf JAMaW2/3 < JMvW but the total assimilation rate f JAMaW2/3 ≥ JMvW (i.e., in a situation where food intake is sufficient to maintain homeostasis but not to grow), then ΔW = 0 as available resources are diverted from the 1-k branch to meet maintenance costs. If the total assimilation rate is insufficient to meet maintenance costs, that is f JAMaW2/3 < JMvW, then the growth rate becomes negative as tissue is metabolized to meet maintenance requirements.where yAV is conversion efficiency of structure to assimilates. The value of k therefore determines the point at which the feeding rate becomes insufficient for growth but does not impact the onset of weight loss. The point at which weight loss results in death is unknown due to a lack of data. Any treatment that induced drastic reductions in feeding would be abandoned as guidelines state that dosing should not cause “death or severe suffering”.[14,15]

Toxicodynamic Model

Finally, the DEBtox toxicodynamic model[18,19,44] was used to link internal toxicant concentration to stress on growth. It is assumed that for any xenobiotic, there exists a “no effect concentration” (NEC) below which it causes no stress to any biological processes. The NEC is a time independent threshold and therefore has no relationship with duration of exposure. Assuming that every molecule of toxicant beyond its NEC creates the same amount of “stress” leads to a “linear with threshold” relationship, which can be modeled quite simply (Figure ):

Figure 2

Stress increases with internal toxicant concentration beyond a threshold, where S is dimensionless stress and is the “tolerance concentration” (Mass(AI) × Mass(BW)–1). Here the NEC determines the point at which stress exceeds zero while CT is the increase in CInt corresponding to an increase in S of 1. This means the gradient of S is 1/CT when CInt exceeds the NEC.

In our approach, stress can be applied to one of three growth parameters (“physiological modes of action”[19,45]) each of which respond differently to stress, these are the maximum assimilation rate, JaAM, the maintenance rate, JvM, or the conversion efficiency, yVA (Figure ). The proportion of resources allocated to the soma, k, could theoretically be affected by a toxicant, but data on reproduction would be required to distinguish this from effects on assimilation and such effects are not well documented.[44]

Figure 3

Plots showing how growth model parameters respond to internal toxicant concentration when stress is applied. A. The maximum assimilation rate JaAM decreases linearly with stress until it reaches zero when S = 1. B. The maintenance rate JvM increases linearly with stress and is doubled when S = 1. C. As costs per unit of tissue synthesis increase linearly with stress, the conversion efficiency yVA approaches zero asymptotically, efficiency is halved when S = 1.

Model Implementation

All models were implemented in Matlab (ver. R2016b). TK and TK–TD growth models were developed with the BYOM[46] flexible model platform (ver. 4.1), several additional functions and scripts were also developed as part of this study. All fitted parameter values were derived using the Nelder Mead simplex algorithm to maximize the likelihood function, given the observed data.[47] Likelihood profiling was also used to check that initial fits were not local optima.[48] TK model parameters were fitted to mean internal pesticide concentration over time, while growth and TD model parameters were fitted to mean body weight over time (Table ).

Table 1

Full List of Model Parametersa

name	explanation	value	dimensions
growth
k	fraction of assimilates for growth and maintenance	0.8*
JaAM	maximum assimilation rate per unit of surface area	fitted to data	g(Assimilates) × g(BW)–2/3 × t–1
J vM	maintenance rate per unit of biomass	fitted to data	g(Assimilates) × g(BW)–1 × t–1
yVA	new biomass per unit assimilates	0.45 (as per Sibly & Calow[51])	g(BW) × g(Assimilates)–1
yAV	yield of assimilates per unit biomass	0.8*	g(Assimilates) × g(BW)–1
toxicokinetics
F	proportion absorbed from gut	calculated from data
ka	absorption rate constant	fitted to data	t–1
ke	elimination rate constant	fitted to data	t–1
toxicodynamics
NEC	no effect concentration	fitted to data	mg(AI) × kg(BW)–1
CT	tolerance concentration	fitted to data	mg(AI) × kg(BW)–1
other parameters
f	scaled feeding rate	calculated from data
I	pesticide ingestion rate	calculated from data	mg(AI) × kg(BW)–1 × t–1
L	volumetric length	W1/3	cm
L∞	ultimate volumetric length	W∞1/3 or	cm
Lm	maximum volumetric length	Wm1/3 or	cm
W∞	ultimate structural body mass	782 (as per Hubert et al.[52]) or L∞3	g
Wm	maximum structural body mass	Lm3	g
state variables
CGut	pesticide concentration in gut (ΔCGut)		mg(AI) × kg(BW)–1 × t–1
CInt	internal pesticide concentration (ΔCGut)		mg(AI) × kg(BW)–1 × t–1
W	structural body mass (ΔW)		g(BW) × t–1

Parameter values marked with a “*” are default values suggested by Jager, Martin, & Zimmer.[26].

Toxicodynamic modeling also required selection of a physiological mode of action. The best fitting physiological mode of action was determined using the Akaike information criterion.[49] As the alternative models were not nested, the likelihood ratio test would be inappropriate however the AIC has no such requirement.[50]

TK Modeling

Parameterization

The percentage of the dose which was absorbed into the body was reported in excretion studies. These data provided the value of the parameter F in the TK model (eq and 2). If the percentage absorbed was very high (>90%), then F was assigned the value 1 as this represents the worst-case scenario. For single dose studies the value of I was zero while the starting gut concentration was the average achieved dose for each treatment. If body weight was recorded at the beginning and end of testing, then ΔW was calculated as the linear growth rate observed in each treatment. If not, then it was assumed that ΔW = 0 as dilution by growth has a minimal effect on the model over the short testing period, typically 48 h.

Ultimate length was calculated as the cube root of 782 cm3 which is the average ultimate volume of male Sprague–Dawley rats with ad libitum food availability,[52] assuming that average wet tissue density is equal to that of water.[40]

Calibration

Having determined the other model inputs directly from experimental data, two free TK model parameters were left to be fitted to data, the absorption and elimination rate constants k and k. The best time course data for internal toxicant concentration (highest number of time points) came from the blood as it can be sampled without terminating animals. In order to ascertain whether pesticide concentration in the blood was a suitable proxy for overall body burden, the Pearson correlation coefficient was used to determine if it was significantly correlated with the concentration in all other sampled tissues. If so, then blood concentration data (in μg× g(blood)–1) were deemed representative of CInt (in mg× kg(BW)–1) and could be used to fit ΔCInt (eq ). Where available, whole carcass concentration could be used as an alternative, however, fewer data points were available.

Data were first inspected to determine whether there were obvious differences in kinetics between males and females. Properties considered were the mean peak concentration (Cmax), the time after dosing at which it was reached (Tmax) and the time taken to eliminate the pesticide from the blood. If clear differences were evident, then models were calibrated separately for males and females.

For some pesticides, the same individual animals were sampled for the whole observation period after dosing, for others only partial time course data were available for each individual. Where complete time course data were available, the model was first fitted to data for each individual. A multiway ANOVA was conducted to determine whether there were significant differences in fitted values of ka and ke associated with sex, dose level (both discrete) or weight (continuous). If more than one radiolabel was used, then this was also included as a factor.

If sex was the only factor to have a significant effect, then the model was calibrated to males and females separately. If both sex and weight were shown to have significant effects, then a Mann–Whitney U test was used to determine whether there was a significant difference in weight between males and females. If so, then the ANOVA was repeated with the data separated by sex. If no effects were found within each sex, then the model was simply fitted to males and females separately. If the rate constants were significantly affected by factors other than sex, then this was noted along with the full TK results (Supporting Information, SI).

Finally, models were calibrated to mean (at each time point) blood concentration observed in the high and low dose groups simultaneously. Where appropriate, this was carried out separately for males and females.

Validation

Toxicokinetics studies sampled pesticide concentration found in different body tissues, including blood, from animals terminated at different time points following a single oral dose. Blood samples were arranged to provide an independent time course data set of internal toxicant concentration. These were then used to assess the model performance by comparing the pesticide concentrations in the blood of the terminated individuals with concentrations in blood simulated by the model with fixed parameters.

Generally, these data covered the same two dose levels as the calibration data set so could not provide validation per se, since the model inputs (i.e., dose levels) were virtually the same. Any differences in achieved doses were generally very small so differences between observed toxicokinetics were primarily due to individual variability. This did however provide an indication of how well the true average response was represented by the calibration data set and therefore the calibrated model.

In most cases blood concentration data sets included only three samples for each time point and responses could be highly variable between individuals. We did not carry out any quantitative assessment of predictions because this would be misleading due to the low sample size combined with strong interindividual variability

TK–TD Growth Modeling

The full TK–TD growth model comprises all the model eqs (eqs , 2 and 3) and simultaneously predicts toxicant concentration in the gut (mg(AI) × kg(BW)–1), internal toxicant concentration (mg(AI) × kg(BW)–1) and body weight (g).

The growth parameter k describes the proportion of resources allocated to maintenance and growth. When only modeling growth, its precise value is not crucial (only its product with the fitted parameter JaAM contributes to the model) so this was fixed at its default value of 0.8,[26] which estimates suggest is reasonable for the species.[53,54] Physiological studies suggest assimilated energy is converted to new tissue by homeotherms with an efficiency between 0.4 and 0.5[51] so the parameter yVA was fixed at 0.45. If strong evidence were provided suggesting different values for either of these parameters, then JaAM and JvM would simply need to be adjusted by the appropriate correction factor. Stress functions would continue to have the same impact and so would not need adjustment.

In repeated dose toxicity tests each animal’s weight and food consumption were recorded at least weekly for some or all of the study period.[6,13−15] In some cases, food consumption was recorded per cage, so the values provided were an average per animal but each weight measurement had a corresponding measurement of food consumption. The achieved toxicant concentration in the diet of each treatment group was also measured. This allowed the growth parameter f (scaled feeding rate) and the TK parameter I (toxicant ingestion rate) to be calculated directly from the data.

The maximum feeding rate at a given food density is assumed to be proportional to surface area. The scaled feeding rate, f, is equal to an individual’s actual feeding rate at a given food density divided by the maximum feeding rate for its size and therefore ranges from 0 to 1.[26]

Measured food consumption per day was converted into surface area specific feeding rate by dividing by the associated body weight raised to the power 2/3. Dividing these values by the maximum feeding rate recorded in the study group (separated by sex) provided scaled f values between 0 and 1 for each individual in each week of the study period. A matrix was then produced containing average weekly feeding rates for each treatment group, these provided the value of f for each treatment in each weekly interval.

Multiplying the achieved toxicant concentration in the diet (mg(AI) × kg(diet)–1) by the mass specific daily feeding rate (kg(diet) × kg(BW)–1 × d–1) provides the ingested dose (mg(AI) × kg(BW)–1 × d–1). Again, a matrix was produced, this time containing the weekly averages of daily ingested dose for each treatment. This provided the values of the toxicant ingestion rate I which fluctuates with feeding rate throughout the study period. All other TK parameters remained fixed at the values determined during TK model calibration (ka and ke were multiplied by 24 to convert them from hourly to daily rates).

The data from 90 day toxicity studies (OECD 408) were intermediate in terms of sample size and dietary dose levels. Thus, these data were more representative than the 28 day studies (OECD 407) while the observed effects on growth were generally larger than in 2 year studies (OECD 414 and 416). For this reason, the 90 day studies were used to calibrate the TK–TD model. Calibration of the growth and TD parameters was conducted separately for males and females.

The two free growth model parameters, JaAM and JvM were fitted to the growth data from the control group and then fixed. Next, the model was run for all treatments but with no stress applied. This step was used to identify the lowest dose group in which observed growth was lower than predicted by the growth model based on feeding rate alone (i.e., the lowest dose group in which chemical stress occurred). The initial value for the estimation of the TD parameter NEC was then set to half the average internal concentration predicted by the TK model for that treatment. Only two treatments (the lowest affected treatment and the top dose) were used for fitting as this allowed more of the data to be used for testing predictions while still providing a wide range of internal concentration predictions.

The TD parameters NEC and CT were then fitted to the lowest affected treatment and the top dose group simultaneously, this was repeated for each physiological mode of action. The fit which produced the lowest AIC value was selected. The physiological mode of action and resulting TD parameter values were then fixed.

For verification, the resulting model was then run for all treatments in the 90 day study producing interpolations to the intermediate dose groups. Then, the model was used to predict growth in the 28 day and two year studies. Only the first 12–14 weeks of growth data from two year studies were used to test model predictions. The reasons for this cutoff point are addressed in detail in the discussion.

Matrices containing weekly averages of feeding rate, f, and toxicant ingestion rate, I, were generated from the data as described previously. While the TK and TD parameters remained fixed, it was necessary to repeat the fitting of growth parameters (JaAM and JvM) to the control group data. This was important so that the effects of feeding rate on growth were predicted relative to the control group of each study rather than to that of a separate study in which conditions (laboratory rat strain, feed, average temperature) may have differed.

All model parameters were then fixed and the model was used to predict effects on growth in all treatments. Predictions were compared to observed data for each treatment at each observed time point. Predictions were considered in terms of animal weight (g) or the proportional effect on body weight relative to the control group (mean weight(treatment) × mean weight(control)–1). Predictions were deemed accurate if they were within one standard deviation of the mean observed value at each time point as this measure takes into account the individual variability within the data. The percentage of predictions that were accurate was reported and any exceptions were noted. Exceptions were used to determine the limitations of the model and to infer the underlying reasons why they arise.

Results

Results are summarized in this section, however, due to the number of studies used to calibrate and test models, it was not practical to include all tables and figures here. These can be found in the SI.

Toxicokinetics

For all the pesticides the concentration reached in the blood was significantly correlated with that in all other tissues sampled (p < 0.01) except for the concentration of fenpropidin in fat which was not significant (p = 0.0574). However, when one outlier (residual >3 s.d. from mean) was removed from the analysis, the correlation was highly significant (p < 0.001).

All the compounds exhibited first order kinetics, producing blood concentration time curves which could be reproduced by fitting of a one compartment TK model. In some cases adjustments to the modeling procedure were required which will be described in turn.

For thiamethoxam, global model fits at both high and low dose levels closely matched the observed data. When model predictions were tested against independent blood concentration data, the predicted curves again closely emulated the observed data.

For four of the compounds, acibenzolar-S-methyl, azoxystrobin, fludioxonil and fenpropidin, the global model fits better represented blood kinetics in the high dose group, with modeled curves not reaching the peak blood concentration (Cmax) observed in the low dose groups. In all these cases this same pattern was observed when model predictions were tested against an independent blood concentration data set. The likely explanation for this phenomenon is that, as pesticide concentration in the gut is increased, absorption rate becomes saturated and reaches a maximum.[55] However, with only two dose levels tested in most toxicokinetics studies, generally differing by a factor of at least 100, it is not possible to estimate the point at which this occurs or to determine whether it is a gradual process or happens suddenly. Lower model accuracy at low internal concentrations has little impact upon eventual predictions of effects, and none at all if below the NEC. As the high dose levels were more relevant to the dietary ingestion rates associated with effects on body weight, these parameter values were accepted.

Male rats administered a high dose (100 mg × kg(BW)–1) of fenpropidin appeared to exhibit a double peak in the concentration reached in the blood. An initial peak was reached 1 h after dosing with a second, lower peak after around 8 h. Double peaks have been attributed to variable absorption in different regions of the gut, enterohepatic recirculation or delayed gastric emptying.[56] Were the dose delivered at a more constant rate in the diet rather than as a single large dose one would expect the overall rates of absorption and elimination to reach equilibrium and so any of these mechanisms would have less impact on blood kinetics. As the intended application of the model was to predict the effects of dietary dosing, it was decided to fit the model for males with the data collected 2 h, 3 h, 4 h, and 6 h after dosing excluded in the high dose group. This allowed the model to fit a single peak in which the observed Cmax, the time at which it was reached (Tmax), and time for total elimination of the dose were matched closely by the model.

Following administration of a high dose, peak concentration of prosulfuron in the blood (μg(AI) × g(Blood)–1) exceeded the body weight normalized dose (mg(AI) × kg(BW)–1). This could not be modeled by our TK model (eqs and 2) which uses blood concentration as a proxy for overall internal concentration. However, blood only accounts for around 7% of body mass[57] so in reality only >7% of the dose needs to be present in the blood at one time for this to occur. Nevertheless, this phenomenon was unusual among the chemicals included in the study.

In order to address this, the relationship between prosulfuron concentration in the blood (CBlood) and overall carcass concentration of terminated animals was investigated. Only male animals were used in the tissue sampling experiments with three animals sampled at each of four time points following a high or low dose. Blood and carcass concentration were strongly and significantly correlated (Pearson’s correlation coefficient, r = 0.99, n = 24, p < 0.0001) . The line of best fit, intercepting the y axis at zero, was derived by finding the least-squares solution to the equation CBlood= XCInt. The gradient, X = 2.4337, was determined as the concentration factor by which prosulfuron concentration in the blood exceeds that in the body as a whole.

A third equation, incorporating that concentration factor but otherwise identical to eq , was then added to the TK model to describe blood concentration over time as follows:This determines that, CBlood = 2.4337 × CInt at any given time point, with ΔC modeled by eq . Blood and whole body internal concentration could then be modeled simultaneously. The fitted model produced curves matching the data well for both variables at the high dose.

In contrast to other compounds the prosulfuron model predicted higher than observed internal concentrations at the low dose level. This was the case for concentration in the whole body as well as in the blood so the concentration factor was not the cause of the discrepancy. In fact, for the low dose group alone the concentration factor was higher than the overall figure. A possible explanation in this case is that the elimination rate becomes saturated beyond a certain internal concentration,[58,59] this would be consistent with the unusually high internal concentrations measured in the high dose groups. The high dose in the prosulfuron toxicokinetics experiments was around 450 mg(AI) × kg(BW)–1, several times higher than the high dose used for other chemicals included in this study (100 mg(AI) × kg(BW)–1). Had the other chemicals been tested at such high doses it is possible that a similar pattern would have been evident. Unfortunately, only two dose levels (differing by a factor of around 900) provide insufficient data to determine the maximum elimination rate and the internal concentration at which it is reached.

With fixed parameters and independent data, the model again predicted higher prosulfuron concentration in the blood than was observed at the low dose. The data were predicted well at the high dose however. For both sexes, Cmax and Tmax were predicted with reasonable accuracy. Elimination of the compound was slower than predicted in females, but the parameters were deemed acceptable and the model was not fitted to males and females separately.

Growth

The fitted DEBkiss growth curve was able to accurately model growth of rats aged around 6–20 weeks. In total, the model was fitted to 34 control group data sets comprising weekly observations of body weight and food consumption rate. Modeled body weight was within 1 standard deviation of the observed mean at all time points in 30 out of 34 cases and at >90% of time points in 32 out of 34 cases. The deviations were most pronounced in two data sets. For the female control group in the 28 day toxicity study of fenpropidin, the modeled body weight was lower than the observed mean by more than one standard deviation at week one only. As a result, only 75% of the modeled weights were within one standard deviation. For the male control group in the 90 day toxicity study of azoxystrobin, the modeled body weight was lower than the mean by more than one standard deviation in weeks two and four. This resulted in only 84.6% of the modeled weights being within one standard deviation of the observed mean. For every data set, all predictions of body weight in the control group were within 10% of the mean at all time points.

Toxicodynamics

With the growth and toxicokinetic parameter values fixed, the toxicodynamic parameters, NEC and CT, were fitted to selected data as described in the Methods section. Predictions were interpolated to other treatments in the 90 day studies and extrapolated to 28 day and 2 year studies.

In terms of body weight ≥75% of predictions were within 1 standard deviation of observed means for 28 out of the 34 study groups. In terms of effect on body weight, ≥ 75% of predictions were within 1 standard deviation of observed means for 30 out of the 34 study groups. A summary of results is shown in Table . For female rats administered thiamethoxam or fludioxonil, TD parameters were fitted to a single treatment group. This was because, in studies of both compounds, body weight reductions beyond those predicted based on feeding rate were only evident in the top dose group. In the case of thiamethoxam, significant body weight reductions were only observed in females dosed with 10 000 mg/kg(diet) of thiamethoxam over 28 days. Females administered fludioxonil in their diet did show significant body weight reductions. However, these were predicted entirely based on reduced feeding rate in all but one treatment, those dosed with 20 000 mg/kg(diet) over 90 days. The consequence was that while the NEC was consistent with data from several treatments, CT was determined using data from only one treatment. However, since toxicant ingestion rate was dynamic, CT was fitted to a range of internal concentrations even within one treatment.

Table 2

Toxicodynamic Parameters Used to Model the Effects of Each Compound on Male and Female Ratsa

The percentage of predictions (in terms of absolute body weight and effect on body weight relative to the control group at each time point) within one standard deviation of the observed mean, are shown. Percentages ≥75% are highlighted in green, those of ≥50% and <75% are shown in blue while those <50% are highlighted in orange. Those marked with a “*” were fitted to only one treatment group. pMoA: best fitting physiological Mode of Action.

Model predictions were also used to investigate the extent to which reductions in body weight could be attributed to reduced feeding or direct toxic action. This was done by comparing experimental data to model simulations in which no stress was applied, producing growth curves predicted based solely on feeding rate. Comparisons between expected body weight modeled with actual feeding rates and observed body weight were conducted separately for every treatment group. The observed and predicted body weights in each treatment were converted to proportions of the control body weight (observed and predicted respectively) at each time point. The proportion of any observed body weight reductions (relative to controls) that were predicted based on feeding rate alone could then be calculated. The remainder was attributed to toxic action. If body weight predicted based purely on actual feeding data was below that observed at a given time point, then any observed weight reduction (relative to controls) was attributed entirely to reduced feeding rate. Likewise, if body weight predicted based on actual feeding rate in a given treatment group was higher than controls, then any observed weight reduction was attributed entirely to toxic action. An example of this process is shown in Figure ; selected results for all compounds are shown in Figure .

Figure 4

Graph showing growth modeled based on feeding rate only (lines) and observed body weights (circles) of male rats. The control group and those dosed with 20 000 mg × kg(diet)–1 fludioxonil are shown. The proportional breakdown of the observed reduction in body weight of treated rats vs controls at the end of testing is represented in a bar chart.

Figure 5

Bar charts showing the proportion of observed weight reductions relative to the control group attributed to reduced feeding rate and/or toxic stress by the growth model. All treatments in which a weight reduction was evident at the end of the analyzed period are included. X-axis labels denote the observation date and dietary dose, in some cases treatments were duplicated between studies. No bar is displayed where there was no reduction in weight. Note that bars are the same size regardless of the magnitude of the observed effect.

For azoxystrobin, prosulfuron, thiamethoxam, and fludioxonil there appeared to be a pattern across the sexes, with reduced growth being driven more by feeding rate in females and by toxicity in males (Figure ). This was most evident in the case of fludioxonil, which was associated with significant body weight reductions in both sexes. While reductions in male body weight were attributed largely to toxicity, the reductions observed in females were predicted based entirely on reduced feeding rate in all but the highest dose group across two studies. A similar pattern was seen for thiamethoxam. Once again toxic effects were only predicted to impact upon female bodyweight in the highest dose group. In this case, however, reductions in female body weight were not observed in most treatments as feeding rate was not affected either.

Discussion

Predicting Growth Under Chemical Stress

Raw data from chronic toxicity studies were used to test DEBtox predictions of sublethal toxic effects in mammals. Weekly measurements of body weight and food consumption as well as precisely measured dietary concentration provide good quality data with which to calibrate models. As the regulatory framework requires several such studies, abundant data are available to test predictions.[6,13−15] As has been noted previously, the way science is generally funded means that corroboration studies are rarely conducted in academia or published in the scientific journals as they are not considered novel.[60]

Our findings showed good agreement between predictions and data. With minimal model fitting, observed effects on body weight were predicted reliably (≥75% of predictions accurate to within one standard deviation of the observed mean in all studies for which data were available) in males and females for four of the six chemicals modeled. This suggests TK predictions were at least proportional to actual internal concentration over time and that the “linear past threshold” TD model (Figure ) is based on reasonable assumptions.

Feeding Rate vs Toxicant Ingestion

An obstacle when analyzing the effects of dietary toxicant exposure is that the ingested dose depends as much on the feeding rate of the study animals as it does on the concentration in the diet. While a high feeding rate will have a positive effect on growth, the corresponding toxicant ingestion will place stress on growth parameters. Understanding of this trade-off is important[61] as ingestion is considered the primary exposure route for terrestrial mammals to pesticides in the field[37,62] although it is argued that other routes should receive greater attention.[63] In the modeling approach used here, both the scaled feeding rate and the pesticide ingestion rate were calculated as dynamic model inputs, directly from the data. This allowed model predictions to separate the competing effects of feeding rate and toxic action on growth. For example, male rats given 2500 mg/kg(diet) thiamethoxam grew larger than those given 1250 mg/kg(diet). This result was correctly predicted by the model as the higher feeding rate of the 2500 mg/kg(diet) group partially counteracted the chemical stress.

Model predictions can therefore provide new insight into the observed data by comparing data to model simulations in which no stress is applied. Such predictions only require the calibrated growth model and data on feeding rate. Thus, even if toxic effects cannot be reliably predicted, growth model predictions can indicate the degree to which observed reductions in body weight were driven by toxicity or reduced feeding relative to controls. Such information is valuable for assessing the risk that a chemical poses to terrestrial mammals in the field. If strong avoidance is observed, then this may increase or decrease the risk posed depending on whether animals would have a choice of food items in the field scenario.[61,62]

For several compounds there appeared to be a pattern across the sexes, with reduced growth being driven more by feeding rate in females but by toxicity in males. This pattern was strongest for fludioxonil and thiamethoxam and in both cases was reflected by a large difference in the values assigned to the NEC for each sex. This would suggest that while females exhibit a higher tolerance for these compounds, at least with respect to growth, they show stronger feeding avoidance. Such inconsistency in the toxicodynamic parameters of males and females may seem surprising however, large differences between reference doses for each sex have long been documented in rats.[64] Moreover, the results were unequivocal with respect to thiamethoxam, as females were unaffected at several dose levels which affected males. In several cases, differences have been noted in the sensitivity of the liver and kidneys to toxicity,[64,65] possibly related to variable enzyme production between the sexes.[66] The liver and kidneys were identified as the target organs of these pesticides in mammals so these results would appear consistent with previous findings (Syngenta, unpublished). Given that such differences in chemical sensitivity can occur between the sexes it should not be surprising either that the models suggested different modes of action for several of the compounds in males and females (Table ). Mode of action in DEBtox refers to abstract processes rather than specific chemical pathways so, this simply implies that the effects on the growth curve differed between the sexes.[44]

Model Limitations

A fundamental limitation to any model of a complex system is the trade-off between realism and simplicity. As user friendliness is a significant consideration for regulatory use,[25] DEBkiss was selected as the simplest possible approach to investigate how raw lab data should be utilized to parametrize and calibrate models with data from dietary toxicity studies. Another reason for prioritising simplicity was so that potential issues could be clearly identified at this early stage for future models to develop, with scope for further elaboration if necessary. Several such issues were highlighted by our results which are discussed in this section.

While our results demonstrate that DEBtox is a useful framework, predictions of growth under chemical stress were not reliable in all cases. The effects of fenpropidin and acibenzolar-S-methyl were predicted accurately in several treatments, but overall model accuracy was lower than for the other chemicals. When considering in vivo effects, any observed deviations from predictions are unlikely to result from measurement error, more often, individual variability is the cause. This represents a hurdle to predictive modeling, even for genetically similar laboratory strains kept in controlled conditions. Indeed, individual variability in growth is still evident in studies using genetically identical springtail (Folsomia candida) clones kept individually and provided unlimited food.[67] As such, we cannot know for sure how treated animals would have grown had they been in the control group but must assume constant growth parameters across treatments. Consequently, models provide predictions of mean body weight over time under specified conditions but do not include uncertainty or variability.

Another consequence of individual variability is that it can be difficult to identify the underlying causes of model inaccuracy. For example, several treatment groups of female rats fed lower dietary doses of acibenzolar-S-methyl grew larger than controls despite feeding at a lower rate. This is clearly a result the model would not predict and could simply be the result of variability in average growth parameters between treatment groups. Alternatively, this could be interpreted as evidence of hormesis, the phenomenon by which lower doses of a chemical have the opposite effect of higher doses on a given end point. If this were the case, then the stress function rather than the growth parameters would require alteration, but we cannot be sure which. In a few treatments the opposite issue arose when the animals fed at a relatively high rate but did not grow as expected. Since feeding rate was calculated based on actual body size rather than predicted size, modeled growth continued at a higher than observed rate. For these treatments, model predictions far exceeded observed growth. This issue could be somewhat resolved by calculating feeding rate relative to the predicted body size over time. However, it is likely that variability in growth parameters also played a role.

Furthermore, growth is not the only modeled property subject to individual variability. Tolerance to a toxicant[68] or the rates at which it is taken up and eliminated[37] may be highly variable among individuals in a population. Intertreatment differences in any one or more of these properties could result in observed effects not being uniformly correlated with internal concentration and are therefore difficult to predict without further knowledge. This was evident for mandipropamid; results were not reported as the data were not suitable to test model predictions. In the 90 day study, males in the highest dose group grew larger than those at a lower dose despite feeding at roughly the same rate. Females meanwhile grew larger than predicted based on their feeding rate in all dose groups. The highest dose groups in the 28 day study were terminated early due to unacceptable reductions to feeding rate, while reduced growth was not observed in the 2 year study so toxicodynamic predictions could not possibly be validated.

Poor model accuracy (<50% predictions accurate to within one standard deviation) only occurred when predicting the effects of 28 day dietary exposure. This is not entirely surprising as there are several factors making these data sets more challenging to model. In 28 day studies sample size was lowest, with only 5 individuals per sex and only four time points observed so naturally, individual and temporal variability would be expected to have a larger impact. Moreover, in the early weeks of dietary studies feeding rate can be highly variable between treatments and over time, as animals react behaviorally to a novel ingredient added to their diet. In all 3 data sets for which predictions were poor, body weight predictions in the highest dose groups were substantially lower than observed. The default value of k could play a role here. A low k value would stop growth when it could still occur, though this would be almost entirely compensated for in the fitting of other growth parameters. Inclusion of a reproduction buffer could also delay the need to metabolize structure,[26] but using up the buffer would still correspond to a reduction in overall body weight. Moreover, starvation was only predicted to occur in the early weeks of testing, around the onset of puberty,[54] so any buffer amassed would be almost negligible. It is quite possible that the dual stresses of reduced feeding and toxicity elicit compensatory physiological or behavioral responses not predicted by the model. Reduced body temperature has been documented as a response to starvation in rats[69] meanwhile chemical stress has been shown to induce reductions in body temperature and activity.[70] Such responses would likely correspond to a reduction in the maintenance rate J, and should be considered in future models of physical and chemical stress.

Individual variability was also evident in the toxicokinetics data. Individual responses varied with regard to the toxicant concentration reached in the blood and the speed with which it was absorbed and eliminated. The low sample size of three individuals per treatment meant that mean observations for each time point, to which the models were fitted, could be heavily influenced by variability on either axis. As an example, even if three individuals exhibited a single peak with very little variability in Cmax, the average data could show a double peak if just one individual peaked later than the other two. Moreover, both peaks would likely be lower than the Cmax of any individual. Even simple quantitative assessment of model fits, such as whether model predictions were in the observed range for toxicokinetic statistics such as Cmax, Tmax or area under the curve, were therefore problematic. Besides, which of these statistics would be more relevant for predicting the internal concentration resulting from prolonged dietary exposure is debatable and none of them they are used by TK–TD models to predict effects. Nevertheless, TK model predictions provide plausible estimates of internal concentration resulting from the recorded time varied ingestion rates, with parameter values fitted to the best available data. It is likely that the predictions are at least proportional to the true values and therefore form a credible basis for the fitting of TD parameters.

While it is not possible to separate variability in TD and growth parameters, future model iterations could incorporate stochasticity in growth parameters by utilizing the wealth of control data available from studies on R. norvegicus.[67] This would also be possible for TK parameters although low sample size would present a challenge.

Implications for DEB Theory

In this study, precise food consumption data were available, rather than simply food availability. This allowed the mean surface area specific feeding rate to be calculated on a weekly basis for each treatment in a study, before being scaled as a proportion of the maximum value in each data set. Consequently, the growth parameters were fitted to controls in each data set to account for variability in feeding rate within and between studies. Previously, it has been assumed that the scaled feeding rate, f, is equal to 1 when food is available ad libitum.[27] For certain purposes this is a reasonable and necessary assumption. When modeling growth in the field for example, detailed data are unlikely to be available and so feeding rate must be estimated as a fraction of ad libitum feeding in laboratory studies. However, this was not satisfactory in this study as variability in feeding rate over time and between treatments was an important driver of effects. Moreover, several treatment groups fed at a higher rate than controls, so it was important not to assume the maximum value of f as the default.

Growth ceases at the ultimate weight, W, as this is the point at which the maximum assimilation rate can only match the maintenance requirements of the organism.[26] Assuming f = 1 means that W is equal to the theoretical maximum weight, W, when food is available ad libitum (see Table ). However, a marked trend in the data was that, even with unrestricted access to food, feeding rate relative to surface area declined as the animals grew. This presents a clear issue in that it places an additional limit on the assimilation rate and ensures that the growth curve reaches a plateau before Wm is reached. In many cases, by week 13–14 of observation, the value of f was well below the maximum and so the theoretical Wm based on fitted parameter values was unrealistically high.

Recent DEBtox studies of toxic effects on rainbow trout have also sought to account for variability in f, however this has focused on differences between study groups rather than within treatment changes over time.[23,24] Much effort has gone into deriving standard DEB parameters for species of interest.[21,53,54] However, lifetime variation in f must be considered for these to be compatible with time varied feeding data rather than the constant food density. For laboratory strains of R.norvegicus this could certainly be addressed; due to their extensive use in regulatory testing there exists a vast database of growth and food consumption in control conditions. If the relationship between feeding rate and body size were described mathematically, then this could be utilized for scaling observed feeding rate such that the maximum feeding rate decreases with size and the resulting value of f remains roughly constant over the lifetime.

Although longer term data were available, predictions of toxic effects on body weight were not reported beyond around 12–14 weeks of dietary exposure for the two year studies. Beyond this point the intervals between body weight and food consumption observations increased. Since observed feeding rate was used as a model input, this reduction in data resolution would be expected to adversely affect predictions. There were also more fundamental reasons behind this cutoff. While sigmoidal curves such as the Von Bertalanffy growth model (to which DEBkiss simplifies[26]) can approximate the growth curve of rats, there are distinct stages where observed growth deviates from such a model. It has previously been reported that the Gompertz function–also sigmoid–matches data closely when fitted to the first 70–105 days of rat growth but that longer term predictions are problematic.[71] When fitted to the full two years of control data the DEBkiss curve also showed systematic errors. For both males and female rats, predicted body weight was lower than observed for roughly the first three months, higher than observed until around month 14 and then lower than observed for the final 10 months.

A possible explanation is that k (the allocation to soma and reproduction), does not remain constant throughout the rat’s lifespan. A logical suggestion is that a greater proportion of energy may be allocated to growth earlier in life with more energy used for sexual maturation during puberty. Indeed, it has been postulated that k may change in humans at puberty.[20] It is also likely that reserve dynamics become more important in adult rats, as continued ad libitum food availability allows animals to develop significant reserve stores after structural growth has ceased. The full DEB model which, unlike DEBkiss, models reserve as well as structure may therefore be better able to address this issue.

Another issue with long-term predictions is that there is insufficient knowledge regarding recovery after inhibition of growth in rats. By default, after a stressor is removed, modeled body weight may theoretically reach W if feeding rate is high enough. On the basis of our analysis, this assumption appeared sound for rats up to around 23 weeks of age. Model predictions agreed well with data in cases where 90 day studies included 4 weeks of recovery for the highest dose group. However, skeletal growth is known to stop in rats at around 26 weeks of age, the underlying processes are complex but appear to be related to age rather than size.[72] Logically, if growth has been suppressed up to a critical age then a full recovery, relative to controls, will not be possible as growth will cease before W, let alone Wm, can be reached. We hypothesize that this occurs because, beyond a certain age, available energy is allocated to processes other than growth, such as maturation and reproduction. This would correspond to a reduction in the parameter k, resulting in a reduced growth rate and, crucially, a lower W, for animals that had experienced stress.[20,26]

Realistic constraint of recovery is essential for long-term predictions to be of use. Otherwise, to match observed data, TD parameters must continually stress growth even when it is no longer possible, and so exaggerate the toxicity of a chemical. Recovery may be limited as a function of the (structural) weight reached by a critical age. However, determining rules by which to accurately decrease the value of k thereafter would likely require significant experimental work. Removing a given stressor at different time points may identify the age at which a full recovery becomes impossible. However, subsequent experiments would still be required to examine how potential for recovery is affected by the level of stress as well as the duration. This is discussed in greater detail in the SI. Such experimental work was beyond the scope of this study and, in any case, for regulatory purposes such long-term predictions are of limited relevance. As pesticides are not applied at constant rates for years at a time, a 21-day exposure scenario is used in ecological risk assessment to determine risk to animals in the field,[62] well below the 12–14 weeks of exposure modeled in this study.

Conclusions

This study shows that DEBtox modeling with DEBkiss can provide an effective and simple to use tool for predicting toxic effects on growth in rats. We show how time varying model inputs for feeding rate and pesticide ingestion rate can be calculated directly from data generated in standard chronic toxicity studies, providing additional insight into data by indicating to what extent body weight is impacted by feeding rate or toxic effects over time.

We also identify several difficulties which future models should aim to overcome. Individual variability presents a significant obstacle to assessing model accuracy. Our models simulate toxicokinetics, toxicodynamics, and growth, all of which may be subject to variability. In most cases, predictions were accurate to within one standard deviation of the observed mean and so provide useful estimates of the mean but not exact projections including variability.

Our findings support DEB theory as an effective basis for predictions of sublethal toxic effects in mammals. However, some issues became apparent regarding its compatibility with chronic toxicity data. Given the extensive use of R. norvegicus in laboratory testing and the resulting wealth of control data, these complications can be addressed. Further analysis of lifetime variation in feeding rate and energy allocation to the soma may improve model accuracy and realistically constrain recovery. Such adjustments would broaden the range of applications for which DEBtox may be used.