Pooled human microsomes were purchased from Tebu-bio ¹ (Peterborough, UK). The microsomes were prepared from a pool of 50 liver samples, mixed gender (20 mg protein ml^⁻1). DEHA and the simple monoester mono-2-ethylhexyl adipate (MEHA) (purity >97.4%) were provided by BASF SE. All chemicals used were of analytical grade or higher; B-nicotinamide adenine dinucleotide phosphate (NADP), purity 97%, Glucose-6-phosphate, 98%-100%, Magnesium chloride, ACS reagent >99%, and Glucose-6-phosphate dehydrogenase (type V from baker’s yeast) were obtained from Sigma Aldrich. Potassium dihydrogen phosphate, analytical grade, and Di-potassium hydrogen phosphate, analytical grade, were obtained from Fisher Scientific.

Samples were analysed by liquid chromatography (Shimadzu Prominence) with tandem mass spectrometry detection (AB Sciex API 3200) using electrospray ionisation. Ion optics, temperatures and gas flows were optimised on our individual system. All analyses used a Synergi Hydro-RP column (150 × 2mm; 4µ; Phenomenex) in conjunction with a methanol:20 mM ammonium acetate (0.1% acetic acid) gradient. Sample injection volume was 2 µl.

Consistent with previous studies with diisononyl-cyclohexane-1, 2-dicarboxylate (Hexamoll^® DINCH) (McNally et al., 2019) and di-(2-propylheptyl) phthalate (DPHP) (McNally et al., 2021; McNally and Loizou 2023) the very high lipophilicity of DEHA resulted in the formation of an insoluble film on the surface of the reaction medium which precluded the measurement of in vitro clearance. Therefore, only the measurement of in vitro clearance of MEHA was possible (Figure 1). In vitro incubations, the determination of in vitro half-life, in vitro intrinsic clearance and the calculation of in vivo clearance were identical to previous studies and are described therein (Pacifici et al., 1988; Soars et al., 2002; Howgate et al., 2006; Barter et al., 2007; McNally et al., 2019; McNally et al., 2021; McNally and Loizou 2023).

Tissue:blood PCs and unbound fractions in plasma were calculated from the log of the octanol to water partition coefficient (Log P_ow) as described in McNally et al. (2019) and McNally et al. (2021). Briefly, the Log P_ow for DEHA and MEHA were calculated using the ACDLogP algorithm (Mannhold et al., 2009) implemented in the ACD/ChemSketch 2019.1.0 ² software (Table 1). The Log P_ows were input into two tissue-composition-based algorithms for the calculation of tissue:blood PCs. The method of Poulin and Haddad (2012), which was developed for the prediction of the tissue distribution of highly lipophilic compounds, defined as chemicals with a Log P_ow > 5.8, was used for DEHA (Table 1). The method of Schmitt (2008), which was developed to predict the tissue distribution of chemicals with Log P_ow < 5.2, was used to predict the PCs of the monoester, MEHA (Table 1). The algorithm of Poulin and Haddad (2012) was implemented as a Microsoft^® Excel Add-in whereas a modified version of the algorithm of Schmitt (2008) was available within the httk: R Package for High-Throughput Toxicokinetics (Pearce et al., 2017). Where the tissue-composition-based algorithms did not provide a tissue:blood PC for a particular compartment, the value from a surrogate organ or tissue with similar blood perfusion rate (i.e., could be lumped within the rapidly or slowly perfused compartments) was assumed. These are presented in italicised text with the surrogate organ or tissue in brackets Table 1.

The fraction unbound (fu) was calculated from log ((1-fu)/fu) with the following equation: f u = 1 10 x + 1 (1)

Where, x = 0.4485 l o g P − 0.4782

When x is the equation for the prediction of fu for a chemical with a predominantly uncharged state at pH 7.4 (Lobell and Sivarajah 2003) (Table 1).

The biological monitoring (BM) data from the human volunteer study of Nehring et al. (2020) were simulated in this investigation. Briefly, four healthy volunteers (2 females, 2 males; aged between 24 and 34 years; bodyweight between 59 and 91 kg), each received a single oral dose of approximately 10 mg of DEHA, dissolved in 1 ml ethanol and diluted with water, administered in a chocolate coated waffle cup. The resulting individual doses of DEHA ranged from 107 to 164 μg/kg bw. (Table 2). The volunteers (A-D) did not have any known occupational exposure to DEHA. The volunteers donated 22, 26, 21, and 20 individual urine samples over a 48-h period with total urine volumes of 2,800, 5,544, 6,576 and 3,606 ml, respectively.

