Development and validation of a physiology-based model for the prediction of pharmacokinetics/toxicokinetics in rabbits

The environmental fates of pharmaceuticals and the effects of crop protection products on non-target species are subjects that are undergoing intense review. Since measuring the concentrations and effects of xenobiotics on all affected species under all conceivable scenarios is not feasible, standard laboratory animals such as rabbits are tested, and the observed adverse effects are translated to focal species for environmental risk assessments. In that respect, mathematical modelling is becoming increasingly important for evaluating the consequences of pesticides in untested scenarios. In particular, physiologically based pharmacokinetic/toxicokinetic (PBPK/TK) modelling is a well-established methodology used to predict tissue concentrations based on the absorption, distribution, metabolism and excretion of drugs and toxicants. In the present work, a rabbit PBPK/TK model is developed and evaluated with data available from the literature. The model predictions include scenarios of both intravenous (i.v.) and oral (p.o.) administration of small and large compounds. The presented rabbit PBPK/TK model predicts the pharmacokinetics (Cmax, AUC) of the tested compounds with an average 1.7-fold error. This result indicates a good predictive capacity of the model, which enables its use for risk assessment modelling and simulations.


Introduction
Environmental and human risk assessments are based on standardized biotests that are conducted on different species. To minimize the use of animals, the 3R principle (i.e., reduction, replacement, and refinement) requires that animal experiments are substituted with appropriate alternative test methods whenever possible [1][2][3][4]. This principle is already a current practice in REACH [5,6], which is a binding European Union regulation regarding the production PLOS ONE | https://doi.org/10.1371/journal.pone.0194294 March 21, 2018 1 / 17 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 data. The discrepancies between observations and the simulation were resolved by optimizing the parameters or retrieving independent information from the literature.

