An existing human PBPK model for DPHP (McNally et al., 2021) was modified further to study the absorption, distribution, metabolism, and elimination of DPHP following single oral doses. The key changes compared with McNally et al. (2021) were the inclusion of a kidney compartment so that an additional dataset of MPHP in urine could be utilised in model calibration, a modification to the uptake of MPHP from the gut to facilitate an improved fit to DPHP and MPHP in blood, and in the absence of Michaelis-Menten constants, a simple modification to the metabolism of MPHP, specified using clearance, to account for a limitation of metabolism, whose mechanism is uncertain, perhaps saturation or inhibition of metabolism, thus allowing a fraction of MPHP created in the liver through metabolism of DPHP to escape into blood. Furthermore, a few simplifications of the model structure were made where data indicated this was appropriate (Figure 2).

Briefly, the model for DPHP described two distinct uptake processes and allowed for a majority fraction to pass directly through the gut and be ultimately eliminated in faeces. The primary uptake process was into blood. The model included a description of absorption from the stomach and gastro-intestinal (GI) tract, simplified to a two-stage intestine compartment (Figure 3). Uptake of DPHP into venous blood from the stomach, metabolism of DPHP to MPHP and uptake of MPHP into venous blood was ascribed to the first gut compartment, and finally uptake of DPHP into venous blood was ascribed to the second-phase gut-compartment. A parameter, Gutlag was included to represent a delay in transport of DPHP through the gut transiting from the first to second gut compartment. An important change compared with McNally et al. (2021) was that a fraction of MPHP absorbed through the first gut compartment was unavailable for first pass metabolism in the liver following uptake into venous blood—this represents incomplete binding. The second uptake mechanism of DPHP was into the lymphatic system. Uptake of DPHP via the lacteals in the intestine and entering venous blood after bypassing the liver was coded. Inclusion of a lymph compartment was based on the assumption that DPHP, 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.

Metabolism of DPHP to MPHP was ascribed to the liver and a section of the gut (Figure 3, gut 1). The model for DPHP additionally encoded the transport process of enterohepatic recirculation. Uptake of DPHP from the liver into bile was modelled as a first order uptake process with a delay (to represent transport from liver to gut) before DPHP appeared in the small intestine where DPHP was available for reabsorption.

The model had a stomach and (the two-phase intestine) draining into the liver and systemically circulated to adipose, kidney, blood (plasma and red blood cell) and slowly and rapidly perfused compartments (Figure 2). Protein binding was described in arterial blood, with only the unbound fraction of DPHP available for distribution to organs and tissues and metabolism.

A sub-model was coded to describe the kinetics of MPHP. As described above, metabolism of DPHP to MPHP was coded in the first intestinal compartment and the liver, therefore models for DPHP and MPHP were connected at these nodes in the model. Metabolism of MPHP was coded in the liver alone. A fraction of MPHP (Escapeli) was coded as unavailable for immediate metabolism: in the absence of in vitro determined Michaelis-Menten constants this is a simple modification to the clearance-based first-order model for metabolism to the represent saturation or inhibition of metabolism. The MPHP sub-model had a stomach and (single-phase) intestine draining into the liver and systemically circulated to adipose, kidney, blood (plasma and red blood cell) and slowly and rapidly perfused compartments (Figure 2). Elimination of MPHP from the kidney was described with a first-order elimination rate; this represents a modification compared to McNally et al. (2021). Furthermore, in contrast to the previously published model, enterohepatic recirculation of MPHP was not considered. As with the DPHP model, binding was described in arterial blood.

To make use of biological monitoring data on two metabolites of MPHP (OH-MPHP and cx-MPHP) it was necessary to include the elimination of these substances into urine. A simplified representation of these downstream metabolites was assessed as being suitable for the aims of modelling. The appearance of OH-MPHP and cx-MPHP in blood were coded as respective fractions of metabolised MPHP. First order elimination constants described the removal of second order metabolites, OH-MPHP and cx-MPHP from blood into the urine. The model did not describe the distribution of these metabolites to organs and tissues.

