In Silico ADME Methods Used in the Evaluation of Natural Products

1. Introduction

Many promising drug candidates fail to reach the market due to poor ADME characteristics (how the drug is absorbed, transported throughout the body, metabolized, and excreted) and safety concerns about toxicity. When these problems are uncovered late in drug development, pharmaceutical companies risk large financial losses since they have already invested heavily in clinical trials and other costly research processes. To address this dilemma, the pharmaceutical industry has significantly changed its strategy in the last decades. Companies are increasingly performing extensive ADMET (absorption, distribution, metabolism, excretion, and toxicity) screening considerably earlier in the drug discovery process. This early screening technique identifies and eliminates problematic compounds before they enter costly development phases, thereby saving money and boosting medication development efficiency [1,2]. The presence of pan-assay interference compounds (PAINS) is an additional significant concern. Deceptive findings in research tests can be generated by these challenging compounds, which appear to be active against a diverse array of targets when they are not. It is a significant waste of time and money to investigate these “frequent hitters” in the context of research and development [2,3].

Compared to synthetic molecules, natural compounds possess unique properties that influence drug discovery. They are more structurally diverse and complex, tend to be larger, contain a greater amount of oxygen (and less nitrogen, sulfur and halogens), more chiral centers, less aromatic rings, and are more water-soluble. This provides them with a distinctive potential as drugs, even when they do not adhere to the conventional principles for drug-like properties, such as Lipinski’s rule of five [4,5]. Nevertheless, the utilization of natural compounds in the drug discovery process is hindered by a number of obstacles. It can be challenging to test natural extracts, identify the active constituents in complex mixtures, obtain sufficient material from nature, and protect discoveries with patents. The pharmaceutical industry has experienced a decline in drug discovery that concentrates on natural products as a result of these obstacles [4,6]. Nevertheless, the discovery and development of natural product-based drugs are being substantially improved by recent developments in biosynthetic engineering technologies, chemical synthesis, as well as “smart screening” techniques or use of robotic equipment in chemical separations and progresses in structural analysis. Complex natural compound scaffolds that were previously considered inaccessible can now be optimized with these technologies. This advancement enables the enrichment of screening libraries with natural products, hybrids derived from natural products, natural product analogs and natural products inspired molecules, as well as superior structure functionalization approaches, including late-stage functionalization, for the optimization of natural product leads [4,6,7,8]. In addition to their potential applications in drug discovery, numerous bioactive natural compounds are essential components of our daily diets. Gaining an understanding of their properties is essential for the advancement of human health, as it facilitates the creation of enhanced foods, such as functional foods and dietary supplements [9].

Many of the challenges applicable in the understanding of the pharmacological or biological properties of natural compounds are also relevant when exploring their ADME properties. For instance, often the available quantities of natural products are limited, and while numerous plant-derived natural products have been isolated and characterized, the amounts available are frequently insufficient for comprehensive ADME testing [10]. Using in silico methods from this point of view has a great advantage as they require no physical sample (not even picograms are necessary once the structural formula is available) or laboratory infrastructure. Besides (and this is an aspect relevant for all products of pharmaceutical interest, irrespective of their origin), the experimental assessment of ADME properties of a substance is costly and time consuming [11], whereas the use of in silico tools is usually very cheap.

Natural compounds are often highly sensitive to environmental factors such as high temperature, moisture, light exposure, oxygen, or pH variations. They may also be volatile or react with other substances and stability issues result in limited shelf-life, making it difficult to develop stable commercial products that maintain their essential properties over time [12]. Besides, many natural compounds may be degraded by stomach acid or undergo extensive metabolism in the liver before reaching their target sites in the body (first-pass metabolism) [13]. This is also a serious issue in assessing their ADME features in animal experiments, whereas using very cheap but valid in silico tools is very useful. An additional challenge for exploring ADME properties in wet lab experiments consists in the low aqueous solubility of natural compounds [14] , which limits the ability to effectively deliver the compound to the biological system examined. Moreover, the growing need to minimize animal use in medical development and research highlights the increasing significance of in silico tools [15].

The prediction of physicochemical properties for understanding ADME characteristics has deep roots in computational chemistry. A pioneering development was the calculated partition coefficient (cLogP), introduced by Hansch, Iwasa, and Fujita in the 1950s and 1960s, which remains one of the most referenced physicochemical parameters. A major breakthrough came with the introduction of the Lipinski’s “Rule of 5,” which was based on the fact that straightforward molecular descriptors – including molecular weight and the number of oxygen, nitrogen, and hydrogen bonds – strongly influence a drug’s absorption and membrane permeation. The field continued to evolve with Leeson and Springthorpe’s criterion (cLogP > 4.5-5), which was later enhanced by new concepts such as lipophilicity ligand efficiency (LLE). Further advances included the development of the solubility forecast index (SFI) and the property forecast index (PFI). These developments were complemented by extensive analyses of clinical drug candidates, which helped refine and update Lipinski’s original rules to better predict drug-like properties [16]. In the last decades substantial efforts have been deployed and marked progress has been made in the field of developing in silico methods and tools for the ADME assessment of new chemical compounds.

In this paper, we have proposed a synthesis of the scientific literature in order to understand the in silico methods and tools available for the prediction of ADME properties, with a focus on those that have been used to evaluate natural compounds. We were interested in understanding which of the available methods and tools have been most frequently used in this field, although such methods and tools have been developed for drug discovery purposes in general and not specifically for natural products. We focus here exclusively on ADME (Absorption, Distribution, Metabolism, and Excretion) rather than ADMET. While both abbreviations are commonly used, our emphasis on ADME is related to its direct relevance to pharmacokinetics. The ’T’ in ADMET primarily addresses safety concerns (a rather distinct field of study) and this limited focus on ADME allows for a clearer focus on the pharmacokinetic aspects of drug development.

2. In silico methods available for ADME predictions

2.1. Quantum mechanics (QM) methods

Early studies of how drugs interact with receptors attempted to use quantum mechanics but these calculations were limited because they required intensive computational resources. Thanks to significant advances in computer speed and new software, quantum mechanics calculations are now used more regularly to study drug-related problems. In the past three decades, such calculations have become much more common in studying ADME properties [17].

Quantum chemistry states that electrons in atoms cannot have any energy level or be located anywhere in space. Their energy and location are determined by the Schrödinger equation. QM tries to solve this equation, often using approximations. The simplest of these approximations are called semiempirical methods. They were originally implemented in computer programs written in FORTRAN and today these semiempirical approaches are widely used and are included in many different computational chemistry software packages [16]. QM is capable of forecasting the molecular properties of a novel compound in the absence of any prior knowledge or parameterization [18]. Semiempirical methods can perform calculations very quickly on modern computers because they use a simplified approach to handle electron interactions and incorporate experimental data, such as ionization energy, bonding energy, and transition energy. These methods achieve their speed by replacing certain electron interaction calculations (e.g., two-electron repulsion integrals) with pre-determined empirical parameters (all QM methods use approximations, but these approximations are most pronounced in semiempirical methods). This simplification makes it possible to perform calculations on very large molecular systems, such as entire proteins, which would be computationally impossible with more detailed methods [16,17,18]. Among the most widely known semiempirical methods available in open-source and proprietary software are the following:

• MNDO (Modified Neglect of Diatomic Overlap) [19] and its improvements, AM1 (Austin Model 1) and PM3 (Parametric Method 3), the latter being an improvement of the former, with higher accuracy [18];
• PM6 and later PM7, which have provided corrections and improvements over AM1 and PM3 [20,21,22];
• the OMx family (orthogonalization-corrected methods - OM1, OM2, OM3), which have been shown to be roughly equivalent or slightly inferior in performance when compared with standard DFT methods (BPE, B3LYP), which will be discussed below [23];
• the density-functional tight-binding (DFTB) method which, with strong theoretical underpinnings, is closely connected to the DFT. However, it is as computationally efficient as empirical tight-binding methods and is used for big molecules, including those present in biological systems [24].

Semiempirical methods (MNDO, PM6) have been used up to date mostly to characterize the chemical stability and reactivity of natural compounds [25,26,27] or to optimize their chemical structures (e.g., for docking purposes) [28,29,30,31], but not to directly make inferences or predictions about specific ADME properties.