Model structure and parameterization
For a comprehensive PBTK model as described here, a stepwise approach is typically used to inform the parameters. The first level retrieves the parameter values from the literature to describe the general properties of the species and is independent of the compound to be modelled. Second, the modelling of a specific compound starts by informing the model with the physicochemical properties of the molecule. This information is readily available or can be calculated independently [19]. Third, the active processes (e.g., transport and metabolism) are parameterized based on the information of in vitro or in vivo assays. Information from other species can be also used when there is no information for the species of interest. Lastly, if exposure data for a compound in the target species is available, this can be used to fine-tune the selected parameters. A schematic workflow of the model establishment is shown in Fig 1. In our model, physiological parameters of rabbit were introduced in the PK-Sim 1 framework in order to represent rabbit physiology (first level). The parameters identified from the literature are shown in supporting information Tables A-J in S1 Appendix. Furthermore, a complete list of the rabbit PBPK model parameters is shown in Table O in S1 Appendix. If parameters of rabbit physiology were not available in the literature, the corresponding parameters of a preexisting PBPK model were used. For this, the model of the taxonomically related mouse was used. Compound related information (second level) and information for respective active processes (third level) were then adopted from online databases and in house experiments and are shown in Tables K-L in S1 Appendix. Finally, fine-tuning (fourth level) was performed by fitting the rabbit model to available experimental data (S1-S10 Figs). The rabbit model is based on the mammalian physiological framework implemented in the PK-Sim 1 software platform [20]. This whole-body physiology-based pharmacokinetic Starting with a generic mammalian PBPK model on top, rabbit specific physiology parameters are added (first level). Next, physicochemistry from a compound (second level), followed by information regarding active processes (third level) are introduced from literature, in-house experiments or transferred from other species. Lastly, PK data are used to identify additional parameter values of the model (fourth level). framework provides a detailed mathematical representation of ADME processes in humans and several common laboratory animals. In essence, the mammalian body is divided into compartments representing the relevant organs or tissues as well as the arterial and venous blood pools that connect the different organs via blood flow (Fig 2A). Organs are further divided into sub-compartments that describe the vascular space (divided into plasma and red blood cells) as well as the avascular space (divided into interstitial and cellular space) (Fig 2B and  2C). The characteristics of the administered compound (small or large molecules such as proteins) determine its distribution by either diffusion through membranes (Fig 2B) or by pores and endosomal clearance (Fig 2C). A substantial amount of experimental information has been used to inform and calibrate the biological and physiological processes included in the model. The model has successfully been used to describe the uptake and distribution of a wide range of compounds in the species included to date [10,[21][22][23][24][25][26].
Based on the current animal models implemented in the PK-Sim 1 software suite, the model for a 2.5 kg rabbit was built by adjusting the values of organ volumes ( Table A in S1 Appendix), haematocrit ( Table B in S1 Appendix), specific blood flow rates (Table C in S1 Appendix), GIT compartment pH (Table D in S1 Appendix), GIT compartment dimensions (Tables E-F in S1 Appendix), GIT compartment transit/emptying times (Table G in S1 Appendix), and the effective surface areas of the GIT-small intestine compartments (Table I in S1 Appendix) according to rabbit physiology values found on literature. For most of the parameters, there was more than one value available from the literature. In most cases, the relative differences among the different sources were not large. The sources that provided information for as many tissues as possible were selected to achieve a consistent representation. Therefore, most organ volumes/specific blood flow rates were taken from [27], GIT-pH values were taken from [28], and GIT-proximal/distal radius and length values were taken from [29]. In the case of gastric emptying time, since the reported values in the literature had large deviations (i.e., from 20 min to 8 hrs), different values were chosen based on the experiment. Similarly, for the large/small intestinal transit rates, the full range of values was considered to set up the model (see Supporting information-GIT-transit emptying times, Tables G-H in S1 Appendix). For the parameter values not found in the literature, the mouse PBPK model parameters were transferred and used, as it is the closest animal species to the rabbit among the species available in the physiology database (e.g., mouse retains a gallbladder unlike the rat).

Validation process
The rabbit model was informed, tested and verified using published data for a range of compounds in which the clearance pathways were known from either the literature or from inhouse data (Table K in S1 Appendix). The choice of compounds took a hierarchical approach beginning with compounds cleared solely by either the kidney or the liver, which are the two main excretory organs, followed by the compounds cleared by both organs (i.v. administration). Subsequently, after assessing the physiology of the kidney and the liver, increasingly more complex clearance processes and administration protocols (p.o. administration) were simulated. For each simulation, the model prediction was first compared with the available data. In the case of deviations, either those parameters are known to vary among species or the parameters with no prior information were optimized. For all simulations, the partition coefficients and the cellular permeabilities were calculated based on the PK-Sim 1 Standard distribution model [22,30]. In short, partition coefficients describe the partition of a compound between water and lipids or water and proteins. They are calculated using a mechanistic formula similar to that described in [31]. The product of organ surface area and permeability controls the rate of permeation across lipid cell membranes from plasma into organs. Calculation of the permeabilities was based on the semi-empirical formula first published in [32] that accounts for transcellular as well as paracellular passive transport for small molecules. For large molecules, the two pores transport model is used [33] without accounting for FcRn binding.

Sensitivity analysis
To evaluate the sensitivity of our model to parameter perturbations, a local sensitivity analysis was performed. At every step of the sensitivity analysis, one parameter j of the model was varied by 10% of its nominal value while the remaining model parameters were kept constant. Next, the sensitivity coefficients (s i,j ) for the AUC or C max outputs were calculated as (Eq 1): where @f i is difference in PK parameter (AUC or C max ) between perturbed and un-perturbed simulation and @x j is the difference between varied and nominal parameter value. PK parameters of the sensitivity analysis were constrained to C max and AUC as the most relevant PK parameters for evaluating drug toxicity/efficiency. The sensitivity analysis was restricted to parameters of the rabbit model while parameters specific for compounds were not considered.

Model validation
First, a rabbit PBPK/TK model was developed, based on the pre-existing mammalian model structure and relative parameters retrieved from the literature or transferred from other mammalian PBPK/TK models. The rabbit PBPK/TK model structure with all considered organs and interactions is shown in Fig 2A. The details of the model, the structure and the parameters are listed in the supporting information. Second, to assess the reliability of the model, six different compounds (i.e., inulin, caffeine, ofloxacin, theophylline, paracetamol, and acyclovir) were modelled. Since inulin is solely filtered by the glomerulus of the kidney but is not secreted or reabsorbed by the tubulus or metabolized, it is well suited to validate the renal glomerular filtration rate (GFR) in the rabbit. To simulate renal excretion, the GFR was set to 100%, a value that is expected for renally cleared molecules. After determining the renal function of the rabbit, caffeine was tested since it is predominantly metabolized by the liver (S2 Fig). The metabolic processes were always localized in the liver as hepatic clearance, and the respective reaction kinetics were optimized based on the available data. Next, more complex compounds such as ofloxacin and theophylline were used (S3 and S4 Figs). To assess the rabbit GIT, oral administration was further added to test the model. The validation of oral administration was based on the compounds paracetamol, theophylline, and acyclovir. In the following sections, the detailed results for the compounds inulin and paracetamol are shown as they provide information for multiple tissues (inulin) and multiple oral formulations (paracetamol).The model results for the additional compounds are described in the supporting information and are only summarized here.  Table 1). For the plasma, skin and lung, the model was adjusted to better agree with the observed data. Here, the specific glomerular filtration rate (GFR specific ) was increased from 0.6 to 0.8 [l/min/kg-kidneyweight] [35] to better explain the observed exposure in the plasma. Similarly, the hydraulic conductivity of the skin was increased in accordance with the inulin simulations performed in rats [21,36]. Finally, to describe the PK profile in the lungs, the interstitial fraction of the lung was increased. These adaptions to the model are all in accordance with the literature. The simulation with the three adjusted parameters indicates a correct distribution of the compound among the different tissues of the rabbit (Fig 3-

Paracetamol case study
For paracetamol, the work of [38] reports three oral formulations: 50 mg of paracetamol administered as an oral solution, a rapidly disintegrating tablet, and a conventional tablet.   Table 1. Comparison of experimentally observed AUC, predicted AUC, and AUC after parameter optimization for available compartments, along with the relative errors, administration details and references. Inulin simulation efficiently described the data adopted from [37] and as such no further optimization was performed.  Fig 4). The simulations after parameter calibration explain the data for a solution and a rapidly disintegrating tablet. For the case of a conventional tablet, the data point towards a slower dissolution. Table 1 summarizes the model performance for a) a prediction based on the rabbit model solely informed by literature based values of rabbit physiology and compound related processes, as well as active transport parameters transferred from other species, and b) the simulations after the appropriate parameters were optimized given the observed data. The endpoint was the area under the curve (AUC) of the tissue concentration-time profiles. The table shows the observed, predicted and optimized AUC values along with the respective prediction errors. The mean fold error of the predicted AUC was 1.7, and with the exception of the tablet formulations, the AUC was always below 2. Similar results were obtained when the C max (maximum concentration) predictions and observations were compared (Table N in S1 Appendix).

