Assessing the pharmacokinetic profile of the CamMedNP natural products database: an in silico approach

Abstract

Background

Drug metabolism and pharmacokinetic (DMPK) assessment has come to occupy a place of interest during the early stages of drug discovery today. Computer-based methods are slowly gaining ground in this area and are often used as initial tools to eliminate compounds likely to present uninteresting pharmacokinetic profiles and unacceptable levels of toxicity from the list of potential drug candidates, hence cutting down the cost of the discovery of a drug.

Results

In the present study, we present an in silico assessment of the DMPK profile of our recently published natural products database of 1,859 unique compounds derived from 224 species of medicinal plants from the Cameroonian forest. In this analysis, we have used 46 computed physico-chemical properties or molecular descriptors to predict the absorption, distribution, metabolism and elimination (ADME) of the compounds. This survey demonstrated that about 50% of the compounds within the Cameroonian medicinal plant and natural products (CamMedNP) database are compliant, having properties which fall within the range of ADME properties of >95% of currently known drugs, while >73% of the compounds have ≤2 violations. Moreover, about 72% of the compounds within the corresponding ‘drug-like’ subset showed compliance.

Conclusions

In addition to the previously verified levels of ‘drug-likeness’ and the diversity and the wide range of measured biological activities, the compounds in the CamMedNP database show interesting DMPK profiles and, hence, could represent an important starting point for hit/lead discovery from medicinal plants in Africa.

Background

Natural products (NPs) play an increasingly important role in drug discovery today [1–5], both serving as drugs and as templates for the design of nature-inspired medicines [3, 6]. In fact, it has been reported that a significant proportion of drugs that undergo clinical trials are either naturally occurring or are derived from NPs [7]. What characterises NPs are their richness in stereogenic centres and coverage of segments of chemical space which are typically not occupied by a majority of synthetic molecules and drugs [8, 9]. In addition, they generally contain more oxygen atoms and less aromatic atoms on average, when compared with ‘drug-like’ molecules [8–11]. It is needless to say that NPs sometimes fail the famous ‘drug-likeness’ test due to the often bulky nature of naturally occurring metabolites [11].

It is also worth mentioning that designing drug-like molecules having interesting pharmacokinetic properties is an important paradigm in drug discovery programs [12, 13]. This entails the search for lead compounds which can be easily orally absorbed, easily transported to their desired site of action, not easily attacked by metabolising enzymes so as to form toxic metabolic products before reaching the targeted site of action and easily eliminated from the body before accumulating in sufficient amounts that may produce adverse side effects. The ensemble of the above properties is often referred to as absorption, distribution, metabolism and elimination (ADME) properties, or better still ADMET or ADME/T or ADMETox (i.e. if toxicity criteria are also taken into consideration).

Computer-based in silico approaches for the prediction of ADMET profiles of drug leads at early stages of drug discovery are increasingly gaining ground [14–16]. This could be explained by the relative cost advantage added to the time factor, when compared to standard experimental approaches for ADMET profiling [17, 18]. On these grounds, several theoretical methods for the determination of ADMET parameters have been developed and implemented in a number of currently available software for drug discovery protocols [19–22], even though the predictions are sometimes disappointing [23]. Such software often make use of quantitative structure-activity relationships [22–24] or knowledge-base methods [25–27]. The goal has been to considerably cut down on the currently very high cost of discovery of a drug [17]. A promising lead is often defined as a compound which combines potency with an admirable ADMET profile. As such, compounds with unfavourably predicted pharmacokinetic profiles are either completely dismissed from the list of potential drug candidates (even if they prove to be highly potent) or the drug metabolism and pharmacokinetics (DMPK) properties are ‘fine tuned’ in order to improve their chances of making it to clinical trials [28]. This explains why the ‘graveyard’ of very highly potent compounds which do not make it to clinical trials keeps filling up, to the extent that the process of drug discovery often presents the challenge of either resorting to new leads or ‘resurrecting’ some buried leads with the view of fine-tuning their ADMET profiles.