Ab initio methods (latin “ab initio”: from the beginning represent a more advanced and accurate approach within QM calculations. Such methods provide a more complete mathematical description of how particles interact within molecules by using a special mathematical function called the Hamiltonian. The Hamiltonian calculates the total energy of a system by adding up two types of energy: the kinetic energy and potential energy for all electrons and nuclei in the molecule. Hartree-Fock (HF) is one of the most basic ab initio methods [16]. It is the foundation of molecular orbital (MO) theory and it is also referred to as the independent particle model or mean-field theory [18]. It uses a self-consistent field (SCF) approach, each electron being treated as moving within an average electrical field created by all the other electrons. However, it doesn’t consider the specific, instantaneous positions of those other electrons. Because of this simplification, the repulsion between electrons is overestimated, because the averaging approach doesn’t account for the fact that electrons naturally avoid each other in real systems [16].

In 1934, Moller and Plesset published a short five-page paper describing how the Hartree-Fock (HF) method might be improved to account for electron-pair correlation using second-order perturbative theory (the interactions between electrons are called electronic correlations in quantum systems). This method is now called Møller-Plesset perturbation theory (MPPT or simply MP). It initially received little attention. Over the past 40 years, however, it has had a major impact on the development of ab initio quantum chemical methods [32,33].

Configuration interaction theory (CI) was the first method used in quantum chemistry to improve on the HF method. E.A. Hylleraas performed the first CI calculations for the helium atom in 1928. Although CI theory started early, it didn’t develop much until the 1950s. This was due to a lack of computing power and the disruptions of World War II. In 1950, S.F. Boys did a small CI study on beryllium, which gave lower energy than the HF method and was computationally cheaper[34]. However, even today both MP and CI methods are still considered too computationally intensive [35]. Full CI is impractical for all but the smallest molecules. Therefore, truncated CI methods, which include only a few selected configurations, are more widely used. For example, CIS only includes configurations in which one electron is “excited” to a higher energy level, while CISD includes both single and double excitations from the H-F reference wave function. [36,37].

Coupled cluster (CC) methods are another improvement, that calculates electronic correlation energy by bulding upon the HF wave function and including increasing levels of electronic excitation (single, double, etc.). Higher-level excitations enhance accuracy but need substantially more computational power. CCSD(T), which contains single, double, and perturbative triple excitations, is regarded as the “gold standard” in computational quantum chemistry [38]. CC methods have very high accuracy but can only be applied to maximum 20 atoms [39], making them rarely suitable for ADME or drug development applications.

A variety of other ab initio methods have also been made available and are implemented in proprietary or open-source sofware, each with their own strengths, limitations and costs (e.g., CASSCF - complete active space self-consistent field, CASPT2 - complete active space perturbation theory, and other multireference methods [40], multi-configurational self-consistent field (MCSCF) [41], quantum Monte Carlo methods [42], explicitly correlated methods e.g .R12/F12 [43], the density matrix renormalization group method (DMRG) [44] , etc.).

Among ab initio methods for the study of large systems, such as those pertinent to ADME, density functional theory (DFT) remains the sole computationally viable approach. DFT’s utility stems from ongoing enhancements in its approximation of complex electron interactions within the Kohn-Sham framework [45]. DFT directly calculates electron density instead of working with the complex many-electron. Due to the incomplete understanding of the theoretical framework for defining the energy function of electron density, Hartree-Fock approaches were deemed more dependable than Density Functional Theory. This, however, has lately altered, and DFT approaches are now gaining significant popularity [16]. DFT methods are not yet fast enough for extensive free energy calculations on systems beyond simple models. Nevertheless, they can serve as a benchmark for quicker methods and give important insights into how energy (enthalpy) and disorder (entropy) contribute to the overall stability of a system. For studying large molecules that interact through weak forces (non-covalent interactions), like proteins binding to drugs (protein/ligand complexes), a modified version of this method, named DFT-D (dispersion-corrected DFT) is now the standard method, but the method is still very slow [46]. Other widely used DFT approaches are those based on generalized gradient approximation (GGA and hybrid-GGA), such as EDF1, EDF2, BP, BP86, BLYP, M06, B97X-D [18], BPE [47] and B3LYP [48]. Among the extensive and expanding collection of DFT functionals, B3LYP has been the most frequently employed functional in the field of biochemical research for the past several decades [49].

Much like semiempirical methods, DFT (the most commonly used ab initio method for biological systems) appears to be the predominant choice for studying natural compounds. In the literature published up to date DFT has been used rather to assess the chemical stability and reactivity of certain compounds [50,51,52,53,54,55] or to optimize ligands before docking procedures [52,56], than for the evaluation proper of ADME properties.

Hybrid QM/MM (quantum mechanics/molecular mechanics) methods are an effective approach that creatively integrates quantum mechanical and classical molecular mechanical calculations. First introduced by Warshel and Levitt in 1976, this method has since been extensively used in computational molecular biology, drug design, nanotechnology, and quantum chemistry. It was also a key factor in the 2013 Nobel Prize in Chemistry [57]. In this approach, quantum mechanics is applied to a small region that requires great accuracy (usually the reactive site), while computationally efficient molecular mechanics modelling is used to simulate the adjacent environment. Since MM computations costs are virtually minimal in comparison to QM, the overall cost is mostly dictated by the size of the QM part. This makes QM/MM approaches especially useful for analyzing large biological systems such as enzymes or protein receptors, where just the active region requires quantum mechanical treatment and the protein framework can be addressed using molecular mechanics [16]. The hybrid QM/MM method enables researchers to simulate large macromolecular systems while maintaining the accuracy of detailed modeling. This approach helps overcome the size limitations of biological systems, which usually have a large or very large number of atoms (quantum mechanics methods are very accurate, but are feasible only for small chemical structures) [58].

The following factors must be addressed in order to utilize QM/MM calculations for biomacromolecules:

a)

Defining the QM region: the accuracy of QM/MM calculations is generally enhanced when the QM region includes larger and less polar functional groups at the frontier between the QM and MM zones;

b)

Handling the QM/MM interface: in biomacromolecular studies, covalent bonds at the QM/MM interface, such as those between protein residues in a receptor, require careful attention. Errors introduced by embedding the QM region into the classical force field of the bulk protein are minimized when the net atomic charges of the boundary atoms are low.

c)

Addressing multiple local minima: the challenge of numerous local minima can be mitigated by performing extensive sampling of protein and ligand conformations and by ensuring the biomacromolecular system is accurately prepared and modeled [59].

While hybrid QM/MM methods are being used to some extent for studying the pharmacodynamics of natural products [60,61,62,63] and can be also employed for ADME assessment applications, we haven’t found any reports of their use for this purpose.

2.2. Molecular docking

Since its emergence in the mid-1970s, docking has played a pivotal role in drug discovery and development, as well as in understanding how chemical compounds interact with their molecular targets. Over time, there has been a substantial rise in studies utilizing molecular docking to uncover the structural features required for effective ligand-receptor binding, alongside advancements in refining docking techniques for greater accuracy [64]. Molecular docking has emerged as the method of choice when researchers have access to the three-dimensional crystal structure of their target protein. It has seen widespread adoption thanks to two key developments: the dramatic increase in computing power and resources, coupled with the exponential growth of structural databases containing both small molecules and proteins. These factors have made docking both more accessible and more powerful as a research tool [65].

Molecular docking is a modeling technique that examines the interactions between a macromolecule (such as a protein, enzyme, or DNA) and a small molecule (ligand, i.e., the potential drug ingredient). Molecular docking primarily has three objectives: virtual screening, binding affinity calculation, and prediction of the specific spatial position, conformation and orientation of the ligand [66]. Sampling in molecular docking is challenging due to the vast conformational space, which includes rotational, translational, and internal flexibility of both the ligand and protein, as well as, in some cases, solvent effects. Exhaustively exploring all possible conformations and orientations within the computational time required for virtual screening is currently unfeasible. As a result, finding efficient sampling methods in docking studies is still an area of active research [65].

Docking studies can be categorized into three types based on molecular flexibility: rigid docking (both molecules fixed), flexible-rigid docking (one molecule flexible, the other fixed), and flexible docking (both molecules flexible) [67]. The first docking methods were based on the idea that molecules have fixed shapes, like a lock and key. This meant that both the receptor (a protein) and the ligand (a drug) were treated as rigid. However, it was later accepted that molecules can change their shapes when they interact (the “induced-fit” model). Because the backbone of a molecule influences the positions of many side chains, allowing both the receptor and ligand to be flexible makes the docking calculation much more complex. However, these flexible docking methods are superior because they more accurately predict not only how molecules fit together but also how strongly they bind to each other [68].