Sensitivity analysis
A sensitivity analysis was performed to further investigate the relationship between the parameter values and the model output for different compounds and administration strategies (i.v. and p.o. administration). Fig 5 illustrates the relative change in the AUC and C max when the ten most sensitive parameters for each are changed. From green to yellow, the square of the sensitivity coefficient (Eq 1) increases denoting a higher sensitivity of the model output (AUC, C max ) to the respective parameter. The parameters listed are related to the distribution, clearance or uptake of the substance. For example, inulin is solely cleared by the GFR in the kidney and therefore, the kidney volume and the GFR (specific) are the most sensitive model parameters regarding inulin's AUC. Similarly, the AUC of caffeine, which is cleared by the liver, remains highly sensitive to the liver volume. On the other hand, C max appears to be more sensitive to variation of parameters that are related to either the route of administration or the distribution of the drugs. For example C max in paracetamol and acyclovir that are administered orally is more sensitive to gastric emptying time whereas in theophylline muscle volume plays a significant role, possibly as relevant volume of distribution. The sensitivity analysis was focused on rabbit-specific parameters and omitted compound specific parameters (e.g. V max , K m of active processes) as it is aimed to understand the impact of choosing alternative literature values of the rabbit physiology. An extended version of the sensitivity maps is shown in S11 and S12 Figs.

