Physiologically-based pharmacokinetic and toxicokinetic models for estimating human exposure to five toxic elements through oral ingestion

Biological monitoring and physiologically-based pharmacokinetic (PBPK) modelling are useful complementary tools in quantifying human exposure to elements in the environment. In this work, we used PBPK models to determine the optimal time for collecting biological samples in a longitudinal study to determine if participants who consumed allotment produce had been exposed to arsenic, cadmium, chromium, nickel or lead. There are a number of PBPK models for these elements published in the literature, which vary in size, complexity and application, given the differences in physiochemical properties of the elements, organs involved in metabolism and exposure pathways affected. We selected PBPK models from the literature to simulate the oral ingestion pathway from consumption of allotment produce. Some models required modification by reducing or removing selected compartments whilst still maintaining their original predictability. The performance of the modified models was evaluated by comparing the predicted urinary and blood elemental levels with experimental data and other model simulations published in the literature. Overall, the model predictions were consistent with literature data (r > 0.7, p < 0.05), and were influential in predicting when samples should be collected. Our results demonstrate the use of mathematical modelling in informing and optimising the design of longitudinal studies.


Introduction
Arsenic (As), cadmium (Cd), chromium (Cr), nickel (Ni) and lead (Pb) are known to occur in urban horticultural sites (allotments) and gardens, as a result of pollution from traffic, general urban activities, industrial emissions and the actions of allotment tenants (Bechet et al., 2016;Kelly et al., 1996;Alloway, 2004). Fruits and vegetables grown in contaminated soil may lead to humans ingesting these elements through the consumption of home grown produce. According to the Environment Agency (2009), the consumption rates for allotment produce by adults in the United Kingdom (UK) range from 0.22 to 2.97 (g −1 fresh weight (fw) kg −1 body weight (bw) day −1 ), depending on the produce categories used in the Contaminated Land Exposure Assessment (CLEA) model. The CLEA model, a tool used in the UK to assess human exposure to soil contaminants, groups allotment produce into six categories (green, root, tuber, herbaceous, shrub and tree) (Environment Agency, 2009). The toxicological effects of these elements are well documented and As, Cd, Cr and Ni are Group 1 inorganic human carcinogens (IARC, 2012). Oral exposure to Pb can cause toxic effects on various body organs and tissues including kidney dysfunction, cardiovascular effects (e.g., increase in blood pressure), gastrointestinal effects (e.g., colic), haematological effects (e.g., anaemia), and brain damage (ATSDR, 2007).
Biological monitoring (biomonitoring) and physiologically-based pharmacokinetic (PBPK) modelling are useful tools in quantifying human exposure to elements in the environment; they are complementary in exposure studies (Clewell et al., 2008;Lyons et al., 2008). Human biomonitoring is the assessment of an individual's exposure to an element or compound through the measurement of a biomarker (e.g., blood, urine) which results from exposure to that element (Coelho et al., 2014). Meanwhile, physiologically-based pharmacokinetic (PBPK) modelling involves predicting the fate of elements in the body (Schmitt and Willmann, 2005), using physiological and biochemical information to describe and quantify the pharmacokinetic processes affecting the absorption, distribution, methylation (biotransformation) and excretion of an element (Wen et al., 1999). PBPK models are based on compartments (e.g., body organs, tissues) and the interconnections among them. The level of model detail relates to the compartments and elements (including associated chemical forms such as metabolites) that are tracked within the model (Krishnan et al., 2010). PBPK models provide a scientific basis for quantitatively estimating risk to human health (Yu, 1999).
In this paper, our aim is to use PBPK models to determine the optimal time for collecting biological samples to determine if participants who consumed allotment produce during a longitudinal biomonitoring study had been exposed to As, Cd, Cr, Ni or Pb. To do this, we sought to obtain applicable published PBPK and physiologically-based toxicokinetic (PBTK) models for these elements from the literature. Some of these models are quite complex in terms of their structure and formulation because of the differences in physiochemical properties of the elements, organs involved in metabolism and exposure pathways simulated by the authors. For example, PBPK models for As published by the El-Masri and Kenyon (2008), Liao et al. (2008), Yu (1999) and Mann et al. (1996) indicate that As undergoes biotransformation into different species within the body, with kidney and liver being the key organs for As methylation. The Cr models published by Finley et al. (1997), O'Flaherty et al. (2001 and Kirman et al. (2013) indicate that Cr species undergo oxidation and reduction in the body. No biological transformations are shown to occur in the models for Cd (Kjellstrom and Nordberg, 1978), Ni (Sunderman et al., 1989) and Pb (Fleming et al., 1999;White et al., 1998;O'Flaherty, 1993). Due to the structural complexity of some of the models, to construct models suitable for our purpose, some of the published models required modification by reducing or removing some compartments whilst maintaining their predictive ability.
The objectives of this part of our research were: (i) to review published PBPK and PBTK models for these elements and select models for use in our study; (ii) to make the necessary modifications to the models to fit our purpose (i.e., to use the models to determine the optimal time for collecting biological samples); (iii) to evaluate their predictive performance by comparing the reduced model predictions with data published in the literature; and (iv) to use the modified models to predict optimal times (following oral exposure to the elements) for collecting biological samples during our biomonitoring study, in order to maximise the detection of the elements in urine and blood samples.
Based on the consumption of allotment produce, we set out to investigate human exposure to As, Cd, Cr, Ni and Pb from soil with low levels of these elements. We identified 30 allotment sites across central Scotland, UK with varying levels of toxic element contamination. Allotment users (36 adults) were recruited to participate in the study. Informed consent was obtained from each participant and ethical approval was granted by the University of Reading Research Ethics Committee (UREC), approval reference UREC 15/21.

