Abstract

Aldose reductase (ALR2), a vital enzyme involved in polyol pathway, has befitted as a novel drug target in antidiabetes drug discovery process. In the present study, the binding mode and pharmacokinetic properties of potential polyphenolic compounds with reported aldose reductase inhibitory activity from the genus Scrophularia have been investigated. The human ALR2 enzyme (PDB ID: 2FZD) acted as the receptor in the current study. Among the compounds investigated, acacetin, a methoxy flavonoid, displayed the stable binding to the active site of ALR2 with least binding energy value. Molecular interaction analysis revealed that acacetin interrupts the proton donation mechanism, necessary for the catalytic activity of ALR2, by forming H-bond with Tyr48 (proton donor). In addition, acacetin also possessed favorable ADME properties and complies with Lipinski’s rule of 5 representing the possible drug-like nature compared to other polyphenols. Interestingly, the biological activity predictions also ranked acacetin with higher probability score for aldose reductase inhibition activity. Moreover, the molecular dynamics simulation of ALR2-acacetin complex was validated for the stability of ligand binding and the refined complex was used for generation of receptor-ligand pharmacophore model. Thus, the molecular insights of receptor-ligand interactions gained from the present study can be utilized for the development of novel aldose reductase inhibitors from Scrophularia.

1. Introduction

Modern lifestyle with poor dieting habits and lack of exercise impose a greater risk of diabetes among millions of people under various age groups around the world. In order to meet the growing number of diabetic complications, the search for new antidiabetic drugs with lesser side effects and increased efficacy has been a requisite of all time. Though the etiology of diabetes depends on various mechanisms, several studies state that the increase in aldose reductase enzyme (ALR2) [EC 1.1.1.21] flux in the polyol pathway acts as the important one [1].

Aldose reductase, a principal enzyme from aldo-keto reductase family, catalyzes the NADPH dependent reduction of excess D-glucose to D-sorbitol in polyol pathway. The sorbitol formed from the pathway has poor membrane penetration and its metabolism leads to various diabetic complications such as neuropathy, nephropathy, retinopathy, and cataract formation [2, 3]. In vitro studies in rat, bovine, and human cataract lenses and in vivo investigations in diabetic mice demonstrated the importance of ALR2 [4, 5]. Therefore, the compounds possessing aldose reductase inhibition property have gained importance in diabetic treatments. Besides the importance of aldose reductase inhibitors (ARIs), due to the inefficient pharmacokinetic properties of several synthetic ARIs such as zopolrestat, ponalrestat, sorbinil, and tolrestat have failed to clear the clinical trials [6]. Thus, there is an urgent need for developing and evaluating new ARI with enhanced pharmacokinetic properties, selectivity, and safety. On the other hand, the growing interest in herbalism and naturceuticals has paved the way for searching potential compounds that lack side effects and toxicity. Plant-derived secondary metabolites, especially polyphenols, exhibited significant anti-ALR2, antioxidant, anti-inflammatory, antiviral, anticancer, and platelet aggregation inhibitory activity [79]. Hence, considering the wide existence of ARIs in plants with incomparable structural diversity, we have explored the genus Scrophularia to identify ARI(s) with higher efficacy. Scrophularia belongs to the Scrophulariaceae family consisting of more than 300 medicinal herbs employed to treat diseases such as diabetes, cancer, cardiac disorders, fever, pain, and other ailments [10]. Polyphenols such as gallic acid, vanillic acid, syringic acid, ferulic acid, and protocatechuic acid, cinnamic acid, previously isolated from Scrophularia frutescens and S. sambucifolia, have been reported to inhibit aldose reductase enzyme activity [11]. The extracts of S. takesimensis, S. kakudensis, S. boreali-koreana, and S. buergeriana displayed the aldose reductase inhibitory activity in rat lens [12]. In order to exploit the enormous amount of the herbal resource from this medicinally important genus for antidiabetic treatment the cutting-edge computational approaches such as molecular docking, in silico ADME profiling (absorption, distribution, metabolism, and excretion), and molecular dynamics can be utilized.

