Design, synthesis, biological evaluation, and molecular modeling studies of pyrazole-benzofuran hybrids as new α-glucosidase inhibitor

In this work, new derivatives of biphenyl pyrazole-benzofuran hybrids were designed, synthesized and evaluated in vitro through enzymatic assay for inhibitory effect against α-glucosidase activity. Newly identified inhibitors were found to be four to eighteen folds more active with IC50 values in the range of 40.6 ± 0.2–164.3 ± 1.8 µM, as compared to the standard drug acarbose (IC50 = 750.0 ± 10.0 μM). Limited Structure-activity relationship was established. A kinetic binding study indicated that most active compound 8e acted as the competitive inhibitors of α-glucosidase with Ki = 38 μM. Molecular docking has also been performed to find the interaction modes responsible for the desired inhibitory activity. As expected, all pharmacophoric features, used in the design of the hybrid, are involved in the interaction with the active site of the enzyme. In addition, molecular dynamic simulations showed compound 8e oriented vertically into the active site from mouth to the bottom and stabilized the enzyme domains by interacting with the interface of domain A and domain B and the back side of the active site while acarbose formed non-binding interaction with the residue belong to the domain A of the enzyme.

www.nature.com/scientificreports/ On the other hand, benzofuran-2-carbohydrazide 7 was prepared by following the reaction scheme. Initially, a mixture of salicylaldehyde 5 with ethyl bromoacetate was refluxed in acetonitrile in the presence of potassium carbonate to afford ethyl benzofuran-2-carboxylate 6. Then, refluxing the ethanolic solution of the ethyl ester 6 and hydrazine hydrate yielded benzofuran-2-carbohydrazide 7. Finally, desired products were obtained by coupling the hydrazide 7 with key aldehydes 4a-o in absolute ethanol and in the presence of glacial acetic acid. The structures of newly synthesized compounds were confirmed by their IR, 1 H NMR, 13 C NMR, MASS and elemental analysis 29,45 . In vitro α-glucosidase inhibitory activity. All of the newly synthesized compounds 8a-o were screened for α-glucosidase inhibition activity. The α-glucosidase enzyme (from Saccharomyces cerevisiae, EC3.2.1.20) was used to evaluate the α-glucosidase inhibitory activity. Compared to human α-glucosidase, they shared similarities in substrate specificity, pH optimum, catalytic residues in the active site and inhibitor sensitivity 46 . Therefore, α-glucosidase from yeast is extensively used for the preliminary screening of α-glucosidase inhibitors 47,48 . Acarbose, a commercially available α-glucosidase inhibitor, was used as the reference drug and the results are expressed as IC 50 values (Table 1). Albeit, all pharmacophoric groups in the designed hybrid seem to have involved in the inhibitory potential but to further elucidate the role of substituents on the aryl rings connected to pyrazole, a wide variety of compounds 8 were synthesized and structure-activity relationship (SAR) has been established based on varying substituents in R 1 and R 2 . Accordingly, SAR was investigated in two categories of 8a-i (Table 1, R2 = H) and 8j-o (Table 1, R2 = CH 3 ). It is worth mentioning that all synthetic compounds showed significant activity with IC 50 values in the range of 40.6 ± 0.2-164.3 ± 1.8 µM, when compared to acarbose (750.0 ± 10.0 µM). Among them, compound 8e (IC 50 = 40.6 ± 0.2 μM), having a nitro group in the R 1 position, was found to be the most potent inhibitor. This is approximately 18-fold more potent than the standard acarbose. The introduction of other electron-withdrawing groups, such as halogen atoms and trifluoromethyl (8f-i), led to reduced inhibitory activity. Halogenated counterparts of compound 8e in the first category of compounds, exhibited inhibitory activity generally depending on the size of halogen substitution. Compound 8f (IC 50 = 79.7 ± 0.5 μM) bearing bromo substituent was found to be the second most active analog in this category and nine-fold more potent than standard. Compounds 8g having a relatively large chlorine substituent and 8h having a small fluorine group www.nature.com/scientificreports/ showed lower activities with IC 50 values of 125.3 ± 1.0 and 164.3 ± 1.8 µM, respectively. It should be noted that change the fluoro group with trifluoromethyl moiety, as in the case of compound 8i (IC 50 = 101.7 ± 0.7 μM), led to a significant improvement in the inhibition potential. Finally, the order of inhibition for compounds 8a-d was 8c > 8b > 8d > 8a (IC 50s = 100.3 ± 0.7, 108.3 ± 0.8, 133.6 ± 1.2, and 141.2 ± 1.4 μM, respectively), which illustrated that inhibitory activity almost affected by the electron-donating property of substituent. For further investigation the SAR, second category of compounds 8j-o possessing methyl group in R 2 were also examined. Compound 8l (IC 50 = 54.7 ± 0.3 μM) containing methoxy group was found to be the most active compound in this series. Also, this compound is the second most potent among the screened compounds and > 13 folds more active than the standard. Replacement of methoxy group in 8l with methyl (8k, IC 50 = 115.6 ± 0.9 μM) or hydroxyl group (8m, IC 50 = 127.0 ± 1.0 μM) resulted in a remarkable decrease the biological activity. The   www.nature.com/scientificreports/ bromo substituted derivative 8n was found as the second most potent molecule among this series of compounds and exhibited approximately eight-fold enhanced activity compared to the standard acarbose. So, as mentioned in the previous sections, the size of the bromine may be responsible for the higher activity of brominated compounds than other halogens. It can be concluded that the in vitro α-glucosidase inhibitory activity mainly depends upon the substituents on R 1 and the nature of substituents at this position affected the efficacy of the methyl group on R 2 . In this regard, the α-glucosidase inhibitory activity of bromine (8f, 8n) and chlorine (8g, 8o) substituted analogs in two series of mentioned compounds revealed that the introduction of methyl groups on R 2 led to lower activity (8f > 8n and 8g > 8o). Similarly, the lack of methyl group in the case of compound 8k resulted in a slight decreased in inhibitory activity as compound 8b. In a different manner, compound 8l was more active than compound 8c. Both compounds have a methoxy group at the R 1 position, but the presence of the methyl group on R 2 in compound 8l led to a two-fold increase in activity. The same is the case of compounds 8d and 8m having hydroxyl group, the methyl substituent on R 2 position was found to confer an increase in inhibitory activity.
Kinetic study. To gain further insight into the mechanism of action of the synthesized compounds against α-glucosidase, a kinetics analysis was performed on the most potent compound 8e. In different concentrations of test compound (0, 10, 25, and 40 µM) and with the incremental concentration (2-10 mM) of substrate, the rate of the enzyme activity was calculated. The type of inhibition and experimental inhibition constant (K i ) value were determined by employing Lineweaver-Burk plots and secondary re-plot of these plots, respectively. The Lineweaver-Burk plot showed that with varying concentrations of compound 8e, V max of enzyme gradually increased without affecting the K m of enzyme (Fig. 3a). This pattern indicates a competitive type of inhibition. The plot of the slope of lines in the Lineweaver-Burk plots (k m ) against the inhibitor concentration gave an approximation of the inhibition constant, K i of 38 μM for compounds 8e (Fig. 3b). www.nature.com/scientificreports/ Cytotoxic activity. Cytotoxicity of the most potent compounds 8e, 8l, 8f and 8n were evaluated against normal 3T3 cell line by using MTT assay. The results of cytotoxic activity were expressed as the IC 50 (µM) and outlined in Table 2. Results revealed that at 150 µM concentration, the selected compounds were noncytotoxic against studied normal cell line.
Homology modeling and molecular docking study. Molecular docking studies were also performed to rationalize the results of biological assays and gain structural insight into the binding of the synthetic deriva-  www.nature.com/scientificreports/ tives against α-glucosidase. Due to the unavailability of the crystallographic structure of α-glucosidase from S. cerevisiae, the 3D structure of α-glucosidase was modeled using MODELER inbuilt in Discovery Studio (DS) package and synthetic compounds were docked against the established homology model. For this purpose, the FASTA format of the primary sequence was downloaded from Uniprot (Access code P53341) and submitted to NCBI BLAST to get a template with a suitable identity for sequence alignment 49 . Isomaltase from S. cerevisiae (PDB ID: 3A4A) with 71.4% identity and 86.7% similarity was selected as the template for modeling 50 . The best model was selected based on the lowest PDF Total Energy (3270.6404) and DOPE Scores (− 73110.257813) and evaluated for further validation. The PROCHECK program was applied to assess the stereochemical quality of the model. The phi/psi Ramachandran plot distributions indicated that 99.6% residues are in the favored and allowed regions and only 0.2% residues lie in the outlier region (Fig. 4). The superposed structure of acarbose (standard inhibitor) and the predicted top-scored conformation of the most potent compound 8e in the active site of a homology model of α-glucosidase was shown in Fig. 5 (left). The detailed binding mode of acarbose showed that it formed hydrogen bonding interactions with residues Asp349, His239, Asp68, Pro309, Glu304, Arg439, Arg212, Glu276 and one hydrophobic interaction with Phe157 ( Fig. 5, right).
The theoretical binding mode of most active compounds 8e, 8l, 8f and 8n were also shown in Fig. 6. Binding mode analysis showed that the following interactions are common among these compounds: (1) Pyrazole moiety interacts with Glu276 through π-anion binding. (2) Phenyl groups attached to the pyrazole moiety, on the other hand, are involved in hydrophobic interactions with Arg439, Ala278 and Tyr71. (3) The carbonyl oxygen of the amide group formed hydrogen bonds with the hydroxyl group of Tyr313, and (4) the planar benzofuran scaffold interacted with Arg312 via π-alkyl interaction. As expected, all pharmacophoric features used in the design of the hybrid, are involved in the interaction with the active site of the enzyme. The most potent compound 8e establishes more interactions with the residues in the binding pocket. The nitro substituent of this compound created a hydrogen bond with Arg439 and also two electrostatic interactions with Asp68 and Tyr71. Besides, the N1-phenyl ring of pyrazole moiety and benzofuran also formed hydrophobic interactions with Leu218 and Phe157, respectively, which leads to a snug fit at the binding site (Fig. 6a). When the interaction mode of compound 8l as the second most potent compound is compared to that of compound 8e, only one hydrogen bond with Asp68 stabilizes diphenyl pyrazole moiety in 8l, while in compound 8e the electrostatic interactions play an important role in the binding of this moiety to the enzyme. Methyl substituent in compound 8l is located in a hydrophobic pocket formed by residues Leu218, His245, His279, Phe300 and Ala278 (Fig. 6b). In the case of compounds 8f and 8n, both have bromine substitution at the R 1 position, the presence of methyl group in R 2 may cause their different orientations in the active site and then the difference in inhibitory activity. The predicted binding mode of compound 8f shows that the NH proton of the amide group forms hydrogen bond interaction with the amide group of Gln350. The hydrogen bond belonging to the carbonyl oxygen of the amide group was not seen in this compound. Moreover, the CH imine group is forming a hydrogen contact with carboxyl oxygen of Asp349 and thus leads to a better fit of this compound in the enzyme's active pocket (Fig. 6c,d). The calculated GOLD Fitness Scores for compounds 8e (70.4927), 8l (65.7091), 8f (64.84.4) and 8n (61.1567) were in good agreement with those results obtained in in vitro assay. www.nature.com/scientificreports/ Molecular dynamic investigation. Molecular dynamic simulation was performed in order to understand the effect of the compound over the enzyme active site. For this purpose, the structural perturbations incurred by the most potent compound (8e) has been investigated through the study of the RMSD, RMSF and its effect on the active site environment in comparison to acarbose as α-glycosidase standard inhibitor and the apoenzyme. Root mean square deviation (RMSD) of the enzymes' backbone was analyzed over 20 ns MD simulation in order to study the stability of the protein-ligand complex trajectories (Fig. 7). The RMSD value of the apo α-glycosidase depicts broad fluctuations throughout simulation time which is higher than the two enzyme www.nature.com/scientificreports/ complexes. The RMSD value increased after about 4 ns and steadily increased up to 16 ns and become more stable for the last 4 ns of simulation time with the value of 2.5 Å. The RMSD value of glycosidase complexed with acarbose was stable until 12 ns and slightly increased through the next 4 ns and become steady for the rest of the simulation time with the RMSD value of 2.1 Å. Although, the mentioned value of α-glycosidase complexed with compound 8e is the same as acarbose bounded state for the first 8 ns. It is observed that compound 8e had higher RMSD than acarbose for its higher number of rings and flexibility, which makes it more deviate from the initial structure for the next 8 ns and finally it decreased and stabilized for the last 2 ns with the same RMSD as acarbose bounded state (2.1 Å). In summary, the RMSD value of the bounded state enzymes deviate from the initial structure of apoenzyme in the early part of the simulation and obviously decreased as a result of α-glycosidase structural rigidity. Thus, the structures at the last 2 ns of the MD equilibrium state used to investigate the structural specificity of the ligand-protein complexes.
The RMSF, which indicates the flexibility of protein structure, refers to the fluctuation of the Cα atom from its average position throughout the simulation time. Figure 8a compares the residue RMSF values of α-glycosidase bound state and unbound state in which, the apoenzyme (yellow color line) had higher RMSF fluctuations compared to the glycosidase bound-states (green and red-colored lines). This observation occurs upon ligand binding to the enzyme, in which residues movement decrease as a result of non-bonding interaction between the ligand and the enzyme. In addition, the structural segments which are affected upon ligand binding have recognized and categorized into four apparent parts, including; B domain loop (residues 139-149), the active site lid, A domain and B domain sides of the active site mouth. Comparing RMSF values shows that the residues of the B domain loop, 139-149, would have significantly lower RMSF value in glycosidase/ acarbose and compound 8e bound-state rather than apoenzyme (Fig. 8a,b). In contrast, the flexibility of the active site lid increased in Figure 9. shows the detailed orientation and ligand atom interactions that occurred more than 30.0% of the simulation time during the equilibrated phase over α-glycosidase complexed with compound 8e (a,b) and acarbose (c,d). Domain A, domain B and the flap region covered the mouth of the active site colored in yellow, blue and pink, respectively. www.nature.com/scientificreports/ enzyme bound-state and is more pronounced through acarbose binding. In order to investigate the reason for the mentioned observation, as noticed in Fig. 8b (glycosidase/ acarbose complex), acarbose interacted with several residues located into the A and B domain side of the active site mouth (the vertical green line). Although compound 8e interacted with the same regions, it formed fewer interactions with the residues of the A domain. So, it can be proposed that the more ligand interaction with A and B domain sides of the active site mouth, the higher the RMSF of the active site lid. Moreover, Fig. 8c,d represent the organization of the α-glycosidase three main domains; A, B, and C and the close-up representation of the active site mouth with the corresponding residues of A and B domain at both sides in which the active side lid and a back-wall helix situated at the front and the back of the mentioned mouth, respectively. Backing to Fig. 8b, compound 8e provides higher interaction with the back-wall of the active site rather than acarbose. Based on the observed result of RMSF plot, although α-glycosidase/ acarbose complex with higher interaction at the entrance region loop covering the active site (310-315) have lower RMSF value than in α-glycosidase-8e complex, the other lid loop covering the active site consists of residues 230-236 51 shows significantly lower RMSF value in α-glycosidase-8e complex rather than α-glycosidase/acarbose complex. Finally, based on RMSF plot, it can be proposed that compound 8e had more interaction over the B domain side and back-wall helix of the active site, while acarbose formed more interaction with the A domain side of the α-glycosidase active site. Figure 9a,b represent the detailed orientation and interactions that occurred more than 30% of the simulation time during the equilibrated phase over α-glycosidase complexed with compound 8e. The interaction analysis depicts compound 8e oriented vertically from the mouth to the bottom of the active site and stabilized the enzyme domains by interacting with Phe311, Tyr313, Arg312 from the A domain side and Phe158, Phe177 and His239 from the B domain side of the active site mouth (Fig. 9a). In addition, polar residues including; Asp214, Asp349 and Arg439 provide polar interactions with compound 8e at the depth part of the active site. In the same way, acarbose disposed vertically and formed non-binding interactions with the Phe311, Asn241, Arg439, Asp68, His245, Asp349, Asp214, which belong to the domain A of the enzyme (Fig. 9c).
Comparing MD simulation of compound 8e and acarbose proposed the long-lasting non-binding interactions with Asp349, Asp214 and Arg439 have a significant role in inhibition activity of the mentioned compounds (Fig. 9b,d). Figure 9a,b represent two important structural moieties in stabilizing compound 8e at the mouth of the active site. The first one is the benzofuran ring, which interacts with Phe311 and Tyr313 through T-shape π-π hydrophobic interactions for 90% and 55% of the simulation time and the next one, pyrazole ring which interacts with Arg312 through π-cation hydrophobic interaction for 96% of the simulation time. Along with the interactions which stabilized compound 8e in front of the A domain side of the active site entrance, the phenyl substituent interacts with Phe157 and His239, which faced at the B domain side of the entrance for about onethird of the simulation time. In addition, the hydrazide moiety can provide H-bond interaction with Asp408 at the back part of the active site (previously known as the back wall side) for about 35% of simulation time. Finally, the nitrophenyl group as a polar moiety interacts with polar catalytic residues Arg439, Asp349 and Asp215 through ion-bridge, H-bond and water-mediated H-bond interactions for a significant amount of simulation time that is similar to the behavior of the NH2 group in the acarbose (Fig. 9d). It is obvious from the MD study that H-bond, hydrophobic interactions and ion-bridge interactions have a critical role in stabilizing compound 8e at different sides of the active site during the simulation time and cause α-glycosidase inhibition activity. This observation may propose the contribution to the higher α-glycosidase inhibition activity.
In addition to the interaction analysis, the Prime/MM-GBSA module was used to estimate the strengths of interactions between the ligand-protein complexes generated by the clustering method. ΔG bind of α-glycosidase/ compound 8e complex and α-glycosidase/acarbose complex was estimated to be − 90.83 and − 62.49 kcal/mol, respectively, revealing stronger binding interaction of compound 8e than acarbose which also supported by experimental assay.