The concentrations of mono-2-ethyl-5-hydroxyhexyl adipate (5OH-MEHA), mono-2-ethyl-5-oxohexyl adipate (5oxo-MEHA), mono-5-carboxy-2-ethylpentyl adipate (5cx-MEPA), and the non-specific hydrolysis product adipic acid (AA) were extracted from the dataset of (Nehring et al., 2020). Concentrations of AA were only available for the first 24 h of study data. The rates of deposition of these metabolites into the urine (mg/h) were calculated based on the concentrations (mg/l), the volume of the urine void (l) and the time between successive voiding events. This rate represents an average rate of deposition since the previous urination event and renders the trends in urine data more clearly compared to concentrations expressed in (mg/l) or concentrations expressed relative to creatinine (Nehring et al., 2020). The derived rate was associated with the mid-point between the two voiding events.

Whilst it is possible to directly estimate the fractions of ingested DEHA eliminated as specific metabolites (5OH-MEHA, 5oxo-MEHA and 5cx-MEPA) from the study of Nehring et al. (2020), for parameterising the PBK model described below it was necessary to 1) estimate the fraction of ingested DEHA that was absorbed and 2) estimate the fractions of absorbed DEHA ultimately eliminated as 5OH-MEHA, 5oxo-MEHA and 5cx-MEPA—these fractions are greater than the fractions of ingested DEHA eliminated in urine as 5OH-MEHA, 5oxo-MEHA and 5cx-MEPA so long as there is incomplete absorption of DEHA.

The fraction of ingested DEHA eliminated as 5OH-MEHA, 5oxo-MEHA and 5cx-MEPA were calculated using Equation 2, where M denotes the molar mass of DEHA and the respective metabolites, D the administered dose (mg) of DEHA and m i denotes the mass of metabolite (mg) deposited in the urine void at time point i . A similar calculation was made for AA but with a correction for background sources. The rate of deposition of AA (mg/h) from background sources was estimated from the final few measurements in the first 24 h of the study (noting that AA was only measured in urine voids within this period). The time points, where the deposition rate of 5cx-MEPA into urine was less than 1% of its maximum value, were used in this calculation, with between 1 and 3 measurement time points used for the four volunteers. The average deposition rate of AA (mg/h) over the selected time points was calculated. A total background deposition of AA (mg) over the 24-h period of the measurements was subsequently calculated. F U E = 1 D ∑ m i M   D E H A M   M e t a b o l i t e (2)

A lower bound for the fraction of DEHA absorbed by each volunteer was estimated as the sum of F U E for 5OH-MEHA, 5oxo-MEHA, 5cx-MEPA and AA. The fractions of absorbed DEHA eliminated as 5OH-MEHA, 5oxo-MEHA, 5cx-MEPA and AA were calculated as the respective F U E divided by the total. These values are provided in Table 2.

A previous human volunteer study described in Loftus et al. (1993), wherein volunteers were administered deuterium labelled DEHA, estimated that approximately 12% of DEHA was eliminated as non-specific metabolites, principally 2-ethylhexanoic acid (EHA), that were not measured in Nehring et al. (2020). Takahashi et al. (1981) reported that a large fraction of DEHA administered to rats was ultimately exhaled as carbon dioxide. Both Takahashi et al. (1991) and Loftus et al. (1993) reported little faecal excretion. Based on these two studies complete absorption of DEHA was assumed for an upper bound case. Lower bounds on the fractions of absorbed DEHA eliminated as 5OH-MEHA, 5oxo-MEHA, 5cx-MEPA and AA were taken as the respective F U E calculated through (2). These values are also provided in Table 2.

An existing human PBK model for DPHP (McNally et al., 2021) was adapted and simplified for studying the absorption, distribution, metabolism, and elimination of DEHA following single oral doses. The DEHA model contained two important simplifications compared to that of DPHP - a single rather than a two-phased gut compartment, and the removal of the coding of enterohepatic recirculation - therefore testing of the PBK model for DEHA, sufficient to verify the coding of the model and its ability to capture the trends seen in the BM data of (Nehring et al., 2020) has been undertaken.