Model selection
We reviewed previously published PBPK/PBTK models for As, Cd, Cr, Ni and Pb in humans, and selected suitable ones based on whether the model was reproducible, relevant to adult humans, included a description of the oral ingestion pathway, and the most recent/most used model. The mathematical equations for each modified model are given in the Supplementary material along with the schematic details of the original models obtained from the literature (Figs. S1-S5). The physiological and chemical-specific parameters were obtained from the literature and are detailed in Tables S1 to S5 of the Supplementary material. In the following sections, we detail the models selected for each toxic element and the modifications made to the original published models.

Arsenic
We identified a number of human PBPK models for inorganic As published in the literature (Mann et al., 1996;El-Masri and Kenyon, 2008;Liao et al., 2008;Yu, 1999). These models are largely similar in structure, and they account for oxidation of trivalent inorganic As (As (III)) to its pentavalent form (As(V)), reduction of As(V) to As(III), and methylation of As(III) to monomethylarsenic acid (MMA) and dimethylarsinic acid (DMA) in the body. The model by Mann et al. (1996) included both oral and inhalation exposure pathways, whereas other models only included the ingestion pathway. Ingestion is the primary route of exposure studied, because exposure to As from allotment land use occurs mainly through oral intake (CL:AIRE, 2014). In addition, a pilot study we carried out identified air As concentrations at the allotments were not significantly elevated enough above background air concentrations to warrant inclusion in our model. The model proposed by Liao et al. (2008) for children was not evaluated using experimental data, whereas the predictive performance of the other models were tested using data from human studies included in the same publications. We adopted the model published by El-Masri and Kenyon (2008) because it is the most recent model (hence the most informed parameter wise) out of those considered. The El-Masri & Kenyon model comprises 9 compartments: gastrointestinal (GI) tract, liver, kidney, blood, muscles, brain, skin, heart and lung. The lung compartment was included in the model because it receives total blood flow, thus mathematically accounting for As reduction that may occur in other body tissues (El-Masri and Kenyon, 2008). We made the following modifications to the model: (i) we considered oxidation and reduction in the lung, liver and kidney only, and assumed that oxidation and reduction occurs in all perfused tissues as previously reported by Mann et al. (1996) and Yu (1999); (ii) we ignored the oxidation and reduction reactions between MMA(III) and MMA(V), DMA(III) and DMA(V) and treated MMA and DMA as single species because in our laboratory analysis we tested for total inorganic As (sum of inorganic As species); and (iii) included biliary excretion of As from the liver as reported by Yu (1999) and Liao et al. (2008). Fig. 1 shows the schematic representation of our resulting modified PBPK model for As.