In a recent paper, we have presented a database of 1,859 compounds derived from the Cameroonian flora, Cameroonian medicinal plant and natural products (CamMedNP), the compounds being predicted to be sufficiently orally available and diverse to be employed in lead discovery programs [29]. Additional arguments in favour of the use of this database are the wide range of the previously observed biological activities of the compounds and the wide range of ailments being treated by traditional medicine with the help of the herbs from which the compounds have been derived [29, 30].

Numerous drugs at a late stage of pharmaceutical development and many more lead compounds fail due to adverse pharmacokinetic properties [18]. It is, therefore, important to incorporate the prediction of the ADME properties into the lead compound selection, by means of molecular descriptors. A molecular descriptor is often defined as a structural or physico-chemical property of a molecule or part of a molecule, for example the logarithm of the n-octanol/water partition coefficient (log P), molar weight (MW) and total polar surface area. A number of relevant molecular properties (descriptors) are often used to help predict the pharmacokinetic behaviour of potential drug leads. In the present study, we have carried out an in silico assessment of the ADMET profile of the CamMedNP database by the use of computed molecular descriptors currently implemented in a wide range of software tools as indicators of the pharmacokinetic properties of a large proportion of currently known drugs.

Methods

Data sources and generation of 3D structures

The plant sources, geographical collection sites, chemical structures of pure compounds and their measured biological activities were retrieved from literature sources and have been previously described [29]. The three-dimensional (3D) structures were generated using the builder module of MOE [31], and energy minimization was subsequently carried out using the MMFF94 [32] until a gradient of 0.01 kcal/mol was reached.

Initial treatment of chemical structures and calculation of ADMET-related descriptors