Briefly, the model for DEHA described two distinct uptake processes and allowed for a fraction to pass directly through the gut and be ultimately eliminated in faeces (Figure 2). The first uptake process was into blood. The model included a description of absorption of DEHA from the stomach and gastro-intestinal (GI) tract. The second important uptake mechanism of DEHA was into the lymphatic system. Uptake of DEHA via the lacteals in the intestine and entering venous blood after bypassing the liver was coded - a delay function (Lymphlag) was coded to account for the relatively slow transportation of DEHA through the lymphatic system. Inclusion of a lymph compartment was based on the assumption that DEHA, like DEHP, binds like lipid to lipoproteins (Griffiths et al., 1988) which are formed in enterocytes and transported in the lymph to enter the venous blood via the thoracic duct (Kessler et al., 2012). The fractions of dose entering venous blood, the lymphatic system and passing straight through the gut summed to unity.

The model for DEHA had stomach and gut compartments draining into the liver and systemically circulated to adipose, blood (plasma and red blood cell) and slowly and rapidly perfused compartments. Protein binding was described in arterial blood, with only the unbound fraction of DEHA available for distribution to organs and tissues and metabolism. Metabolism of DEHA to MEHA was ascribed to the liver and the gut.

A sub-model was coded to describe the kinetics of MEHA. As described above, metabolism of DEHA to MEHA was coded in the gut and the liver, therefore models for DEHA and MEHA were connected at these nodes in the model. Metabolism of MEHA was coded in the liver alone. The MEHA sub-model had a stomach (Loizou and Spendiff 2004) and intestine draining into the liver and systemically circulated to adipose, blood (plasma and red blood cell) and slowly and rapidly perfused compartments. As with the DEHA model, binding was described in arterial blood.

To make use of biological monitoring data from the human volunteer study of Nehring et al. (2020) it was necessary to describe further metabolism of MEHA in the PBK model. As indicated in earlier discussion, the metabolic pathways of DEHA are non-trivial and downstream metabolites of DEHA may be ultimately eliminated in both urine and exhaled breath. A simplified representation of metabolism of MEHA was assessed as being suitable for the aims of the study with the focus on two immediate metabolites of MEHA, 5OH-MEHA and 5cx-MEPA, eliminated in urine. Amounts of 5OH-MEHA and 5cx-MEPA produced were expressed as fractions of metabolised MEHA, and these were eliminated from blood into urine at a rate proportional to the amount in blood. First order elimination constants described the removal of these respective fractions from blood and into urine. The kinetics of these second order metabolites were thus described using four parameters in all; the model did not describe the distribution of these metabolites to organs and tissues. Whilst neither data on 5oxo-MEHA nor background corrected AA were ultimately used (Table 2), these data were indirectly utilised in forming appropriate ranges for fractions of MEHA metabolised to and eliminated from blood as 5OH-MEHA and 5cx-MEPA.

The model code is available in Supplementary Material.

Probability distributions for uncertainty and sensitivity analysis of the final PBK model are listed in Table 3. Anatomical and physiological parameter distributions were obtained from the freely available web-based application PopGen (McNally et al., 2014). A population of 10,000 individuals comprising 50% Caucasian males and 50% Caucasian females was generated. The range of ages, heights and body weights supplied as input to PopGen were chosen to encompass the characteristics of the volunteers who participated in the human volunteer study (Nehring et al., 2020). Parameter ranges for organ masses and blood flows were modelled by normal or log-normal distributions as appropriate with parameters estimated from the sample and truncated at the 5th and 95th percentiles.

Uniform distributions were ascribed to the various delay terms and uptake and elimination rates. The upper and lower bounds in Table 3 were refined during the model development process. The tabulated values are therefore based upon expert judgement and represent conservative yet credible bounding estimates.

McNally et al. (2021) describe an interactive approach for development and testing of the human PBK model for DPHP using techniques for uncertainty and sensitivity analysis to study the behaviour of the model and the key parameters that drove variability in the model outputs. The principal techniques used for model evaluation were Latin Hypercube Sampling (LHS), to evaluate the qualitative behaviour of the model, and a two-phased sensitivity analysis consisting of elementary effects screening and a variance-based sensitivity analysis to identify the important uncertain parameters in the model to be refined in calibration.

As described previously the DEHA model contained two important simplifications compared to that of DPHP - a single rather than a two-phased gut compartment, and the removal of the coding of enterohepatic recirculation - therefore testing of the PBK model for DEHA, sufficient to verify the coding of the model and its ability to capture the trends seen in the BM data of (Nehring et al., 2020) has been undertaken.

In the first phase of analysis a 200-point Latin Hypercube Design (LHD) was used to draw a sample of parameter sets that efficiently explored the parameter space defined by the parameter distributions given in Table 3. For each of these design points the PBK model was run and data from four outputs—concentrations of DEHA and MEHA (mg/l) in blood, and rates of deposition of 5OH-MEHA and 5cx-MEPA into urine (mg/h) were extracted. The concentration-time profiles over the design points were used to visually assess the bounding behaviour of the model (model form coupled with parameter distributions) and assess whether the model was broadly consistent with trends in the volunteer BM data.