Cadmium
There are a number of PBTK models for Cd published in the literature (Nordberg and Kjellstrom, 1979;Choudhury et al., 2001;Amzal et al., 2009;Ju et al., 2012;Fransson et al., 2014) based on the original PBTK model for Cd published by Kjellstrom and Nordberg (1978) (KN model). The KN model consists of eight compartments, describing Cd uptake from the GI tract and the lungs, distribution of absorbed Cd to three blood compartments, liver, kidney and 'other tissues', and Cd elimination through urine and faeces. The KN model was tested using data from studies involving Cd exposure through inhalation (Nordberg and Kjellstrom, 1979). The model has been used in other studies (Ruiz et al., 2010;Ju et al., 2012;Fransson et al., 2014) to simulate Cd exposure through oral ingestion. It was modified by Ruiz et al. (2010) and used to accurately predict urinary Cd concentrations following oral ingestion using data from the National Health and Nutrition Examination Survey (NHANES). A study by Choudhury et al. (2001) also used a modified version of the KN model to predict urinary Cd concentrations consistent with the NHANES data.
We adopted the KN model in this study because it is the most used model and it is also the basis for other published Cd PBTK models.
Given that our study involves exposure through oral ingestion of food, we excluded the inhalation pathway from the model. In addition, the direct transfer of unabsorbed Cd to faeces was added to the modified model to account for unabsorbed Cd. The schematic representation of the modified PBTK model for Cd is given in Fig. 2.

Chromium
Two key models for Cr were identified in the literature. A physiologically-based model for the ingestion of Cr(III) and Cr(VI) by humans was developed by O'Flaherty et al. (2001). The compartments included in the model are the GI tract, blood (plasma and red blood cells), liver, kidney, bone (trabecular and cortical), well-and poorly-perfused tissues with pulmonary absorption of inhaled Cr and excretion pathways via faeces and urine. The model was calibrated using blood and urine Cr concentration data from controlled studies in which human volunteers drank solutions containing Cr(III) and Cr(VI) (O'Flaherty et al., 2001). Another PBPK model for humans orally exposed to Cr was developed by Kirman et al. (2013). This model includes compartments for the GI tract (stomach and intestines), blood (systemic and portal plasma and red blood cells), liver, kidney, bone and 'other tissues'. Model performance was evaluated using toxicokinetic data for Cr in human tissues and excreta obtained from the literature (Kirman et al., 2013). The model was observed to provide a good description of Cr toxicokinetics in humans. One key difference between these two models is that Kirman et al. (2013) included detailed Cr toxicokinetics in the GI tract compartments, while O'Flaherty et al. (2001) modelled the GI tract as a single compartment. We adopted the model by Kirman et al. (2013) since it is the most recent model, and modified it by ignoring the detailed competing toxicokinetic processes of Cr in the stomach and intestines, because there are published absorption rates for Cr in the GI tract (Kirman et al., 2013;Sasso and Schlosser, 2015). Therefore, the 'stomach', 'small intestines' and 'large intestines' were lumped into a single compartment (GI tract). Fig. 3 shows the structure of our modified PBPK model for Cr.

Nickel
There are limited PBPK models for Ni compounds in the literature relevant to humans. Models describing the deposition, retention and (A) Oral absorption of As(III), As(V), MMA, DMA was accounted for in the GI tract. eF is faecal excretion rate of As species (day −1 ), eB is biliary excretion rate of As species (day −1 ), eU is urinary excretion rate of As species (day −1 ). K12 to K92 refer to the calculated transfer rates (day −1 ) of As species between compartments. (B) Shows the oxidation/reduction of inorganic arsenic in all tissues and methylation of As(III) to MMA and DMA in kidney and liver (Liao et al., 2008). K ox and K red are metabolic constants (day −1 ) for oxidation and reduction, respectively, V max (μmol day −1 ) and K m (μmol L −1 ) are metabolic constants for methylation of As (III), and i refers to the respective compartment (liver, kidney). Kjellström & Nordberg (1978). Blood1 to Blood3 refer to 'plasma other', 'red blood cells' and 'plasma metallothionein', respectively; C5 to C19 and CX refer to the parameters describing the transfer of Cd between compartments as originally defined by Kjellström & Nordberg (1978). clearance of inhaled Ni (in the lung) have been developed (e.g., Hsieh et al., 1999;Yu et al., 2001). However, given that inhalation route is not the dominant pathway for allotment land use, these models were not reviewed. The only PBTK model in the literature for human exposure to Ni through the oral pathway is that of Sunderman et al. (1989), following controlled oral exposure to Ni ions. They developed the model based on two experiments in which Ni levels in serum, urinary and faecal excretions were monitored after 8 human subjects were given an oral dose of Ni (as NiSO 4 ) in either water (experiment 1) or food (experiment 2). The model comprises four compartments (gut, serum, urine and tissues) with the parameters based on model-data fitting to the two experiments. This allowed rate constants for alimentary Ni absorption from the gut, Ni transfer rate constants from serum to tissues and urine, and from tissues to serum to be determined. Given that our study involves exposure through oral ingestion of food, the model parameters provided by Sunderman et al. (1989) from experiment 2 were considered more relevant to our study than parameters from experiment 1. However, Sunderman et al. (1989) did not determine the rate transfer from tissues to serum in experiment 2, which they indicated was due to the small mass of Ni absorbed from the gut into subsequent compartments. Therefore, we used the Ni rate transfer from tissues to serum from experiment 1 in our simulations. The unabsorbed fraction of Ni was accounted for by adding a 'faeces' compartment to the model. A schematic representation of the modified PBTK model for Ni is given in Fig. 4.