Docking is not regarded mainly as a tool particularly apt for ADME purposes, and there have been literature reviews that have not even considered ADME as a potential application of docking [64,67]. Molecular docking can, however, be used to evaluate the interaction between natural compounds (as well as synthetic or semisynthetic ones) with a variety of proteins involved in different steps of the ADME processes. Thus, a variety of natural compounds have already been evaluated through molecular docking and reviewed for their ability to interact with various transporter proteins relevant for the absorption, excretion or cell efflux [69,70,71]. The effects of 75 flavonoids on the P-glycoprotein have also been evaluated through in vitro and in vivo models, after which the authors have used molecular docking to better understand the nature of the interactions between the five flavonoids that demonstrated significant inhibition and the transporter protein [72]. In addition to in vitro uptake tests, molecular docking was used to examine the selectivity of 10 natural-origin inhibitors of the organic cation transporter 1 (OTC1), and benzoylpaeoniflorin was found to be the most selective [73].

The interactions between natural products and proteins involved in the drug distribution step, such as albumin, α1-acid glycoprotein or hemoglobin can also be investigated in silico using molecular docking. This is how, for instance, the interactions between enantiomers of stipuol (a polyacetylene from Panax notoginseng (Burkill) F.H.Chen) and human serum albumin have been evaluated, the authors finding that the natural compounds bind at subdomain III of albumin [74].

In a similar way, the interactions between natural compounds and metabolizing enzymes such as the cytochrome P-450 family of proteins can be evaluated using molecular docking and a review of this approach for drugs in general (not necessarily natural products) has also been published [75]. Virtual screening models based on machine learning have been used to identify selective CYP1B1 inhibitors (different libraries of compounds, including one of natural products), and the compounds thus identified have then been evaluated for their selectivity against other isoforms by using molecular docking [76]. Inhibitors of other CYP450 isoforms have also been evaluated using molecular docking [77,78,79]. Drug transporters in the kidney are particularly relevant for drug excretion and the interactions of numerous natural compounds with such transporters have been explored through molecular docking [69].

2.3. Pharmacophore modeling

The notion of „pharmacophore” was pioneered by Ehrlich in 1909, who described it as “a molecular framework that carries (Gr. phoros) the essential features responsible for a drug’s (Gr. pharmacon) biological activity” Over the past century, the concept has remained, but its scope and significance have been remarkably broadened. Currently (according to IUPAC) a pharmacophore is understood as the complex combination of spatial and electronic features that are essential for optimal interactions with a specific biological target, ultimately initiating or inhibiting a biological response [80]. The pharmacophore model relies on the interaction between ligands and receptors. These interactions encompass all information pertaining to the structural, spatial, and chemical features that underlie specific pharmacological effects. The interaction mostly entails non-covalent bonding, including hydrogen bonding, pi-pi stacking, and ion-dipole interactions, among others [81].

Pharmacophore models are often developed using chemical features found in common among known ligands for a specific target. These are known as ligand-based pharmacophores. Obviously, when only one or several ligands are known for a specific target, or when the available ligands are devoid of structural diversity, ligand-based pharmacophores cannot be used. In such cases, the so-called (protein) structure-based pharmacophore (SBP) models are used [82]. SBP modeling has become increasingly prominent in recent years, driven by the significant growth in high-resolution protein structures. As of January 2025, the Protein Data Bank (PDB) houses approximately 123,000 three-dimensional structures of biological macromolecules, primarily proteins, providing an unprecedented foundation for understanding these molecular therapeutic targets [83]. If experimental structural data for a particular protein are not available, computational methods such as homology modeling can provide an alternative approach to generate a 3D model. Machine learning techniques have also proven effective in predicting protein structures, AlphaFold2 being probably the most widely known example [84].

Developing an effective pharmacophore model demands a deep knowledge of the target protein’s structure and binding characteristics. Proteins commonly feature multiple binding sites, and within each site, different ligands can adopt various binding orientations by interacting with distinct regions of the binding pocket. This complexity is further increased by protein subfamily variations, such as those seen among Src kinases, where subtle structural differences between related proteins can significantly impact pharmacophore selectivity [85]. The initial step in structure-based pharmacophore modeling involves selecting and preparing the target protein structure (a step similar in approach and importance for other computational approaches, too, such as molecular docking). This is followed by identifying potential binding sites on the protein. A detailed analysis is then conducted to determine the complementary chemical properties and spatial arrangement of the amino acids at the binding site. Based on this information, pharmacophore features are generated and refined using specialized tools available within the selected software. Finally, the critical pharmacophore features responsible for the biological activity are identified [86]. A detailed description of these steps can be found in the referenced publications and interested readers are encouraged to consult these sources for further information [83,84].

A pharmacophore model has been used recently in the study of a number of natural compounds (mainly flavonoids) to explore their ability to inhibit organic anion transporter 3 (OAT3) [87], organic cation transporter 1 [88] and organic cation transporter 2 [89]. Although the focus on those transporters was rather related to safety endpoint (prevention of nephrotoxicity or hepatotoxicity), the same approach can be equally relevant for ADME purposes, as biological transporters are involved in ADME processes and drug interactions related to the excretion of certain active substances [90]. Pharmacophore models were also used to evaluate the effects of a number of flavonoids on urate transporter 1 (URAT1) [91].

Pharmacophore models have also been used to explore the effects of certain natural products on specific CYP450 fractions, such as those of saponins [92] or flavonoids [93] on CYP3A4 or those of a variety of natural compounds on CYP1A1 [94,95].

2.4. Quantitative Structure-Activity Relationship (QSAR) models

QSAR modeling is a computer-based method that establishes empirical correlations between chemical descriptor values derived from molecular structures and experimentally obtained properties or bioactivities. These models are then used to predict or design new chemicals with specific properties [96]. In other words, the chemical structure of a substance can be described by a series of numeric features named generically “molecular descriptors” (e.g., molecular weight, number of double bonds, percentage of oxygen atoms, and others much more sophisticated) and QSAR establishes mathematical relationships between these numeric values and certain biological activities of interest (for instance, the inhibition or not of a certain CYP450 isoform, the percentage of oral bioavailability etc).

The use of QSAR has undergone significant changes in recent decades, which include: (1) using more detailed descriptions of molecules, from simple one-dimensional (1D) representations to more complex multi-dimensional (nD) ones; and (2) employing more sophisticated methods, including machine learning, to correlate molecular structure with biological properties. Initially, QSAR studies were limited to small groups of closely related chemicals and relied on basic regression analysis. Now, QSAR can handle much larger and more diverse datasets, using a wide range of machine learning techniques for both modeling and virtual screening [97].

Based on the assumed relationship between the chemical descriptors and the biological effect modelled, QSAR models can be classified into linear (e.g., based on multiple linear regression or partial least squares) and non-linear (e.g., support vector machines or artificial neural networks) [98,99].

Based on the mode in which descriptors are computed and models are built, QSAR models with different dimensions are distinguished: 2D-QSAR, 3D-QSAR, and more recently 4D-, 5D-, 6D-QSAR models have been proposed. The 2D-QSAR models are based on descriptors computed based on 2D representation of ligand structures (but ignoring completely the 3D shape and spatial arrangement of atoms within the ligand) [100]. 3D-QSAR correlates biological target properties with molecular descriptors derived from 3D chemical structures; it uses probe-based sampling within a molecular lattice to compute the 3D molecular descriptors [101]. In a paper published in 1997, Hopfinger et al. proposed for the first time the concept of 4D-QSAR [102]. The fourth dimension, known as ’ensemble sampling,’ consists in a set of molecular configurations, describing ligands in different molecular and spatial forms, such as conformations, orientations, stereoisomers, and protonation states [100,103]. First proposed by Vedani and Dobler in 2002 [104], the 5D-QSAR approach addresses the induced-fit phenomenon in the ligand-protein interactions by considering multiple conformations that a ligand can adopt when binding to its target. This flavor of QSAR, which builds upon 4D-QSAR, creates virtual representations of these different binding scenarios to improve the predictivity of the model(s) [105]. 6D-QSAR builds upon previous QSAR versions by adding a sixth dimension: solvation. It allows the simultaneous evaluation of multiple solvation models to better understand how solvents influence molecular interactions. This can be carried out explicitly by mapping solvent properties onto the molecular surface or implicitly by adjusting the contribution of solvation terms (like ligand desolvation and solvent elimination) for different models within the analysis. During the simulation, weights associated to different conformations are adjusted through an evolution process [106].

Mechanistic QSAR is an approach that employs parameters/descriptors that have a consistent and unambiguous meaning, aiming to provide insights into the underlying mechanisms of the predicted biological effects. The alternative, named sometimes non-mechanistic QSAR, prioritizes the development of equations with predictive value, even if the parameters used are not readily interpretable. This approach focuses on statistical correlations between molecular descriptors and biological activity, but it does not necessarily offer mechanisti insights into the relationship between the descriptors and the biological effects (e.g., Comparative Molecular Field Analysis - CoMFA) [107,108].