In recent days, the identification of novel lead molecules by modern computational approaches based on the three-dimensional structure of the therapeutic target has been increasing. Particularly, molecular docking simulations employed in the computer-aided drug discovery aid in the prediction of promising small molecules based on the binding structure of the ligand to its receptor [13]. In most cases, the docking protocols predict the optimal conformation of the ligands in the biding pocket of the receptor. This approach can be used to model the protein small molecule interaction at atomic level, which renders information about the behavior of small molecules in the binding site of the target protein [14]. One of the primary goals in drug designing process is the identification of small molecules exhibiting high binding affinity towards the target along with reasonable pharmacokinetic properties such as ADME [15]. This approach has resulted in the early determination of drug-likeliness of the lead molecules during the drug discovery process. According to Smith et al. [16], the utilization of in silico and in vitro methods for the investigation of ADME profile of new lead molecules is extensively employed in drug discovery process along with molecular modeling approaches. Furthermore, the evaluation of ligand orientation and the stability of receptor-ligand complex are befitted as an important factor in modern computer-aided drug designing process that can be validated by molecular dynamics simulation [17]. Moreover, the refined structure resulting from the molecular dynamics simulation can be utilized for the generation of pharmacophore model of the complex. The generated pharmacophore can be used for the identification of perquisite chemical features required for the optimal binding of the ligand to the receptor with high efficacy [18]. Thus, owing to the immense need of efficient ARI(s), for the first time this study has attempted to explore the genus Scrophularia for ARIs and investigated the binding mechanism of the polyphenolic compounds and modeled a receptor-ligand pharmacophore using computational chemistry approaches.

2. Materials and Methods

2.1. Molecular Docking Simulation
2.1.1. Small Molecule Preparation

