As described in (McNally, et al, 2019) in vitro intrinsic clearance for MPHP, CL _{in vitro} (ml min⁻¹ mg⁻¹ microsomal protein) in human hepatic microsomes was calculated using the half-life (T _½ ) derived from the decay constant (k) using the following equations (Obach et al., 1997): i n v i t r o T 1 / 2 = ln ( 2 ) k (1) C L i n v i t r o = ln ( 2 ) i n v i t r o T 1 / 2 × m l i n c u b a t i o n m g m i c r o s o m e s (2) Where, ml incubation is the volume (ml) of the incubation medium and mg microsomes is the mass (mg) of microsomes in the incubation medium.

The intrinsic hepatic clearance CL_{int_H} (L h⁻¹) was calculated using the following formula adapted from (Obach, 1999): C L int _ H = C L i n v i t r o × M P Y × V l i × 60 (3) Where, MPY is the microsomal protein yield per g liver (mg g⁻¹), Vli is mass of the liver (g) and the 60 converts from minutes to hours.

Whole liver plasma clearance CL _H (L h⁻¹) was calculated assuming the well-stirred model of hepatic clearance taking into account the unbound fraction in plasma, fu and the red blood cells to plasma ratio, C_RBC/C_P, using the following equation (Yang et al., 2007): C L H = Q H · f u · C L int _ H Q H + f u · C L int _ H / ( C R B C / C P ) (4) Where, Q _H (L h⁻¹) is the blood flow to the liver as a proportion of cardiac output.

The intrinsic gut clearance CL_{int_gut} was calculated similarly as described for hepatic clearance but substituting MPY _gut and Vgu for MPY and Vli, respectively, in Eq. 4. The resulting calculated CL_{int_gut} was used in place of CL_{int_H} for calculation of CL_gut.

The tissue:blood PCs and unbound fractions in plasma were calculated from the logarithm of the octanol–water partition coefficient, Log P_ow as described in McNally, et al. (2019). Briefly, the Log P_ow for DPHP and MPHP were calculated using the ACDLogP algorithm (Mannhold et al., 2009) implemented in the ACD/ChemSketch 2014 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 DPHP (Table 1). The method of (Schmitt, 2008) , which was developed to predict the tissue distribution of chemicals with Log P_ow < 5.17, was used to predict the PCs of the monoester, MPHP (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 partition coefficient for a particular compartment, the value from a surrogate organ or tissue 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 (5) Where, x = 0 . 4485 logP - 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 proportion of MPHP metabolised to cx- and OH-MPHP, represented by FracMetab (FracMetabcx to cx- MPHP and FracMetabOH to OH-MPHP) (Table 2) for each volunteer was estimated by expressing all the biological monitoring (BM) data (MPHP, OH-MPHP, cx-MPHP, oxo-MPHP) in moles and dividing the amount of cx- and OH-MPHP each by the sum total of all metabolites (Table 2).

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, aged between 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).

A human PBPK model was developed to study the fate of DPHP following single oral doses. The initial model structure was based upon the PBPK model for the plasticizer DINCH described in (McNally, et al, 2019), with a minor modification to account for a large proportion of administered oral dose that is not absorbed, but eliminated by faecal excretion (Klein et al., 2018). The simulation of urinary excretion of metabolites as amount per hour was preferred as they were less variable than when expressed against creatinine to adjust for urinary dilution (Lermen et al., 2019). The model included a description of absorption from the stomach and gastro-intestinal (GI) tract and a simple model of the lymphatic system describing uptake of DPHP via the lacteals in the intestine and entering venous blood after bypassing the liver. Inclusion of a lymph compartment was based on the assumption that DPHP like di(2-ethylhexyl) phthalate (DEHP) binds like lipid to lipoproteins (Griffiths et al., 1988) which are formed in enterocytes and transported in the lymph of the thoracic duct (Kessler et al., 2012). The dose that entered the lymphatic system via the GI tract was coded as a fraction of the administered dose, with the complementary proportion entering the liver via the portal vein. The model described the metabolism of DPHP to MPHP in both liver and gut; therefore, both DPHP and MPHP entered systemic circulation via uptake from the gut. Enterohepatic circulation was also described as a possible explanation of the second small peak in urinary metabolite concentration observed in several datasets. A sub-model was added to describe the kinetics of MPHP, with the two models connected through the gut and liver compartments. The model for DPHP differed from the sub-model only with the presence of a lymphatic component and the MPHP sub-model describing urinary excretion of metabolites. Both models had a stomach and GI tract draining into the liver and systemically circulated to adipose, blood (plasma and red blood cell) and slowly and rapidly perfused compartments. First order elimination constants described the removal of second order metabolites OH-MPHP and cx-MPHP from blood into the urine. A representation of the model structure is given in Figure 2. The model code is available in Supplementary Materials.