A global QSAR model has been defined as a model characterized by a dataset that is highly structurally diverse, which in turn reflects a variety of mechanistic actions among the chemical compounds of the dataset. In such a model statistical methods are employed to attempt to identify structure/activity patterns that are independent of the mechanistic differences [109]. It has been stated that despite the fact that statistics might indicate that these models are robust, they tend to exhibit relatively low predictivity. Local models, of the contrary, are characterized by structural or mechanistic similarity [109]. The heterogenous definition of global QSAR models is somewhat problematic, and it could gain clarity by distinguishing between structurally global models (models built on a variety of structurally diverse compounds) and mechanistically global models (models built on chemical compounds acting through a variety of mechanisms and on distinct biological targets). There are theoretical reasons to believe that the latter have poorer performances than the former.

Sometimes QSAR models are qualified by the nature of the chemical descriptors used in their construction, for instance “topological QSAR” models (based on topological descriptors) [110,111] or “quantum mechanical QSAR” models (based on quantum mechanical descriptors) [112]. The “hologram QSAR” technique involves the decomposition of each molecule in the dataset into a molecular hologram. This hologram is primarily composed of linear, branched, and overlapping fragments, which are then arranged in an array of fixed length (ranging from 53 to 401 bins). These bin values serve as X variables in QSAR modeling, capturing both the composition and topology of molecules. Several factors affect how the hologram is generated and how the resulting HQSAR models perform statistically, including: the length of the hologram, the size of fragments, and how fragments are distinguished (based on atom types, bonds, connections, hydrogen atoms, chirality, and groups of hydrogen bond donors or acceptors) [113].

QSAR models can also be classified based on the categorical or continuous nature of the predicted outcome, distinguishing between classification and regression models. Classification models use chemical information of the training data set to make predictions about an outcome category, most often using a binary category (e.g., active/inactive, orally bioavailable/non-bioavailable etc). Often such a classification involves an arbitrary conversion of a continuous variable in two categories at a certain threshold. It is argued that in such a case, “the predictive reliability of the resultant QSAR model will increase at the expense of the resolution of the prediction” [114]. Regression QSAR models make predictions about the outcome variable on a continuous scale (predicting, for instance, a numerical value of ki, IC50 or of a percentage for the oral bioavailability and not a mere category) [115].

Biological effects that can be predicted using QSAR models can be vastly diverse, from the effects on a specific biological target to the prediction of a specific adverse effect or a specific ADME feature. Among the in silico methods for ADME predictions, probably QSAR models have been used the most extensively, and generally the same models can be applied for both natural and synthetic organic chemical compounds, although often the specific use of such models for natural compounds has not been explicitly mentioned. In the case of intestinal absorption prediction, the performance of QSAR models evaluated in a large benchmark study was only slightly inferior to that of in vitro methods (83% of QSAR predictions and 87% of in vitro method predictions fell within 2-fold of observed values [116]. Similarly, an earlier study indicated that a pair of computer models achieved reliability comparable to the Caco-2 and 2/4/A1 cell lines, with one model predicting the absorption of a collection of 65 medicines nearly as well as in vivo absorption experiments conducted in rats [117].

Kansy et al. introduced the parallel artificial membrane permeation assay (PAMPA) in 1998 [118]. PAMPA is composed of lipophilic filters that are coated with lecithin in an organic solvent solution. Transcellular permeation is assessed by a rapid in vitro assay system that is suitable for high-throughput screening. Akamatsu et al. (2009) constructed an in silico prediction model of human oral absorption for potentially transported compounds by developing a QSAR model that describes correlations between chemical structure and PAMPA permeability [119]. This study has not mentioned explicitly the use of the QSAR models for natural compounds, but similar QSAR models based on PAMPA have been developed and one has been used to predict the gastro-intestinal absorption of flavonoids from Silybum marianum (L.) Gaertn. [120].

Oral clearance refers to the rate at which a substance/drug ingredient is removed from the body following oral administration. Boik and Newman (2008) developed three QSAR classification models based on data from human oral clearance, rodent LD50, and in-vitro cytotoxicity studies. Then they used these models to analyze over 115,000 natural compounds, resulting in the prediction that hundreds of these substances exhibit low systemic toxicity, low to moderate oral clearance, and manageable cytotoxicity levels [121]. There are numerous other QSAR/QSPR models that have been developed to predict oral absorption of drugs, although they have not necessarily been applied to natural compound datasets (Table 1).

QSAR models have also been developed to predict plasma protein binding [145,146,147,148] or the unbound fraction [149], the volume of distribution in humans, rodents or other mammals [150,151,152], the blood-brain barrier (BBB) permeability [153,154,155], inhibition of specific CYP450 fractions [156,157,158,159], elimination half-life (149), rodent (150,151) and human clearance (152), effect of substances on key transporters such as P-glycoprotein [160,161,162], organic anion transporters [163,164,165], organic cation transporters [166] etc. Some of the models referenced here have been specifically developed for natural products (for instance, Li et al. (2018) explored the quantitative relationship between the structural characteristics of flavonoids and their ability to inhibit CYP3A4 [158]), but most of the published models have been developed for medicines or organic substances irrespective of their natural or synthetic origin.

Certain molecular descriptors can be computed/estimated using molecular dynamic simulations, case in which QSAR and MD simulations are integrated in a single final model. For instance, Iyer et al. (2007) have used MD to simulate the behaviour of different ligands across a dimyristoylphosphatidylcholine (DMPC) monolayer membrane and compute a number of descriptors used then to build QSAR models predicting the intestinal oral absorption of the drugs included in the models. The test and train compound datasets included both natural and synthetic substances [132]. Other relevant descriptors can be computed using the DFT theory, case in which QSAR and DFT are integrated in a single model (see, for instance, ref. [130]). Another hybrid approach has been recently (year 2024) reported for 48 compounds, mostly natural (including alkaloids, flavonoids and coumarins). The authors built a QSPR model to predict the apparent permeability coefficient (Papp) (measured on Caco-2 cells), then used MD simulation to evaluate the molecular mechanisms involved in the absorption differences observed for ligustrazine and EGCG, as well as molecular docking to understand the interactions between the compounds and the P-glycoprotein [167].

2.5. Molecular dynamics (MD) simulations

X-ray crystallography provides static snapshots of molecules because they are locked in a crystal. This led to an initial perception that large molecules like enzymes and receptors were rigid. However, combining multiple snapshots often reveals that these molecules exist in different states, indicating they are actually dynamic and fluctuate between these states [168]. Molecular dynamics simulations use physics-based modeling to predict how individual atoms move within proteins and other molecular systems over time. These highly precise simulations can track atomic positions down to the femtosecond level, allowing researchers to acquire knowledge about important biological processes like protein folding and ligand binding. The simulations provide valuable insights into how biomolecules respond at the atomic scale to various changes like mutations or chemical modifications. This computational approach has become a powerful tool for understanding molecular behavior and interactions [169]. The method was for the first time applied to proteins at the end of the 1970s [170]. The applications of MD simulations in the fields of molecular biology and drug discovery are vast, but we are here focused on their potential use in making ADME predictions. It seems likely that the most important such applications consist in improving molecular docking results for the interactions between different ligands and the proteins of interest, relevant for pharmacokinetics purposes.

While molecular docking can predict how a protein and a ligand interact, the initial prediction of how the ligand is positioned and oriented within the protein’s binding site is often approximate. To get a more accurate picture of the interaction, further refinement is necessary. MD simulations can refine the docking results by minimizing steric clashes (where atoms are too close together) and adjusting the initial binding position of the ligand [171].

To enhance the effectiveness of virtual screening, MD simulations can be combined with binding free energy calculations to evaluate the binding energies of small molecules with a target (thus understanding how strongly a ligand binds to the target of interest). These energies can then be used to rank candidates and refine the scoring of generated poses. Methods like MM-PBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) and MM-GBSA (Molecular Mechanics Generalized Born Surface Area) enable more precise calculations of the binding free energy between the ligand and receptor. Entropy changes are estimated through normal mode analysis, while enthalpy changes for the entire system are statistically derived from electrostatic, van der Waals, and solvation energies, all based on MD simulation trajectories. This approach allows for the accurate determination of binding free energies [171]. Whereas MM-PBSA and MM-GBSA are not themselves MD simulations, they use the MD-generated data to estimate more or less accurately the free energy and thus to better understand the ligand-target interaction.

