Inhibition of microsomal prostaglandin E2 synthase-1 (mPGES-1) is promising for designing novel nonsteroidal anti-inflammatory drugs, as they lack side-effects associated with inhibition of cyclooxygenase enzymes. Azole compounds are nitrogen-containing heterocycles and have a wide use in medicine and are considered as promising compounds in medicinal chemistry. Various computer-aided drug design strategies are incorporated in this study. Structure-based virtual screening was performed employing various docking programs. Receiver Operator Characteristic (ROC) curves were used to evaluate the selectivity of each program. Furthermore, scoring power of Autodock4 and Autodock Vina was assessed by Pearson’s correlation coefficients. Pharmacophore models were generated and Güner-Henry score of the best model was calculated as 0.89. Binding modes of the final 10 azole compounds were analyzed and further investigation of the best binding (-8.38 kcal/mol) compound was performed using molecular dynamics simulation, revealing that furazan1224 (ZINC001142847306) occupied the binding site of the substrate, prostaglandin H2 (PGH2) and remained stable for 100 ns. Continuous hydrogen bonds with amino acids in the active site supported the stability of furazan1224 throughout the trajectory. Pharmacokinetic profile showed that furazan1224 lacks the risks of inhibiting cytochrome P450 3A4 enzyme and central nervous system-related side-effects.
Research Article
In Silico Discovery of Potential Azole-Containing mPGES-1 Inhibitors By Virtual Screening, Pharmacophore Modeling and Molecular Dynamics Simulations
https://doi.org/10.21203/rs.3.rs-1228363/v1
This work is licensed under a CC BY 4.0 License
You are reading this latest preprint version
mPGES-1
anti-inflammatory drug
virtual screening
pharmacophore modeling
molecular dynamics
Nonsteroidal anti-inflammatory drugs (NSAIDs) are very popular today as inflammatory responses are very common conditions. Compared to corticosteroid-based therapy, NSAIDs have significantly lower risks and side effects [1–3]. However, cardiovascular and gastrointestinal (GI) complications may develop due to NSAID use [4–6]. Since most NSAIDs are COX (cyclooxygenase-1, cyclooxygenase-2) inhibitors, new targets have been explored for inhibition. Prostaglandin E2 (PGE2) is one of the main inflammatory mediators produced in the arachidonic acid cascade and synthesized by cytosolic prostaglandin E2 synthase (cPGES) and microsomal prostaglandin E2 synthases (mPGES-1 and mPGES-2). mPGES-2 and cPGES enzymes are not inducible, whereas mPGES-1 is known to release PGE2 consequent to an inflammatory signal. Inhibition of mPGES-1 is therefore promising for designing novel NSAIDs without the side-effects which have been linked to COX-inhibition. Moreover, mPGES-1 has been reported to play important roles in carcinogenesis, angiogenesis, migration, invasion and apoptosis. There are also studies demonstrating that, increased amount of PGE2 is linked with many pathological conditions like gastric ulcer, GI cancer, colon, prostate, lung cancer, glioma and neuroblastoma, as well as non-carcinogenic pathologies like Alzheimer’s disease, Parkinson’s disease and Multiple Sclerosis [7].
Being a group of nitrogen-containing heterocycles, azole compounds have a wide variety of use. They are effective against mycosis for both humans and animals; therefore, they are used in many anti-fungal medications including veterinary medicines [8]. Being electron-rich heterocycles, imidazole and benzimidazole easily accept and donate protons. They also form various weak physical bonds. These features allow them bind with therapeutic targets which brings about pharmacological activity [9–11]. Besides having medicinal potential against many biological targets [12, 13], azole derivatives are also considered promising for mPGES-1 inhibition, e.g., 2-aminothiazole derivatives were found to have inhibitory activity against PGE2 production [14]. Some triazole and oxadiazole derivatives have also been reported to inhibit mPGES-1 in cell-free assay models [15, 16].
Virtual screening (VS) is a cheminformatics tool to discover new molecules, usually expected to bind with a certain target, e.g. enzyme, from a large compound set [17–19]. The main goal of VS is to identify compounds with a specific therapeutic potential very fast so as to filter down for biological assays. VS can be performed with two approaches; ligand-based and structure-based. Ligand-based methods exploit the similarity of compounds with the known bioactive compound. Structure-based methods rely on the cavity of the target receptor and the interactions a compound can form within. In the structure-based virtual screening (SBVS), docking stands as a very important tool to predict the affinity of a molecule to another, namely, a macromolecule (target). Although molecular docking applications are very useful at identifying hit compounds, they have weaknesses as well as strengths. For example, a scoring function (an essential component of docking applications) might be very good at discarding true negatives (inactive molecules), but it might not be successful at ranking binding affinities, hence the need for incorporating various scoring functions. Pharmacophore models are another computer-aided drug design (CADD) tools that are created according to the knowledge derived from structure-activity relationships and structural properties of active analogs, therefore they are better at distinguishing inactive ligands. A pharmacophore can be determined by direct method which uses receptor-ligand complexes or indirect method which makes use of the information that a collection of ligands known to interact with a certain receptor [20]. Molecular dynamics (MD) simulations offer a visualization of stable binding of the ligand in the active site of the macromolecule. This allows researchers to estimate ligand binding free energies [21]. In the present study, a CADD protocol incorporating VS, pharmacophore modelling and MD was followed. Binding mode of the best candidate, obtained after filtering ~1.000.000 azole compounds using various methods, were closely investigated by means of MD simulations.
2.1 Library construction and filtering
A compound library was constructed with 1,011,637 azole compounds retrieved from ZINC database [22]. Hydrogens were added and three-dimensional structures were generated with Babel [23]. The library was “cleaned” by using SMARTS (SMiles ARbitary Target Specification) patterns implemented in Discovery Studio Client 2016 [24]. The patterns that were used for filtering are provided in the Supplementary Material Table S1. Next, the entire library was filtered according to Lipinski’s rule of five [25]. Briefly, the rule states that the molecular weight of an oral drug should be lower than 500Da, the number of hydrogen bond acceptors should be at most 10, the number of hydrogen bond donors should be at most 5, and log P should be lower than 5, allowing only one violation. Therefore, filtered ligands were decreased to 971,040.
2.2. Design of the workflow and evaluation of selectivity of the scoring functions
Design of the workflow in VS is critical as capabilities of computational tools are different. For the purpose of discarding inactive compounds early, a scoring function that is good at distinguishing between actives and inactives (selectivity) would be the best option for the first stage of screening. rDock is known to eliminate inactive molecules at an early stage [26]. However, scoring functions are target-dependent. For this reason, 3 different docking programs were evaluated for their use at a specific stage. Receiver Operator Characteristic (ROC) curve for each program was obtained and Area Under the Curve (AUC) was calculated. Toward this end, 50 known active compounds were downloaded from cHEMBL database and 36 decoys for each ligand [27], 1800 decoys were retrieved from DUD-E database [28] (A database of useful decoys: Enhanced, http://dude.docking.org/). rDock was selected as the first screening tool according to the AUC profile. Moreover, Autodock Vina was used for the second screening. Next, pharmacophore-based screening was performed to identify potential mPGES-1 inhibitors. Finally, Autodock4 was used to rank compounds according to their binding affinities before MD simulations. The workflow followed in this study is demonstrated in Figure 1.
2.2.1 Screening with rDock
Prior to screening with rDock [26], the cavity definition was made by reference molecule method as described in the manual. Toward this end, structure of a known inhibitor [29] with a high in vitro inhibition activity (0.241 µM) was used (Figure 2). Crystal structure of mPGES-1 (PDB ID: 4AL0) was downloaded from PDB Database [30] and prepared in Discovery Studio Client 2016.
The inhibitor was docked using rDock to the crystal structure of mPGES-1 enzyme in the presence of glutathione to ensure a minimized ligand-receptor interaction energy. Number of Genetic Algorithm (GA) runs were kept at 50. The binding pose with the lowest energy was used for cavity mapping. Furthermore, parameter file was prepared where the following parameters were set: Radius: 6.0 Å; Small sphere: 1.5 Å; Minimum volume: 100 Å (default); Maximum cavities: 1 (default); Grid step: 0.5 (default); RECEPTOR_FLEX: 3.0 Å (default).
Considering the size of the chemical space (971,040 compounds), High-Throughput Virtual Screening (HTVS) protocol for rDock was followed. Toward this end, the first docking with all compounds were performed with 5 runs. Later, ligands exhibiting binding score of -10 or below were subjected to further docking with 10 runs. Next, for ligands which scored below -20 were filtered for further docking simulations with Autodock Vina. Filtering at every step was performed according to SCORE.INTER, which refers to the ligand’s binding energy.
2.2.2. Further structure-based virtual screening with Autodock Vina
According to ROC profile, Autodock Vina ranks the second best at distinguishing active and inactive compounds. Furthermore, the scoring function of Vina is known for its speed and accuracy [31]. For this reason, screening was pursued with Autodock Vina. Filtered ligands according to their scores in rDock were further subjected to docking with Autodock Vina. Prior to this, tautomers and protonation states were generated using Prepare Ligands in Discovery Studio Client 2016. Ligands which raised problems of generating 3D coordinates, bad valences or being duplicates were eliminated. Consequently, 44,366 compounds were subjected to molecular docking simulations with Vina. A grid box was prepared with center coordinates of x: 11.749, y: -12.53, z: -12.147 and dimensions of x: 60 Å, y: 60 Å and z: 60 Å. The center was chosen at the sulfur atom of GSH.
2.2.3. Pharmacophore generation
Remaining compounds in the library were screened with pharmacophore to recover more drug-like compounds. In order to generate a pharmacophore model, a training set comprising known active mPGES-1 inhibitors was prepared. Toward this end, 20 molecules were retrieved from the literature and cHEMBL database. The entire set was docked to the crystal structure of mPGES-1 enzyme with Autodock Vina to determine their binding affinities (Table 1).
Next, the training set was subjected to Common feature pharmacophore generation protocol in Discovery Studio Client 2016. 10 hypotheses were generated by HipHop algorithm. Conformation generation was selected BEST, with an energy threshold of 20 kcal/mol and maximum conformation of 250. Maximum interference distance was set to 2.97 and the number of maximum pharmacophores to be generated was set to 10.
2.2.4. Validation of the pharmacophore model
Each generated pharmacophore model was validated by using a test set including 380 decoys and 20 compounds that are known to be active mPGES-1 inhibitors (the training set). Decoys were downloaded from DUD-E database (A database of useful decoys: Enhanced, http://dude.docking.org/). The quality of the pharmacophore often refers to its power to differentiate between inactive compounds and active compounds. Towards this end, the test set was screened against each pharmacophore and several statistical parameters such as percentage yield of actives, ratio of actives and enrichment factor (EF) were estimated. Güner-Henry score [32], which determines the quality of a pharmacophore model, was also calculated.
$$\%A=\frac{Ha}{\text{A}} \times 100\%{(Eq.1)}$$
(Eq. 2)
(Eq. 3)
\(\)\(GH=\frac{Ha\left(3Ha+Ht\right)}{4\text{H}\text{t} \times \text{A}}\times \left(1-\frac{Ht-Ha}{D-A}\right)\) (Eq. 4)
Where %A is the ratio of actives in the entire hit list, %Ya is the fraction of active hits relative to total hits, Ht refers to total hits acquired during database search, Ha refers to active compounds in the hit list, D is the total number of compounds, and A is the number of active compounds in the test set.
2.2.5. Screening the ligand library with pharmacophore model
After validating each pharmacophore model, hypothesis 01 was selected to further screen the ligand library. Upon filtering with this model, 48 compounds remained in the library. 10 compounds having the lowest Autodock Vina scores and Fit Value ≥3 were selected for further investigation.
2.3. Optimization with quantum mechanical methods
Selected compounds were subjected to geometry optimization using Gaussian [33]. M062X method and 6-31G(d,p) basis set were selected for the calculations. All calculations were performed on Linux platform.
2.4. Evaluation of scoring power and docking with Autodock4
Experimental data (IC50, nM) of 50 known active compounds (2.2. Design of the workflow and evaluation of selectivity of the scoring functions) were obtained from cHEMBL database plotted against the highest score estimated for each compound. Pearson’s correlation coefficient was calculated for each docking program and Autodock4 was found to be better at ranking binding affinities of ligands. Optimized structures were subjected to docking using Autodock4 [34]. Toward this end, the crystal structure of mPGES-1 enzyme was cleaned and prepared. Grid box was defined as described in section 2.2.2 and Autogrid4 was employed for grid calculation. Conformational space was searched by Lamarckian Genetic Algorithm [35], which was implemented in AutoDock4. Maximum number of energy evaluations was set to 2,500,000 and 100 runs were performed for each ligand. All calculations were performed on Linux platform.
2.5. Molecular dynamics simulation
The azole derivative with the highest affinity in complex with mPGES-1 enzyme was selected for MD simulations. mPGES-1 enzyme was defined with CHARMM36 force-field while GSH and ligand were parameterized by using CGenFF server (https://cgenff.umaryland.edu). The complex (enzyme with GSH and ligand) was solvated in TIP3P water [36] using a periodic cubic box and was minimized by using steepest descent algorithm with GROMACS 5.1.4 package [37]. Minimization was complete as the maximum force of 10 kJ/mol was reached. Next, the system was subjected to NVT equilibration for 100 ps with leap-frog integrator, V-rescale thermostat [38] and T=300 K. This was followed by NPT equilibration to equilibrate volume for 100 ps, also with leap-frog integrator, V-rescale thermostat and T=300 K. Finally, the equilibrated complex was submitted to MD simulations for 100 ns. For both stages of equilibration and MD simulation, LINCS algorithm [39] was used for constraining bond lengths and Particle Mesh Ewald (PME) [40] was employed for long-range electrostatic interactions (with a cut-off of 1.2 nm). V-rescale algorithm was selected for temperature coupling and Parrinello-Rahman barostat [41] for pressure coupling. The simulation time-step was set to 2 fs and data were recorded every 10 ps. Clustering of conformations was performed based on root-mean-square deviation using gmx cluster tool. GROMOS clustering algorithm [42] and a cut-off of 0.16 nm was selected. Optimal cut-off value was chosen according to the plot provided in the Supplementary Material Figure S1. Central complex structure was generated from the most populated cluster.
3.1. ROC curves
ROC (Receiver operating characteristic) curve is a recognized method to evaluate the ability of a test to distinguish between two sets of data [43]. In the field of drug design, this approach is used to discriminate “active” and “inactive” compounds in virtual screening. Eliminating inactive compounds from a large set of compounds at an early stage in VS is critical. Therefore, decision on the right scoring function is beneficial for a successful screening. ROC curves were generated for three different docking programs; rDock, Autodock Vina and Autodock4 (Figure 3). Sensitivity was plotted against 1-specificity for each threshold of the score estimated by the scoring function. Sensitivity can be described as the ratio of true positives to the sum of true positive and false negative, and specificity can be described as the ratio of true negatives to the sum of true negatives and false positives. AUC (area under the curve) values were calculated for each curve. An AUC value of 0.5 indicates a test that selects randomly, whereas an AUC value of 1 (highest possible value) indicates an ideal test for discriminating two different populations. According to AUC values calculated for each program, rDock is better at discriminating actives and inactives than the other two.
3.2. Evaluation of scoring power
A scoring function with a high ability to determine the correct binding pose does not ensure its performance on ranking the binding affinity [44]. For this reason, Pearson’s correlation coefficients (rP) were calculated using experimental binding affinities (IC50) to evaluate the scoring power of Autodock4 and Autodock Vina (Table 2).
|
Docking Program |
rP |
p |
|---|---|---|
|
Autodock4 |
0.584 |
< 0.05 |
|
Autodock Vina |
0.513 |
< 0.05 |
| p: Significance level | ||
A value of rP close to 1 indicates a very strong correlation between experimental data and the docking scores, whereas rP being close to 0 shows a poor correlation. Both programs showed a fair correlation while Autodock4 was found to be superior at ranking binding affinities for mPGES-1 enzyme.
3.3. Common feature pharmacophore model
10 best pharmacophore hypotheses were generated after Common Feature Generation was employed for the training set, scores ranging between 78.870 – 65.864 (Table 3). Hypotheses can be grouped as HDAA (01-07) and HAD (08-10).
|
Hypothesis |
Featuresa |
Rankb |
Direct Hitc |
Partial Hitd |
Max Fit |
|---|---|---|---|---|---|
|
01 |
HDAA |
78.870 |
1111111111 |
0000000000 |
4 |
|
02 |
HDAA |
77.058 |
1111111111 |
0000000000 |
4 |
|
03 |
HDAA |
76.993 |
1111111111 |
0000000000 |
4 |
|
04 |
HDAA |
74.604 |
1111111111 |
0000000000 |
4 |
|
05 |
HDAA |
74.569 |
1111111111 |
0000000000 |
4 |
|
06 |
HDAA |
71.187 |
1111111111 |
0000000000 |
4 |
|
07 |
HDAA |
70.856 |
1111111111 |
0000000000 |
4 |
|
08 |
HDA |
68.675 |
1111111111 |
0000000000 |
3 |
|
09 |
HDA |
66.399 |
1111111111 |
0000000000 |
3 |
|
10 |
HDA |
65.864 |
1111111111 |
0000000000 |
3 |
| a H depicts hydrophobic group; A, hydrogen bond acceptor; D, hydrogen bond donor. | |||||
| b Score of the training set compounds’ ranking that fit the hypothesis. | |||||
| c Indicates whether a molecule in the training set mapped each feature in the hypothesis. “1” depicts being mapped and “0” depicts not. | |||||
| d Indicates a particular molecule in the training set mapped all but one feature in the hypothesis. Both Direct Hit and Partial Hit show numeration of molecules from right to left. | |||||
3.4. Pharmacophore validation
Each generated hypothesis was investigated by using certain statistical parameters such as the number of total hits (Ht), the number of active hits (Ha), enrichment factor (EF), % yield of actives and % ratio of actives. Güner-Henry scoring method was used to determine the discriminating power of a pharmacophore model. Table 4 shows these parameters for each hypothesis.
|
Hypothesis |
Da |
Ab |
Htc |
Had |
%Yae |
%Af |
EFg |
GHh |
|---|---|---|---|---|---|---|---|---|
|
01 |
400 |
20 |
23 |
20 |
86.95 |
100 |
17.39 |
0.89 |
|
02 |
400 |
20 |
39 |
16 |
41.02 |
80 |
8.20 |
0.48 |
|
03 |
400 |
20 |
28 |
14 |
50 |
70 |
10 |
0.53 |
|
04 |
400 |
20 |
49 |
19 |
38.77 |
95 |
7.75 |
0.49 |
|
05 |
400 |
20 |
28 |
16 |
57.14 |
80 |
11.42 |
0.72 |
|
06 |
400 |
20 |
34 |
17 |
50 |
85 |
10 |
0.56 |
|
07 |
400 |
20 |
56 |
20 |
35.71 |
100 |
7.14 |
0.46 |
|
08 |
400 |
20 |
39 |
15 |
38.46 |
75 |
7.69 |
0.44 |
|
09 |
400 |
20 |
63 |
18 |
28.57 |
90 |
5.71 |
0.39 |
|
10 |
400 |
20 |
36 |
17 |
47.22 |
85 |
9.44 |
0.54 |
| a D is the total number of compounds in the test set. | ||||||||
| b A is the number of active compounds. | ||||||||
| c Ht is the number of total hits. | ||||||||
| d Ha is the number of hits for active compounds. | ||||||||
| e %Ya is the yield of actives, refers to the selectivity. | ||||||||
f %A is the ratio of actives retrieved in the hit list.
g EF is the enrichment factor.
h GH refers to the Güner-Henry score.
According to Table 4, hypotheses 01 and 07 identify all of the 20 active compounds. However, as Ht also has a high value, hypothesis 07 presents poor selectivity. A GH score higher than 0.7 indicates a very good and reliable model [45, 46]. Hypothesis 05 exhibit a sufficient GH score, nevertheless, 16 active compounds mapped to the model out of 20 lowered the selectivity. Hypothesis 01, displaying an excellent GH score of 0.89 suggests that this is a reliable model for identifying active compounds in a database. Therefore, hypothesis 01 was chosen as the best pharmacophore model. Figure 4 shows the hypothesis model 01 with its features: Hydrophobic point (H), hydrogen bond acceptor (A) and hydrogen bond donor (D).
3.5. Screening the ligand library with pharmacophore model
After validating each pharmacophore model, hypothesis 01 was selected to screen the ligand library. 10 compounds with fit values ≥ 3 and having the lowest Autodock Vina scores were selected for Autodock4 docking.
3.6. Autodock4 results
Remaining ten compounds after pharmacophore screening were subjected to geometry optimization and docking with Autodock4. Binding energies and calculated inhibition constants are given in Table 5.
3.7. Docking and binding mode analysis
Detailed interaction types and distances are provided in the Supplementary Material Table S2. The most commonly observed binding mode seems to take advantage of surrounding hydrogen bond donors in the catalytic channel. A more detailed look to the binding mode and interaction types reveals that this binding mode, which the final ten compounds share (Supplementary Material, Figure S2), is in the substrate-binding region. This is supported by the mechanism suggested by Sjögren et al, which requires GSH, Arg126 and Ser127 in close contact with the substrate, PGH2 [47]. Additionally, Jegerschold et al proposed that GSH and Arg126 are directly involved in the mechanism [48]. Therefore, it is important to note these close interactions between the cofactor GSH, Arg126, Ser127 and the inhibitor. Since furazan1224 forms all the critical interactions for the enzyme reaction and binds with the lowest energy (-8.38 kcal/mol), this inhibitor was used to demonstrate interactions in the cavity of mPGES-1 (Figure 5). The reader is referred to Figure S3-S11 for interactions of other inhibitors. Furazan1224 forms hydrogen bonds with Arg126, Ser127, Leu121 and Arg70. Guanidine groups of Arginine residues are mostly responsible for the interactions. Furthermore, aromatic rings of the ligand make π-cation interactions with the guanidine groups of Arg126 and Arg67. Ala123 forms π-alkyl interactions with furazan1224.
3.8. Molecular dynamics simulation
The strongest-binding azole compound (furazan1224) was chosen for 100 ns molecular dynamics simulation as the other azole derivatives share the same binding mode but differ in affinity. Root mean square deviation of the ligand is shown in Figure 6, which reveals the stability of furazan1224 throughout the trajectory. RMSD of the protein and GSH are shown in Figure S12-S13. The central complex structure of the most populated cluster was defined as to represent the dominant conformation during the simulation (Figure S14). The central structure has the smallest average RMSD from all other structures in the same cluster. An RMSD cut-off of 0.16 nm yielded 19 clusters with only 2 one-member clusters. Output files of clusters calculated with various cut-off values are provided in the Supplementary Material for interested readers. Figure S15 displays new amino acids contributing to the stability of furazan1224 during the simulation (Lys120, Arg122) while some amino acids observed in the docking complex are absent (Arg70, Ala123, Ser127). However, contacts with Arg67, Leu121 and Arg126 remain throughout the trajectory. Stability of the ligand is supported by hydrogen bonding profile involving certain amino acids in the binding site (Figure 7, 9, 10). Being the cofactor of the enzyme, GSH plays a key role in the enzymatic reaction in concord with Arg126 and Ser127. Arg126 is very important as it has a key role in the final step of the enzymatic reaction which gives PGE2 [47] and is in close contact with GSH and furazan1224 in the docking complex (Figure 5). Amide oxygen of the ligand (O2) makes very close and stable contact with the nearest hydrogen of guanidine group of Arg126 (Figure 7). It is often expected to see shifts in hydrogen bond distances when there is another hydrogen bond donor nearby, due to rotation of single bonds. However, electron delocalization constraints this movement, giving a very stable hydrogen bond between HH11 and O2 (Figure 8). Therefore, the distance between HH12 and O2 remains at a certain range throughout the trajectory (Figure 7).
Another close contact of furazan1224 was observed with Leu121. Two hydrogen bonds are observed which seem consistent throughout the trajectory (Figure 9). However, H1 of furazan1224 makes more consistent hydrogen bond with amide oxygen of Leu121, with respect to H4. 3D representation of these contacts is given in Figure S16.
Arg122 is observed to contribute for stabilization of furazan1224 throughout the simulation although it is absent in the interaction map of the docking complex (Figure 5). Figure 10 shows stable hydrogen bond interaction between amide oxygen (O) of Arg122 and the hydrogen (H4) belonging to the furazan ring of furazan1224. Neighboring hydrogen (H1) does not make very close contact with oxygen of Arg122, but forms continuous hydrogen bond with oxygen (O) of Leu121 (Figure 9). These stable hydrogen bonds and steady distances between the ligand and amino acids show that furazan1224 is very stable during the simulation.
3.9. Pharmacokinetics
Since the compound library was constructed according to Lipinski’s rule of five, all molecules in the final set comply to the rule. However, certain pharmacokinetic properties were not considered during library construction and were calculated by using SwissADME (http://www.swissadme.ch). The results are given in Table 6. P-glycoprotein (P-gp) is an ATP-dependent transmembrane protein that is a member of ABC transporter family. It has a crucial role in the absorption and penetration of drugs through blood-brain barrier (BBB) [49]. Furazan1224 was predicted as a P-gp substrate, however its GI absorption was predicted high and it does not appear BBB permeable. For drugs designed for peripheral action, when central nervous system (CNS) is not their primary site of action, BBB permeability can cause undesirable CNS-related side-effects [50]. The estimated pharmacokinetic profile of the selected compounds showed that none of them has BBB permeability.
Cytochrome P450 (CYP) isoforms are very important for drug metabolism in humans. Inhibition of these CYPs may lead to serious clinical consequences like drug-drug interactions. Especially, inhibition of CYP3A4 can lead to serious adverse effects as nearly half of the drugs in the market are metabolized by this enzyme [51]. Of the compounds selected by pharmacophore analysis and molecular docking simulations, only furazan1199 was estimated to inhibit CYP3A4 isozyme. The best inhibitor (furazan1224) identified by virtual screening appears not to inhibit CYP3A4 showing that it does not have a potential to elevate blood levels of drugs metabolized by this enzyme.
|
ID |
GI absorption1 |
BBB permeant2 |
P-gp substrate3 |
CYP1A2 inhibitor4 |
CYP2C19 inhibitor5 |
CYP2C9 inhibitor6 |
CYP2D6 inhibitor7 |
CYP3A4 inhibitor8 |
|---|---|---|---|---|---|---|---|---|
|
furazan1224 |
High |
No |
Yes |
No |
Yes |
No |
No |
No |
|
furazan3510 |
High |
No |
Yes |
No |
No |
No |
No |
No |
|
furazan1199 |
High |
No |
Yes |
Yes |
No |
No |
Yes |
Yes |
|
furazan909 |
High |
No |
Yes |
No |
No |
No |
No |
No |
|
furazan856 |
High |
No |
No |
No |
No |
No |
No |
No |
|
furazan1014 |
High |
No |
Yes |
Yes |
No |
No |
No |
No |
|
oxadiazole1769 |
High |
No |
Yes |
No |
No |
No |
No |
No |
|
triazole1047 |
Low |
No |
No |
No |
No |
No |
No |
No |
|
oxadiazole1006 |
High |
No |
No |
No |
No |
No |
No |
No |
|
oxadiazole1011 |
High |
No |
No |
No |
No |
No |
No |
No |
| 1 Gastrointestinal absorption | ||||||||
| 2 Blood-brain barrier permeability | ||||||||
| 3 P-glycoprotein substrate | ||||||||
| 4 Cytochrome P450 1A2 inhibitor | ||||||||
| 5 Cytochrome P450 2C19 inhibitor | ||||||||
6 Cytochrome P450 2C9 inhibitor
7 Cytochrome P450 2D6 inhibitor
8 Cytochrome P450 3A4 inhibitor
Various CADD methods were pursued to identify novel azole compounds that would inhibit human mPGES-1 enzyme. Docking studies provided filtering down a large molecule set to only the best-binding compounds. Choice of scoring function at the former and succeeding stages of VS was made according to ROC results and calculated Pearson’s correlation coefficients. Pharmacophore models were successfully constructed. The best hypothesis (Güner-Henry score: 0.89) was used to identify compounds that would possess similar features with known mPGES-1 inhibitors. Further docking studies and binding mode analysis showed that furazan1224, similar to nine other compounds, occupied the binding site of the substrate PGH2, and took advantage of both hydrogen bond donors and acceptors in the mPGES-1 cavity. Molecular dynamics simulations also demonstrated that furazan1224 occupied the binding site of PGH2, and remained stable for 100 ns. Hydrogen bonds analysis supports this as stable hydrogen bonds interactions were observed with amino acids in the binding site, including Leu121 and Arg126. Calculated pharmacokinetic properties show that only furazan1199 among the final azole compounds inhibits the CYP3A4 isozyme, which is a critical enzyme for metabolizing most drugs. Furazan1224 is promising as a potential mPGES-1 inhibitor regarding its stability in the enzyme cavity and considerably safe pharmacokinetic profile.
Funding
Author L. Ozalp has received research support from the Council of Higher Education (YÖK) of Turkey (PhD scholarship of 100/2000 program) and TUBİTAK 2211/C National PhD Scholarship Program in the Priority Fields in Science and Technology.
Competing Interests
The authors have no relevant financial or non-financial interests to disclose.
Availability of data and material
Not applicable
Code availability
Not applicable
Authors' contributions
Lalehan Ozalp: Research, analysis, writing
İlkay Küçükgüzel: Writing, revision, editing
Ayşe Ogan: Revision, editing
Acknowledgments
We thank Ammar Elmezayen and Gizem Nur Duran for guidance and helpful discussions. We are also grateful to Prof. Kemal Yelekçi for Discovery Studio Client 2016 package. Molecular dynamics simulation calculations reported in this paper were performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources). L. Ozalp acknowledges the Council of Higher Education (YÖK) of Turkey for the PhD scholarship of 100/2000 program and TUBİTAK 2211/C National PhD Scholarship Program in the Priority Fields in Science and Technology.
- Cushing H (1964) The basophil adenomas of the pituitary body and their clinical manifestations (pituitary basophilism). J Neurosurg 21:318–347
- Williams TJ, RE O’Hehir D, Czarny M, Horne, Bowes G (1988) Acute Myopathy in Severe Acute Asthma Treated with Intravenously. Am Rev Respir Dis 137:460–463
- Hasselgren P-O (1999) Glucocorticoids and muscle catabolism. Curr Opin Clin Nutr Metab Care 2:201–205
- Mukherjee D, Nissen SE, Topol EJ (2001) Risk of cardiovascular events associated with selective COX-2 inhibitors. Jama 286:954–959
- Praveen Rao PN, Knaus EE (2008) Evolution of nonsteroidal anti-inflammatory drugs (NSAIDs): Cyclooxygenase (COX) inhibition and beyond. J Pharm Pharm Sci 11:81–110
- Serhan CN, Ward PA, Gilroy DW (2010) Fundamentals of inflammation. Cambridge University Press
- Bülbül B, İ Küçükgüzel (2019) Microsomal prostaglandin E2 synthase-1 as a new macromolecular drug target in the prevention of inflammation and cancer. Anticancer Agents Med Chem 19: 1205–1222
- Bhagat J, Singh N, Nishimura N, Shimada Y (2021) A comprehensive review on environmental toxicity of azole compounds to fish. Chemosphere 262:128335
- Wright JB (1951) The Chemistry of the Benzimidazoles. Chem Rev 48:397–541
- Bhatnagar A, Sharma PK, Kumar N (2011) A review on “Imidazoles”: Their chemistry and pharmacological potentials. Int J PharmTech Res 3:268–282
- Gaba M, Mohan C (2016) Development of drugs based on imidazole and benzimidazole bioactive heterocycles: Recent advances and future directions. Springer US)
- Sharma PC, Sinhmar A, Sharma A, Rajak H, Pathak DP (2013) Medicinal significance of benzothiazole scaffold: an insight view. J Enzyme Inhib Med Chem 28:240–266
- Keri RS, Patil MR, Patil SA, Budagumpi S (2015) A comprehensive review in current developments of benzothiazole-based molecules in medicinal chemistry. Eur J Med Chem 89:207–251
- Smith B, Chang H-H, Medda F, Gokhale V, Dietrich J, Davis A, Meuillet EJ, Hulme C (2012) Synthesis and biological activity of 2-aminothiazoles as novel inhibitors of PGE2 production in cells. Bioorg Med Chem Lett 22:3567–3570
- He S, Li C, Liu Y, Lai L (2013) Discovery of highly potent microsomal prostaglandin e2 synthase 1 inhibitors using the active conformation structural model and virtual screen. J Med Chem 56:3296–3309
- Erensoy G, Ding K, Zhan C-G, Elmezayen AD, Yelekçi K, Duracık M, Özakpınar ÖB, İ Küçükgüzel (2020) Synthesis, in silico studies and cytotoxicity evaluation of novel 1, 3, 4-oxadiazole derivatives designed as potential mPGES-1 inhibitors. J Res Pharm 24: 436–451
- Shoichet BK (2006) Virtual screening of chemical libraries. Nature 432:862–865
- Bologa CG, Ursu O, Oprea TI (2019) How to Prepare a Compound Collection Prior to Virtual Screening. In: Bioinformatics and Drug Discovery. Humana Press, New York, pp 119–138
- Gimeno A, Ojeda-Montes MJ, Tomás-Hernández S, Cereto-Massagué A, Beltrán-Debón R, Mulero M, Pujadas G, Garcia-Vallvé S (2019) The light and dark sides of virtual screening: What is there to know? Int J Mol Sci 20
- Khedkar SA, Malde AK, Coutinho EC, Srivastava S (2007) Pharmacophore modeling in drug discovery and development: an overview. Med Chem (Los Angeles) 3:187–197
- Ozalp L, S Sağ Erdem B, Yüce-Dursun Ö, Mutlu, Özbil M (2018) Computational insight into the phthalocyanine-DNA binding via docking and molecular dynamics simulations. Comput Biol Chem 77
- Irwin JJ, Shoichet BK (2005) ZINC– a free database of commercially available compounds for virtual screening. J Chem Inf Model 45:177–182
- O’Boyle NM, Banck M, James CA, Morley C, Vandermeersch T, GR Hutchison (2011) Open Babel: An open chemical toolbox. J Cheminform 3: 33
- BIOVIA (2016) Discovery Studio Modeling Environment
- Lipinski CA, Lombardo F, Dominy BW, Feeney PJ (1997) Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev 23:3–25
- Ruiz-Carmona S, Alvarez-Garcia D, Foloppe N, Garmendia-Doval AB, Juhos S, Schmidtke P, Barril X, Hubbard RE, Morley SD (2014) rDock: A Fast, Versatile and Open Source Program for Docking Ligands to Proteins and Nucleic Acids. PLoS Comput Biol 10:1–7
- Huang N, Shoichet BK, Irwin JJ (2006) Benchmarking sets for molecular docking. J Med Chem 49:6789–6801
- Mysinger MM, Carchia M, Irwin JJ, Shoichet BK (2012) Directory of useful decoys, enhanced (DUD-E): better ligands and decoys for better benchmarking. J Med Chem 55:6582–6594
- Partridge KM, Antonysamy S, Bhattachar SN, Chandrasekhar S, Fisher MJ, Fretland A, Gooding K, Harvey A, Hughes NE, Kuklish SL, Luz JG, Manninen PR, McGee JE, Mudra DR, Navarro A, Norman BH, Quimby SJ, Schiffler MA, Sloan AV, Warshawsky AM, Weller JM, York JS, Yu XP (2017) Discovery and characterization of [(cyclopentyl)ethyl]benzoic acid inhibitors of microsomal prostaglandin E synthase-1. Bioorganic Med Chem Lett 27:1478–1483
- Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, Shindyalov IN, PE Bourne (2000) The protein data bank. Nucleic Acids Res 28: 235–242
- Trott O, Olson AJ (2010) Software news and update autodock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 31:455–461
- Guner O, Clement O, Kurogi Y (2004) Pharmacophore modeling and three dimensional database searching for drug design using catalyst: recent advances. Curr Med Chem 11:2991–3005
- Frisch M, Trucks GW, Schlegel HB, Scuseria GE, Scalmani RMACJRG, Barone V, Mennucci B, Petersson GA (2009) Gaussian 09. Gaussian Inc.), Wallingford
- Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS (2009) AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility. J Comput Chem 30:2785–2791
- Morris GM, Goodsell DS, Halliday RS, Huey R, Hart WE, Belew RK, Olson AJ (1998) Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J Comput Chem 19:1639–1662
- Jorgensen W, Chandrasekhar J, Madura J, Impey R, Klein M (1983) Comparison of Simple Potential Functions for Simulating Liquid Water. J Chem Phys 79:926–935
- Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, Lindah E (2015) Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2:19–25
- Bussi G, Donadio D, Parrinello M (2007) Canonical sampling through velocity rescaling. J Chem Phys 126:14101
- Hess B, Bekker H, Berendsen HJC, Fraaije JGEM (1997) LINCS: A linear constraint solver for molecular simulations. J Comput Chem 18:1463–1472
- Darden TA, York D, Pedersen LJ (1993) Particle mesh Ewald: An Nṡlog(N) method for Ewald sums in large systems. J Chem Phys 98:10089–10092
- Parrinello M, Rahman A (1981) Polymorphic transitions in single crystals: A new molecular dynamics method. J Appl Phys 52:7182–7190
- Daura X, Gademann K, Jaun B, Seebach D, van Gunsteren WF, AE Mark (1999) Peptide Folding: When Simulation Meets Experiment. Angew Chemie Int Ed 38: 236–240
- Triballeau N, Acher F, Brabet I, Pin J-P, Bertrand H-O (2005) Virtual screening workflow development guided by the “receiver operating characteristic” curve approach. Application to high-throughput docking on metabotropic glutamate receptor subtype 4. J Med Chem 48:2534–2547
- Wang Z, Sun H, Yao X, Li D, Xu L, Li Y, Tian S, Hou T (2016) Comprehensive evaluation of ten docking programs on a diverse set of protein–ligand complexes: the prediction accuracy of sampling power and scoring power. Phys Chem Chem Phys 18:12964–12975
- Jang C, Yadav DK, Subedi L, Venkatesan R, Venkanna A, Afzal S, Lee E, Yoo J, Ji E, Kim SY (2018) Identification of novel acetylcholinesterase inhibitors designed by pharmacophore-based virtual screening, molecular docking and bioassay. Sci Rep 8:1–21
- Zeb A, Son M, Yoon S, Kim JH, Park SJ, Lee KW (2019) Computational simulations identified two candidate inhibitors of Cdk5/p25 to abrogate tau-associated neurological disorders. Comput Struct Biotechnol J 17:579–590
- Sjögren T, J Nord, M Ek, P Johansson, G Liu and S Geschwindner (2013). Crystal structure of microsomal prostaglandin E₂ synthase provides insight into diversity in the MAPEG superfamily. Proc Natl Acad Sci U S A 110: 3806–3811.
- Jegerschold C, Pawelzik SC, Purhonen P, Bhakat P, Gheorghe K, Gyobu N, Mitsuoka K, Morgenstern R, Jakobsson PJ, Hebert H (2008) Structural basis for induced formation of the inflammatory mediator prostaglandin E2. Proc Natl Acad Sci U S A 105:11110–11115
- Lin JH, Yamazaki M (2003) Role of P-glycoprotein in pharmacokinetics. Clin Pharmacokinet 42:59–98
- Van De Waterbeemd H, Gifford E (2003) ADMET in silico modelling: towards prediction paradise? Nat Rev Drug Discov 2:192–204
- Manikandan P, Nagini S (2018) Cytochrome P450 structure, function and clinical significance: a review. Curr Drug Targets 19:38–54
Due to technical limitations, table 1 & 5 is only available as a download in the Supplemental Files section.