Lead
There are a number of Pb kinetic models for humans in the literature. The most cited in the literature reviewed are the Leggett model (Leggett, 1993), the Integrated Exposure Uptake Biokinetic (IEUBK) model (White et al., 1998) and the O'Flaherty model (Flaherty, 1991(Flaherty, , 1993;Fleming et al., 1999). The IEUBK model deals with Pb kinetics in children only whilst our study involved adults. The Leggett model does not account for physiological factors in detail because it utilises 'transfer rates' that are age-specific and regards those older than 25 years old as one age-class, which does not account for physiological variabilities that can be associated with varying body weights (Bailey et al., 2004;Prakash et al., 2013). Therefore, both the IEUBK and the Leggett models were not considered suitable for this study. The O'Flaherty model was chosen for our work because it is physiologically-based and is only applicable to adults. Full details of the construction of this model are provided elsewhere (O'Flaherty, 1991;O'Flaherty, 1993;Fleming et al., 1999). Briefly, the model includes 9 body compartments, namely: the GI tract, blood, liver, kidney, bone (cortical and trabecular), other tissues (well-and poorly-perfused) and the lungs. This model has been evaluated against data from human subjects exposed to Pb through oral and inhalation pathways (O'Flaherty, 1993). We excluded the inhalation component of this model because from our preliminary investigation, air Pb concentrations at allotments were not elevated enough above background air concentrations. In addition, we simplified the model by neglecting the detailed Pb kinetics in the bone associated with bone growth, since human skeletal growth is considered to be complete by the age of 25 (O'Flaherty, 1993) and all of our study participants were aged 30 years and above. Fig. 5 provides a schematic representation of the modified PBPK model for Pb.

Model solutions and analysis
Each model was solved using the SimBiology application in MATLAB (version R2016a, MathWorks ® ) with the mathematical equations and parameters detailed in the Supplementary Material. Given the variation in parameter values, the stiff solver ode15s was chosen to provide accurate numerical simulations with the relative tolerance set at the default value (0.001).