The 1,859 low-energy 3D chemical structures in the CamMedNP library were saved in mol2 format and initially treated with LigPrep [33], distributed by Schrodinger, Inc. (New York, USA). This implementation was carried out with the graphical user interface of the Maestro software package (New York, USA) [34], using the OPLS force field [35–37]. Protonation states at biologically relevant pH were correctly assigned (group I metals in simple salts were disconnected, strong acids were deprotonated and strong bases protonated, while topological duplicates and explicit hydrogens were added). All molecular modelling was carried out on a Linux workstation (San Francisco, USA) with a 3.5 GHz Intel Core2 Duo processor (Santa Clara, USA). A set of the ADMET-related properties (a total of 46 molecular descriptors) were calculated using the QikProp program (New York, USA) [21] running in normal mode. QikProp generates physically relevant descriptors and uses them to perform ADMET predictions. An overall ADME-compliance score, drug-likeness parameter (indicated by #stars), was used to assess the pharmacokinetic profiles of the compounds within the CamMedNP library. The #stars parameter indicates the number of property descriptors computed by QikProp, which falls outside the optimum range of values for 95% of known drugs. The methods implemented were developed by Jorgensen et al. [38–40]. Among the calculated descriptors are the total solvent-accessible molecular surface, S _mol in Ų (probe radius 1.4 Å; range for 95% of drugs is 300 to 1,000 Ų); the hydrophobic portion of the solvent-accessible molecular surface, S _mol,hfob in Ų (probe radius 1.4 Å; range for 95% of drugs is 0 to 750 Ų); the total volume of molecule enclosed by solvent-accessible molecular surface, V _mol in ų (probe radius 1.4 Å; range for 95% of drugs is 500 to 2,000 ų); the logarithm of aqueous solubility, logS _wat (range for 95% of drugs is −6.0 to 0.5) [36, 38]; the logarithm of predicted binding constant to human serum albumin, logK _HSA (range for 95% of drugs is −1.5 to 1.2) [41]; the logarithm of predicted blood/brain barrier partition coefficient, log B/B (range for 95% of drugs is −3.0 to 1.0) [42–44]; the predicted apparent Caco-2 cell membrane permeability (BIP _Caco-2) in Boehringer-Ingelheim scale, in nm/s (range for 95% of drugs is <5 low, >100 high) [45–47]; the predicted apparent Madin-Darby canine kidney (MDCK) cell permeability in nm s⁻¹ (<25 poor, >500 great) [46]; the index of cohesion interaction in solids, Ind_coh, calculated from the number of hydrogen bond acceptors (HBA), hydrogen bond donors (HBD) and the surface area accessible to the solvent (S _mol) by the relation $\mathrm{In}{\mathrm{d}}_{\mathrm{coh}}=\mathrm{HBA}\times \sqrt{\mathrm{HBD}}/S_{\mathrm{mol}}$ (0.0 to 0.05 for 95% of drugs) [40]; the globularity descriptor, Glob = (4πr²)/S _mol, where r is the radius of the sphere whose volume is equal to the molecular volume (0.75 to 0.95 for 95% of drugs); the predicted polarizability, QP_polrz (13.0 to 70.0 for 95% of drugs); the predicted IC₅₀ value for blockage of HERG K⁺ channels, logHERG (concern <−5) [48, 49]; the predicted skin permeability, logK _p (−8.0 to −1.0 for 95% of drugs) [50, 51]; and the number of likely metabolic reactions, #metab (range for 95% of drugs is 0 to 15).

Results and discussion

Overall DMPK compliance of the CamMedNP library

The 24 most relevant molecular descriptors calculated by QikProp are used to determine the #star parameter [52]. A plot of the #stars parameter (on the x-axis) against the corresponding counts (on the y-axis) in the CamMedNP is shown within the same set of axes with those of the ‘drug-like’, ‘lead-like’ and ‘fragment-like’ standard subsets, Figure 1. The criteria for the respective standard subsets were defined as MW < 500, log P < 5, HBD ≤ 5, HBA ≤ 10 [14]; 150 ≤ MW ≤ 350, log P ≤ 4, HBD ≤ 3, HBA ≤ 6 [53–55] and MW ≤ 250, −2 ≤ log P ≤ 3, HBD < 3, HBA < 6, NRB < 3 [56]. QikProp was unable to compute the ADMET descriptors for 25 compounds out of the total library due to limitations that were not clear to us. Of the remaining 1,834 compounds, 48.04% showed #star = 0, while 74.21% had #star ≤ 2. Among the 1,122 compounds of the drug-like subset, 79.12% had pharmacokinetic descriptors within the acceptable range for 95% of known drugs, while 97.33% showed #stars ≤ 2. The lead-like and fragment-like subsets were, respectively, 81.15% and 55.56% compliant for all of the 24 most relevant computed descriptors. The mean values for 19 selected computed descriptors have been shown in Table 1 for all four compound libraries, while the percentage compliances for 14 selected ADMET-related descriptors are shown in Table 2. The mean values and percentage compliances indicate a high probability of finding drug leads within the CamMedNP compound library.

Figure 1

Distribution curves for #stars within the CamMedNP library and subsets. Blue = CamMedNP library, red = drug-like subset, green = lead-like subset and violet = fragment-like subset.

Table 1 Average pharmacokinetic property distributions of total CamMedNP library in comparison with various subsets

Table 2 Percentage compliances of selected ADMET-related descriptors of total CamMedNP library in comparison with various subsets

Bioavailability prediction

The bioavailability of a compound depends on the processes of absorption and liver first-pass metabolism [57]. The absorption, in turn, depends on the solubility and permeability of the compound, as well as on the interactions with transporters and metabolizing enzymes in the gut wall. The computed parameters used to assess oral absorption are the predicted aqueous solubility, logS _wat, the conformation-independent predicted aqueous solubility, CI logS _wat, the predicted qualitative human oral absorption, the predicted % human oral absorption and compliance to Jorgensen's famous ‘Rule of Three’ (ro3). The solubility calculation procedure implemented depends on the similarity property space between the given molecule and its most similar analogue within the experimental training set used to develop the model implemented in QikProp, i.e. if the similarity is <0.9, then the QikProp predicted value is taken; otherwise, the predicted property, P _pred, is adjusted such that

$$P_{\mathrm{pred}}=S{P}_{\exp }+\left(1-S\right)P_{\mathit{QP}}$$

(1)

where S is the similarity and P _exp and P _QP are, respectively, the experimental and QikProp predictions for the most similar molecule within the training set. In Equation 1, if S = 1, then the predicted property is equal to the measured experimental property of the training set compound. According to Jorgensen's ro3, if a compound complies to all or some of the rules (logS _wat > −5.7, BIP _Caco-2 > 22 nm/s and number of primary metabolites < seven), then it is more likely to be orally available. The distribution curves for two of the three determinants for the ro3 (logS _wat and BIP _Caco-2) are shown in Figure 2A,B. In general, 47.22% of the CamMedNP library was compliant to the ro3, while the respective percentage compliances for the various subsets were 72.28%, 92.11% and 100% for the drug-like, lead-like and fragment-like libraries. Among the individual computed parameters, the most remarkable was logS _wat, which was met by 75.74% of the compounds within the CamMedNP library, while this property shows a Gaussian distribution for the drug-like and lead-like subsets. Only 37.94% of the compounds fell within the respected range for the BIP _Caco-2 criterion. The predicted apparent Caco-2 cell permeability, BIP _Caco-2 (in nm s⁻¹), models the permeability of the gut-blood barrier (for non-active transport), even though this parameter is not often correctly predicted computationally [58]. The histograms of the predicted qualitative human oral absorption parameter (in the scale 1 = low, 2 = medium and 3 = high) are shown in Figure 3. It was observed that 52.45% of the compounds in CamMedNP were predicted to have high human oral absorption. The predicted % human oral absorption (on 0 to 100% scale) shows a similar trend, with 41.06% of the compounds being predicted to be absorbed at 100%, and 57.96% of the compounds predicted to be absorbed at >90%.

Figure 2

Distribution curves for compliance to Jorgensen's ‘Rule of Three’. (A) calculated logS _wat against count. (B) Predicted BIP _Caco-2 against count. Blue = CamMedNP library, red = drug-like subset, green = lead-like subset and violet = fragment-like subset.

Figure 3

Histograms showing the distribution of human oral absorption predictions.

The size of a molecule, as well as its capacity to make hydrogen bonds, its overall lipophilicity, its shape and flexibility are important properties to consider when determining permeability. Molecular flexibility has been seen as a parameter which is dependent on the number of rotatable bonds (NRB), a property which influences the bioavailability in rats [58]. The distribution of the NRB for this dataset has been previously discussed [29] and revealed that the compounds within the CamMedNP library show some degree of conformational flexibility, the peak value for the NRB being between 1 and 2, while the average value is 5.31 (Table 1).

Prediction of blood-brain barrier penetration

Too polar drugs do not cross the BBB. The blood/brain partition coefficients (log B/B) were computed and used as a predictor for access to the central nervous system (CNS). The predicted CNS activity was computed on a −2 (inactive) to +2 (active) scale and showed that only 1.85% of the compounds in the CamMedNP could be active in the CNS (predicted CNS activity >1). A distribution of the log B/B (Figure 4) shows a right-slanted Gaussian-shaped curve with a peak at −0.5 log B/B units (the same for all the standard subsets), with >88% of the compounds in the CamMedNP falling within the recommended range for the predicted brain/blood partition coefficient (−3.0 to 1.2). The MDCK monolayers are widely used to make oral absorption estimates, the reason being that these cells also express transporter proteins, but only express very low levels of metabolizing enzymes [58]. They are also used as an additional criterion to predict BBB penetration. Thus, our calculated apparent MDCK cell permeability could be considered to be a good mimic for the BBB (for non-active transport). It was estimated that only about 50% of the compounds had apparent MDCK cell permeabilities falling within the recommended range of 25 to 500 nm s⁻¹ for 95% of known drugs. This situation was not greatly improved in the drug-like and lead-like subsets (58% and 57%, respectively).

Figure 4

Plot of the physico-chemical descriptor used to predict BBB penetration. Predicted log B/B against count. The x-axis label is the lower limit of the binned data, e.g. 0 is equivalent to 0.0 to 1.0. Blue = CamMedNP library, red = drug-like subset, green = lead-like subset and violet = fragment-like subset.

Prediction of dermal penetration

This factor is important for drugs administered through the skin. The distribution of computed skin permeability parameter, logK _p, showed smooth Gaussian-shaped graphs centred at −2.5 logK _p units for all the four datasets (Figure 5), with approximately 91% of the compounds in the CamMedNP database falling within the recommended range for >95% of known drugs. The predicted maximum transdermal transport rates, J _m (in μ cm⁻² h⁻¹), were computed from the aqueous solubility (S _wat), the MW and skin permeability (K _p) using the relation (2):

$$J_{\mathrm{m}}=K_{\mathrm{p}}\times \mathrm{MW}\times {\mathrm{S}}_{\mathrm{wat}}$$

(2)

Figure 5

Distribution curves for the predicted skin penetration parameter. Blue = CamMedNP library, red = drug-like subset, green = lead-like subset and violet = fragment-like subset.

This parameter showed variations from 0 to 1,603 μ cm⁻² h⁻¹, with only about 1.39% of the compounds in CamMedNP having the predicted value of J _m > 100 μ cm⁻² h⁻¹.

Prediction of plasma-protein binding

The efficiency of a drug may be affected by the degree to which it binds to the proteins within the blood plasma. It is noteworthy that the binding of drugs to the plasma proteins (like human serum albumin, lipoprotein, glycoprotein, α, β and γ globulins) greatly reduces the quantity of the drug in the general blood circulation, and hence, the less bound a drug is, the more efficiently it can traverse cell membranes or diffuse. The predicted plasma-protein binding has been estimated by the prediction of binding to human serum albumin; the logK _HSA parameter recommended range is −1.5 to 1.5 for 95% of known drugs. Figure 6 shows the variation of this calculated parameter within the CamMedNP dataset, as well as for the standard subsets. This equally gave smooth Gaussian-shaped curves centred on −0.5 logK _HSA units for all the four datasets. In addition, our calculations revealed that >85% of the compounds within the CamMedNP library are compliant to this parameter, indicating that a majority of the compounds are likely to circulate freely within the blood stream and, hence, have access to the target site.

Figure 6

Distribution curves for predicted plasma-protein binding. Blue = CamMedNP library, red = drug-like subset, green = lead-like subset and violet = fragment-like subset.

Metabolism prediction

An estimated number of possible metabolic reactions has also been predicted by QikProp and used to determine whether the molecules can easily gain access to the target site after entering the blood stream. The average estimated number of possible metabolic reactions for the CamMedNP library was between five and six, while those of the standard subsets drop sequentially by one step in a progressive manner (Table 1). Even though some of the compounds are likely to undergo as many as up to 26 metabolic reactions due to the complexity of some of the plant secondary metabolites within the database (Figure 7), about 80% of the compounds are predicted to undergo the recommended number of metabolic steps (one to eight reactions), with the situation improving to around 90% and approximately 97% in the drug-like and lead-like subsets, respectively. From Figure 7, it can be observed that, except for the fragment-like subsets which peaks at two predicted metabolic reactions, the peak values for the number of predicted metabolic reactions were at three for all of the datasets.

Figure 7

Distribution of the predicted number of metabolic reactions for compounds in the CamMedNP. Blue = CamMedNP library, red = drug-like subset, green = lead-like subset and violet = fragment-like subset.

Prediction of blockage of human ether-a-go-go-related gene potassium channel

Human ether-a-go-go-related gene (HERG) encodes a potassium ion (K⁺) channel that is implicated in the fatal arrhythmia known as torsade de pointes or the long QT syndrome [59]. The HERG K⁺ channel, which is best known for its contribution to the electrical activity of the heart which coordinates the heart's beating, appears to be the molecular target responsible for the cardiac toxicity of a wide range of therapeutic drugs [60]. HERG has also been associated with modulating the functions of some cells of the nervous system and with establishing and maintaining cancer-like features in leukemic cells [61]. Thus, HERG K⁺ channel blockers are potentially toxic, and the predicted IC₅₀ values often provide reasonable predictions for cardiac toxicity of drugs in the early stages of drug discovery [62]. In this work, the estimated or predicted IC₅₀ values for blockage of this channel have been used to model the process in silico. The recommended range for the predicted log IC₅₀ values for blockage of the HERG K⁺ channels (logHERG) is >−5. A distribution curve for the variation of the predicted logHERG is shown in Figure 8, which is left-slanted Gaussian-shaped curve with a peak at −5.5 logHERG units for both the total library and the drug-like subset, meanwhile the lead-like library rather peaks at −4.5 units. It was observed that, in general, this parameter is predicted to fall within the recommended range for about 55% of the compounds within the CamMedNP database, approximately 62% for the drug-like subset and around 73% for the lead-like subset.

Figure 8

A plot of the predicted logHERG values for the CamMedNP and standard subsets. Blue = CamMedNP library, red = drug-like subset, green = lead-like subset and violet = fragment-like subset.

Usefulness of the CamMedNP library

The usefulness of the CamMedNP database in lead generation has been exemplified with the docking and pharmacophore-based screening for potential inhibitors of a validated anti-malarial drug target in our laboratory, and the results will be published in a subsequent paper. It is important to mention that virtual screening results could provide insight and direct natural products chemists to search for theoretically active principles with attractive ADMET profiles, which have been previously isolated, but not tested for activity against specified drug targets (if samples are absent). This ‘resurrection’ process could prove to be a better procedure for lead search than the random screening, which is a common practice in our Cameroonian laboratories. CamMedNP is constantly being updated; meanwhile, a MySQL platform (Cupertino, USA) to facilitate the searching of this database and ordering of compound samples is under development within our group and will also be published subsequently. However, 3D structures of the compounds, as well as their physico-chemical properties that were used to evaluate the DMPK profile, can be freely downloaded as additional files accompanying this publication (see Additional file 1, Additional file 2, Additional file 3, Additional file 4). In addition, information about compound sample availability can be obtained on request from the authors of this paper or from the pan-African Natural Products Library (p-ANAPL) project [63, 64].

Conclusion

Modern drug discovery programs usually involve the search for small molecule leads with attractive pharmacokinetic profiles. The presence of such within the CamMedNP library is of major importance and, therefore, renders the database attractive, in addition to the already-known properties (drug-like, lead-like fragment-like and diverse). This is an indication that the 3D structures of naturally occurring compounds within the CamMedNP could be a good starting point for docking, neural networking and pharmacophore-based virtual screening campaigns, thus rendering the CamMedNP as a useful asset for the drug discovery community. 3D structures of the compounds, as well as their physico-chemical properties that were used to evaluate the DMPK profile of the CamMedNP library, can be freely downloaded (for non-commercial use) as additional files which accompany this publication (see Additional file 1, Additional file 2, Additional file 3, Additional file 4). The physical samples for testing are available at the various research laboratories in Cameroon in varying quantities. Questions regarding the availability of the compound samples could be addressed directly to the authors of this paper. Otherwise, the samples could be obtainable from the p-ANAPL consortium, which has a mandate to collect samples of NPs from the entire continent of Africa and make them available for biological screening. This network is being set up under the auspices of the Network for Analytical and Bioassay Services in Africa [63, 64].

Authors' information

WS and SMNE are professors of Medicinal Chemistry with an interest in CADD, while SMNE also focuses on organic synthesis and on natural product leads from the Cameroonian medicinal plants. LMM and JAM are natural products chemists actively involved in the isolation and characterization of secondary metabolites from the Cameroonian medicinal plants. LLL holds a PhD in Environmental Science and manages the Chemical and Bioactivity Information Centre (CBIC) with a focus on developing databases for information from medicinal herbs in Africa. PNJ is a retired research officer of Lhasa Ltd., who currently leads the CBIC branch in Leeds, UK. FNK is a PhD student working on CADD under the joint supervision of LCOO and EM.