Sensitivity analysis using elementary effects screening was subsequently applied to determine the subset of sensitive parameters to take forward to calibration. A total of 52 parameters were varied with seven elementary effects per model parameter computed, leading to a design of 371 runs of the PBK model. The ranges for each parameter in elementary effects screening were derived from the 2.5th and 97.5^th percentiles for the respective probability distributions (Table 3). The Morris test was applied to the model outputs of: DEHA and MEHA concentrations in venous blood (mg/l) at 0.5- and 5-h following ingestion of DEHA; and for rates of deposition of 5OH-MEHA and 5cx-MEPA into the bladder (mg/h) at 1-, 3-, 5- and 10-h following ingestion of DEHA. Euclidean distance from the origin was computed from the Morris Test output for each parameter, with parameter rankings at each time point based upon this measure. The results were normalised at each time point such that a value of unity corresponded to the most important parameter at a given time point.

The initial filter for further consideration of a parameter to be taken forward into calibration was a normalised Euclidean distance in excess of 0.1 for at least one of the twelve measures. A final subset of sensitive parameters was obtained following a further phase of review.

Calibration is the process of tuning a subset of model parameters such that the discrepancy between model predictions and comparable measurement data is minimised. This is achieved through the specification of an error model that links predictions to measurements. A Bayesian approach to calibration was followed (McNally et al., 2012) as this allows the uncertainty in the concentration-response predicted by the PBK model, which is a function of a subset of sensitive model parameters, to be explicitly quantified.

A Bayesian approach requires the specification of a joint prior distribution for the parameters under study. It is necessary to distinguish between two classes of parameters: global parameters which are common to all individuals, which are appropriate for various constants and physicochemical properties such as partition coefficients etc; local parameters, which vary between individuals are suitable for accounting for variability in the physiology and modelling the participant specific uptake of DEHA etc. These two classes of model parameters are denoted by the vectors θ and ω j respectively, where the subscript j = 1 … 4, denotes the participant. A prior distribution for each global parameter was specified through the distributions provided in Table 3. A prior distribution for each individual, four copies in all, was specified for each of the local parameters. These distributions are also provided in Table 3. A median and 95% credible interval for global and local parameters is provided in Tables 4 (global) and Table 5 (locals) respectively.

The second facet of model specification is the statistical error model. The final calibration model utilised HBM data from the four volunteers with two specific outputs formally compared within the error model. The rates of deposition of 5OH-MEHA and 5cx-MEPA (mg/h) into the urine (RUrine OH and RUrine cx) were computed from the raw data of Nehring et al. (2020) as described earlier, and compared with equivalent predictions extracted from the PBK model through Equations 3, 4.

The terms R U r i n e O H i j and R U r i n e c x i j denote measurement i (at time t i ) for individual j (for j in 1:4) for the two respective model outputs, whereas μ O H _ U θ , ω j i j and μ c x _ U θ , ω j i j , denote the predictions from the PBK model corresponding to parameters ( θ , ω j ). Normal distributions, truncated at zero were assumed for both these relationships, where σ O H _ U and σ c x _ U denote the respective error standard deviations. R U r i n e O H i j ∼   N μ O H _ U θ , ω j i j , σ O H _ U 0 , ∞ (3) R U r i n e c x i j ∼   N μ c x _ U θ , ω j i j , σ c x _ U 0 , ∞ (4)

Weakly informative, half-normal prior distributions with standard deviations of 1 were assumed for the standard deviation parameters in Equations 2, 3.

Inference for the model parameters was made using Markov chain Monte Carlo (MCMC) implemented in MCSim (see Software). Inference for model parameters in the final calibration model was made using thermodynamic integration (TI) as described in Bois et al. (2020). A single chain of 150,000 iterations was run with every 10th retained.

The PBK model was written in the GNU MCSim ³ language and run using the RStudio Version 1.3.1093 ⁴ . The DiceDesign package of R ⁵ . was used for generating Latin Hypercube designs. GSA of model outputs through elementary effects screening was conducted using the Sensitivity package of R. The reshape2 package of R was used for reshaping of data for plotting and other processing of results. MCMC was undertaken using the thermodynamic integration (TI) option within GNU MCSim. All plots were created using R and the ggplot2⁶ package.