A Deep Neural Network -- Mechanistic Hybrid Model to Predict
  Pharmacokinetics in Rat
An important aspect in the development of small molecules as drugs or
agro-chemicals is their systemic availability after intravenous and oral
administration. The prediction of the systemic availability from the chemical
structure of a potential candidate is highly desirable, as it allows to focus
the drug or agrochemical development on compounds with a favorable kinetic
profile. However, such pre-dictions are challenging as the availability is the
result of the complex interplay between molecular properties, biology and
physiology and training data is rare. In this work we improve the hybrid model
developed earlier [1]. We reduce the median fold change error for the total
oral exposure from 2.85 to 2.35 and for intravenous administration from 1.95 to
1.62. This is achieved by training on a larger data set, improving the neural
network architecture as well as the parametrization of mechanistic model.
Further, we extend our approach to predict additional endpoints and to handle
different covariates, like sex and dosage form. In contrast to a pure machine
learning model, our model is able to predict new end points on which it has not
been trained. We demonstrate this feature by predicting the exposure over the
first 24h, while the model has only been trained on the total exposure.
Florian Führer, Andrea Gruber, Holger Diedam, Andreas H. Göller, Stephan Menz, Sebastian Schneckener
# A Deep Neural Network - Mechanistic Hybrid Model to Predict Pharmacokinetics in Rat

###### Abstract

An important aspect in the development of small molecules as drugs or agrochemicals is their systemic availability after intravenous and oral administration. The prediction of the systemic availability from the chemical structure of a potential candidate is highly desirable, as it allows to focus the drug or agrochemical development on compounds with a favorable kinetic profile. However, such predictions are challenging as the availability is the result of the complex interplay between molecular properties, biology and physiology and training data is rare.

In this work we improve the hybrid model developed earlier [1]. We reduce the median fold change error for the total oral exposure from **2.85** to **2.35** and for intravenous administration from **1.95** to **1.62**. This is achieved by training on a larger data set, improving the neural network architecture as well as the parametrization of mechanistic model. Further, we extend our approach to predict additional endpoints and to handle different covariates, like sex and dosage form. In contrast to a pure machine learning model, our model is able to predict new end points on which it has not been trained. We demonstrate this feature by
predicting the exposure over the first 24h, while the model has only been trained on the total exposure.

Keywords:Hybrid modelling, Deep Learning, Property prediction, PBPK modelling, Drug design, Bioavailability, Pharmacokinetics

## 1 Introduction

Drug discovery is about the optimization of the interaction of molecules with biological targets to achieve the desired therapeutic effect, while reducing toxical effects. The same is true for developing compounds for applications in agriculture, with a large focus on the reduction of toxical effects in mammals. Both development processes are long and risky. Numerous methods and tools have therefore been established to support decisions in the search for best performing candidates. Experimental characterization of compounds can take up considerable resources and time. The targeted, early identification of favorable properties and consequently informed selection of compounds can significantly reduce development cycles and the associated costs. Selection criteria include both, pharmacological and toxicological effects, as well as pharmacokinetics (PK)1, in particular the availability of the compound in the body.