Conclusion
With aim of developing a novel class of α-glucosidase inhibitors, new series of biphenyl pyrazole-benzofuran hybrids derivatives were designed, synthesized and evaluated for their α-glucosidase inhibition. All screened compounds displayed multifold enhanced inhibitory strength in the range of 40.6 ± 0.2-164.3 ± 1.8 µM when compared to acarbose (IC 50 = 750.0 ± 10.0 µM). Among them, compound 8e, having a nitro group in R 1 position, was found to be the most potent inhibitor. This is approximately 18-fold more potent than the standard acarbose. Also, the kinetic analysis revealed that compound 8e compete with the substrate for binding to the binding site of the enzyme. Limited SAR studies indicated that the in vitro α-glucosidase inhibitory activity mainly depends upon the substituents on R 1 and the nature of substituents at this position affected the efficacy of the methyl group on R 2 . Binding mode analysis showed that almost all structural features such as pyrazole ring, Phenyl groups attached to the pyrazole moiety, amide linkage and benzofuran scaffold are contributing to binding affinity through hydrogen bonding, hydrophobic and electrostatic interactions. In addition, MD simulations showed compound 8e oriented vertically into the active site from mouth to the bottom and stabilized the enzyme domains by interacting with the interface of domain A and domain B and the back side of the active site, while acarbose formed non-binding interaction with the residue belong to the domain A of the enzyme. Moreover, carbonyl hydrazide linker, pyrazole and its related substituents provide such a strategic point with the ability to interact with various parts of the active site, which has the binding and catalytic role for α-glycosidase activity.
Taken together, the above results suggest that newly synthesized hybrids could be promising hits for the further development of α-glucosidase inhibitors for the treatment of diabetes patients. General procedure for the preparation of 1,3-disubstituted-4-pyrazole carbaldehydes 4a-n. To a solution of 4-substituted phenylhydrazine hydrochloride 2 (20 mmol) in ethanol (15 mL), substituted acetophenone 1 (20 mmol) and catalytic amounts of sulfuric acid was added and then the mixture was refluxed for 8-12 h. After reaction completion, the mixture was cooled to room temperature and poured on crushed ice to afford hydrazone intermediate 3a-n. The resulting solid was filtered and recrystallized from ethanol. A solution of hydrazone 3a-n (20 mmol) in DMF (5 mL) was added drop-wise to an ice-cold solution of DMF (15 mL) and phosphorus oxychloride (60 mmol) and the resulting mixture refluxed at 60-70 °C for 5-8 h. After the completion of the reaction, the reaction mixture allowed to cool, poured into ice-cold water and then neutralized with saturated aqueous sodium hydroxide solution. Further, the solid precipitated was filtered, washed with excess cold water and recrystallized from ethanol to afford aldehydes 4a-n.

