M. C.
Wenlock
InSilicoLynx Ltd, BioHub at Alderley Park, Mereside, Alderley Park, Cheshire, SK10 4TG, UK. E-mail: mark.wenlock@insilicolynx.com
First published on 5th February 2018
Assessing the oral drug suitability of compounds as early as possible within drug discovery is an important objective. This study describes a methodology that attempts to simplify the evaluation of compounds based on their in vivo quantity levels within a mammalian body, represented using a mathematical model that imposes a time limitation on oral absorption and assumes non-instantaneous drug distribution between plasma and tissue. This simplification results in two new oral drug suitability parameters that can quantitatively relate oral dose to in vivo exposure for compounds with vastly different tendencies in terms of absorption into, and elimination from, the body. Consequently, the complexities associated with evaluating a compound's oral drug suitability are simplified to an assessment of these two new parameters. Application of this methodology at the virtual design stage is discussed, along with functionality that accounts for uncertainty related to a compound's distribution kinetics and errors associated to in silico QSAR predictions for the required input data.
The octanol–water partitioning system from which a logP is derived, and on which a predicted log
P is based, is intended as a model system for compound transfer between aqueous and lipid phases and is arguably too simplistic for a whole mammalian system. This study addresses this point by proposing alternative oral drug suitability parameters for consideration during virtual drug design.
The oral drug suitability parameters are derived from rate equations for determining in vivo quantity levels associated to a relevant mathematical model of a mammalian system. The focus of this study is a mathematical model for a human system in the form of a three-compartment pharmacokinetic (PK) model, which imposes a time limitation on the absorption of an orally dosed drug from the small intestines, assumes non-instantaneous drug distribution between plasma and tissue, and only allows for elimination of the absorbed drug from one compartment.11,12 Application of the proposed PK model in a repeat-dosing simulation can approximate steady-state in vivo data for a quantity (e.g., maximum plasma concentration) as a function of the oral dose, resulting in a dose–quantity curve (Fig. 1).
Each dose–quantity curve has two compound-specific features, which are treated as a pair of oral drug suitability parameters: the first is the dose–quantity intercept (DQI) for the regression line fitted through the curve's linear PK region; the second is the highest dose that still leads to linear PK (or near-linear PK) for the particular quantity (highest linear dose [HLD]). This study describes the methodology for calculating a compound's DQI and HLD, and demonstrates how to interpret the results to assess a compound's oral drug suitability. Importantly, the PK model is not intended to be physiologically accurate, rather it is a minimum model system that permits consideration of key factors related to a drug's in vivo exposure. Hence, the DQI and HLD contain information on the absorption, distribution and elimination properties of a compound based on this model system, which can facilitate the ranking of compounds. However, the calculated DQI and HLD values may differ from an in vivo measurement.
This study considers the application of the aforementioned PK model to 15 known oral drugs. The size of the set was restricted to 15 because of the limited number of compounds (within the literature and the European Bioinformatics Institute's ChEMBL database13) with experimentally derived values for pKa, aqueous solubility at pH 7.4 (and room temperature, solubilitypH7.4), apparent human Caco2 membrane permeability (A to B) at pH 6.5 (Papp,Caco2,pH6.5), human volume of distribution at steady state (Vss), human in vivo plasma clearance (Cl), human in vivo efficacy and, optionally, human plasma protein binding (PPB). However, it is envisaged that the DQI and HLD should be used at the virtual design stage and that such compound information would be generated using in silico QSAR models. Such methods have associated prediction errors, which can be large, so the impact of these errors on the DQI and HLD values is discussed.
Drug | Charge type | pKa | SolubilitypH7.4 (M) | P app,Caco,pH6.5 (cm s−1) | V ss (L kg−1) | Cl (mL min−1 kg−1) | PPB (% bound) |
---|---|---|---|---|---|---|---|
Chlorpromazine | Monobase | 9.24 (ref. 14 and 15) | 3.72 × 10−4 (ref. 14, 15, 21 and 22) | 2.00 × 10−5 (ref. 26) | 10.00 (ref. 27) | 16.00 (ref. 29) | 98.64 (ref. 22) |
Diazepam | Neutral | 1.57 × 10−4 (ref. 22 and 23) | 6.03 × 10−5 (ref. 26) | 1.00 (ref. 27) | 0.38 (ref. 29) | 98.61 (ref. 22) | |
Diclofenac | Monoacid | 4.26 (ref. 15–18) | 6.53 × 10−3 (ref. 15, 16, 18, 23 and 24) | 5.37 × 10−5 (ref. 26) | 0.19 (ref. 27 and 28) | 3.83 (ref. 28 and 29) | 99.68 (ref. 22) |
Furosemide | Monoacid | 3.94 (ref. 14, 17 and 19) | 1.78 × 10−1 (ref. 14) | 2.82 × 10−7 (ref. 26) | 0.11 (ref. 27 and 28) | 2.19 (ref. 28 and 29) | 98.37 (ref. 22) |
Haloperidol | Monobase | 8.63 (ref. 17 and 18) | 7.56 × 10−5 (ref. 18, 21 and 22) | 1.58 × 10−5 (ref. 26) | 17.00 (ref. 27) | 7.80 (ref. 29) | 85.48 (ref. 22) |
Imipramine | Monobase | 9.50 (ref. 14 and 17) | 9.85 × 10−3 (ref. 14) | 2.75 × 10−5 (ref. 26) | 12.00 (ref. 27) | 13.00 (ref. 29) | 87.87 (ref. 22) |
Indomethacin | Monoacid | 4.31 (ref. 15, 17 and 20) | 1.44 × 10−3 (ref. 15, 20 and 25) | 2.95 × 10−5 (ref. 26) | 0.16 (ref. 27 and 28) | 1.61 (ref. 28 and 29) | |
Ketoprofen | Monoacid | 4.25 (ref. 15–17 and 19) | 5.55 × 10−1 (ref. 15 and 16) | 4.68 × 10−5 (ref. 26) | 0.13 (ref. 27) | 1.60 (ref. 29) | 98.73 (ref. 22) |
Naproxen | Monoacid | 4.23 (ref. 15–19) | 9.87 × 10−2 (ref. 15, 16, 18 and 23) | 4.68 × 10−5 (ref. 26) | 0.09 (ref. 28) | 0.07 (ref. 28) | |
Nifedipine | Neutral | 5.23 × 10−5 (ref. 18, 21, 22 and 24) | 3.24 × 10−5 (ref. 26) | 0.79 (ref. 27) | 7.30 (ref. 29) | 96.17 (ref. 22) | |
Phenytoin | Monoacid | 8.06 (ref. 15–18) | 8.51 × 10−5 (ref. 15, 16, 18, 21, 22 and 24) | 3.98 × 10−5 (ref. 26) | 1.40 (ref. 28) | 0.57 (ref. 28) | 87.11 (ref. 22) |
Pindolol | Monobase | 9.54 (ref. 15) | 2.62 × 10−2 (ref. 15) | 3.89 × 10−5 (ref. 26) | 1.20 (ref. 27) | 7.70 (ref. 29) | |
Prazosin | Monobase | 6.50 (ref. 17) | 1.48 × 10−5 (ref. 22) | 4.07 × 10−6 (ref. 26) | 0.73 (ref. 27) | 4.70 (ref. 29) | 96.00 (ref. 22) |
Trimethoprim | Monobase | 7.26 (ref. 17) | 2.17 × 10−3 (ref. 23) | 8.71 × 10−5 (ref. 26) | 1.50 (ref. 27) | 2.10 (ref. 29) | 66.10 (ref. 22) |
Warfarin | Monoacid | 4.97 (ref. 14, 15, 17 and 18) | 3.35 × 10−3 (ref. 14, 15 and 18) | 4.37 × 10−5 (ref. 26) | 0.12 (ref. 27 and 28) | 0.05 (ref. 28 and 29) | 99.31 (ref. 22) |
The application of the rate equations, along with the input data required to calculate compound compartment levels at different time points, constituted a simulation. Simulations were run at a series of doses and included multiple distribution kinetics scenarios. An open-source Python (software) library was written to perform and analyse these simulations.31
To establish steady-state conditions, simulations considered repeat dosing. Using the principle of superposition, compound compartment levels were considered for each repeated dose from time of dose to end of the simulation, and then all levels at each time point were summed for the repeated doses. Compartment A represented a hypothetical uniform cylindrical intestinal segment, filled with a constant intestinal fluid volume (Vintestinal, L) of 0.08 L,32 into which the compound (in solid form) was delivered as a bolus; for this study, the intestinal fluid was represented by aqueous pH 6.5 buffer. The uniform mixing and dissolution process within compartment A was treated as instantaneous and only solubilised compound could undergo the one-way transfer into compartment B; this transfer was treated as either a zero-order or a first-order process, characterised by the absorption rate constant k1 (min−1), similar to the approach used by Dressman et al.12 The model restricted this process to an absorption window, during which compound could only be removed from compartment A by said absorption; any remaining compound at the end of the absorption window was discarded, and the model reverted to a classic two-compartment PK model.33 For clarity, the absorption window represented the transit of this uniform cylindrical slug of intestinal fluid (containing the dosed drug) through the small intestines, during which its shape and volume remained unchanged.
Compartment B (the central compartment) represented the body spaces into which a compound could distribute extremely rapidly (i.e., plasma and well-perfused tissues, including the major eliminating organs). The volume of compartment B was compound specific and represented the compound's initial dilution volume (Vcentral, L). Compartment C (the peripheral compartment) represented the body spaces into which a compound distributed more slowly. Compound transfer from compartment B to compartment C was treated as a first-order process, characterised by the rate constant k2 (min−1), and the reverse process was characterised by the rate constant k3 (min−1). It was assumed that, once absorbed, compound could only be removed permanently from the body from compartment B via a first-order process, characterised by the rate constant k4 (min−1).
The solubilitypH6.5 was calculated by multiplying a compound's solubilitypH7.4 by the ratio of its fraction neutral at pH 7.4 to its fraction neutral at pH 6.5; the fraction neutral values were calculated using the method described by Wenlock.37
k1 = (SA·Peff,human,pH6.5·60)/Vintestinal | (1) |
P eff,human,pH6.5 was estimated based on a compound's Papp,Caco2,pH6.5 using a similar approach to that previously described.37,41 Critically, this approach focuses on the establishment of a (log10–log10) linear regression equation for a set of compounds between their neutral species human membrane permeability at pH 6.5 (Pm,human,neutral,pH6.5, cm s−1) and their neutral species human Caco2 membrane permeability (A to B) at pH 6.5 (Pm,Caco2,neutral,pH6.5, cm s−1). For this study, a larger set of 32 compounds was used to establish the following regression equation:
log10(Pm,human,neutral,pH6.5) = 0.916·log10(Pm,Caco2,neutral,pH6.5) + 1.579 | (2) |
Fh = 1 − Clb/Qh | (3) |
For each scenarios, values for k2, k3 and k4 were derived using the values for Vcentral, Vss, Vterminal and Cl, in conjunction with rearrangement of standard mathematical equations associated to a two-compartment model.33 Details of these equations can be found in the ESI.†
Five simulations were performed for each compound and each dose, based on the five scenarios.
The quantities calculated for each simulation were collated and used to create a (log10–log10) dose–quantity curve for each distribution kinetics scenario.
Fig. 1 highlights a typical dose–quantity curve resulting from the application of the PK model. It is characterised by a linear region with a slope of 1 for doses below a certain limit. As the dose increases above this limit, compartment A becomes increasingly saturated for the duration of the absorption window. This results in the dose–quantity curve bending as the dose increases, eventually plateauing. For the linear region, the corresponding regression equation takes the form of:
log10(quantity) = 1·log10(dose) + DQI | (4) |
Rearrangement of eqn (4) leads to eqn (5) or (6):
quantity = dose·10DQI | (5) |
![]() | (6) |
However, these equations only apply up to a log10(dose) equal to the HLD; this equates to the extrapolated log10(dose), using eqn (4), for the log10(max_quantity) observed for a particular distribution kinetics scenario (Fig. 1). The HLD value arrived at using this method is not strictly the HLD, as the dose–quantity curve has already begun to deviate from linearity at this log10(dose) (Fig. 1); rather, it is an approximation using a standardised approach.
Importantly, the DQI and HLD provide a simple way to understand a compound's dose–quantity curve. A plot of DQI against HLD provides an alternative scale against which to evaluate different compounds.
Other PK models based on this linear differential equation approach are possible – rate equations for some alternative PK models can be found in the ESI.† The simplest of these is a one-compartment intravenous (bolus dosing) model that only considers a first-order rate constant for elimination.43 This model can be extended to include an oral absorption step – the approached used by Wenlock and Page.37,42 Such a model is applicable for compounds where oral absorption cannot be assumed to be instantaneous; it can also be refined to account for time-limited oral absorption, and details of this can be found in the ESI.† A limitation of one-compartment models is the assumption of instantaneous drug distribution between plasma and tissue; this is overcome in the present PK model by use of a second compartment that permits distribution kinetics to be modelled. For reference, details of a two-compartment intravenous model can be found in the ESI.†
It is intended that DQI and HLD be calculated at the virtual drug design stage using in silico technologies37,44 to provide the input data, but, for simplicity, this study only considers experimentally derived input data. With respect to a compound's distribution kinetics, it is unlikely that there will be any insight into these at the virtual design stage. To account for this, the model considers five hypothetical distribution kinetics scenarios for each compound, intended to cover a range of possibilities that encompass the true situation. These scenarios depend on a range of hypothetical values for a compound's Vcentral and Vterminal. Vcentral ranges from 3.0 L, which is the approximate value for the plasma volume in human, to a value that is 50% of the compound's Vss; a mid-point value is also considered. Vterminal ranges from 1.1 times the Vss to 2.0 times the Vss; a mid-point value is also considered. For comparison, Rowland and Tozer33 provide human details for aspirin, salicylic acid and gentamicin C: their Vcentral values can be calculated as 5.3, 6.8 and 14.0 L, respectively, and their Vterminal values as 1.1-, 1.3- and 5.1-fold greater than Vss, respectively. The Vterminal of gentamicin C is considered an extreme and a value 2.0-fold greater than Vss deemed appropriate. To illustrate the effect of these different distribution kinetics scenarios, Fig. 3 shows the dose versus Css,central,max curves generated by the PK model for two of these scenarios for diazepam (further details can be found in the ESI†).
It can be seen from Fig. 3 that differing distribution kinetics scenarios can lead to significant variation in values for DQI (i.e., 2.1-fold) and HLD (i.e., 1.3-fold).
The polygons shown in Fig. 4 highlight the DQI–HLD space occupied by a compound. The extent of the area overlap of the polygons for two different compounds reflects their similarity. Importantly, at the virtual design stage it is not possible to define precise values for a compound's DQI and HLD. Instead, it is possible to determine an area of DQI–HLD space, defined by a non-self-intersecting closed polygon, where the true value lies. Provided that the two polygons do not overlap, the compounds can be argued to have distinct DQI and HLD values.
This method can be expanded to incorporate errors in the estimation of the input data. A Gaussian distribution can be assumed, based on the original value and an estimation of the associated standard deviation. Additional simulations using randomly selected values from such a Gaussian distribution can then be considered. To illustrate, DQI and HLD were calculated for the 15 compounds, using 50 additional simulations to account for an error of 0.3 on a logarithmic base 10 scale for the values of Vss and Cl. To simplify the plot in Fig. 5, DQI–HLD polygons are shown for only six compounds: chlorpromazine (monobase), diazepam (neutral), ketoprofen (monoacid), nifedipine (neutral), trimethoprim (monobase) and warfarin (monoacid) (further details can be found in the ESI†). It is clear from Fig. 5 that the polygons are very irregular, but again, provided that two polygons do not overlap, the compounds can be argued to have distinct DQI and HLD values. This approach also applies when the input data are sourced from in silico QSAR models, where any prediction error is expressed as a standard deviation.45
The areas covered by the different polygons vary, with that for chlorpromazine being the largest and that for warfarin the smallest. The calculated areas are summarised in Table 2, along with the polygons' second moment of area and centroid values in the DQI and HLD dimensions. The second moment of area indicates the sensitivity of the compound to errors in the input data with respect to the DQI and HLD dimensions; the larger the value, the greater the sensitivity. With respect to warfarin, the DQI value is more sensitive to the combined errors in the values of Vss and Cl than that for HLD; the reverse is true for prazosin. It is useful to think of the HLD dimension reflecting the ability of a compound to be absorbed into the body, while the DQI dimension reflects the extent of elimination of a compound from the body.
Drug | Area | Second moment of area | Centroid | ||
---|---|---|---|---|---|
HLD | DQI | HLD | DQI | ||
Chlorpromazine | 1.72 | 15.30 | 11.88 | 2.61 | −2.93 |
Diazepam | 0.33 | 0.31 | 0.89 | 1.63 | −0.94 |
Diclofenac | 0.73 | 1.82 | 3.04 | 2.03 | −1.56 |
Furosemide | 0.09 | 0.36 | 0.74 | 2.93 | −2.03 |
Haloperidol | 1.58 | 13.27 | 5.88 | 1.90 | −2.83 |
Imipramine | 1.00 | 9.68 | 13.96 | 3.74 | −3.06 |
Indomethacin | 0.20 | 0.21 | 0.47 | 1.52 | −1.02 |
Ketoprofen | 0.11 | 0.11 | 1.63 | 3.88 | −1.00 |
Naproxen | 0.05 | 0.01 | 0.70 | 3.72 | −0.45 |
Nifedipine | 0.78 | 3.54 | 0.51 | 0.77 | −2.07 |
Phenytoin | 0.46 | 0.62 | 0.60 | 1.11 | −1.13 |
Pindolol | 0.33 | 1.32 | 5.18 | 3.98 | −1.95 |
Prazosin | 0.80 | 4.35 | 0.17 | 0.39 | −2.30 |
Trimethoprim | 0.49 | 0.99 | 5.28 | 3.26 | −1.39 |
Warfarin | 0.05 | 0.01 | 0.29 | 2.38 | −0.30 |
To simplify the information represented by the polygons in Fig. 5, it is proposed that the centroid coordinate could be used to standardise comparisons. Fig. 6 shows a DQI–HLD plot for the total level Css,central,max centroid values for the 15 compounds considered.
log10(max_quantity) = HLD + DQI | (7) |
These observations are well established, but evaluation of the properties of oral compounds (i.e., pKa, solubilitypH7.4, Papp,Caco2,pH6.5, Vss, Cl and, optionally, PPB) can be vastly simplified to a quantitative comparison within the two dimensions of DQI and HLD. The magnitude of a compound's DQI and HDI should be considered in conjunction with the desired log10(quantity).
Human in vivo efficacy (mg L−1) data for the 15 compounds considered have been sourced from Schulz et al. (further details can be found in the ESI†).46 Comparison to predicted logP (ref. 13) indicates a very weak, non-significant linear relationship, but comparison to Cl shows a significant linear relationship (Fig. 7). This is not unexpected, as a compound's Cl will heavily influence in vivo levels – with respect to an intravenous one-compartment model, the steady-state concentration upon repeat dosing is, in theory, inversely proportional to Cl.43
Fig. 8 shows two plots between the 15 compounds' human in vivo efficacy and their centroid DQI for (a) total level Css,central,max and (b) total level AUCss,central quantities (further details can be found in the ESI†). Both show significant linear relationships, such that those that displaying human in vivo efficacy at higher levels tend to have higher centroid DQI and vice versa. The statistics for both plots are slightly better than for the plot in Fig. 7. It can be argued that these centroid DQI values contain slightly more information than Cl for this small data set of 15 compounds. This can be understood from the perspective that oral absorption and drug distribution will have influence over in vivo levels in addition to that of Cl. It can be inferred from these plots that compounds with lower potency require a higher centroid DQI (total level Css,central,max or AUCss,central). This can be understood from the perspective of a fixed dose, such that a less potent compound requires more of the oral dose to be in the central compartment at steady state to drive the therapeutic effect, whereas a more potent compound requires less of the oral dose to be in the central compartment to drive a similar effect.
![]() | ||
Fig. 8 Human in vivo log10(efficacy) versus the total level centroid DQI for (a) Css,central,max and (b) AUCss,central. |
Clearly, an appreciation of potency is a key factor in understanding the plots in Fig. 8. Still, the centroid DQI values for the total level Css,central,max or AUCss,central quantities are sufficient to explain approximately half of the information in the human in vivo efficacy values for the 15 compounds.
Relatedly, DQI and HLD values can be proposed as novel compound descriptors for use in conjunction with QSAR modelling methods (and other descriptors) to model other in vivo quantity endpoints. At the virtual design stage, such an approach would involve use of in silico predictions for input data for the model to derive DQI and HLD values, which could then be used as descriptors in further in silico models.
Although this study focuses on determining DQI and HLD for the total and free level Css,central,max and AUCss,central, other in vivo quantities could also be assessed using this model. These include total and free levels for the average concentration (Css,central,average) and the minimum concentration (Css,central,min) at steady state in compartment B. In addition, the corresponding levels in compartment C (peripheral compartment) can be considered, including steady-state total and free levels for the maximum concentration (Css,peripheral,max), the average concentration (Css,peripheral,average), the minimum concentration (Css,peripheral,min) and the area under the curve AUCss,peripheral.
Importantly, a twice-daily dosing scenario is considered; changing the dosing interval will lead to different DQI and HLD values for the same compound. Changing the simulation length can also affect these values, in particular if it is shortened. From a virtual drug design perspective, consideration of 14 repeat doses using a once- or twice-daily dosing scenario is recommended.
With respect to the use of 20 different dose simulations to define a dose–quantity curve, this can be reduced to two in theory. The DQI can be determined from a simulation using a very low dose (e.g., 0.000001 mg), where linear PK can be assumed. The HLD can be determined from a simulation using a very high dose (e.g., 10000.0 mg), where compartment A can be assumed to be saturated throughout the absorption window and the corresponding log10(quantity) is at its maximum (i.e., log10(max_quantity)). Relatedly, the standardised method used to determine the HLD only provides an approximation – other methods can be used, including fitting the dose–quantity curve to a power function of the form:
log10(quantity) = log10(max_quantity) − log10(1 + 10HLD−log10(dose)) | (8) |
The application of eqn (8) using a least-fit method benefits from the use of as many data points as possible to define the transition region of the dose–quantity curve from linear to non-linear PK.
The DQI and HLD calculations can be made in different mammalian species by adjusting the PK model settings accordingly and by using species-specific in silico predictions for Vss, Cl and PPB (and assuming the use of Papp,Caco2,pH6.5 for other species).
Finally, an extension of this study, which is beyond the scope of the present work due to its size and complexity, would be to use in silico QSAR models for pKa, solubilitypH7.4, Papp,Caco2,pH6.5, Vss, Cl and (optionally) PPB to predict the input data for a larger set of compounds, calculate their DQI values, and assess how well they relate to in vivo efficacy and toxicity data.
Application of this methodology to 15 known oral drugs demonstrates how different compounds that have vastly different tendencies in terms of absorption into, and elimination from, the body can be compared on the same scale. In spite of the computational complexities associated with gaining such insight, DQI and HLD provide a direct relationship between a compound's dose and in vivo exposure. Evaluation of a compound's oral drug suitability is simply dependent on the ability to match a compound's dose to the required in vivo exposure.
Footnote |
† Electronic supplementary information (ESI) available: (1) A file containing the raw data for Table 1; (2) a file containing details of the mathematical equations associated to the PK models discussed; (3) a file containing the raw data for eqn (2); (4) a file containing the summary of raw data for Fig. 3 and 4a; (5) a file containing the summary of raw data for Fig. 4b; (6) a file containing the plot data for Fig. 4a and b, 5, 7 and 8a and b; (7) a file containing the command line code for the insilicolynxdqi Python library to replicate the raw data; (8) a file containing the required input data formatted for use with the insilicolynxdqi Python library. See DOI: 10.1039/c7md00586e |
This journal is © The Royal Society of Chemistry 2018 |