A set of pharmacologically important polyphenolic compounds present in the genus Scrophularia were collected from the literature survey. The two-dimensional structures of the compounds were sketched using ArgusLab V. 4.0.1 molecular graphics tool (http://www.arguslab.com/arguslab.com/ArgusLab.html) with the structural coordinates retrieved from ChemSpider database (http://www.chemspider.com/) and the ligand geometry was optimized using Universal Force Field (UFF) according to Rashmi et al. [19]. The optimized geometries were converted to three-dimensional PDB format and employed for docking analysis.

2.1.2. Macromolecule Preparation

Three-dimensional (3D) atomic coordinates of human aldose reductase enzyme conjugated with tolrestat and NADP+ (PDB ID: 2FZD) determined by X-ray crystallography with 1.08 Å resolution were retrieved from the protein data bank (PDB) (http://www.rcsb.org/pdb). The respective ligands and water molecules were removed and subsequent structure optimization was performed using CHARMM force field by performing 500 steps of steepest descent energy minimization and employed as the receptor for the docking analysis.

2.1.3. Molecular Docking Protocol

Docking simulations were carried out using AutoDock 4.2 [20]. AutoDock calculates the grid map based energy of the rigid protein-flexible ligand complex. For computing the grid maps, protein coordinate files were prepared by adding polar hydrogens, partial charges, and solvation parameters. Atomic affinity grid maps were calculated for all atom types in the protein and the ligands. Grid maps were centered on the catalytic region with 60 Å × 60 Å × 60 Å dimensions and 0.375 Å spacing. All the dockings were initiated with a random population of 50 binding orientations and completed after 1.5 × 10−6 energy evaluations. Point mutation and crossover rates were set at 2% and 80%, respectively. Local search probability performing on an individual was set at 6%. Due to the nondeterministic nature of genetic search algorithms, 50 trials of each docking were performed. Results with clusters varying by less than 1.5 Å in positional root mean square deviation (RMSD) were generated and sorted based on the most favorable free energy for binding. All the docked conformations were visualized using Python Molecular Viewer 1.4.3 (PMV).

2.1.4. Validation of Docking Protocol

The molecular docking protocol was validated by redocking approach. RMSD value was calculated between the atomic coordinates of crystal structure (2FZD) retrieved from PDB and the AutoDock generated docked pose of ALR-2 and tolrestat according to Vianna and de Azevedo Jr. [21]. Both crystal structure and the docked complex were superimposed for RMSD calculation and the observations were performed by using Python Molecular Viewer 1.4.3 (PMV).

2.2. ADME Analysis

The ADME properties of the ligands were predicted using QikProp [22]. It predicts the physiochemical descriptors and pharmaceutically important properties of ligands by comparing the query molecule with those of 95% of known available drugs. Prior to the prediction, the entire known and searched compounds were neutralized and the program was run on the normal mode. The physiochemical properties such as percentage human oral absorption (QP %), QPlogHERG (Predicted IC50 value for blockage of HERG K+ channels), QPPCaco (Predicted Caco-2 cell permeability in nm/sec.), and QPPMDCK (Predicted apparent MDCK cell permeability in nm/sec.). In addition, Lipinski’s rule of 5, an important factor for rational drug design, was also calculated for the small molecules.

2.3. Prediction of Biological Activity Using PASS

The PASS online tool evaluates the general biological potential of drug-like molecules based on the structure of the compound. The statistical probability score for aldose reductase inhibition property of all the polyphenols included in this study was predicted using PASS (http://www.pharmaexpert.ru/passonline/index.php) according to Tripathi et al. [15].

2.4. Molecular Dynamics Simulation

The molecular dynamics validation of ALR2 apo-form (native structure) and ALR2-acacetin complex was performed using Gromacs (Ver 4.5.4) software [23] by applying leap-frog integration steps to solve the equations of motion. The Amber03 force field was employed to simulate the structures. Appropriate number of sodium (Na+) and chloride (Cl) counterions were added to neutralize the simulation system and steepest decent algorithm was employed for the energy minimization with 1000 steps. All the bond constraints were described using LINCS algorithm [24] and the system was equilibrated at 250 K for 100 ps with the solute atoms fixed. After equilibration of the systems, 25 ns production MD simulation with a time step of 2 ps at constant temperature (300 K), pressure, (1 atm) and number of particles, without any position restraints, was performed. The representative structures were selected and analyzed using the standard Gromacs package tools.

2.5. Receptor-Ligand Pharmacophore Modeling

The refined ALR2-acacetin complex was employed for the generation of pharmacophore model using Receptor-Ligand Pharmacophore Generation Protocol in Discovery Studio v3.0 (DS), Accelrys, San Diego, USA, with default parameters. Based on the receptor-ligand interactions, pharmacophore models were generated by the program using predefined chemical descriptors, such as hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic center (HY), ring aromatic (RA), negative ionizable (NI), and positive ionizable (PI). After model generation, the top ranked pharmacophores were returned based on the measure of selectivity score. The score was estimated by the Genetic Function Approximation (GFA) model with a training set consisting of 3285 pharmacophore models searched against the CapDiverse database in DS [18].

3. Results and Discussion

3.1. Molecular Docking Analysis of ARIs

Docking simulations of the secondary metabolites with anti-ALR2 activity in Scrophularia were performed to gain insight into the binding mode of the screened small molecules to the crystallographically determined active site of ALR2 crystal structure [25, 26]. In order to determine the accuracy of the docking protocol, prior to docking simulation of the polyphenolic compounds the docking protocol was validated by redocking method. The results of redocking indicated the accuracy of the docking protocol with the RMSD value of 0.53 Å between the crystal structure and the docked pose of ALR2 with tolrestat (Figure 1). Since the best docking simulation is the outcome of the binding pose closer to the crystal structure, the lesser RMSD value obtained during the validation denotes the successful docking protocol that can be implemented for further binding analysis of the ligands [21].

Docking results with the free energy of binding, number of hydrogen bonds formed, residues participating in H-bond formation, and electrostatic interaction along with predicted and experimental IC50 values are listed in Table 1. All the 50 docked conformations of the inhibitors were clustered according to the RMSD criterion of 1.5 Å. To enhance the accuracy, the docked poses were visually screened according to Cosconati et al. [27]. However, it has to be concerned that the scoring functions used in docking algorithms provide only the approximate estimate of binding energy; therefore the results must be further validated with other physiochemical properties and biological activities [28]. Consequently, the free energy of binding resulting from the AutoDock calculations has been correlated with the reported experimental IC50 values. The docking energy displayed strong positive correlation with the IC50 values () illustrating the reliability of the receptor-ligand interaction. AutoDock predicts the IC50 value solely based on the binding energy of the complex [29]. However, the Pearson correlation analysis revealed the existence of strong positive association () between the predicted IC50 values generated by the AutoDock and experimental IC50. Overall the aforementioned results demonstrate the fidelity of the docking approach employed to investigate the binding mode of ARIs in the current endeavor.

3.1.1. Interaction Analysis of Acacetin with ALR2

Among all the compounds, acacetin (−10.08 KJ·mol−1) displayed promising binding with the least reported experimental IC50 value (4.76 μM). Acacetin (5,7-dihydroxy-4′-methoxyflavone) is an O-methylated flavonoid with anticancer, antioxidant, anti-inflammatory, and antialdose reductase property [30]. This potent hit with stable binding mostly exhibited polar interaction with the hydrophobic residues located in a deep cleft by forming four hydrogen bonds in the C-terminal end of the β barrel (Figure 2). In detail, it anchored to the active site of ALR2 by forming hydrogen bonds with 5- and 7-hydroxyl oxygen of the benzopyran ring with Tyr48, Lys77, and Cys298, respectively. The activity of ALR2 depends on the Tyr48 residue (proton donor) during the catalytic mechanism and it acts as a general acid catalyst to assist the hydride transfer from NADPH to the substrate’s carbonyl carbon [31]. Since, acacetin forms H-bond with Tyr48 it can possibly disrupt the proton donation of Tyr48. Our results were also in accordance with Rastelli et al. [32], which suggest that the hydrogen bond formation with Tyr48 is the minimal requisite for any compound to consider as an ARI. Acacetin interacts not only with the proton donor but also with Lys77 and Asp43 residues that are responsible for the salt bridge formation necessary for the proton transfer (Table 1). Subsequently, the hydrogen bond of acacetin with Cys298 could destabilize the NADP+ binding, since Cys298 is the key residue that contributes to the NADP+ stabilization [32, 33]. One more hydrogen bond was observed between 4-hydroxyl hydrogen of the pyran ring with Ile26. This interaction might stabilize the docked complex and it has been noted that the involvement of Ile26 is specific for acacetin alone while being compared with all other ligands used in this study. Further, 4′-methoxy ring of acacetin attached to the benzopyran ring with single free rotational bond stacked into the hydrophilic outer region covering the active site and electrostatic interaction has been observed with Ser210, Pro211, Leu212, Ser214, and Pro215 residues. Consequently, the aromatic substituent is inserted into the hydrophobic pocket of the active site composed of Gly18, Trp20, Lys21, Ile26, Asp43, and Tyr209. Interestingly, acacetin interacts with Trp20 which is primarily required for the tight binding of any ARI to the active site of ALR2. Additionally, it interacts with Tyr209 that balances the positively charged aromatic system of NADPH to NADP+ upon oxidation during the enzyme activity.

3.1.2. Interaction Analysis of Cinnamic Acid Derivatives

Ferulic acid a notable hydroxycinnamic acid derivative showed good interaction with binding energy of −7.21 KJ·mol−1 with ALR2. It formed three hydrogen bonds with Tyr19, Asn160, and Lys77 residues. The phenol ring of ferulic acid is exactly positioned into the active site of ALR2 and interacts with Gly18, Tyr48, His110, Ser159, Tyr209, Ser210, and Lys262 (Figure 3(a)). Notable σ-π interaction has been observed between the phenol ring of ferulic acid and Trp20. However, due to the lack of hydroxyl groups in the six-carbon ring for H-bond formation the binding of transcinnamic acid appears to be less stable than ferulic acid. Yet, the phenol ring of transcinnamic acid managed to fold into the active site of ALR2 and displayed prominent electrostatic interaction with Tyr48, Trp20, Asp43, His110, Gln183, Ser159, and Tyr209 (Figure 3(b)).

3.1.3. Interaction Analysis of Hydroxybenzoic Acid Derivatives

Protocatechuic acid and gallic acid are the hydroxybenzoic acid derivatives which exhibited similar interaction and H-bond formation with ALR2 because of their structural similarity. In both acids, 3-hydroxyl oxygen and hydrogen interacted with Thr19 and Asp43, respectively. Interestingly, the 4-hydroxyl oxygen of gallic acid formed H-bond with Lys77, while the same functional group in protocatechuic acid formed H-bonds with both Tyr48 and Lys77, respectively (Figures 4(a) and 4(b)).

3.1.4. Interaction Analysis of Methoxybenzoic Acid Derivatives

Similarly, methoxybenzoic acid derivatives, vanillic acid, and syringic acid displayed H-bond interaction with Asp43 and Lys77 (Figures 5(a) and 5(b)). Nearly, all the polyphenols investigated in this study interacted with the major active site residues such as Trp20, Asp43, Tyr48, Lys77, and Tyr209. In addition, ferulic acid and cinnamic acid displayed remarkable interaction with His110 which has also been reported as a vital residue for ALR2 activity. Hence, the docking results of polyphenols illustrate that the phenol ring in all the polyphenols with ionized hydroxyl group and carboxyl groups is inevitable for the stable binding to ALR2. Further, the docking pattern revealed that the polyphenols shared a similar mode of binding into the active site of ALR2.

Based on the molecular docking results, it can be evident that the acacetin displays the stable binding towards ALR2 compared to the other polyphenols employed in this study.

3.2. ADME Analysis

The results of ADME profiles predicted by QuikProp are listed in Table 2. All the ADME descriptor values for the polyphenols were in permissible range except for IC50 value of the HERG K+ with exception for acacetin. Interestingly, the pharmacokinetic analysis revealed that the top docking hit acacetin possesses the most acceptable range of values for all physiochemical descriptors studied. In detail, the higher human oral absorption percentage (QP = 87.784%) was exhibited by acacetin. The QP% values of the other polyphenolic compounds were in the range of 41.486 to 79.474. Similarly, acacetin also displayed the prominent permissible predicted IC50 range for HERG K+ channel (QPlogHERG K+ = −5.201), Caco-2 cell permeability (QPPCaco = 384.043), and MDCK cell permeability (QPPMDCK = 175.838) compared to other compounds. The Caco-2 cells are considered as the model for gut-blood barrier and its permeability governs the drug metabolism and transportation to the membranes; likewise the MDCK cells represent the blood-brain barrier [34].

Along with the ADME properties the investigation of Lipinski’s rule of five has been widely accepted as physicochemical descriptors for predicting the drug-likeliness of a compound. The results of Lipinski’s rule of five descriptors along with the number of rotational bonds for the polyphenols were shown in Table 3. The rule of five proposed by Lipinski is essential for the estimation of drug-likeliness that determines the oral bioactivity of a compound or drug [15]. According to Lipinski’s rule [35], a compound with high drug-likeliness should satisfy the following criteria: the compound should not have more than 5 hydrogen bond donors and 10 H-bond acceptors, the molecular weight should not exceed 500 Daltons, and the octanol-water partition coefficient () value should be less than 5. In our results all the phytochemicals comply with Lipinski’s rule of five. Moreover acacetin comprises the maximum value of 2.999. In general, the value is a critical factor for the determination of drug absorption within the human body [34].

3.3. Prediction of Biological Activity Using PASS

In order to investigate the biological activity of the promising ARIs based on the structural formula PASS analysis was employed. In the present study the compounds were investigated for the aldose reductase inhibitory activity and the results are listed in Table 4. PASS predicts the pharmacological influences, action mechanism, and toxicity exhibited by a compound during its interaction with biological molecules [15]. In PASS the structures of the ligands are described by universal descriptors called multilevel neighborhoods of atoms (MNA) that represent several structure-property relationships, with reference to biological activity, mutagenicity, carcinogenicity, and drug likeness [36]. Thus, based on the MNA descriptor statistics for active and inactive compounds from the original training set, the probability of the compound being active () and inactive () was calculated [37]. Based on the probability value the compounds were classified to be active or inactive for a particular class of action (aldose reductase inhibition in this case). A compound with most probable activity should have the value close to 1 and value close to 0 [15]. Most interestingly, acacetin displayed the highest value (0.341) and the least value (0.003) among the polyphenols. Thus the biological activity prediction also suggested that acacetin could serve as a better ligand molecule with aldose reductase inhibition activity than other phytochemicals under investigation.

Overall, the molecular docking simulations along with the reported experimental IC50 value, ADME analysis, and PASS predictions demonstrated the higher affinity and favorable drug-like property of acacetin towards aldose reductase than other polyphenolic compounds involved in this study. Therefore to validate stability of acacetin with ALR2, molecular dynamics based evaluation has been performed.

3.4. Molecular Dynamics Simulation of Apo-Form and ALR2-Acacetin Complex

The stability and fluctuation of the residues present in the ALR2 apo- and complex structures were validated based on the root mean square deviation (RMSD) and root mean square fluctuation (RMSF), respectively, resulting from 25 ns production run. The ALR2 apo-form and the ALR2 complex reached the equilibrium state around 2 ns, respectively. The RMSD result suggested that the docked complex was stable from 7 ns throughout the entire simulation period (Figure 6). Furthermore the local mobility of the secondary structures analyzed by RMSF with respect to the Cα atoms was shown in Figure 7 along with the conformational changes observed in apo-form and docked complex. In RMSF plot, five distinct peaks were observed and the corresponding structural variation has been analyzed. Peak 1 corresponded to the amino acid residue Asn50 which is proximal to the proton donor Tyr48 of ALR2. The loop region formed by Asn50 in apo-form has been converted to the part of the helix upon acacetin binding. Similarly, Peak 2 represented the formation of short helical turn by Val67 in the docked complex whereas in apo-form it formed a loop region. Furthermore Peaks 3 and 4 denoted the residues Ser127 and Leu175, respectively, in loop regions with slight displacement in the apo- and docked forms. Peak 5 corresponded to Asp224 which forms a helical turn in the native structure and relaxed to loop upon the docking of acacetin. Thus, the molecular dynamics simulation of the docking result illustrates that the binding position of the acacetin did not affect the stability of the ALR2 native structure and also revealed the conformational changes occurring during the acacetin binding.

3.5. Receptor-Ligand Pharmacophore Model

The receptor-ligand pharmacophore modeling helps to determine the important chemical features that instigate the vital interactions between the receptor and ligand along with the determination of excluded volume spheres corresponding to the receptor structure [18]. In this study totally 10 pharmacophore models for ALR2-acacetin complex were generated with the 5 to 4 features set with the selectivity scores ranging from 9.9832 to 7.5549. The pharmacophore model with the highest selectivity score of 9.9832 has been illustrated in Figure 8 with five features such as HBA = 1, HBD = 2, RA = 1, and PI = 1. Interestingly, the pharmacophore model consisted of chemical features that interacted with the vital residues of ALR2 as shown in Figure 9. In detail, the HBA and HBD features oriented towards key residues such as Asp43, Cys44, Ala45, Tyr48, Lys77, His110, Gln183, Gly18, Tyr209, Ser210, Ile260, and Lys262. Subsequently the features RA, HY, and PI were located in such a way that interacts with important amino acids such as Pro211, Leu212, Gly213, Ser214, and Val297. Our results are in concordance with the previous reports suggesting the pharmacophore features required for the ARIs necessary for the active site binding of aldose reductase [31, 32, 38]. Thus, the occurrence of chemical features vital for the ligand interaction with key residues implicate that the pharmacophore model based on the ALR2-acacetin complex may offer an efficient approach to identify novel aldose reductase inhibitors.

4. Conclusion

Overall the current study has attempted to investigate the binding modes of polyphenolic compounds from Scrophularia towards aldose reductase enzyme and to identify a potential ligand molecule with efficient binding using computational chemistry approaches. The molecular docking results revealed the higher binding affinity of acacetin to the active site of ALR2 withfavorable free energy of binding and H-bonds. In addition acacetin forms H-bond with Tyr48 that could aid in the disruption of proton donation necessary for the catalytic mechanism. Furthermore, acacetin also encompassed the optimal ADME properties and higher biological activity than the other phytochemicals. In addition molecular dynamics simulation has been carried out to validate the stability of acacetin binding towards ALR2. Additionally, the refined ALR2-acacetin complex has been used to generate a receptor-ligand pharmacophore model that can be utilized for the identification of novel ARIs.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

Abinaya Manivannan and Prabhakaran Soundararajan were supported by a scholarship from the Brain Korea 21 plus (BK21 Plus) Program, Ministry of Education, Science and Technology, Korea.