The model code is available in Supplementary Materials.

Anatomical, physiological, physicochemical, and metabolic parameters were very similar to the previous model for DPHP (McNally et al., 2021) with small differences arising due to changes to model structure. These parameters are provided in Tables 1–3. Sources, derivation and calculation of parameters are described in Supplementary Materials.

The BM data described in (Klein et al., 2018) were kindly provided by Dr. Rainer Otter of BASF, SE. Briefly, DPHP was administered orally to six healthy male volunteers (designated A to F), aged be-tween 30 and 64 years, weighing between 74 and 108 kg. A single dose of 738 ± 56 μg/kg BW DPHP was administered as an emulsion of 7% (w/v) in an aqueous saccharose solution (70% w/v) between 45 and 140 min after breakfast. The resulting respective doses for the six individuals were between 0.639 and 0.783 mg/kg body weight (Table 2).

Klein et al. (2018) report on the trends in blood specimens of volunteers over a 24-h period following ingestion of DPHP with samples collected at 0.25-, 0.5-, 0.75-, 1-, 2-, 3-, 4-, 5-, 6-, 7-, 8-, 10- and 24-h following exposure. The total concentrations of DPHP and MPHP (nM) were extracted from the dataset. BM data were converted to units of mg/L or µg/L prior to use in calibration. Klein et al. (2018) also reported on trends of metabolites of DPHP in urine specimens over a 46-h period following ingestion, with samples collected at 1-, 2-, 3-, 4-, 5-, 6-, 8-, 10-, 12-, 14-, 18-, 22-, 26-, 30-, 34-, 38-, 42- and 46 h. The concentrations of MPHP, cx-MPHP and OH-MPHP (mg/L) were extracted from the dataset. The rates of deposition of metabolite into the bladder (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 (Nehring et al., 2020). The derived rate was associated with the mid-point between the two voiding events.

Following a review of data, the BM data from volunteers C and E were inconsistent with data from the other four volunteers; for these two volunteers a second large spike at 15 and 20 h of second order metabolites (cx-MPHP and OH-MPHP) was observed in urine voids (Figure 4). This feature of the data would be consistent with a second significant uptake of DPHP from the gut associated with the intake of food/drink. An exclusion of these data was regarded by the authors as preferable to further modifications to model structure based on weak evidence. Only data from volunteers A, B, D and F were taken forward into calibration. We comment further on exclusions of data in the discussion.

Probability distributions for uncertainty and sensitivity analysis of the final PBPK 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 of 100% Caucasian males 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 (Klein et al., 2018). 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 PBPK 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. For this refined model a less intensive testing process was required. The key goal of this analysis was to determine whether the change to encode partial binding of MPHP in the gut was sufficient to allow the model to better fit the BM data on concentrations of MPHP in venous blood, whilst not sacrificing the quality of fit of predictions of the deposition of OH-MPHP and cx-MPHP in urine. Furthermore, the study of MPHP in urine, a new measure for use in calibration, was necessary to ensure model predictions were broadly credible. Sensitivity analysis of the concentration of MPHP in venous blood and the rate of deposition of MPHP in the bladder (mg/h) was subsequently conducted to assess whether further sensitive parameters should be taken forward to calibration.

Uncertainty analysis was conducted though rejection sampling. A 500-point maxi-min Latin Hypercube Design (LHD) based upon the probability distributions given in Table 3 was created and the PBPK model was run for each of these design points. All simulations where predictions of the peak concentration of MPHP in blood (mg/L) was less than 0.025 or greater than 0.3 mg/L and where predictions of MPHP in blood (mg/L) at 48 h were more than 0.025 mg/L were rejected. These were generous constraints based upon the trends in BM from the volunteers. The concentration response trends for MPHP in blood (mg/L), and the rate of deposition of MPHP and OH-MPHP into urine (mg/h) were studied for the retained sample to assess the credibility of the model form.

Sensitivity analysis of MPHP in blood (mg/L) and rate of deposition of MPHP into the bladder was conducted using elementary effects screening (Morris test). A total of 61 parameters were varied with five elementary effects per input computed, leading to a design of 325 runs of the PBPK model. The Morris test was applied to model output at 0.5- and 3-h following ingestion, which were broadly representative of the prior to peak concentration and post-peak concentration periods.

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 was followed (McNally et al., 2012) for calibration in this work since this allows the uncertainty in parameters (and thus on the concentration response predictions from the PBPK model) 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 (appropriate for various constants and physicochemical properties such as partition coefficients etc.,); and local parameters, which vary between individuals (suitable for accounting for variability in the physiology and modelling the participant specific uptake of DPHP etc.,). These two classes are denoted by the vectors θ ; ω 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% 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. As described earlier, the final calibration model utilised data from four of the six individuals with data on five specific outputs formally compared within the error model. Concentrations of DPHP and MPHP (CBlood DPHP and CBlood MPHP) (mg/L) and the rates of deposition of MPHP, OH-MPHP and cx-MPHP (mg/h) into the urine (RUrine MPHP, RUrine OH and RUrine cx) were computed from the raw data of (Klein et al., 2018), as described earlier, and compared with equivalent predictions extracted from the PBPK model through Eqs. 1–5.

The terms C B l o o d D P H P i j , C B l o o d M P H P i j , R U r i n e M P H P i j , R U r i n e O H i j ad 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 five respective model outputs, whereas μ D P H P _ B θ , ω j i j , μ M P H P _ B θ , ω j i j , μ M P H P _ U θ , ω j i j , μ O H _ U θ , ω j i j and μ c x _ U θ , ω j i j , denote the predictions from the PBPK model corresponding to parameters ( θ , ω j ). Normal distributions, truncated at zero were assumed for all five relationships and where σ D P H P _ B , σ M P H P _ B , σ M P H P _ U , σ O H _ U and σ c x _ U denote the respective error standard deviations C B l o o d D P H P i j ∼ N μ D P H P _ B θ , ω j i j , σ D P H P _ B 0 , ∞ (1) C B l o o d M P H P i j ∼ N μ M P H P _ B θ , ω j i j , σ M P H P _ B 0 , ∞ (2) R U r i n e M P H P i j ∼   N μ M P H P _ U θ , ω j i j , σ M P H P _ U 0 , ∞ (3) R U r i n e O H i j ∼   N μ O H _ U θ , ω j i j , σ O H _ U 0 , ∞ (4) R U r i n e c x i j ∼   N μ c x _ U θ , ω j i j , σ c x _ U 0 , ∞ (5)

Weakly informative, half-normal prior distributions with standard deviations of 1 were assumed for the five standard deviation parameters in Eqs. 1–5.

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 PBPK model was written in the GNU MCSim language (version 6.1.0) ¹ and run under Windows 10 Pro. Files for running MCSim under windows, tools and instructions for installation are available from Github ² . Scripts were written for calling the model from RStudio (RStudio Team 2016). R version 4.0.2 and the ggplot2, DiceDesign, Sensitivity and reshape2 packages were used in analysis (Dupuy et al., 2015; Pujol and Iooss 2015; R Development Core Team 2008; Wickham 2007; 2016).

Extensive results from uncertainty and sensitivity analysis of our earlier PBPK model for DPHP are presented in McNally et al. (2021). We therefore limit our attention to the novel results from this analysis; our efforts to better fit the trends of MPHP in venous blood from the human volunteer study of Klein et al. (2018) and the change to model form so that the rate of deposition of MPHP into the urine could be simulated.

The results from uncertainty analysis are shown in Figure 5. Panel Figure 5A shows the concentration response profiles for MPHP in venous blood (mg/L) for the full 500 samples generated using a computationally efficient Latin hypercube design. The profiles indicate a wide range of behaviour for this output that is consistent with the prior distributions for model parameters. After applying the rejection criteria, specified in methods, 248 simulations were retained, and this subset of concentration response profiles is shown in Figure 5B; this analysis indicated that the model could simulate concentration-time profiles for MPHP in venous blood that were broadly consistent with the BM data of Klein et al. (2018); Figures 5C, D show predictions from the retained sample for the rates of deposition of MPHP and OH-MPHP into the urine: this analysis aims to study whether the modifications necessary to better fit the data trends in venous blood compromised the ability of the model to fit trends in urine data. These results (Figures 5C, D) showed considerable variation in the profiles from the different runs in the retained sample, with a subset of the simulations appearing to be broadly consistent with trends in BM data observed in urine data (Klein et al., 2018).

Detailed results from elementary effects screening are not provided, however briefly this work aimed to explore whether additional parameters should be taken forward into calibration. The analysis was limited to concentrations of MPHP in venous blood (mg/L) at 0.5- and 3- hours and the rate of deposition of MPHP in the bladder (mg/h) also at 0.5- and 3-h following ingestion of DPHP. The proportion of MPHP unavailable for metabolism due to partial binding (escapeFracgu), the fraction of MPHP escaping from the liver due to saturation or inhibition of metabolism (escapeFracli), the mass of the kidney (as a fraction of body weight) (Vki), the blood:kidney partion coefficient for MPHP (PkiM) and the elimination rate of MPHP from kidney tissue (K1_MPHP) were identified as additional parameters to tune in model calibration. In total, the five model outputs (Figures 6–9) used in calibrations showed some sensitivity to 34 parameters (Tables 4 and 5).

Summary statistics based upon the retained sample (posterior median and a 95% credible interval) for the global and local (volunteer specific) parameters are provided in Tables 4, 5 respectively. The fit of the calibrated model is shown in Figures 6–9 for individuals A, B, D and F respectively. The five panels in each figure correspond to A) DPHP in venous blood (mg/L); B) MPHP in venous blood (mg/L); C) deposition of MPHP into urine (mg/h); D) deposition of OH-MPHP in urine (mg/h); E) deposition of cx-MPHP in urine (mg/h). The central estimates indicated in plots correspond to the posterior mode parameter set, the single best fitting parameter set over the 20 measures (5 outputs for each of 4 individuals) used for calibration. The shaded regions represent pointwise 95% credible intervals for the respective curves. This interval was derived by running each retained sample drawn from the posterior through the PBPK model and storing the output from each model output from 0 to 48 h in 0.05- hour increments. Output at each time point was ordered with the 2.5th and 97.5th percentiles read off; the plotted 2.5% and 97.5% bounds are a smooth interpolation of these series of pointwise values.