All ligand-protein interactions that can be studied for ADME purposes using molecular docking can also be better understood using MD simulations, and the applications of the former are also applications of the latter. But MD simulations have also been used for certain processes relevant for ADME but less investigated with molecular docking. The most important such application is the study of membrane permeability. In 1994, Marrink and Berendsen performed the first atomistic simulation of water permeability in a lipid bilayer. Since then, hundreds of coarse-grained and atomistic simulations have been carried out to study the passive transport of permeants across different lipid bilayers and synthetic membranes [172]. Molecular simulations are useful for studying both synthetic and natural substances. Two studies specifically investigated how menthol enhances the skin’s ability to absorb substances [173,174]. These studies showed that menthol readily penetrates the skin’s protective outer layer (lipid bilayer), making its interior more fluid. However, at high concentrations, menthol can damage this layer. Borneol, another natural compound, at low concentrations helps substances penetrate the skin, but at high concentrations creates pores and unusual structures called reverse micelles, a mechanism similar to how some antimicrobial peptides disrupt cell membranes (the “carpet mechanism”) [172,175,176].

2.6. Physiologically-Based Pharmacokinetics (PBPK) Modeling

PBPK models simulate how a drug and its metabolites change in concentration over time in the blood and target organs. Building PBPK models, however, requires more effort and data compared to simpler predictive methods. They necessitate estimating a larger number of parameters and demand extensive knowledge of the body’s physiology and the specific drug’s properties [177]. The model’s parameters are obtained either from non-clinical experimental data, clinical trials or calculated based on appropriately validated formulas. Broadly, these parameters can be categorized into two distinct groups: drug-specific parameters and organism-specific parameters. Although a PBPK model is usually very complex, with hundreds of ordinary differential equations, a large number of parameters, particularly organism-specific ones, are already known (e.g., organ blood flows, tissue volumes, membrane permeabilities etc); therefore, for a new medicine, often less than five drug-parameters are needed (e.g., clearance, volume of distribution, binding parameters etc) [178].

PBPK models estimate drug kinetics in one or several organs, necessitating a submodel for each organ or tissue considered and integrate then the organ submodels in a whole-body, comprehensive model. Drug transport occurs via blood circulation as determined by anatomical structures. PBPK strive to replicate authentic, quantifiable physiological and/or pharmacokinetic processes instead of more abstract rate constants and volumes. Parameters used in PBPK models may be ’intrinsic’ covariates, such as weight, body surface area or glomerular filtration rate, rather than covariates identified post-hoc to improve model fit (such as, for instance, the patient age, in the absence of a clear mechanism or explanation of why it would improve fit) [179]. A PBPK model for a specific medicine can be validated using both in vitro data and in vivo oral pharmacokinetic data from rats. To predict the human oral pharmacokinetic profile, the physiological and in vitro parameters from the rat model can be substituted with human-specific parameters, following a successful prediction of the rat’s in vivo oral pharmacokinetics from in vitro data [180].

Compartmental models, widely used in the traditional pharmacokinetic research, are built using a “top-down” method, where all the model’s information is derived from a specific set of pharmacokinetic and covariate data. Conversely, PBPK models tend to use a “bottom-up” approach, where the model’s information is derived from in vitro data and pre-existing knowledge of physiological and pharmacological mechanisms involved in the drug pathways through the body, the model parameters being thus derived from mechanistic principles (the focus is not so much on merely fitting the data as to providing a mechanistic understanding of the drug behavior) [179].

PBPK models offer valuable predictions for various drug development purposes:

• At the clinical trial design stage, they can predict how drug formulation and food intake will influence pharmacokinetics, guiding initial human studies and forecast drug-drug interactions mediated by enzymes or transporters, informing inclusion/exclusion criteria, dose selection, and potentially waiving unnecessary clinical interaction studies or studies where enrolling subjects is anticipated to be difficult;
• They can predict appropriate dosing regimens for different pediatric subsets, from newborns to adolescents, by enabling informed selection of sampling timepoints and proposing suitable doses.
• They can predict exposure to the drug in patients with impaired renal or hepatic function, guiding organ impairment studies or supporting decisions to waive such studies.
• They can estimate the drug disposition in the mother and fetus, aiding in optimizing the therapeutic benefit-risk ratio during pregnancy.
• They can predict pH-mediated drug-drug interactions in patients receiving proton pump inhibitors or antacids, guiding formulation development and efforts to minimize food-drug interactions [181].

It is not feasible to cover technical details on building PBPK within the limited space of this paper, but they are comprehensively addressed in the literature, with several high-quality tutorials providing valuable guidance on different aspects of the modelling process [178,182,183].

A limited number of studies have used PBPK modelling up to date to approach ADME properties for natural products. For instance, a recent study reported its use to estimate the biodistribution of Centella asiatica administered orally in a rat model [184]. Another study used a PBPK model to simulate the drug-drug interaction between aconitine and the glycyrrhizic acid, using parameters derived from in vitro experiments. The model simulations were validated against experimental observations to evaluate their accuracy. The analysis focused on understanding how CYP3A and P-glycoprotein contribute to the reduction in aconitine toxicity induced by the glycyrrhizic acid [185]. Similarly, PBPK models were used to explore the interactions between silybin A and losartan [186], between schisandrin A, schisantherin A and cyclophosphamide [187], the drug-drug interactions involving bergamottin (from grapefruit juice), curcumin (from turmeric) and hyperforin (from St. John’s wort) on the one hand, and anticancer drugs on the other [188], or the interaction between rifampicin and retrorsine (a pyrolizzidine alkaloid occurring occasionally as a contaminant in herbal teas) [189]. A PBPK model was used to show that cichoriin is capable of reaching high concentration levels in the lung (when administered intravenously) and thus, that it is a potential candidate for development against COVID-19 [190] and another to predict the plasma concentration variation for schaftoside, subsequent to the ingestion of a capsule containing total flavonoids extracted from Desmodium styracifolium [191].

3. Comprehensive ADME Tools

Whereas many of the approaches discussed above involve the building and validation of a single model or a family of models that are then used “inhouse” for the virtual screening of natural or synthetic compounds, there are a number of software tools that have been developed specifically for in silico ADME predictions. They are either available as server services or as standalone products, for free or for a fee. Although a larger number of commercially or freely available software is available for ADME predictions, our search of PubMed using appropriate keywords (“natural” OR “herbal” plus the name of each software product, and when not very specific, also the name of the manufacturing company) has shown that for natural products only a part of the available software has been used, and only a few have been applied in multiple papers in the field of natural products: SwissADME webserver[192] (> 250 papers), pkCSM [193] (> 60 papers), ProTox-II/Protox 3.0 [194,195] (> 45 papers), QikProp [196] (part of the Schrödinger Suite – > 40 papers), ADMETSAR [197] (> 39 papers), ADMETlab 2.0/3.0 [198,199] (> 30 papers), ADMET Predictor [200] (> 12 papers), FAF-Drugs4 [201] (6 papers), ACD/Percepta [202] (4 papers), ADMET module of Biovia Discovery Studio (2 paper [71]). Admetboost [203] and Interpretable-ADMET [204] have been relatively recently published and made available for free, but no paper has apparently used them as yet in the assessment of natural products. We have limited the discussion to general-purpose ADMET tools, but specialized tools for specific ADMET are also available; for instance, a review of software applications used to specifically predict CYP450-mediated drug metabolism is available elsewhere [205].

In what follows we present in more details the most widely used such free application (SwissADME) and the most widely used commercial application (QikProp). We will then provide concise overviews of other commonly utilized tools, including pkCSM, Protox 3.0, ADMETSAR, and ADMETlab.

3.1. SwissADME