Discussion
PBPK/TK modelling is progressively gaining acceptance in risk assessment as it allows for the simulation of tissue-specific compound concentrations and as it is a suitable framework to translate in vitro experiments to in vivo exposures when considering ADME processes for relevant species [43][44][45]. In this work, a PBPK/TK model was developed for rabbits, and its predictive capacity was tested based on data from the current literature.
There are a number of reports, opinions and reviews that stress the significance of PBPK/ TK modelling in risk assessment and propose the best modelling practices [46][47][48][49][50][51][52][53]. Although there may be deviations in model development strategies based on the questions of interest, the first critical step in PBPK modelling is the review of existing datasets for the relevant species. The results from the rabbit physiology literature search are shown in Tables A-J in S1 Appendix, and the results include information about rabbit organ volumes, blood flow rates, GIT-compartments length/pH/radius/transit times and effective surface areas enhancement factors that represent the various folds, villi and microvilli that ultimately determine the absorption area in the GIT. After informing the model parameter describing the rabbit physiology, the model was validated by testing its predictive capacity on PK data from selected compounds. The validation was conducted to gain more knowledge by either re-estimating the relevant parameters and observing how the predictions improved or by evaluating the sensitivities of the physiological parameters to certain administration scenarios. The rabbit physiology was validated by compounds of gradually increasing complexity regarding clearance, so that any model inconsistency between the simulation and observations could be traced back to the clearance mechanism. For the compounds tested in this work, there were no evidence for different metabolic/transport pathways between rabbit and other species. However, in a PBPK setting, inter-species differences regarding active processes (e.g. an extra transporter) can be considered by re-evaluating compound-related parameters.
The first compound tested was inulin, which is a "gold standard" for measuring GFR [54] and thereby renal function. In our PBPK/TK model, the filtrated amount depends on the fraction unbound (fu) of the compound, the weight of the kidney, and GFR specific , which is the volume of the compound filtered per time per weight of the kidney (i.e., l/min/kg organ). These parameters in the rabbit PBPK/TK model were retrieved from the literature and included the deviation of the filtering capacity (GFR specific ) between individual animals [35]. Fig 3 shows the prediction of our model in dotted lines for six relevant rabbit tissues namely, plasma, lung, skin, bone, heart, and muscle. Interestingly, our model has an exceptionally high predictive capacity even before the experimental data are used to further optimize the model parameters. As Table 1 indicates, the % error of the PBPK/TK model predictions for this experiment is less than 70% (the mean value for all tissues is approximately 35%). To further inform the rabbit PBPK/TK model, additional data from the literature were considered, and the skin's hydraulic conductivity was recalibrated in accordance with the values observed in rats [36]. Finally, the interstitial fraction of the lungs (i.e., the interstitial volume of the lung over whole tissue volume) was re-estimated to better explain the data. In our panel of compounds, only concentration levels of inulin were observed in other tissues than blood plasma. This data was utilised to adjust parameters of the skin and the lung. As more data in those organs with other compounds become available, those parameters need to be re-evaluated.
In pharmacokinetics/toxicokinetics, there is no a priori threshold on an acceptable model error, but PBPK/TK models are generally accepted and considered useful when the prediction error on AUC or C max is in the 2-fold or 3-fold range [55][56][57]. In the rabbit model presented in this study, a mean fold error of 1.7 was calculated for the AUC among all compounds tested and simulated (Table 1). This result illustrates the high predictive power of our PBPK/TK model in i.v. and p.o. administration scenarios. Specifically, the validation included the kidney function through inulin i.v. administration (Fig 3, S1 Fig), liver function through caffeine i.v. administration (S2 Fig), and simultaneous kidney and liver function through ofloxacin and theophylline i.v. administration (S3-S5 Figs). Additionally, the model was tested on an oral administration of paracetamol through solution, rapidly disintegrating tablet, and conventional tablet formulations. The model predictions for the three formulations are shown in dotted lines in Fig 4. The predictions for the solution and the rapidly disintegrating tablet are in close accordance with the observed data, whereas the 'conventional' tablet prediction suggests a significantly different dissolution time (DT). This result is reasonable considering the available experimental data. For the case of solutions, the model assumes an instantaneous availability of the drug in the stomach compartment of the GIT, whereas for both formulations of tablet, a Weibull function was used to empirically describe the tablet dissolution process. For the case of the rapidly disintegrating tablet, the information for the dissolution time found in the published literature (i.e., 0.25 min) [38] was included. However, the dissolution time for the 'conventional' tablet formulation of paracetamol was not reported. In practice, dissolution times are experimentally accessible ex vivo and could be used in modelling, as done for the rapidly disintegrating tablet. However, the parameters were fit to the data in this study.
As a final step in the PBPK/TK model development, a local sensitivity analysis was performed. A sensitivity analysis is described as a crucial step in the best practices for modelling in risk assessment [47] and is recently used increasingly in numerous mechanism-based modelling efforts in PK/PD [58,59]. A sensitivity analysis not only evaluates the robustness of the PBPK/TK models but in a PBPK/TK setting where several literature values are available for certain physiologically-based parameters, it reveals whether the choice of alternative values would significantly impact the observed behaviour. The results of the sensitivity analysis output parameters are shown in Fig 5 and S11 Fig for AUC, and Fig 5 and S12 Fig for C max . As determined by the change in AUC, the most sensitive parameter was the liver volume followed by the kidney and venous blood volumes. This result is expected since most of the compounds tested were partially or totally cleared through metabolism in the liver. Therefore, a relative change in the liver mass or volume would highly affect the amount of drug cleared and as such the AUC. Similarly, for the compounds cleared renally, either through GFR (inulin) or tubular secretion and GFR (acyclovir), the parameters highly affecting the AUC are kidney volume and filtration capacity (GFR specific ). The parameters such as gastric emptying time, or GITlength and the effective surface area are sensitive only in scenarios involving oral administration (i.e., paracetamol PO, acyclovir PO). In contrast, gastric emptying time appears to be the most sensitive parameter regarding C max variation on compounds involving oral administration such as paracetamol and acyclovir. Furthermore, the muscle volume is the largest organ of the rabbit and the most important regarding distribution of drugs in the body hence it appears to be one of the most sensitive parameters. Overall, our analysis reveals sensitivities that are consistent with the biology and clearance pathways of each compound.
The workflow for using the rabbit PBTK model in risk assessment starts with identifying the physicochemical properties of a new active ingredient and its formulation. This step is followed by the use of in vitro assays to understand the active transport and metabolism. Ideally, before informing an initial PBTK model, TK data are generated in common species such as mice or rats in addition to the in vitro experiments. The translation of the toxicokinetics to the relevant species, such as rabbits and other wild animals, is then an exchange of the physiologyrelated parameters to the relevant physiology. The method can be used to predict the toxicokinetics of a new substance for different tissues and organs in rabbits. Based on our work, an average fold error of 1.7 can be expected for a compound of typical complexity.
In the present study, we demonstrated the development of a rabbit PBTK model based on data from the available literature. Specifically, a generic rabbit model based on well verified mammalian physiology was further informed by data from the available rabbit-related literature. The model yielded good TK predictions for different compounds with different exposure and clearance routes. By comparing predictions to measured toxicokinetic data, further species-specific parameters could be refined. Overall, this work presents a rabbit PBPK/TK model and its successful validation, and can additionally serve as a blueprint for the development of PBPK/TK models for additional species. The green dots are data adopted from the work of [42]. To describe the data, the gastric emptying time was increased to 4.91 hr from 0.5 hr, which is inside the bounds observed in the literature ( Table F in