General procedure for the preparation of ethyl benzofuran-2-carboxylate 6.
A mixture of salicylaldehyde (20 mmol), ethyl bromoacetate (20 mmol) and K 2 CO 3 (40 mmol) in acetonitrile (10 mL) heated under reflux for 4 h. After completion of the reaction (monitored by TLC), the reaction mixture was allowed to cool to room temperature and poured into crushed ice. After extracting the product with ethyl acetate (2 × 25 mL), the organic layer was washed using brine solution (2 × 20 mL) and dried over anhydrous MgSO 4 . The solvent was evaporated under vacuum to afford the product as an oil.
General procedure for the preparation of benzofuran-2-carbohydrazide 7. benzofuran-2-carboxylate 6 (20 mmol) and hydrazine hydrate (30 mmol) in EtOH (10 mL) heat under reflux overnight. Upon cooling, the product precipitated was filtered, washed with cold water and recrystallized in methanol to afford the pure product.
General procedure for the preparation of pyrazole-benzofuran hybrids 8a-n. To a mixture of the appropriate pyrazole aldehydes 4 (1 mmol) and a catalytic amount of glacial acetic acid (3-4 drops) in absolute ethanol (10 mL), were added hydrazide 7 (1.1 mmol) and refluxed for 12-18 h. As the reaction was completed, the mixture was allowed to cool to room temperature, the precipitate was filtered off and crystallized from ethanol to give the pure final derivatives 8a-n (see related spectra in Supplementary data).    13 13 13          www.nature.com/scientificreports/ Cytotoxicity evaluation of the most potent compounds on 3T3 cell line. The cytotoxicity of the compounds 8e, 8l, 8f and 8n was determined using the 3-(4,5 Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assay according to previously described methods 44 .

N'-((3-(4-Methoxyphenyl)-1-p-tolyl-1H-pyrazol-4-yl)methylene)benzofuran-2-carbohydrazide
Homology modeling and docking study. The primary sequence of Saccharomyces cerevisiae α-glucosidase downloaded from UniProtKB database (Uni-ProtKB, http:// www. unipr ot. org/) with accession number P53341. Hence, a search was executed to identify a protein with a high sequence similarity using NCBI BLAST server (https:// blast. ncbi. nlm. nih. gov/ Blast. cgi). The crystallographic structure of Saccharomyces cerevisiae (PDB ID: 3A4A) was chosen as a template and subjected to sequence alignment using Align sequence to template protocol in Discovery Studio v4.1 (DS) (Accelrys, San Diego, CA) (http:// accel rys. com). The 3D structure of α-glucosidase for S cerevisiae was predicted using the modeler module of the DS4.1. The number of models was set to 10, and the optimization level was changed to high. DOPE score and profile-3D for preliminary evaluation of developed model were carried out using verify protein protocol in DS. Ramachandran plot with PROCHECK program (http:// servi cesn. mbi. ucla. edu/ PROCH ECK) was also applied to verify the quality of the obtained homology model 52,53 . With the modeled structure, the docking of selected compounds was carried out. The predicted model was subjected to protein preparation using the prepare protein protocol in DS software. In this step, the complex was typed with CHARMm force field, hydrogen atoms were added, all water molecules were removed and pH of protein was adjusted to almost neutral, 7.4. The 3D structure of the most active compounds 8e, 8l, 8f and 8n were built using the Marvin v15.4.6 program (https:// chema xon. com/ produ cts/ marvin) and then were transferred into DS. Ligand structures were typed with CHARMm force field and partial charges were calculated by Momany-Rone option. Subsequently, resulting structures were minimized with Smart Minimizer, which performs 1000 steps of steepest descent with a RMS gradient tolerance of 3, followed by conjugate gradient minimization. Molecular docking approach has been carried out using genetic algorithm based docking program (GOLD v5.3.0) (https:// www. ccdc. cam. ac. uk/ solut ions/ csd-disco very/ compo nents/ gold/), which considers complete flexibility of side chains of amino acids at the active site. A 9 A˚ radius sphere was defined as a docking pocket. Acarbose and most potent compounds were docked into the active site of the protein and generated binding modes of ligands were ranked based on GOLD score fitness function. A top-score binding pose was selected for analyzing the interactions between the enzyme and the inhibitors using Discovery Studio Visualizer v4.1 (http:// accel rys. com/ produ cts/ disco very-studio) 54 .

Molecular dynamic simulation. Molecular dynamic (MD) simulation of this study was performed by
using the Desmond v5.3 module (https:// www. schro dinger. com/ produ cts/ desmo nd) implemented in Maestro interface (from Schrödinger 2018-4 suite) 55 . The appropriate pose for MD simulation procedure of the compounds was obtained by the docking method.
In order to build the system for MD simulation, the protein-ligand complexes were solvated with SPC explicit water molecules and placed in the center of an orthorhombic box of appropriate size in the Periodic Boundary Condition. Sufficient counter-ions and a 0.15 M solution of NaCl were also utilized to neutralize the system and to simulate the real cellular ionic concentrations, respectively. The MD protocol involved minimization, pre-production, and finally, production MD simulation steps. In the minimization procedure, the entire system was allowed to relax for 2500 steps by the steepest descent approach. Then the temperature of the system was raised from 0 to 300 K with a small force constant on the enzyme in order to restrict any drastic changes. MD simulations were performed via NPT (constant number of atoms, constant pressure i.e. 1.01325 bar and constant temperature i.e. 300 K) ensemble. The Nose-Hoover chain method was used as the default thermostat with 1.0 ps interval and Martyna-Tobias-Klein as the default barostat with 2.0 ps interval by applying isotropic coupling style. Long-range electrostatic forces were calculated based on Particle-mesh-based Ewald approach with the he cut-off radius for columbic forces set to 9.0 Å. Finally, the system was ubjected to produce MD simulations for 20 ns for the protein-ligand complex. During the simulation, every 1000 ps of the actual frame was stored. The dynamic behavior and structural changes of the systems were analyzed by the calculation of the root mean square deviation (RMSD) and RMSF. Subsequently, the energy-minimized structure calculated from the equilibrated trajectory system was evaluated to investigate of each ligand-protein complex interaction. where ΔG Bind is the calculated relative free energy which includes both ligand and receptor strain energy. E Complex is the MM-GBSA energy of the minimized complex, and E Ligand is the MM-GBSA energy of the ligand after removing it from the complex and allowing it to relax. E Receptor is the MM-GBSA energy of relaxed protein after separating it from the ligand. The MM-GBSA calculation was performed based on the clustering method for energy calculation.