The SwissADME web application (available http://www.swissadme.ch/) is a free and user-friendly platform used to predict ADME properties of potential drug candidates. No computer-aided drug design (CADD) expertise is required. Its key features include:

• Simple Submission: single or multiple molecules can be uploaded for analysis with ease.
• Clear Results: straightforward visualizations are available for the predictions made for each molecule.
• Data Sharing: the results for individual molecules can be saved and shared, or can be used to generate comprehensive, interactive graphs for deeper analysis.
• Integrated CADD Tools: the application seamlessly connects to the SwissDrugDesign workspace, granting the user access to a suite of other useful applications developed by the SIB Swiss Institute of Bioinformatics, to identify promising drug candidates based on their similarity to known active compounds (ligand-based virtual screening), to predict potential biological targets for a specific ligand, to perform ligand-protein docking, to create bioisosteric replacement of functional groups, or to perform molecular mechanics (e.g., to analyze the 3D structure and energy of a molecule) [192].

The following properties and values are predicted by SwissADME:

3.1.1. Biovailability radar

A bioavailability radar map is generated by the server application for a swift evaluation of a molecule drug-likeness. For this purpose, six physicochemical features are used: size, flexibility, polarity, lipophilicity, solubility, and saturation [192].

3.1.2. Physicochemical properties

A small number of physicochemical properties are computed by the application using OpenBabel (molecular weight, molecular refractivity, counts of heavy and aromatic atoms, counts of H-bonds donors and acceptors, as well as the polar surface area, PSA, the latter being known to be often a good predictor of oral absorption or of the blodd-brain barrier permeability) [192].

3.1.3. Lipophilicity

The n-octanol - water partition coefficient (log Po/w) is a classic measure of lipophilicity, playing a crucial role in the pharmacokinetic behavior of drugs and drug candidates, being regarded as “the single most informative and successful physicochemical property in medicinal chemistry” [206]. SwissADME offers five predictive models to calculate this parameter and a “consensus log Po/w” is derived by taking the arithmetic mean of the predictions generated by these five approaches, but the five individual results are also returned [192].

3.1.4. Water solubility

Three methods are used by the application to predict water solubility, presenting the results in a variety of formats: logarithmic scale (the decimal logarithm of the molar solubility - log S -, allowing comparisons of values across a wide range), molar and mass concentrations (mol/l and mg/ml), and qualitative classes (e.g., “soluble”, ”moderately soluble”, “poorly soluble” etc) [192]

3.1.5. The skin permeability coefficient

Skin permeability is crucial in a variety of applications, such as the risk evaluation of hazardous compounds, the design of skin cream products, and transdermal medication delivery [207]. The assessment of skin permeability often involves measuring the skin permeation coefficient (Kp) of a chemical through the stratum corneum, which is the top layer of the epidermis [208]. SwissADME uses a model adapted from Potts and Guy (1992), based primarily on the molecule lipophilicity and its size [192,209]. For this model, the authors have claimed an r²=0.67, but in a benchmark study on 30 molecules the correlation coefficient was only 0.544, whereas QikProp predictions achieved a correlation coefficient slightly higher (0.612) [210].

3.1.6. Passive HIA and BBB permeability

Both these PK variables are predicted by SwissADME using the BOILED-Egg model, which is not based on regression, but rather on a classification algorithm, providing results in a binary form (GI absorption high vs. low, and permeant vs. non-permeant, respectively) [211].

3.1.7. Pgp substrate property

P-glycoprotein (ABCB1) is an ATP-dependent efflux pump that transports a large variety of structurally diverse hydrophobic and amphipathic compounds. These substances include pharmaceuticals, peptides and lipid-type molecules. The pump plays a critical role in the efflux of these compounds from cells due to its broad substrate specificity [212]. SwissADME uses a machine-learning algorithm based on support-vector machines (SVM) to classify the evaluated compounds as substrates of Pgp or non-substrates [213].

3.1.8. CYP fraction substrate

The different fractions of CYP450 enzyme system metabolizes over two-thirds of xenobiotics, therefore being able to predict whether a potential drug or toxic will be a substrate of one or several of those fractions is important and useful [205]. Similarly to the Pgp substrate predictions, SwissADME uses several classification models to predict whether a substance is or not a potential inhibitor of five CYP450 fractions: CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4 [192].

3.1.9. Drug likeness and medicinal chemistry properties

Besides the ADME properties mentioned above, the web server application also computes several drug-likeness scores widely used in the field of drug discovery (Lipinski [214], Ghose [215], Veber [216], Egan [217], Muegge [218]), as well as the Abbot bioavailability score [219]. Such scores are relevant in the understanding ADME properties of a substance, because compounds with appropriate scores of drug-likenss/druggability also have desirable ADME characteristics [11]. Several medicinal chemistry alerts, filters or scores are also implemented: PAINS alerts [220], structural alerts (e.g., potentially mutagenic groups, other groups with potential negative effects on PK or very reactive groups) [221], a lead-like filter [222], and a synthetic accessibility score (ranging between 1 and 10, from very easy to very difficult).

3.2. QikProp

Schrödinger’s proprietary software, QikProp, is widely used to generate ADME predictions for various drug candidates, including natural products [196], but its use seems to become outpaced (at least in the field of natural compounds) by free webservers such as SwissADME or pkCSM. The software estimates a total of 46 molecular descriptors relevant to understanding ADME properties. One of the descriptors calculated by QikProp is a “drug-likeness parameter” represented by a “star” rating system. This parameter serves as an overall ADME compliance score, where the number of “stars” assigned corresponds to the number of descriptors that that are outside the range of values that are considered optimal for 95% of known medications. The higher the number of “stars,” the greater the deviation from the typical drug-like profile, potentially indicating poorer ADME characteristics for the compound [196]. Jorgensen and Duffy [223,224,225] are credited as the creators of this approach [196]. This global ADME score generated by QikProp is used to help prioritize and select promising drug candidates for further development, based on their predicted ADME suitability. QikProp predictions are calculation-based, as opposed to fragment-based. Fragment-based methods may be problematic when they fail to identify specific components of a structure or encounter unfamiliar fragment interactions. In contrast, QikProp will calculate properties based on the entire molecule. This allows QikProp to make its predictions on scaffolds that are new and unfamiliar [226].

Among the calculated descriptors are the following [196]:

3.2.1. Total Solvent-Accessible Molecular Surface (Smol or SASA, solvent-accessible surface area)

It is expressed in ų and is calculated using a probe radius of 1.4 Å (the approximate radius of a water molecule). Approximately 95% of pharmaceutical compounds exhibit S mol values between 300–1000 Ų. It quantifies the surface area of a molecule accessible to a solvent (usually water) and reflects the molecule’s hydrogen bonding capacity. It is an important property useful and used in predicting solubility and permeability of a potential drug [227,228].

3.2.2. Hydrophobic Portion of the Solvent-Accessible Molecular Surface (Smol, hfob)

It is expressed in ų and is calculated using a probe radius of 1.4 Å. Approximately 95% of pharmaceutical compounds exhibit S mol values between 0 and 750 Ų. Also known as FOSA, it describes the nonpolar interactions of the ligand, being useful in estimating the permeability capability of a small molecule [229].

3.2.3. Total Volume of Molecule Enclosed by Solvent-Accessible Molecular Surface (Vmol, hfob).

It is a measure used to quantify the volume occupied by a molecule within its solvent-accessible surface. It is expressed in ų and computed using a probe radius of 1.4 Å. Approximately 95% of pharmaceutical compounds exhibit V mol, hfob values between 500 and 2000 ų [196,230].

3.2.4. Logarithm of Aqueous Solubility (log Swat)

It is the logarithm to base 10 of a molecule’s solubility, measured in mol/L. In drug development, the ability to dissolve in water (aqueous solubility) is a key physical and chemical property that needs to be optimized, particularly for aspects ’A’ and ’D’ in ADME-Tox. This solubility, along with the capacity to penetrate membranes (membrane permeability), plays a major role in determining how well a drug is absorbed into the body when taken orally (oral bioavailability) [231]. Approximately 95% of pharmaceutical compounds exhibit log S wat values between -6.0 and 0.5 [196,230].

3.2.5. Predicted octanol / water partition coefficient (QPlogPo/w).

Recommended values for this parameter (see als section 3.1.3), as satisfied by 95% of pharmaceutical compounds, are between –2.0 and – 6.5 [232].

3.2.6. Logarithm of Predicted Binding Constant to Human Serum Albumin (log K HSA or logKhsa)

Approximately 95% of pharmaceutical compounds exhibit logKhsa values between -1.5 and 1.2. It is used to predict a compound’s binding affinity to plasma proteins such globulins, human serum albumin, lipoproteins, and glycoproteins [229].

3.2.7. Logarithm of Predicted Blood/Brain Barrier Partition Coefficient (log B/B).

In pharmacokinetics, BBB distribution is quantitatively evaluated using the blood-brain partition coefficient, mathematically represented as the logarithm of the concentration ratio between brain tissue and blood (and referred to as log B/B) [233]. A positive logB/B signifies a greater concentration of the compound in brain tissue compared to blood, suggesting effective brain penetration, whereas a negative logB/B indicates limited brain accessibility. Approximately 95% of pharmaceutical compounds exhibit log B/B values between -3.0 and 1.0 [229] (this means that drugs with excellent blood-brain barrier penetration can reach levels about 10 times higher in the brain than in the blood, whereas for drugs with low blood-brain barrier penetration, the concentration in the brain is only 1/1000th of the concentration in the blood).

3.2.8. Predicted Apparent Caco-2 Cell Membrane Permeability (BIP caco-2)

Permeability assays are an essential part of ADME profiling, as adequate permeability across the intestinal membranes can improve oral absorption, which is a highly sought after property in most drug discovery initiatives. The Caco-2 cell permeability model is useful because it is human-derived and expresses transporters common to the human gastrointestinal tract [234]. The model can predict how well a drug can pass through the gut wall into the bloodstream. This parameter reflects (mainly) non-active transport across gut barrier, and is influenced by factors such as the drug’s solubility and its ability to cross cell membranes. However, the absorption of drugs can also be affected by interactions with transporters and metabolizing enzymes present in the gut wall [235]. Approximately 95% of pharmaceutical compounds exhibit BIP caco-2 values between less than 5 (low) and higher than 100 (high). It has been stated that if BIP caco-2 is > 22 nm/s, there is a higher probability that a drug will be orally bioavailable [196]. However, this descriptor, even if experimentally measured (computational methods may have larger margins of error), is only of a limited use in predicting oral bioavailability. Broadly speaking, it seems that once the permeation rate exceeds a certain threshold—approximately 100 nm/s—it is no longer the primary factor limiting oral bioavailability. It has been estimated that compounds with permeation rates below 100 nm/s exhibit an average oral bioavailability of 13%, whereas those with rates above 100 nm/s have an average of 35% [236]. In vitro apparent permeability (Papp) measurements done using Caco-2 or MDCK cells (discussed below) cannot be directly translated to in situ effective permeability (Peff) values due to inherent methodological and physiological differences (but mathematical models can use those values in predicting relatively well the human absorption properties) [237].

3.2.9. Predicted Apparent Madin-Darby Canine Kidney (MDCK) Cell Permeability (QPMDCK)

Similar to the BIP Caco-2, MDCK cell permeability serves as a valuable tool for predicting oral drug absorption. This is attributed to the elevated expression of transporter proteins and the scarcity of metabolizing enzymes in these cells, making them an effective model for assessing drug permeability [238]. QPMDCK is measured in nm s⁻¹ and about 95% of pharmaceutical compounds exhibit QPMDCK values between less than 25 (poor) to over 500 (great) [196,239].

3.2.10. Index of Cohesion Interaction in Solids (Indcoh)

Calculated from the number of hydrogen bond acceptors (HBA), donors (HBD), and the surface area accessible to the solvent (SASA, Smol) by the relation Indcoh = HBA × HBD^1/2/Smol. Approximately 95% of pharmaceutical compounds exhibit Indcoh values between 0.00 to 0.05 [196,240]. The descriptor in question appears to have gotten little attention in the current body of scientific literature. There have been few in-depth debates or studies of its significance, applicability, or implications in a variety of fields. This scarcity indicates a potential research gap, which may allow for additional analysis and exploration to better understand the descriptor’s utility and relevance.

3.2.11. Globularity Descriptor (Glob)

Calculated as (4πr²) / S mol, where r is the radius of the sphere whose volume is equal to the molecular volume. Approximately 95% of pharmaceutical compounds exhibit Glob values between 0.75 and 0.95 [241]. The 4πr² term represents the surface area of a perfect sphere with the same volume as the molecule. The denominator represents the actual surface area of the molecule, considering its real shape and any irregularities. This descriptor quantifies the molecular shape’s conformational compactness relative to an ideal spherical configuration, i.e., is a measure of the molecule bulkiness [242]. A value close to 1 indicates a spherical and compact molecule, while lower values suggest a more elongated, irregular, or flexible shape. Flexible molecules typically have lower globularity values as they can easily change their shape and become more elongated [243].

3.2.12. Predicted Polarizability (QPpolrz)

Polarizability quantifies how easily the cloud of electrons around an atom or molecule can be distorted or displaced from its equilibrium configuration when subjected to an external electric field. For electric fields significantly weaker than the intrinsic electrostatic interactions within the atom or molecule, the energy associated with this distortion increases with the square of the field’s strength. While polarizability can be fully described using a tensor (more precisely, a 3D tensor), in cases where the charge distribution around the molecule is spherically symmetric (or approximated so), a single number (a scalar value) can effectively represent polarizability. In many instances, using an average polarizability is sufficient for calculations [244]. In pharmaceutical research, polarizability has been used as a key descriptor for a variety of QSPR and QSAR models, inter alia showing a moderate relationship with the n-octanol/water partition coefficient logarithm (log P, explaining approximately 21% of the variance), and a more substantial connection to aqueous solubility, (accounting for roughly 44% of the observed variability) [244]. Other aspects of pharmacokinetic relevance seem also to be related to polarizability, e.g., the van der Waals interactions between different drugs or metabolites and body albumin or lipids, are favored by the polarizability of those drugs or molecules [245]. The way chemicals bind with body fluids or cells, their movement, and other interactions are influenced by the polarizability of both the chemical and the system. Despite its importance, most QSAR studies up to date have not utilized this descriptor [246]. Whereas it seems to have great potential in offering insight in the pharmacokinetics of a molecule, QPpolrz is not a PK parameter as such, but rather an important descriptor that can be used in predicting various PK properties. Approximately 95% of pharmaceutical compounds exhibit polarizability values between 13.0 to 70.0 [196].

3.2.13. Predicted Skin Permeability (log Kp)

A benchmarking study comparing experimental Kp values for 30 compounds found that QikProp predictions achieved a correlation coefficient of 0.612 with the experimental data. In contrast, the simpler Potts and Guy model, which uses linear regression based solely on logP and molecular weight (MW), had a lower correlation coefficient of 0.544 [210]. While QikProp demonstrated better predictive accuracy compared to the Potts and Guy model, a correlation coefficient of 0.61 is still relatively modest and it indicates that its predictions are not highly accurate. Approximately 95% of pharmaceutical compounds exhibit log Kp values between -8.0 and -1.0, this being therefore considered “the recommended range” [196,247].

3.2.14. Number of Likely Metabolic Reactions (#metab)

The number of metabolic steps a molecule undergoes before reaching its target site influences its ability to access that site after entering the bloodstream. Fewer metabolic steps generally mean easier access to the target site [248]. Molecules that experience numerous metabolic processes may be destroyed or altered prior to arriving at the target site, hence diminishing their likelihood of efficacy. Drug developers are interested in identifying or creating compounds that exhibit metabolic stability and can effectively reach their target locations. Most pharmaceutical compounds (95%) undergo between 1–8 metabolic reactions [196].

3.3. pkCSM

Similarly to SwissADME, this is a freely available web server, using machine learning algorithms to model the relationships between chemical structure and different ADME values. It incorporates two distinct types of predictive models:

• Quantitative models: 14 regression models are available to forecast numerical values for various pharmacokinetic and toxicity properties.
• Qualitative models: 16 classification models are available to predict the likelihood of a specific outcome falling into one of two categories.

When published in 2015, the authors claimed equality or superiority over other available models published in the literature [193].

The 30 models available in pkCSM are grouped under several headings:

a)

Absorption: numeric predictions are available for water solubility, Caco2 permeability, HIA (%), skin permeability (log Kp), and binary categorical predictions are made for P-glycoprotein substrates (yes/no), P-glycoprotein I inhibitors and I P-glycoprotein II inhibitors (inhibitors of P-glycoprotein and P-glycoprotein transport, respectively).

b)