Evaluation of model performance
We sought to evaluate the predictive performance of the modified models by comparing model simulations with data presented in the literature. Firstly, the relationship (goodness of fit) between results from each mathematical model and the corresponding time course of literature data were examined using linear regression (Pearson's correlation test) using a statistical level of significance of p < 0.05.
Secondly, the predictive accuracy of the models was also assessed using the root mean square error (RMSE) as calculated by Ju et al. (2012) and Walther and Moore (2005) Here, C s and C lit refer to the simulated data and literature data values, respectively at time point i for a total of n number of data points. We expressed maximum C s as a percentage of C lit to determine the magnitude of over-or under-prediction of the models. RMSE explicitly  Kirman et al. (2013). RBC refers to red blood cells. All compartments contain Cr(III) and Cr(VI). Reduction of Cr(VI) to Cr(III) occurs in the GI tract, blood, liver, kidney, and other tissues. Cr is excreted in urine mainly as Cr(III) due to rapid reduction of Cr(VI) to Cr(III) in the body (O'Flaherty et al., 2001).  Sunderman et al. (1989). K1, Kf, K12 and K21 refer to the transfer rates of Ni between compartments and eU is the rate constant for urinary elimination. The absorbed fraction of Ni dose in the gut (A gut ) was 0.011, as determined by Sunderman et al. (1989). The rate constant for faecal excretion of unabsorbed Ni in dose was expressed as ρK1, and ρ was calculated as (1-A gut )/ A gut .
shows how much the model predictions deviate from the corresponding literature data; hence RMSE gives an absolute measure of fit. Statistical analyses and processing of data were carried out using R statistical software version 3.3.0 (R Core Team, 2016) and Microsoft Office Excel 2010.

Sensitivity analysis
Parametric sensitivity analysis was performed using sensitivity coefficients (SC) to determine which parameters were more sensitive to change. Values of SC were calculated using the following expression (Choudhury et al., 2001;Evans and Andersen, 2000).
Here, δm is the change in model output (m) resulting from the change (δp) in an input parameter value (p). An SC of zero implies that there is no change in model output regardless of the parameter value used. Positive and negative SC values indicate an increase in model output with a given increase in parameter value, and a decrease in model output with a given increase in parameter value, respectively. High SC values indicate a model is highly sensitivity to a given parameter value (Choudhury et al., 2001). Sensitivity analysis was performed in MA-TLAB (using full normalisation option); input values were varied by ± 50%.

Results and discussion
3.1. Evaluation of model performance 3.1.1. Arsenic The modified As model was used to simulate cumulative urinary As metabolites based on single and multiple oral doses of As in the form of As(V) and As(III), as reported in the literature (Mann et al., 1996;El-Masri and Kenyon, 2008;Buchet et al., 1981). Selected results from the simulations are plotted in Fig. 6. It has been reported that inorganic As undergoes oxidation and reduction in body tissues and methylation to MMA and DMA in the liver and kidney (Liao et al., 2008;El-Masri and Kenyon, 2008). Generally, studies have found DMA to be the major metabolite in urine following exposure to inorganic As (Yu, 1999;Buchet et al., 1981;Hughes, 2006;Hwang et al., 2002), as biotransformation of inorganic As to MMA and DMA occurs rapidly making DMA the dominant metabolic species after approximately 2 days. Our model simulations (Fig. 6) agree with this observation, with the modified As model able to reproduce the literature data well (r > 0.9 for total As in urine). The RMSE values for total inorganic As were 4.7 × 10 −6 and 0.4 for Fig. 6(A) and (B), respectively. Further evaluation of the model was carried out using ingestion of 6.67 μmol (500 μg) As (Buchet et al., 1981) and 1.33 μmol (100 μg) As (El-Masri and Kenyon, 2008). Model predictions of cumulative urinary As (total) were within 8% of the published data in the literature. Analytical experimental procedures for determining As species in urine are often complex and costly. Therefore, the PBPK model can be used to estimate internal dose and urinary concentration of speciated As, thus providing a proxy for the analysis of speciated As in urine.
The SC for all model input parameters were generated to determine the influence of parameter variation to the model output (urinary As). The modified As model is most sensitive to urinary excretion of As (eU) with an SC value > 0.9. The model also showed moderate sensitivity to As methylation constants in the kidney (V max and K m ) and the oral absorption constant (K a ), with SC values varying between 0.5 and 0.7. Therefore, parameters eU, V max , K m and K a for the modified As model require more attention during simulations and model fitting with measured data. Our findings compare well with published literature. For example, the sensitivity analysis performed by Yu (1999) for their PBPK model for inorganic As determined that kidney V max coefficients and urinary elimination constants were key input parameters that affected the model output.

Cadmium
We compared data published by Ju et al. (2012) regarding Cd ingestion with the concentrations of Cd in blood and urine predictd by our modified model. Simulations were performed with all parameters fixed to the values given in the original KN model and only the bioaccessibility values were varied to correspond to values used by Ju et al. (2012). The creatinine values reported by Ju et al. (2012) were used in the conversion of the predicted urinary Cd concentrations. The predicted results and literature data are presented in Fig. 7, which shows that predicted concentrations are in good agreement with literature data, with low RMSE values between 0.01 and 0.02 μg g −1 creatinine (r > 0.9).
We also used the dietary intake of Cd reported by Berglund et al. (1994) in their study investigating intestinal absorption of dietary Cd in women subjects (20-50 years of age) to predict Cd concentrations in urine and blood following daily dietary exposure. Our simulations mimicked five long-term exposure scenarios lasting between 10 and 50 years which are compared in Table 1, against the results Cd measured in urine and blood by Berglund et al. (1994). Although no corresponding data (with respect to duration after exposure) was given by the authors (this would allow direct comparison with our simulated results) the model predictions fit within the reported data ranges.
Sensitivity analysis conducted for the modified Cd model input parameters indicated that rates for oral absorption (C5), systemic absorption (C6 and 1-C7) and urinary excretion (C7) of Cd were the input WP is well-perfused tissues, PP is poorly-perfused tissues, IR gi is the oral intake rate of Pb (mg day −1 ), A gi is the Pb absorption coefficient from the GI tract (unitless), the unabsorbed fraction is represented by (1-A gi ), eB is the biliary excretion rate of Pb (day −1 ), eU is the urinary excretion rate of Pb (day −1 ), K12 to K62 refer to the calculated Pb transfer rates (day −1 ) between compartments.

Chromium
The predictive performance of the Cr modified model was tested using data published by Paustenbach et al. (1996), Kirman et al. (2013) and Kerger et al. (1996). The results are plotted in Fig. 8.
In Fig. 8(A), the model predicts a faster urinary excretion of Cr than the measured data in the literature. However, the cumulative urinary Cr is consistent with literature data. The predicted results show a reasonable fit to the literature data following oral ingestion of Cr(III) and Cr(VI) in single and multiple doses.
It has been observed that human absorption of Cr(VI) in the GI tract can vary between individuals and also in the same individual at different times, with suggestions that physiological fluids in the GI tract  (Mann et al., 1996), (B) following repeated daily oral ingestion of 1.67 μmol As for 5 days (Buchet et al., 1981). (iAs is the sum of As(III) and As(V); m, model; d, experimental data). Fig. 7. Comparison of the predicted Cd levels in urine and blood with data reported by Ju et al. (2012). M and F refer to male and female non-smoking human subjects, respectively. Six scenarios of Cd bioaccessibility as used by Ju et al. (2012) were simulated, I-V represent Cd bioaccessibility values (dimensionless) of 0.021, 0.032, 0.044, 0.057 and 0.094, respectively. VI represents the scenario of using the original fraction of Cd absorbed to gastrointestinal tract and systematic circulation (0.048) adopted by Kjellström and Nordberg (1978). such as gastric juice and diet constituents (for instance orange juice) can lead to poor intestinal absorption of Cr(VI), because of their capacities to reduce Cr(VI) to Cr(III) (De Flora et al., 1997;Kerger et al., 1996;Sasso and Schlosser, 2015). Intra-individual variability in Cr(VI) absorption due to differences in this reduction capacity has been noted in studies involving human subjects (Finley et al., 1997;Paustenbach et al., 1996;O'Flaherty et al., 2001). Absorption values of 0.25 day −1 for Cr(III) and 2.5 day −1 for Cr(VI) were specified by O'Flaherty et al. (2001). In addition, we used the mean rate constants for absorption of Cr(III) (4.6 × 10 −6 L h-cm −1 ) and Cr(VI) (3.2 × 10 −4 L h-cm −1 ) in the small intestines given by Kirman et al. (2013) to estimate corresponding intestinal absorption rates of 0.05 day −1 for Cr(III) and 3.5 day −1 for Cr (VI). This suggests that oral absorption values for Cr should be carefully Table 1 Modified model predicted urine and blood Cd concentrations and data reported by Berglund et al. (1994).

Mixed diet
High fibre All concentrations are expressed in μg L −1 . The median daily dietary intakes of Cd reported by Berglund et al. (1994) (10 μg day −1 for mixed diet and 13 μg day −1 for high fibre diet) were used in the simulations. Predicted Cd mass in urine were converted into concentrations using the daily urine volumes reported by Berglund et al. (1994). Likewise, Cd loads in blood were converted to concentrations using blood volume of 5.2 L (Ju et al., 2012). a Berglund et al. (1994). Fig. 8. Comparing model predicted cumulative urinary excretion of Cr and the measured urinary Cr (mean values) in human subjects. (A) for humans exposed to a single dose of 5 mg of Cr(III) as reported by Kerger et al. (1996), (B) for humans exposed to 0.4 mg of Cr(III) per day for 3 days (Kirman et al., 2013), (C) for humans exposed to a single dose of 5 mg of Cr(VI) as reported by Kerger et al. (1996), and (D) for a human volunteer exposed to 4 mg Cr(VI) per day for 17 days . selected when fitting a model using experimental data. In our simulations we used a range of GI tract absorption rates, 0.05-0.25 day −1 for Cr(III) and 1-2.5 day −1 for Cr(VI) to fit model predictions with the literature data (r > 0.8, RMSE ≤0.01 mg). Sensitivity analysis identified the oral absorption rate as the most sensitive parameter (SC > 0.9) in the modified Cr model. The second most sensitive parameter in the Cr modified model was the urinary elimination rate of Cr, which had an SC value of 0.4. Our findings are consistent with those published by Kirman et al. (2013) who also identified a number of sensitive parameters to their Cr model, including parameters associated with Cr absorption in the GI tract and urinary excretion.

Nickel
The predictive performance of the Ni modified model was tested using data obtained from Sunderman et al. (1989) and Nielsen et al. (1999). The model results are presented in Fig. 9, which indicates that the predicted urinary excretion of Ni match closely with data from Sunderman et al. (1989) thus showing the capability of the model to reproduce literature data (r > 0.9 with low RMSE values ranging from 0.02 to 0.9 μg). The mass fraction of Ni dose absorbed from the gut in the experiment by Sunderman et al. (1989) was 0.7 ± 0.4% for Ni dose ingested in food. In a similar study involving controlled ingestion of Ni dose in food, Nielsen et al. (1999) reported a median value of Ni oral absorption of 2.95 ± 1.32%, which is notably higher than oral absorption values reported by Sunderman et al. (1989). This may explain why measurements of urinary Ni published by Nielsen et al. (1999) were higher than those from the Sunderman et al. (1989) study (Fig. 9).
A range of oral absorption values of Ni in the gut were also noted in the literature. Studies by Horak andSunderman (1973), McNeely et al. (1972) and McNeely et al. (1971) reported Ni oral absorption values between 1 and 1.6%. Overall, these oral absorption values indicate that faecal excretion is a major route for elimination of Ni from the human body.
The most sensitive parameters for the Ni modified model were related to oral absorption of Ni (A gut and K1), which had SC values > 0.9. The model also showed some sensitive to the urinary elimination rate of Ni which had an SC value of 0.7.

Lead
Evaluation of the modified Pb model was carried out using experimental data of Rabinowitz et al. (1976), who studied the steady-state kinetics of Pb in five healthy men (subjects A to E) with stable isotope tracers. Since the transport of Pb throughout the body is governed by its concentration in the plasma (O'Flaherty, 1993;Fleming et al., 1999;Leggett 1993), Pb concentration in whole blood was calculated from the model-predicted plasma Pb concentrations using the expression (MacMillan et al., 2015). . ( Here CB is Pb concentration in whole blood, HCT is the haematocrit fraction of whole blood (0.45), CPLASMA is Pb concentration in blood plasma, G is the ratio of unbound erythrocyte Pb concentration to plasma Pb concentration (1.2), BIND is the Pb binding capacity of erythrocytes (0.437 mg Pb L −1 cell), and KBIND is the binding constant of erythrocytes (3.72 × 10 −4 mg Pb L −1 cell). This allowed for comparison to be made between model-predicted plasma concentrations and Pb measurements in whole blood. Simulated results for the five subjects are presented in Fig. 10, which shows maximum Pb concentrations similar to the observations made by Rabinowitz et al. (1976). The predicted Pb concentrations in blood were highly correlated to the literature data (r > 0.9 for subjects A, B, D, E, and r = 0.7 for subject C). Notably, there were fewer data points for subject C, which could be the reason for the reduced correlation. In addition, the literature data for subject C shows a constant blood concentration beyond day 2, while the model simulated declining concentrations (Fig. 10). Individual variabilities (such as differences in absorption rates and stomach clearance rates) cannot be ruled out as potential contributing factors to this variation. Despite the reduced correlation for subject C, the predicted maximum Pb in blood (which is relevant in understanding exposure) was consisted with experimental data. Predicted maximum Pb concentrations were within 10% of the literature data, apart from subject E which recorded 38% above the literature data. Similar observations were made by O'Flaherty (1993) and Morisawa et al. (2001) who tested their Pb models using the same data by Rabinowitz et al. (1976); although data for subject C was not used in either of these publications. Parametric sensitivity analysis indicated that our modified model for Pb was most sensitive to the oral absorption coefficient (A gi ) with an SC value > 0.9. The model also indicated some sensitivity to Pb transfer rate from liver to blood (K12) with an SC value of 0.6, and to Pb transfer from blood to liver (K21) with an SC value of 0.4.

Application of the models in planning our biomonitoring study
The modified models for As, Cd, Cr, Ni and Pb were used to predict optimal times for the collection of biological samples during our biomonitoring study. The aim was to identify the best time points, following oral ingestion of allotment produce, at which to collect urine and blood samples (the biomarkers of exposure) given their expected low element concentrations. Model results would thus maximise the potential of detecting these elements in biological samples.
The simulated element doses (mg day −1 ) were calculated based on the following considerations: a) Average produce concentrations recorded during our preliminary investigation. For Cd, we assumed a value equivalent to the limit of detection (0.03 mg kg −1 ) because all samples recorded Cd concentrations below the limit of detection; b) A produce consumption rate of 3.34 (g −1 fw kg −1 bw day −1 ) for a 'high end' consumer, derived from the data used in the CLEA model (as illustrated in Table 2). This consumption rate relates to CLEA age classes 17 and 18, which correspond to adults aged between 16 and 65 (age class 17) and 65-75 (age class 18), respectively. Our study participants are aged over 30 years old. c) An average adult body weight of 70 kg (Brown et al., 1997); and d) Short durations of exposure (1, 3 and 7 days) were simulated to mimic hypothetical minimal (worst-case) exposure scenarios.
The predicted optimal sample collections times are presented in Table 3. Model predictions indicate that at very low levels of exposure, detection of Cd and Cr in urine might not be achievable, unless the participants have higher intake than simulated. We relied on a detection limit of 0.01 μg L −1 using inductively coupled plasma mass spectrometry (ICP-MS). Higher intake would increase the 'window' of optimal sampling times. On the contrary, detection of urinary As (total) and Ni, and blood Pb would not prove difficult. These optimal times informed the sampling frequencies adopted in our biomonitoring study.