Footnote 1: Even though the name Pharmacokinetics implies that the field is only concerned with pharmaceutical substances, the field is concerned with all types of xenobiotic substances, see [https://en.wikipedia.org/wiki/Pharmacokinetics](https://en.wikipedia.org/wiki/Pharmacokinetics)

In this multi-parameter optimization of physicochemical properties, efficacy, safety and PK, many compounds are usually tested in different high-throughput assays to generate a basic understanding of a compound's characteristics. However, as PK is determined by the complex non-linear interplay of compound properties and physiology, using these assays to test and optimize all aspects and parameters relevant for PK is usually not possible. Therefore, animal studies remain an important contribution to understanding the PK characteristics of a potential drug candidate. However, animal studies are usually performed later in research for selected compounds that are already optimized with respect to the early accessible assays. This approach helps to keep the number of animal experiments low, but unfortunately often struggles from eventual limitations in further PK optimization.

The most import quantity in PK is the blood plasma concentration \(C\) as a function of time after a an oral (per os, PO) or intravenous (IV) administration. In this publication we are mainly interested in a few key parameters characterizing the concentration-time curve, such as the maximal concentration (\(C_{max}\)) and the exposure between two time points \(t_{1}\) and \(t_{2}\):

\[AUC_{t_{1},t_{2}}=\int_{t_{1}}^{t_{2}}dt\,C(t). \tag{1}\]

Most important is the total exposure, i.e. the exposure between the time of administration and infinity, here simply denoted as \(AUC\), sometimes also the exposure during the first 24h after administration AUC\({}_{24\text{h}}\) is considered. For pharmaceutical compounds,
oral drug delivery plays an important role. It represents the most common administration route and is convenient for patients and physicians leading to high patient compliance. The extent to which the systemic exposure of a drug after PO administration (AUC\({}_{\text{PO}}\)) differs from the exposure after intravenous (IV) administration (AUC\({}_{\text{IV}}\)) is quantified by the oral bioavailability \(F\) defined as

\[F=\frac{AUC_{PO}}{AUC_{IV}}\cdot\frac{D_{IV}}{D_{PO}}, \tag{2}\]

where \(D_{PO}\) denotes the oral dose and \(D_{IV}\) the IV dose. We would like to stress that AUC\({}_{\text{PO}}\), AUC\({}_{\text{IV}}\) and hence \(F\) depends on the compound as well as on the dose and the formulation. In general the \(AUC\) is a non-linear function of the dose, depending on e.g. the metabolic capacity or the availability of binding proteins. In addition, the AUC\({}_{\text{PO}}\) can have non-linear dose dependencies related to the oral absorption, e.g. due to a limited solubility or a transport mechanism. Typically, non-linearities in the dose dependence of AUC\({}_{\text{IV}}\) can be neglected, this allows to extrapolate the AUC\({}_{\text{IV}}\) from \(D_{IV}\) to \(D_{PO}\) hence it is sufficient to consider \(D_{PO}=D_{IV}\).

Furthermore, the oral dosage form can affect the PK, i.e. whether the compound is administered as solution, suspension or tablet. The AUC\({}_{\text{PO}}\) is typically lower for administration of a suspension or tablet than for a solution, since for a suspension or a tablet the particles first have to be released from the formulation and dissolve in the gastrointestinal tract (GIT). Note, that for tablets special so-called enabling formulations exist, which can increase the dissolution rate hence increase AUC\({}_{\text{PO}}\) and \(F\). Those formulations are usually not used in the early phases of drug and agrochemical development, hence we do not consider them in this publication.

While in pharmaceutical development a compound's systemic availability is desired to be as high as possible, for agrochemicals a compounds' it should be as low as possible to minimize the risks associated with safety and toxicities.

As the determination of the PK-parameters \(AUC\), \(C_{max}\) and \(F\) requires performing in-vivo studies in animals or even humans, they can not be used as a selection criterion for the early screening and optimization phases due to effort, cost and animal welfare considerations. Therefore, being able to predict them as early as possible, preferably directly from a compound's chemical structure, would reduce the risk during lead identification and optimization phases. Equally important, focusing only on the most promising compounds reduces the number of animal experiments.

There have been several attempts to predict PK-parameters from chemical structure [2; 3]. However, most of them are purely data-driven, hence do not exploit the available mechanistic knowledge about the different processes determining PK. In the present work, we combine Deep Learning with a mechanistic model to predict \(AUC_{PO}\), AUC\({}_{\text{IV}}\) and C\({}_{\text{max,PO}}\) in rats from the chemical structures only. Predictions for \(F\) can be calculated using (2) from the prediction for AUC\({}_{\text{PO}}\) and AUC\({}_{\text{IV}}\).

Our approach builds on the recent progress in applying Deep Learning to molecule property predictions [4; 5; 6; 7]. But, in contrast to these works, our data set is rather small with only a few thousand compounds. To compensate for this, we combine Deep Learning for property prediction with physiological based pharmacokinetic (PBPK) models.
PBPK models are well established mechanistic models, describing the kinetics of compounds in physiological environments [8, 9]. Doing so we benefit from our knowledge about rat physiology and the interplay of different processes, and make much more efficient use of the available data to learn relevant molecule characteristics.

## 2 Methods and materials

In this section we first describe our hybrid model and give a brief overview on PBPK models. Then we describe the training procedure for our model. Finally, we give an overview over the data used to train our hybrid model

### Hybrid modelling

To predict pharmacokinetics in rats we combine Deep Learning for molecular property prediction with PBPK models. A PBPK model is a system of ordinary differential equations (ODE) describing the PK processes a compound is undergoing within an organism. The processes are usually referred to as ADME processes, which stands for administration, distribution, metabolism and excretion. PBPK models are compartment models in which organs are represented by the compartments and the processes

Figure 1: Overview over our hybrid model structure consisting of a graph convolutional neural network for predicting a set of molecule properties. These molecule properties are the free parameters of a physiological model of rats predicting the pharmacokinetics. In practice, we approximate the PBPK model by a surrogate neural network.
are parametrized by physicochemical and other molecule properties. For readers not familiar with PBPK models we provide a brief overview in A.

In our approach, depicted in figure 1, the molecular properties predicted by a neural network, here called property net, correspond to parameters of the PBPK model, e.g. the solubility and the amount of substance cleared in the liver (hepatic clearance). The other parameters of the PBPK models are compound independent and describe the physiology of the organism, e.g., the organ volumes, blood flows, or specify the drug administration, e.g., administration route (IV/PO), dose and formulation.

The clear split between parameters describing the physiology and the molecule in our model highlights the expected advantages of our hybrid model in terms of required data and generalization. In our model the physiological parameters are fixed by the choice of target organism, so the neural network bridges between the molecular structure and the physiology. Furthermore, as certain aspects of the problem, e.g. the dose dependency, are mechanistically modelled, our model is able to exploit the extrapolation capabilities of the PBPK model, for example by generalizing to dosages outside the training range or predicting properties it has not been trained on, e.g. concentrations in different tissues.

Furthermore, we can exploit the flexibility of the property net to compensate for a misspecified or inaccurate mechanistic model. Two examples we consider in this publication, are the differences between male and female rats and different formulations used for the oral administration. We model those cases by a mechanistic model for male rats and solution as formulation, but let the molecule properties depend on these covariates by passing them to the (last layers) of the property net. By doing so, the property net learns to adapt its outputs such that even though we are using a wrong mechanistic model, we are still able to obtain an accurate model for female rats and suspensions.

### Physiologically based pharmacokinetic models

For the mechanistic model, we use the generic rat PBPK model available in the Open Systems Pharmacology (OSP) Suite [11] and add a generic hepatic (metabolism) and renal clearance (glomerular filtration in kidney) as well as a generic global P-glycoprotein (P-gp)-like active transport, which causes a flow from the inside to the outside of cells. For our purpose of predicting bioavailability in the early phases of drug and agrochemical development, it is sufficient to fix the physiological parameters to those of a typical rat, by using the OSP default values. The compound properties used as input for the PBPK model are listed in table 1. For oral administrations we assume a solution as formulation and account for differences between solution and suspension as well as differences between male and female rats by passing formulation and sex to the property net as described in section 2.1.

### Neural network architecture

Compared to the model developed in [1] we here replace the SMILES string representation of molecules and the corresponding 1D convolutional architecture with a graph
convolutional network (GCN) architecture directly acting on the graph representation of molecules. As SMILES representations are generated by a depth-first traversal through the molecular graph with an arbitrary starting node, they are not unique and connected sub-graphs are neither represented as contiguous sub-strings nor are they represented in same way when occurring in different molecules. In contrast, graph convolutional layers explicitly respect permutation invariance of the graph nodes and the connectivity of the graph. Hence, they do not suffer from the non-uniqueness and connectivity issues of SMILES based architectures. We therefore expect a GCN architecture to be superior to a SMILES based architecture. Indeed, we even have difficulties finding a set of hyperparameters for the SMILES based architecture, which

\begin{table}
\begin{tabular}{|p{56.9pt}|p{85.4pt}|p{85.4pt}|p{85.4pt}|} \hline
**Parameter** & **Short description** & **Pretraining data source** & **Distribution for surrogate model training** \\ \hline Hepatic clearance & Clearance rate in liver & In-vitro data (from hepatocyte assay) & Log-normal \\ \hline V\({}_{\text{max}}\) & P-gp like transporter & In-vitro data (Caco-assay) & Mixture of point mass at 0 and a half-normal \\ \hline GFR fraction & Fraction filtered in kidney & No pretraining & Uniform \\ \hline Fraction unbound & Fraction unbound in plasma & Predicted by independent DNN & Truncated normal \\ \hline Lipophilicity & Membrane affinity (log(MA)) & Predicted by independent DNN & Normal \\ \hline Effective molecular weight & Surrogate for molecule size & Calculated molecular weight excluding halogens & Half-normal \\ \hline Stomach solubility & Solubility in stomach & Predicted using Henderson-Hasselbach equation, DNN for solubility and pKa & Log-normal \\ \hline Small intestine & Solubility in small intestine & Predicted using Henderson-Hasselbach equation, DNN for solubility and pKa & Log-normal \\ \hline Large intestine & Solubility in large intestine & Predicted using Henderson-Hasselbach equation, DNN for solubility and pKa & Log-normal \\ \hline Small intestine permeation & Absorption rate in small intestine & Predicted from predicted log(MA) and molecular weight & Log-normal \\ \hline Large intestine permeation & Absorption rate in large intestine & Predicted from predicted log(MA) and molecular weight & Log-normal \\ \hline \end{tabular}
\end{table}
Table 1: Overview over the compound properties used as input of the PBPK model and whether predicted or in-vitro observed values are used for pertaining. Details on the used prediction models can be found in [10].
result in stable pre-training and acceptable accuracy, while we easily found those hyperparameters for the GCN.

We use the GCN-architecture proposed in [12] implemented in deepchem [13]. This GCN-architecture uses differentiable operations inspired by those used to calculate circular fingerprints and equips them with learnable weights.

### Model training

#### 2.4.1 Surrogate

For end-to-end training of the hybrid model, we need to back-propagate through the PBPK model. Even though this is possible for small ODE systems, it is computationally prohibitive for our model with about 300 stiff ODEs. Therefore, we replace the PBPK model by a surrogate neural network approximating the PBPK model. Here, we use a fully connected neural network. We train the surrogate model on 2.4M simulations with random model parameters, sampled using latin hyper-cube sampling, and test it on additional 0.72M randomly sampled simulations. Each parameter is distributed according to a simple parametric distribution, e.g. normal or log-normal, roughly matching the distribution of values in our training data. The functional form of the distributions are summarized in table 1. We increase the variances of the distributions by 50 % to avoid the predictions of the property net leaving the training range of the surrogate. As can be seen from figure 2, the surrogate is able to reproduce the PBPK model accurately, with a fold change error of about 1.04 for the two PO endpoints and 1.02 for the IV endpoint. The error of the surrogate is negligible compared to the expected error caused by the high biological variability of about 50% of the in-vivo data.

To be able to reliably back-propagate through the surrogate, good point wise approximations of the PBPK model are not sufficient. Also, the surrogate gradients and ideally higher order derivatives need to be good approximations of the PBPK models derivatives, i.e. the surrogate needs to reproduce the response of the PBPK model to changes in the molecule properties. We confirm this qualitatively by randomly generating a set of points in the PBPK input space and then vary each molecule property individually, while holding the others fixed, some examples are shown in figure 3. Overall, we find good agreement between the curves predicted by the surrogate and the PBPK model.

Figure 2: Overall accuracy of surrogate neural evaluated on a hold out test set of simulations. A median fold change error of 2-4% is small to the expected biological variability of the data of 50%.
#### Training strategy

To overcome the small data set of about 7000 compounds, we pre-train the property net on molecule properties of about 140k compounds. As this is only a pretraining step, we use predicted properties where available and measured values otherwise. Table 1 indicates for which properties predicted and for which measured values are used.

Figure 3: Some examples showing the full PBPK model and the surrogate model as a function of a single model parameter, while keeping the others fixed. In these examples the dependence on the hepatic clearance (top left) and dose (bottom right) is very accurately described by the surrogate. In the GFR example (top right) the surrogate is able to reproduce the shape of the PBPK, but shows a constant offset of about 20%, which is acceptable given the variability of the PK-data. The solubility example (bottom left) shows an offset of similar size, but is able to qualitatively reproduce the step-like behavior seen in the PBPK model. The small oscillations of about few % do not introduce major problems during training of the hybrid model.
Details on the pretraining data set are given in section 2.5. The pretrained property net is then trained end-to-end, as part of the hybrid model, on the in-vivo data to predict the target PK-parameters.

To constrain the model parameters predicted by the property net to physiological values, e.g. a non-negative clearance, and to the range of the surrogate training data, we add a penalty term to the loss function \(L_{total}\)

\[L_{total}=\sum_{i}L\left(y_{i},\widehat{y_{i}}\right)+\lambda\sum_{i,j}max \left(p_{j,i}-p_{max,j},0\right)^{2}+max\left(p_{min,j}-p_{j,i},0\right)^{2}. \tag{3}\]

The first sum is over all data points \(i\), the second over all data points \(i\) and all molecule properties indexed by \(j\). \(p_{i,j}\) denotes all predicted molecule properties for all data points. The penalty is zero as long as \(p_{min,j}\leq p_{i,j}\leq p_{max,j}\) and positive otherwise. When during training the penalty term is not decreasing for several epochs, we increase the weight \(\lambda\) until a pre-defined tolerance, here \(10^{-8}\), is reached. Empirically, we find that this is sufficient to constrain the PBPK model parameters to the viable range, see section B.

### Data

We retrieved all in-vivo data taken after PO or IV administration in Wistar rats from the Bayer data warehouse. After filtering out pro-drugs, salts, molecules heavier than \(1500\,g/mol\) and non-standard formulations we are left with 7192 compounds, with in total 5731 \(AUC_{PO}\), 6183 \(C_{max,PO}\) and 6408 \(AUC_{IV}\) measurements. In contrast to [1], in addition to male rats and PO administrations as solution we consider also female rats and suspensions. Furthermore, in contrast to our previous work [1], we do not restrict the dose, therefore, our data covers a dose range from \(0.0024\,mg/kg\) to \(1000\,mg/kg\). Overall, there is more data available for low dosages than for high dosages, see figure 5. The large dose range reflects the fact that our data set includes relative low dose (\(\sim 1\,mg/kg\)) data mainly taken from male rats at Bayer Pharmaceuticals, as well as high dose data (\(\gtrsim 10\,mg/kg\)) mainly taken from female rats at Bayer

Figure 4: Number of data points for different sub-sets of the data set. Since the standard test for pharmaceuticals is on male rats using, for PO, a solution, most of our compounds are tested on male rats.
CropScience. Furthermore, the compounds from both divisions are expected to have different properties. This increases the diversity in our data set. We expect that this results in a better generalization of the model. To challenge the capabilities of our hybrid model to generalize to new observables we also collect available AUC\({}_{\rm 24h}\) data.

For pretraining, we use about 100k compounds from the Bayer data warehouse. We use our internally available models to predict solubility, pKa values, lipophilicity and plasma protein binding in human for all compounds. For the hepatic clearance and membrane permeation no model is available, so we use all available in-vitro measurements resulting in an additional 40k compounds for pretraining. Note that usually no urine data is collected, so neither data nor a model is available for the GFR fraction, hence the GFR fraction is only trained end-to-end. In total, we use about 140k compounds for pre-training.

## 3 Results

In this section we validate the predictive performance of our model and compare it to a standard GCN. We further challenge the generalization capabilities of our hybrid model by using it to predict the AUC\({}_{\rm 24h}\), a quantity the model has not been trained on.

### Model performance

We optimize the hyperparameters of our hybrid model using the HORD algorithm [15]. The model architecture is optimized on the pretraining set while training hyperparameters are optimized on the training set. We validate the best model on our 20% hold out test set in figure 6. For evaluation, we use the median fold change error defined as:

\[mfce=\ exp\left(median\left|log\left(y\right)-log\left(\hat{y}\right)\right| \right), \tag{4}\]

such that for a perfect model \(mfce=1\). A fold change error of 2 to 3 is considered sufficient to inform compound selection [16, 17, 18, 19, 20, 21, 22]. For all targets, except for the AUC\({}_{\rm PO}\) for male rats and suspension, we reach this goal. For AUC\({}_{\rm IV}\) and \(F\) our model

Figure 5: Distribution of used dose in PO (left) and IV (right) measurements. High doses are typically tested only in PO experiments, hence they span much large dose range then the iv experiments.
even achieves \(mfce<2\). Compared to previous work [1], which uses a slightly different test set, the \(mfce\) of AUC predictions for male rats for PO (solution) and IV has improved from 2.85 to 2.35 (PO) and from 1.95 to 1.62 (IV). Also the \(mfce\) of \(F\) predictions improved from 1.83 to 1.62. Additionally, we observe a more stable training and easier to tune hyperparameters. \(\mathrm{C_{max,PO}}\) after oral administration, which has not been considered in [1], can be predicted with a slightly higher accuracy than the AUC. Even though there is less data available for female rats than for male rats, predictions for female rats after an IV administration or a PO administration using solution can be made with similar accuracy as for male rats. We observe that PO predictions for suspensions are less accurate than predictions for solutions. This is expected, given that the dissolution of a suspension adds complexity to the dynamics in the GIT. Likewise, predictions for IV are more accurate than predictions for PO. As described in section 2.1 we predict suspension by using a mechanistic model for solutions with adapted molecule properties, we expect that using a mechanistic model for suspension, the prediction accuracy for suspension can be improved.

Figure 6: Model predictions vs observed values for our three training tasks \(\mathrm{AUC_{PO}}\) (top left), \(\mathrm{AUC_{IV}}\) (top right) and \(\mathrm{C_{max,PO}}\) (bottom left). For comparison with earlier work [1; 14] we also show the derived predictions for \(F\) (bottom). The predictions for \(\mathrm{AUC_{IV}}\) are more accurate than the PO predictions, this is expected since the processes involved in a PO administration are more complex than those involved in an IV administration. We observe that despite having less data for female rats than for male rats the predictions have a similar accuracy, while predictions for suspensions are a bit less accurate than those for solutions.
### Advantages of hybrid models

Figure 7 shows for comparison the performance of a standard GCN, having the same architecture as the property net except for the output layer's size. The predictions of the GCN are for all three endpoints worse than those of the hybrid model. While the performance drop of the AUC\({}_{\text{IV}}\) predictions is moderate, the performance drop of AUC\({}_{\text{PO}}\) and C\({}_{\text{max,PO}}\) is of practical relevance as the standard GCN does not reach \(mfce<3\).

In addition to the improved performance of the hybrid model, compared to a pure Deep Learning model, we can expect that the hybrid model is able to extrapolate and predict target parameters on which it has not been trained. For a first assessment of the extrapolation capability of our hybrid model we use the AUC\({}_{\text{24h}}\). Figure 8 shows the AUC\({}_{\text{24h}}\) predictions of our hybrid model compared to the observed values. The accuracy of the AUC\({}_{\text{24h}}\) predictions are comparable to the endpoints the model was trained on. Note that the AUC\({}_{\text{24h}}\) are predicted using the full PBPK model instead of the surrogate. The high predictive accuracy reconfirms that our surrogate is an accurate approximation of the PBPK model. This is further confirmed by the predictions of our training targets using the PBPK model instead of the surrogate model in section C.

Figure 8: Predicted AUC\({}_{\text{24h}}\) compared to the observed values. The AUC\({}_{\text{24h}}\) are predicted from the molecule properties predicted by the property net, using the exact PBPK model. Despite not being trained on the AUC\({}_{\text{24h}}\) our hybrid model achieves an accuracy comparable to the total \(AUC\).

Figure 7: Model accuracy for a pure GCN model for the same three predictions tasks AUC\({}_{\text{PO}}\) (top), AUC\({}_{\text{IV}}\) (mid) and C\({}_{\text{max,PO}}\) (bottom) as for the hybrid model. The accuracy of all 3 end-points is higher for the hybrid model than for the pure deep learning model, see figure 2.
## 4 Summary and conclusion

In this work we present a hybrid model to predict the pharmacokinetics of pharmaceutical and agrochemical compounds bioavailability in rats directly from chemical structure. As predicting in-vivo targets is challenging due to the complex non-linear interplay of many processes as well as the low amount and high variability of data. We tackled these challenges by combining expert knowledge about rat physiology and processes affecting pharmacokinetics with Deep Learning for molecular property predictions.

The work in [1] was extended by using a GCN for the property net and changing the parametrization of the mechanistic model which improves the interface between the neural network and the mechanistic model. Additionally, we increased the available training data set by less restrictive filtering of the data and the inclusion of additional training endpoints, such as \(\mathrm{C_{max,PO}}\) or predictions for female rats and suspensions. Furthermore, we used a larger internal data set for pretraining, which is expected to be better correlated to bioavailability than the previously used public data from the TOX21 challenge [23]. That lead to an improved performance of the model compared to [1].

For all except one end-point our model has an \(mfce<3\). Our \(AUC_{IV}\) and \(F\) predictions even have an mfce \(<2\). This is expected to be accurate enough to inform decisions during early phases of drug discovery [16, 17, 18, 19, 20, 21, 22]. Furthermore, our model enables the selection and prioritization of compounds which are directly optimized with respect to their pharmacokinetic profile [24, 25, 26].

Our prediction accuracy is competitive to the \(F\) prediction accuracy in [14], when their model is only trained on in-vivo data from the chemical series it is applied to and superior otherwise. Our model shows a similar performance to the different models in [27], which seem to have a slightly higher accuracy, but to achieve this, the chemical structure as well as in-vitro parameters are required, whereas our approach does not rely on in-vitro measurements. We like to stress that such a comparison should not be over interpreted, as different data-sets are used for training and validating the different models.

Additionally, our approach is able to handle different covariates like sex and formulation by predicting effective molecule properties, which are sex and formulation dependent. We expect that prediction the accuracy can be improved by accounting for these covariates in the mechanistic model. The inclusion of further covariates like body weight is therefore likely to result in even better predictions. However, the incorporation of more covariates is limited in practice as the typically available covariates do not fully specify the physiology of an individual. To account for the residual variability and to possibly improve predictions further probabilistic models which estimate population parameter distributions are required. To do so one could build on the recent progress in building deep generative models [28, 29, 30].

Incorporating more covariates - either deterministic or probabilistic - requires more input parameters to the mechanistic model, which complicates the use of a surrogate neural network for the PBPK model. Training and using a surrogate worked well for the small number of inputs and outputs considered in this work, but becomes harder if their number increases. In such cases the use of a full PBPK model can become
superior. However, this is currently computationally infeasible. Recently, there has been increasing interest in combining differential equations and Deep Learning [31, 32, 33] and, consequently, in tools to train these models [34, 35, 36], such that one can expect using the full PBPK model will become feasible in the near future. A complementary approach is to use simpler and hence computationally cheaper PK models, such as compartmental models or a reduced version of the full PBPK model.

Using PBPK models directly would also alleviate the need for constraining the molecule properties to the validity range of the surrogate. Constraining the molecule properties by introducing a penalty term in the loss worked in our case, but still complicated model training. Using a PBPK model directly would also enable to train the hybrid model on concentration-time profiles, which would be highly desirable, since more accurate predictions of concentration time profiles would allow a much more detailed description of the pharmacokinetics of a compound.

Successful application of PK models not only depends on the prediction accuracy, but also on the possibility to estimate uncertainty on the prediction. Such, that decisions based on predictions are only made for molecules for which the model is expected to be accurate. Ref. [27] assesses prediction uncertainty. For none of the tested approaches the uncertainty estimates are fully satisfactory. But, approaches with statistical correct Bayesian or Frequentist epistemic uncertainty provide better uncertainty estimates. We expect that in both cases well calibrated uncertainties can be provided by computationally expensive ensembles techniques [37, 38, 39].

Our model has the potential to reduce cost, development time and animal experiments in drug and agrochemical research by focusing the development on the most promising candidates and being able to directly optimize a compounds PK. Furthermore, our approach can be used to predict human PK [40], therefore directly optimizing for clinical use.

## Declaration of Competing Interest

The authors declare no competing financial interest.

## Acknowledgments

### Appendix A Physiologically based pharmacokinetic models

Physiologically based pharmacokinetic (PBPK) models are ordinary differential equation models describing how a substance, e.g. a drug, is absorbed, distributed, metabolized, and excreted in an organism. For the reader not familiar with PBPK models we provide a brief overview over the basic concepts, building blocks and equations forming a PBPK model. For more details we refer to [41].

In PBPK models physiological organs and tissues are represented by compartments. The transport of substance via the blood is modeled by balance equations of
the form

\[\frac{dC_{i}}{dt}=\frac{Q_{i}}{V_{i}}\left(C_{art}-\frac{C_{i}}{P_{i}}\right), \tag{1}\]

where \(C_{i}\) denotes the compound concentration in the compartment \(i\), \(V_{i}\) its volume, \(Q_{i}\) the blood flow, \(P_{i}\) the partition coefficient between blood and tissue, and \(C_{art}\) the compound concentration in arterial blood, which is governed by

\[\frac{dC_{art}}{dt}=-\sum_{i}\frac{Q_{i}}{V_{i}}\left(C_{art}-\frac{C_{i}}{P_{i }}\right)\!, \tag{2}\]

To describe dissolution, absorption, metabolism and excretion, as well as additional distribution mechanism the equations 1 and 2 need to be extended. For example, dissolution and absorption in a single GIT compartment is described by the following equations:

\[\frac{dC_{g}}{dt} =\frac{Q_{g}}{V_{g}}\left(C_{art}-\frac{C_{g}}{P_{g}}\right)+K_{a} C_{lum}, \tag{3}\] \[\frac{dC_{lum}}{dt} =-K_{a}C_{lum}+\frac{dC_{dis}}{dt},\] (4) \[\frac{dC_{dis}}{dt} =K\left(C_{0}-C_{dis}\right)^{2/3}\left(C_{s}-C_{lum}\right), \tag{5}\]

Equation 3 describes concentration in the GIT tissue \(C_{g}\), which is sourced by a linear absorption process from the GIT lumen. Equation 4 is describes the compound concentration in the GIT lumen \(C_{lum}\), which is sourced by the dissolved compound \(C_{dis}\). Equation 5 is the Noyse-Whitney equation describing the dissolution of the compound in the GIT lumen, with \(K\) being a compound dependent constant, \(C_{0}\) is the total amount of compound administered divided by the administered volume and \(C_{s}\) is the solubility, i.e. the compound concentration the GIT lumen at (thermal) equilibrium. Metabolism is described by the Michaelis-Menten-Kinetics, which for \(C\ll K_{m}\) can be linearized:

\[\frac{dC}{dt}=-V_{max}\frac{C}{K_{M}+C}=-\frac{V_{max}}{K_{M}}C+O\left(\left( \frac{C}{K_{m}}\right)^{2}\right), \tag{6}\]

The constants \(V_{max}\) and \(K_{M}\) depend on the compound and the metabolizing enzyme and control the speed and saturation of metabolism. We assume a single generic metabolizing enzyme, hence in our hybrid model hepatic clearance is fully characterized by the rate \(\frac{V_{max}}{K_{M}}\).

An active P-gp like transport via membrane proteins, assuming a constant protein concentration, follows also a Michaelis-Menten-Kinetics

\[\frac{dC_{1}}{dt} =-V_{max}\frac{C_{1}}{K_{M}+C_{1}} \tag{7}\] \[\frac{dC_{2}}{dt} =V_{max}\frac{C_{1}}{K_{M}+C_{1}}. \tag{8}\]
As for the metabolism, the constants \(V_{max}\) and \(K_{M}\) control the speed and saturation of the transport are compound and are transport protein dependent. For our purpose it is sufficient to set \(K_{M}=1\frac{\mu\text{mol}}{\text{L}}\), i.e. use the OSP default value, hence the transport is parametrized by its maximal velocity \(V_{max}\).

## Appendix B Validation of property constraints

In figure B1 the distribution of predicted molecule properties of the test set are shown together with the maximal and minimal values in the surrogate training data set. All predicted molecule properties lie within in the surrogates training range, confirming the effectiveness of the penalized loss described in section 2.4.2. Note that for \(V_{max}\) and \(FU\) we used heavy tailed distributions for generating the surrogate training data, resulting in the large range shown in figure B1. For the \(FU\) this results in unphysiological values \(>1\), for which the equations of the PBPK model are still defined. But
in practice the property net does not predict a \(FU>1\). Furthermore, to increase the flexibility of our clearance model we increased the maximal allowed value for the GFR fraction from 1 to 5.25.

## Appendix C A posteriori surrogate validation

We can validate the surrogate model a posteriori by predicting the training targets of our hybrid model using the PBPK model instead of the surrogate. Figure C2 shows the predictions obtained using the PBPK model vs those obtained using the surrogate. The accuracy is not as good as expected from the analysis in section 2.4.1, but still accurate enough to be used in the hybrid model, the \(mfce\) of the surrogate (\(1.2-1.4\)) is clearly better than the \(mfce\) of the hybrid model (\(mfce\gtrsim 1.6\)). Additionally, figure C3 shows the predictions using the full PBPK vs the observed values. These predictions are almost as accurate as those using the surrogate model. A maximal difference of 0.24 in the \(mfce\) can be observed, and no additional features are visible. This highlights again the accuracy of the used surrogate model.