The baseline model was subsequently refined using an iterative model development process to better represent the trends in BM (blood and urine) data from the human volunteer study reported in (Klein et al., 2018). Techniques for uncertainty and sensitivity analysis (described in the statistical analysis section) were deployed, at each iteration of model development, to establish the bounding behaviour of the model and the key uncertain parameters that the model outputs under study were sensitive to. The improvements and a brief justification are described in the points below:

1. The model was adapted to account for a majority fraction of DPHP passing through the GI tract without being absorbed, as suggested by BM data, ranging from 75% (Wittassek and Angerer 2008; Leng et al., 2014) to 94% (Klein et al., 2018). FracDOSEHep and FracDOSELymph modelled the fractions of the administered dose entering hepatic and lymphatic circulation respectively, with the complementary fraction (1 – FracDOSEHep – FracDOSELymph) passing directly in faeces.

Baseline estimates of organ and tissue masses and regional blood flows were taken from Brown et al. (1997) and ICRP (2002). The mass of the lymphatic system was obtained from Offman et al. (2016).

Tissue: blood partition coefficients were estimated using algorithms as described previously.

The bio-transformation of DPHP to MPHP in the liver was described by an intrinsic clearance term determined in vitro and scaled to in-vivo: the half-life of DPHP was estimated using the PBPK model. The in-vivo intrinsic clearance of DPHP in the gut was calculated using the in vitro hepatic clearance scaled to in vivo using gut microsomal protein yield and gut volume. The bio-transformation of MPHP to second order metabolites (OH-MPHP and cx-MPHP) in the liver was of the same form as the expression for DPHP, however an estimate of the half-life for MPHP was determined experimentally (as described above). A single term for metabolism of MPHP was coded in the PBPK model with the rates of removal of two direct metabolites, (OH-MPHP and cx-MPHP), from plasma assumed to be proportional to the rate of metabolism of MPHP. A urinary elimination constant was estimated for OH-MPHP and cx-MPHP.

Baseline values for parameters for which there was no prior knowledge such as FracDOSEHep, FracDOSELymph and the various delay terms and uptake and elimination rates were determined during the model development and testing process to provide a reasonable (but not optimised) fit to BM data.

Baseline (default) values are given in Table 3.

As described above, uncertainty analysis was conducted throughout the model development process in order to efficiently establish the bounding behaviour of the model (i.e. the variations in model outputs under study that were consistent with the current version of the model, and parameter value uncertainty defined through probability distributions). A 200 point maxi-min Latin Hypercube Design (LHD) was created based upon the probability distributions ascribed to model parameters and the PBPK model was run for each of these design points; the behaviour of the final model was studied based upon the probability distributions given in Table 3.

The development process followed here was broadly similar to that of (McNally et al., 2019). However whereas (McNally et al, 2019) monitored only three outputs from their PBPK model for DINCH, eight outputs from the model for DPHP were monitored in order to study the absorption, uptake, metabolism and excretion of DPHP, these were: amount of DPHP (mg) in the bowel compartment (i.e. DPHP excreted in faeces), amount of DPHP (mg) in the lymph compartment, concentrations of DPHP and MPHP in venous blood, masses of DPHP and MPHP in the plasma compartment (i.e. bound to proteins within plasma and hence unavailable for metabolism), rates of deposition of OH-MPHP and cx-MPHP in urine (mg/hour). Furthermore mass balance of DPHP and MPHP were monitored to ensure that mass balance was retained for all the tested parameter variations. The differing units of the outputs under study reflect the different aspects of model outputs that the uncertainty analysis was designed to study. This phase of work is only briefly reported on in results.

Sensitivity analysis was conducted throughout the model development process in order to study the key model output sensitivities for each version of the model under development. A two-phased GSA was implemented (McNally et al., 2011; Loizou et al., 2015) comprising of elementary effects screening (Morris Test) followed by a variance-based approach. Results from sensitivity analysis for the final model were obtained using the probability distributions given in Table 3.

A total of 59 parameters were varied in elementary effects screening, with five elementary effects per input computed, leading to a design of 300 runs of the PBPK model. The model outputs studied are described below.