The fits to BM data shown in Figures 6–9 demonstrate the PBPK model was able to simulate the diverse trends of DPHP and downstream metabolites in the blood and urine specimens from the four volunteers. A very high quality of fit to the deposition rates of OH-MPHP and cx-MPHP in urine was previously achieved in (McNally et al., 2021) and this has not been improved upon with the updated model. However, the fits to DPHP and MPHP in blood represent very significant improvements compared to those of (McNally et al., 2021). The shape of the data series on MPHP was generally captured for the four volunteers, however the model generally struggled to simulate the very rapid reduction in concentrations following the peak—the discrepancy was greatest for individual A (Figure 6), The new calibration measure considered in this work, rate of deposition of MPHP in urine (panel C of Figures 6–9), was generally well fitted, although as with MPHP in blood, the model struggled to capture the very sharp rate of decline following the peak.

The two parameters that were most critical to the improvement relative to (McNally et al., 2021) were the proportion of MPHP binding in blood within the gut and unavailable for first pass metabolism in the liver (Escape_gu), and the proportion of MPHP escaping from the liver due to saturation or inhibition of metabolism (Escape_li). In the absence of these parameters the model was unable to approximate the trends of MPHP in blood. The uniform priors U(0,1) for these parameters substantially narrowed following calibration with estimates of 0.59 (0.31, 0.83) and 0.59 (0.45, 0.75) for Escape_gu and Escape_li respectively, suggesting that a majority fraction of MPHP evades first pass metabolism in the liver, and binds with plasma, however it is subsequently rapidly removed.