Distribution: numeric predictions are available for VDss (human), fraction unbound (human), BBB permeability and CNS permeability. BBB permeability predictions are derived from in vivo measurements conducted in animal models. In contrast, CNS permeability predictions are based on data obtained through direct carotid artery injection, eliminating the influence of systemic distribution effects [193].

c)

Metabolism: only binary categorical predictions are available under this heading. The application predicts whether a substance is a substrate for the CYP2D6 or CYP3A4 fractions, and whether it is an inhibitor of one or all of the CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4 fractions.

d)

Excretion: total clearance is predicted numerically, whereas the property of a substance to be a renal OCT2 substrate is predicted categorically (yes/no).

Toxicity – a number of ten toxicity endpoints are predicted under this heading, but those are not really ADME features, but safety features, and we will not enter into details here. There is also an option to make all available predictions in a single click.

3.4. ProTox-II

ProTox-II was initially launched as a web-based platform for predicting the toxicity of chemical substances [249], and it has been used in at least 45 scholarly articles focusing on the ADME evaluation of natural compounds. However, in the meantime its successor Protox-3 has been developed and made available [195]. Protox-3 makes a relatively large number of predictions relevant from the toxicity standpoint, but with respect to ADME, it only features a number of six models predicting whether the compund is metabolized by the CYP fractions CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4, and CYP2E1. The prediction is based on classification models (active/inactive), and the predicted class is accompanied by the probability of that prediction, thus offering some additional quantiative information on the extent of the trust one should have in that prediction. The paper describing the platform [249] has not reported on the calibration validation of those predictors, therefore it is not clear to what extent the probability output of those models are well aligned with true class probability (in other words, the authors have not reported on how accurate are those probabilities).

3.5. admetSAR

Similarly to ProTox, admetSAR has also recently reached a third version [197]. In our experience, the access to the website and the speed of the calculations was slower than for the other web-servers mentioned above (but the experience might be different for users from other countries). It also groups its different predictions under several headings:

a)

Physicochemical properties - molecular weight, nAtom (number of atoms), nHet (number of heteroatoms), nRing (number of rings), nRot (number of rotable bounds), HBA (hydrogen bound acceptors), HBD (hydroen bound donors), TPSA, SlogP, „application domain” (applicability domain, i.e., the range of chemical space or feature space where a predictive model is reliable and valid, in other words, a measure about how trustful should be one in the validity of the predictions performed).

b)

Absorption – logS, logP, pka, acidic pKa, basic pKa, Caco-2 (two predictions, one for the permeability value, the other to classify a substance as of high or low permeability), HIA, MDCK, F50%, F30%, F20% (oral bioavailability, a compound will be classified as belonging to one of the three and not belonging to the other two).

c)

Distribution – BBB; inhibition of OATP1B1, OATP1B3, OATP2B1, OCT1, OCT2, BCRP, BSEP, MATE1, and Pgp; property of being a Pgp substrate; plasma protein binding ratio and VDss (only the last two are numeric values, all others are based on categorical models).

d)

Metabolism – inhibition of different CYP fractions (CYP1A2, CYP3A4, CYP2B6, CYP2C9, CYP2C19, CYP2D6, CYP1A2, CYP3A4, CYP2B6, CYP2C9, CYP2C19, and CYP2D6), human liver microsomal stability (HLM), rat liver microsomal stability (RLM), and being a UGT substrate (all are binary predictions).

e)

Excretion – plasma clearance (CLp), renal clearance (CLr), half-life (t1/2), and mean retention time (MRT).

The application leverages the CLMGraph framework, employing a deep-learning approach that combines regression and classification for predictions. The average AUC (Area Under the Curve) across 90 classification endpoints, as reported by the model developers, is 0.870, indicating high predictive accuracy for most endpoints. However, a few endpoints, notably renal clearance, demonstrate lower predictive performance. Over 82% of regression endpoints exhibit moderate or strong predictive power, with Pearson correlation coefficients exceeding 0.70. Certain parameters, though, such as half-life and mean residence time, are affected by a limited predictive accuracy [197].

3.6. ADMETlab

First launched in 2018, ADMETlab has reached its third version in 2024 [198]. It has been developed based on 400,000 data points of high quality and makes predictions for 119 endpoints (an increase of 31 additional endpoints as compared with the previous version). ADMTlab is based on what is in essence a QSAR model built with multi-task deep message passing neural networks (DMPNN) based on a variety of chemical descriptors.

The ADMET 3.0 platform encompasses a comprehensive set of 119 endpoints, categorized as follows:

a)

Physicochemical Properties (21 endpoints): molecular weight, van der Waals volume, density, nHA (number of hydrogen acceptors), nHD (number of hydrogen donors), nRot (number of rotatable bonds), nRing (number of rings), MaxRing (number of atoms in the largest ring), nHet (number of heteroatoms), fChar (formal charge), nRig (number of rigid bonds), Flexibility (ratio of rotatable bonds and rigid bonds), Stereo Centers, TPSA, logS, logD7.4, logP, melting point, boiling point, pKa (acidic), pKa (basic).

b)

Medical Chemistry Properties (20 endpoints): QED (Quantitative Estimate of Drug-likeness, a metric to assess how likely a compound is to have the necessary features to become a drug), SAscore (synthetic accessibility score, estimating the easiness or difficulty of synthesizing the chemical compound), GASA (Graph Attention-based assessment of Synthetic Accessibility), Fsp³ (the ratio between the count of sp³ hybridized carbon atoms and all atom carbons in the molecule), MCE-18 (Medicinal Chemistry Evolution, a descriptor evaluating the novelty of organic molecules through their tetrahedral carbon atoms), NPscore (Natural product likeness score), Lipinski Rule, Pfizer Rule, GSK Rule, GoldenTriangle, PAINS, Alarm_NMR Rule (used to predict thiol reactive compounds), BMS Rule, Chelating Rule, Colloidal aggregators, FLuc inhibitors (predict inhibitors of firefly luciferase), Blue fluorescence, Green fluorescence, Reactive compounds, Promiscuous compounds.

c)

Absorption Properties (9 endpoints): Caco-2 Permeability, MDCK Permeability, PAMPA, Pgp inhibitor, Pgp substrate, HIA, F20%, F30%, and F50%.

d)

Distribution Properties (9 endpoints): PPB (Plasma Protein Binding), VDss, BBB, Fu (fraction unbound in the plasm), OATP1B1 inhibitor, OATP1B3 inhibitor, BCRP inhibitor, MRP1 inhibitor, BSEP inhibitor.

e)

Metabolism Properties (14 endpoints): the property of being an inhibitor or substrate of CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4, CYP2B6, CYP2C8, as well as HLM stability.

f)

Excretion Properties (2 endpoints): CLplasma (plasma clearance) and T1/2.

36 toxicity properties and 8 toxicophore rules are additionally predicted for safety evaluation properties.

A synthetic comparison of the four main free ADME web-server applications is shown in Table 2.

4. Conclusions

The application of in silico ADME-PK methods in pharmaceutical research became widespread following Lipinski’s rule of 5, though it has historically been viewed as an auxiliary component of disciplines like drug metabolism, pharmacokinetics, and computational chemistry. Most practitioners in this field come from backgrounds in computer science, chemistry, or drug metabolism and PK rather than having specialized education in in silico ADME-PK techniques. While this lack of formal training might seem limiting, it has actually fostered a collaborative environment where experts from diverse fields - including in vitro research, statistics, analytical chemistry, pharmacokinetics, structural biology, medicinal chemistry, and machine learning - can work together to drive innovation [250].

Initially met with some hesitation, the application of in silico ADME models and tools to natural compounds, particularly those of herbal origin, has witnessed a significant surge in recent years. While this field is relatively young, having emerged within the last three decades, its growth has been remarkable. The number of research publications utilizing in silico approaches for assessing natural compounds has surpassed 3000 in our estimation, with a substantial majority appearing within the last five years. This rapid growth is evident in the publication trends:

• Pre-2020: Fewer than 100 publications containing the phrase “in silico” in the title or abstract.
• 2021: Approximately 200 such publications.
• 2023-2024: Over 270 such publications annually.

Similarly, if only seven papers indexed by PubMed seems to have included “natural” and “docking” in their title or abstract before 2010, in 2024 alone the number of this kind of papers has exceeded 300. This exponential growth underscores the increasing recognition and adoption of in silico methods in natural product research.

As seen in this review, there are multiple ways of approaching ADME predictions using computational methods, each with their own challenges, costs and limitations. Quantum mechanics methods have been less used, but it is likely that their applications will also increase more and more, as already evidenced by the trend in numbers in the last ten years. The usefulness of molecular docking in speeding drug development remains a topic of debate, despite significant refinement and

improvement of this method over the years. It is widely known that the predictive accuracy of docking often falls short in complex biological systems and it is expected to improve the currently available docking tools to make them more reliable. Even so, molecular docking continues and will probably continue to play a valuable role, particularly in the field of natural products, where resources for wet lab experiments are scarcer. One would expect that pharmacophore modeling will also improve moving beyond the conventional approach based on key interaction points to more sophisticated approaches, for instance integrating quantum mechanics calculations or and machine learning tools. This will allow a better understanding of pharmacophoric features particular to natural compounds and better predictions not only of pharmacodynamic interactions, but also of their ADME profiles.

The increasing computational power and refinement of force fields is likely to allow more accurate predictions of protein-ligand interactions, membrane permeability, and metabolic transformations for natural compounds in the future. The extension of these techniques and improvement with the help of machine learning approaches, will probably contribute even more to a better understanding of protein interactions with natural compounds. They might even be successful where traditional QSAR methods may fall short due to the structural complexity, novelty and diversity. However, QSAR methods, particularly in implementations based on deep learning approaches seem to be the de facto workhorse for ADME predictions, particularly in the case of the ADME(T) prediction servers. An impressive number of such servers have become available for free in the last decade and they seem to be particularly valued by the natural compound researchers, who tend to use them more and more in their research.

Finally, the continued advancement of -omics technologies, coupled with the generation of additional wet-lab experimental data, will allow a refinement and improvement of predictions of natural compound ADME properties using the PBPK modelling. Most likely all of these methods will continue to evolve and improve, and they will allow a better understanding of the ADME profiles and druggability not only for natural compounds, but also for various chemical modifications of such compounds, intended specifically to make them more appropriate for drug development purposes. The future will most likely be brighter, but for the computational chemist, the drug developer, as well as for the phytochemist, these are exciting times to be alive.