Concentrations of DPHP and MPHP in venous blood at 0.5, 3 and 12 h and 1, 3 and 12 h following ingestion, respectively were studied using elementary effects screening. The three output times studied were broadly of representative of the following periods in the concentration-time courses: prior to peak concentration (of DPHP and MPHP); post peak concentration; and returning to baseline (zero) concentrations. Rates of deposition of OH-MPHP and cx-MPHP in urine were studied at 3-, 12- and 20-h following ingestion for both model outputs with these times corresponding to the periods where peak concentration in urine was reached; when the first harmonic (due to `hepatic recirculation) was predicted (12 h); and returning to baseline. Finally, the concentrations of DPHP and MPHP in plasma were studied. Rather than studying model output at specific time points, instead the peak concentrations of DPHP and MPHP in plasma; the times that corresponded to these peak concentrations; and the rate of change of DPHP and MPHP in plasma over the hour following the peak concentrations, were extracted from each of the 300 model runs. These measures were chosen since they proved to be more useful metrics for understanding the BM data of Klein et al. (2018), and in particular the more rapid clearance from plasma than would be expected given the predictions of logP (and from this an estimate of protein binding) (McNally et al., 2019). This phase of sensitivity analysis is only briefly described in results.

All parameters that were within 0.2 of the maximum value of µ* (one of the two measures computed in elementary effects screening) for any of the 18 metrics studied (three metrics for each of the six outputs described above) were retained and studied in the second phase of analysis using the variance-based analysis.

In the second phase of analysis 31 retained parameters were studied using the extended Fourier Amplitude Sensitivity Test (eFAST). In this analysis 1,000 runs per retained parameter were conducted, leading to 31,000 simulations of the PBPK model. In this more computationally expensive phase of sensitivity analysis, the rates of deposition of OH-MPHP and cx-MPHP in urine and concentrations of DPHP and MPHP in plasma were studied.

Calibration of a subset of sensitive model parameters using the BM data of (Klein, et al., 2018) was attempted. A Bayesian approach was followed (McNally et al., 2012). This requires the specification of a joint prior distribution for the parameters under study, which is refined through a comparison of PBPK model predictions and measurements within a statistical model. The resulting (refined) parameter space that is consistent with the prior specification and measurements is the posterior distribution.

The final calibration model utilised data from five of the six individuals (data from individual E were unusual and thus excluded) from the BM study of (Klein et al., 2018) with data on four specific outputs formally compared within the calibration model. Concentrations of DPHP and MPHP (CBlood DPHP and CBlood MPHP) and the rates of deposition of OH-MPHP and cx-MPHP (mg/hour) into the bladder (RUrine OH and RUrine cx) were computed from the raw data (Klein et al., 2018). The latter measure i.e. the concentration (mg/L), the volume of the void (ml) and the times between voids (hours)) – (see for example Nehring et al. (2020) represents the underlying trends in concentration response data in a more precise manner than expressing the metabolite concentration relative to creatinine concentration. These measurements were compared to corresponding predictions from the PBPK model using the statistical models depicted in Eqs. 7–10.

The terms RUrineOHij,RUrinecxij,CBloodDPHPij and CBloodMPHPij denote measurement i (at time ti) for individual j (for j in 1:5) for the four respective model outputs, whereas μ O H ( θ , ω j ) i j , μ CX ( θ , ω j ) i j , μ D P H P ( θ , ω j ) i j and μ M P H P ( θ , ω j ) i j denote the predictions from the PBPK model corresponding to parameters ( θ , ω j ) . The vectors θ and ω j denote the global parameters common to all individuals (suitable for partition coefficients, fractions metabolised to OH-MPHP and cx-MPHP etc.) and participant specific parameters (suitable for accounting for variability in the physiology and modelling the participant specific uptake of DPHP etc.). σ O H , σ C X , σ D P H P and σ M P H P   denote the respective error standard deviations. Normal distributions, truncated at zero were assumed for all four relationships. R U r i n e O H i j ∼ N ( μ O H ( θ , ω j ) i j , σ O H ) [ 0 , ∞ ] (7) R U r i n e c x i j ∼ N ( μ c x ( θ , ω j ) i j , σ c x ) [ 0 , ∞ ] (8) C B l o o d D P H P i j ∼ N μ D P H P ( θ , ω j ) i j , σ D P H P [ 0 , ∞ ] (9) C B l o o d M P H P i j ∼ N ( μ M P H P ( θ , ω j ) i j , σ M P H P ) [ 0 , ∞ ] (10)

Prior distributions for global and local parameters in the PBPK model were taken from Table 3 (for the sensitive parameters that were studied). Non-informative gamma (0.01, 0.01) prior distributions were assumed for the four standard deviation parameters.

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 thermo-dynamic integration (TI) as described in Bois et al. (2020). A single chain of 1,000,000 iterations was run with every 10th retained.

The PBPK model was written in the R language and run using the RStudio and RVis software applications during the development of the PBPK model. PBPK models were solved using the deSolve package of R. The DiceDesign package of R was used for generating Latin Hypercube designs. GSA of model outputs (through elementary effects screening and eFAST) were 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.

The PBPK model was rewritten in the GNU MCSim language prior to calibration. MCMC was undertaken using the TI option within GNU MCSim

All plots were created using R and the gg2plot package.

In Figure 4 a small subset of results from uncertainty analysis of the final model are shown: each curve shows the predictions of a particular model output corresponding to a design point. These plots indicate that the model structure remains stable over a wide range parameter values. Figures 4A–C show the mass of DPHP within the plasma, lymph and bowel compartments and were used to study the range of behaviours of specific aspects of the model that could be achieved based upon the form of the PBPK model and through variations in the uncertain parameters. Through this uncertainty analysis of the mass of DPHP in plasma (Figure 4A) the effect of protein binding on the mass retained and subsequent elimination from plasma could be studied. Specifically, in this exploratory phase of work we had a particular interest in attempting to replicate the unusual findings from (Klein et al., 2016), who reported peak concentrations of DPHP in the blood of volunteers which occurred later than the corresponding peak of the metabolite MPHP in blood and also later than peak concentrations of second order metabolites in urine. Through uncertainty analysis of DPHP in the lymph, variability in the uptake, retention and elimination of DPHP in the lymph compartment could be studied. The uncertainty analysis of DPHP in the bowel allowed the study of variability in the fraction of DPHP that was initially unabsorbed and the further removal of DPHP following elimination in bile. Figure 4D shows predictions of the rate of elimination of OH-MPHP in urine (mg/hr), one of the chosen metrics for calibration, and was studied to establish how well the unique trends in urine voids from the six volunteers could be captured by the final model. Other checks on a range of outputs or functions of model outputs (masses, rates and concentrations) were also undertaken in this phase of modelling to ensure behaviour of the model appeared reasonable over the range of parameter space specified through probability distributions.

Figure 5 shows the results from elementary effects screening (Morris Test) applied to the mass of DPHP in plasma. A high μ* indicates a factor with an important overall influence on model output; a high σ indicates either a factor interacting with other factors or a factor whose effects are non-linear. The magnitude of μ* and σ for each model parameter is relative, i.e., a parameter has a low μ* relative to the parameter with the highest μ*. Results from this technique are usually obtained at specific time points, i.e. sensitivity analysis of a given model output at say, 1 h following dosing. However sensitivity analysis can be applied to any chosen model outputs calculated from each model run specified through the design. The results shown Figure 5 correspond to some unusual measures calculated from model output that were chosen to study particular aspects of model behaviour: sensitivity analysis of the peak mass of DPHP in plasma (mg) (A), the time when peak concentration was reached (hours after dosing) (B) and the rate of change of DPHP in plasma (mg/hour) in the hour following peak concentration (C). The parameters with lower overall importance are clustered toward zero of both axes. Unfortunately, we could not prevent the overlapping of some parameter labels in this region (Figure 5).

The parameters ranked as most important by the Morris test were analysed by eFAST. The period from the start of the simulation to 20 h showed the most variance in blood concentrations of DPHP and MPHP. The most important parameters were reduced further in number and ranked as follows: FBMPHP, FBDPHP, FracDoseLymph, FracDoseHep, Lymphlag, PbaM, QCC, K1Lymph, VspdC, VliC, QguC, and VbldC. The fractions of MPHP and DPHP bound to plasma proteins were significantly more important than the other parameters over this period.

The period from the start of the simulation to 25 h showed the most variance in urinary excretion of OH-MPHP and cx-MPHP. The most important parameters were ranked as follows: FracDoseHep, FracMetabMOH, FracMetabcx, BW, K1_MOH, K1_cx, QCC, GIPERM1, PguM, QguC, PliM, and PbaM. The first four parameters were significantly more important for variance in urinary excretion than the other parameters over this period.

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 Table 3 and Table 4 respectively. The fit of the calibrated model is shown in Figures 6A–D, Figures 7A–D and Figures 8A–D for three of the five participants. The trends for the individual shown in Figure 8 are broadly representative of the two individuals whose data are not shown. The central estimates indicated in plots correspond to the posterior mode whereas the shaded regions represent 95% intervals for the respective curves. This is a pointwise credible interval which was derived through running the PBPK model for each retained parameter set, ordering the predictions by magnitude at each time point and reading off the 2.5^th and 97.5^th percentiles.

The BM data from each of the study participants showed unique trends. The blood DPHP data indicated differences in lag prior to uptake, rates of uptake and fractions absorbed from the gut. In the blood data of some participants, the first order metabolite MPHP was measurable in blood prior to the appearance of parent chemical (Figure 6A and Figure 7A). Uptake of DPHP appeared to be multi-phased in some participants. The PBPK model contained nine parameters that were tuned to the BM data from each participant for fitting the hepatic uptake of DPHP and a further three parameters governing uptake into the lymphatic system and subsequent deposition into blood at the thoracic duct (Table 3). The results in Figures 6–8 demonstrate that the BM data on DPHP in blood were captured by the calibrated model. The differences in the calibrated individual specific parameters (Table 3) reflect the large differences in the trends of DPHP in blood for the study participants.

In contrast, the rate of deposition of second order metabolites in urine was mainly governed through global (non-volunteer specific) parameters (Table 4). Reasonable fits were obtained for the urinary excretion of OH-MPHP and cx-MPHP with the model successfully fitting earlier peaks in the urine compared to blood. Enterohepatic recirculation of DPHP was an important mechanism which is observed in the harmonics seen in the OH-MPHP and cx-MPHP time courses at intervals of approximately 8 h post peak concentration (Figures 6C,D and Figures 8C,D). The fit to the data shown in Figures 7C,D was poorer and appears to be consistent with a second absorption event (not explicitly captured in the model). Data from another volunteer, not used in final calibration, showed even stronger evidence of a second absorption event. As enterohepatic recirculation was modelled using global parameters, the differences seen in the simulations of urine data for these three participants appear to arise as a consequence of differences in uptake of DPHP.

The fits to MPHP in blood were generally quite poor; particularly when MPHP spiked in blood specimens shortly after DPHP was consumed by study participants. This highlights a deficiency in the model, which does not impact upon the ability of the model to accurately predict deposition of second order metabolites in urine. We address this deficiency further in the discussion section.

The bound fractions of DPHP and MPHP were 0.991 (0.989, 0.995) and 0.975 (0.956, 0.986), which represent very high binding in blood, although somewhat smaller than the very high fractions predicted by algorithms. There are no known direct measurements of protein binding of DPHP or MPHP in blood although there are estimates of the “free” area-under-the curve (AUC) concentrations of MPHP as a proportion of total DPHP concentration which are approximately 66% (Klein et al., 2018). The half -life of DPHP, which could not be estimated in vitro incubations, was estimated in the model to be very short in the liver at 4.38 min (3.05, 9.10), and approximately a factor of 10 greater than in the gut 56.80 min (46.96, 59.90).

The vast majority of DPHP was unabsorbed from the gut. Total uptake ranged from 16.3% (13.1%, 19.3%) to 2% (1.5–2.6%), dominated by hepatic uptake at around a factor of five greater than lymphatic uptake. The fractions of absorbed DPHP excreted as OH-MPHP and cx-MPHP were estimated as 0.24 (0.21, 0.31) and 0.011 (0.01, 0.014).

A comparison of simulations of entry of DPHP through the hepatic (black line) and lymphatic (blue line) routes is shown in Figure 9. These simulations were based upon optimised parameters for individual A where with lymphatic fraction was set to zero (black lines) or the hepatic fraction set to zero (blue lines). Thus, the independent effects of absorption through the two routes at otherwise credible values for model parameters could be studied. The key biological difference between the two routes of entry is that DPHP absorbed through the hepatic route is subject to first pass metabolism, primarily in the liver, whereas the lymphatic fraction by-passes first pass metabolism. Despite the hepatic fraction being a factor of five greater, the peak plasma concentration from the two routes was similar (Figure 9A), which indicates that a very large fraction of DPHP absorbed through the hepatic route is intercepted by the liver and never enters systemic circulation. Uptake of DPHP into the systemic circulation from either hepatic or lymphatic routes is almost completely bound to proteins in plasma and is thus neither available for distribution to organs and tissues nor for metabolism. Our simulations indicate that DPHP entering the systemic blood circulation via the lymphatic route is almost entirely held in plasma and is important for understanding trends in plasma. It does however represent a small mass of DPHP, which has a negligible impact on trends in MPHP in blood and second-order metabolites in urine. In particular, due to the high binding in plasma and the consequent reduced rate of metabolism only a shallow peak with a long tail can be seen in urinary metabolite simulations.