Conclusions
Overall, our five modified models for oral ingestion of As, Cd, Cr, Ni and Pb have shown good agreement with literature data. As such, they can be considered in the analysis of human exposure to the selected elements when only the oral ingestion pathway requires investigation, thus simplifying the overall model structures, in some cases to a large degree. Model parameters can also be adjusted to reflect the actual exposure characteristics being investigated.
Our work has demonstrated that mathematical modelling can be a useful tool in planning the collection of human studies during, in this case, a biomonitoring study. It indicates the period during which the likelihood of sample detection is highest, thus informing the best periods to collect samplesa time consuming and costly task in itself. We envisage the use of this technique could be more widely employed in environmental studies and provides a relatively inexpensive means for optimising sample collection in conjunction with the known biology and chemistry of the respective compounds. Fig. 10. Comparison of our modified PBPK model predicted Pb concentration in blood with measured Pb in whole blood from the Rabinowitz et al. (1976) study. Duration of study refers to time after the beginning of controlled ingestion of Pb isotope. Ingestion periods for Subjects A (204 μg day −1 ), B (185 μg day −1 ) and D (105 μg day −1 ) were 104, 124 and 82 days, respectively. Subject C ingested 68 μg day −1 for 1 day, while subject E ingested 99 μg day −1 for the first 8 days and again from days 42-52. We assumed individuals consuming all CLEA produce categories. c Vegetables (Green, Root, Tuber) and Fruits (Herbaceous, Shrub, Tree). d Units (g fw kg −1 bw day −1 ), fw (fresh weight, produce), bw (body weight). e Allotment-related consumption rate is the product of 'consumption rate' and 'homegrown fraction'.