Exploring Novel N-Myristoyltransferase Inhibitors: A Molecular Dynamics Simulation Approach.

N-Myristoyltransferase (NMT) is a cytosolic monomeric enzyme involved in the allocation of the myristoyl group to the aminoterminal of glycine in several viral and eukaryotic cellular proteins. NMT has been validated as a potential drug target against kinetoplastid for parasitic protozoa. A multistep virtual screening protocol based on the pharmacophore modeling, molecular docking, and molecular dynamics simulation was carried out. Initially, Maybridge database was virtually screened via a validated pharmacophore model. The effective pharmacophore models were accompanied with exclusion volumes to improve their receiver operating characteristic curve to identify potential NMT inhibitors. The hits identified as actives based on the 3D-pharmacophore model were evaluated by molecular docking studies. In stepwise screening, six compounds were shortlisted for the dynamic simulation to get insights into their binding mode. In conclusion, this study provides fundamental information about the architecture of the binding site and some crucial residues that may provide insights into the development of new antiparasitic agents.

Introduction

Protein modification is of utmost importance in a cellular dilemma, as it is essential for signal transduction, protein trafficking, and a variety of other cellular functions. The fatty acylation of proteins, covalent addition of fatty acid to the proteins, is one of the essential post-translational modification. The process primarily mediated by the enzymes N-myristoyltransferase (NMT) and palmitoyl-acyl transferase (PAT)[1,2] is fundamental for a variety of biological processes, including cell signaling, stabilization of protein–membrane interactions, and vesicular targeting.[2,3]

The enzyme NMT removes initial methionine and facilitates the transfer of myristate moiety to terminal glycine (Figure S1). The NMTs appear to be ubiquitous in eukaryotes. Fungi, protozoa, and insects possess a single NMT, essential for their survival while mammals possess various NMTs.[4] The overall structure of NMT share isomorphism among the orthologues. The essentiality of NMT as an antiviral drug target was established as early as 1985.[5] However, the highly conserved orthologous nature of the Myristoyl binding pocket hampered any practical efforts focused at the development of its inhibitors because of toxicity to humans. The discovery of specificity in the peptide binding sequence of Saccharomyces cerevisiae and humans can be regarded as a cornerstone in the development of NMT inhibitors.[6−8]

NMT has been extensively investigated as a potential drug target in various parasitic protozoa including Trypanosoma brucei, Leishmania major, Leishmania donovani, and Plasmodium falciparum.[9−12] Myristoylation is conclusive for growth and cellular proliferation of Eukaryota.[13] The essentiality of NMT has also been validated in mammal dwelling stages of Trypanosoma cruzi, the vector of Chagas disease.[14] The study has also reported constitutive expression of NMT in all of the stages of T. cruzi life cycle. The perspective is further supported by studies on Toxoplasma gondii, demonstrating the adverse effects on virulence following ablation of a single PAT isoform.[15]

Several hits to lead optimization campaigns have identified specific and potent NMT inhibitors. For instance, in 2014, Brand and fellows have developed a series of brain permeable pyrazole sulphonamide class of inhibitors against T. Brucei.[16] The discovery is based on the optimization of brain permeability of a previously reported pyrazole sulphonamide compound DDD646.[17] The blood–brain permeability of the series was further improved by capping the sulphonamide group and reducing the polar surface area. The most active compound of series DDD100097 demonstrated partial efficacy in stage II mouse model of HAT.[16] Recent efforts based on the hybridization of the previously reported series of inhibitors, resulted in the identification of a novel highly selective series of tbNMT inhibitors.[18] The principal lead compound of the series, 54, presented higher than 1000-fold selectivity for T. brucei NMT (IC50 = 5 nM) as compared to human NMT (IC50 = 230 μM). In another hit to lead campaign, the authors have reported NMT inhibitors with a 40-fold increase in the activity via hybridization of binding modes of two reported inhibitors.[19] The compounds also display modest selectivity against human NMT. A comprehensive review of reported NMT inhibitors can be found in the references of Bell, 2015; Zhao & Ma, 2014.[20,21]

In the present study, we have carried out multistage virtual screening (VS) to identify novel N-myristoyl inhibitors as antiprotozoal agents. In this connection, multimodel pharmacophore modeling was carried out, followed by docking simulation. The obtained hits were clustered using MACCS-based clustering approach. The shortlisted hits were subjected to a short production run of 80 ns. The binding energies of the compounds were extracted from the molecular dynamics (MD) simulation using MM-GBSA. We further evaluated the druggability and ADME properties of the compounds to demonstrate the effectiveness of selected leads.

Results and Discussion

VS has become an indispensable technique in early-stage drug discovery to identify bioactive lead against a particular macromolecule target in a cost-effective and time-efficient manner.[22] VS filters out hundreds and thousands of nonbinders from the available chemical space, reducing the cost associated with the synthesis and bioscreening. In the present study, we aim to implement VS strategies to identify novel NMT inhibitors. In this connection, multistage VS was carried out to identify potential leads from Maybridge database.

Multimodel Pharmacophore Modeling

The prime objective of the pharmacophore modeling process is to develop an “abstract of chemical features” to identify as many active compounds from the dataset of active and as few inactive from the reported inactive as possible. Further, the model should be able to differentiate between the active and inactive compounds with modest accuracy.

The quality of a structure-based pharmacophore model highly depends on the fidelity of the crystal structure and accuracy of the annotation scheme used for the recognition of molecular features.[23,24] In the present study, we used the “unified” annotation scheme implemented in MOE to develop pharmacophore models.

At the time of writing (December 2018), PDB holds evidence for L. major NMT (lmNMT) crystal structures with cognate ligands belonging to nine different classes (Table S1). Brannigan et al. have reported significant diversity in the binding of lmNMT with selective ligands, belonging to four different chemically distinct scaffolds.[25] During a hit expansion campaign, Spinks and group have pointed out that the thiazolidinone derivatives adopt two binding modes with lmNMT.[26] Moreover, the observed binding modes also differ from the previously reported binding modes. Taken together, the evidence suggests that a single receptor cannot entail the diversity in the binding pattern of the reported lmNMT inhibitors.

To overcome this difficulty, we have developed pharmacophore models using 12 different crystal structures (Table S1). For each complex, a structure-based pharmacophore hypothesis was generated in the pharmacophore module in MOE. The initially developed hypotheses were refined to optimize predictive performance. Initially, the active and inactive datasets were screened, and the feature itself or its size was altered to include the active selectively and exclude inactive compounds. The features defining the interactions with key residues such as Phe90, Tyr217, Ser330, and Asn396 were retained.

In the next step, exclusion volumes were added to remove the inactive compounds that sterically clash with the protein. Furthermore, when required the occupied volumes were also added to confer specificity. The models were further evaluated with an aim to select the best models for VS. Among the generated pharmacophore models, the models with optimum performance in the theoretical validations were employed for searching NMT inhibitors. Table presents the results of various statistical metrics calculated in this regard.

Table 1

Details of the Developed Pharmacophore Models: PDB Entries, Number of Active, Inactive, and Decoys Retrieve from the Dataset, Number of Maybridge Compounds Screened by Pharmacophore, and the Values of Sensitivity, Specificity, Precision, and NPV

model no.	1	2	3	4	5	6	7	8	9	10	11	12
PDB ID	2WSA	4A2Z	4A30	4A31	4A32	4A33	4A28	5AG4	5AG5	5AG6	5AG7	5AGE
active recognized, % (number of actives recognized out of 60 actives)	16 (11)	46 (32)	41 (29)	53 (37)	46 (32)	11 (8)	43 (30)	31 (22)	41 (29)	31 (22)	13 (9)	39 (27)
inactive recognized, % (number of inactives recognized out of 30 inactives)	0 (0)	10 (3)	6.7 (2)	3.3 (1)	6.7 (2)	0 (0)	10 (3)	6.7 (2)	3.3 (1)	3.3 (1)	3.3 (1)	0 (0)
decoys recognized, % (number of decoys recognized out of 6264 decoys)	0 (1)	1.4 (88)	1.5 (95)	1.8 (114)	3 (188)	2 (125)	4.2 (267)	5.7 (361)	6.8 (430)	3.3 (208)	4.5 (287)	2.1 (131)
number of VHs from the Maybridge screening collection dataset	21	881	1198	728	1781	150	2610	6363	3583	1890	1782	1947
specificity	1.0	0.9	0.9	1.0	0.9	1.0	0.9	0.9	1.0	1.0	1.0	1.0
sensitivity	0.16	0.46	0.41	0.53	0.46	0.11	0.43	0.31	0.41	0.31	0.13	0.39
precision	1	0.91	0.94	0.97	0.94	1.00	0.91	0.92	0.97	0.96	0.9	1
NPV	0.34	0.42	0.41	0.47	0.42	0.33	0.40	0.37	0.41	0.38	0.32	0.32

The specificity evaluates the ability of a binary classifier to differentiate between the members of active and inactive datasets. In our case, all of the models were found to be highly specific (0.9–1). The sensitivity of a model depicts the percentage of active compounds retrieved by the model from the active database. As mentioned in the Table , the pharmacophore models, 2–5, 7, and 9–10 presented decent sensitivity values (0.41–0.53). In a binary classification system, there exists a compromise between the sensitivity and specificity.[27] Therefore, a specific pharmacophore is not expected to retrieve all of the members of the active dataset.

The precision value describes the ability of a pharmacophore model to fetch true positive and to reject false positive outcome. A precision score closer to 1 depicts the bias towards the actual positive outcomes. A higher precision: negative predictive value (NPV) ratio pronounces the ability of the model to identify a positive compound accurately and to reject specific negative instance. In conclusion, all of the developed pharmacophores present significant bias toward the members of the active dataset. However, for VS the most restrictive models with nonoverlapping hit spectra were selected. As demonstrated earlier, using a combination of multiple restrictive models representing different binding modes rather than a single broad model restricts the retrieval of false positive hits.[28]

The hits obtained from all of the models were mapped to evaluate the hit spectra. In case of models with overlapping hit spectra, only the most restrictive model was retained. As a result, three pharmacophore models 2, 4, and 10 (PDB IDs; 4A2Z, 4A31, and 5AG6) were selected for VS of Maybridge database. The finally selected models were able to identify all of the active compounds. Altogether, the three models recognized three inactives out of 30 compounds in the dataset. Further, each of the finally selected pharmacophore models had a decoy hit-rate of less than 5%, which pronounce the specificity of our models. Herein, we will restrict our discussion to the finally selected pharmacophore models.

The pharmacophore model 2 was developed using the protein–ligand interaction pattern of PDB 4A2Z. The finally selected pharmacophore model consisted of four features: two aromatic, a donor, and an acceptor (Figure a,b). The fundamental pharmacophore interactions of the finally selected pharmacophore included interactions with Phe88 and Phe90 and Tyr217 as well as two hydrogen bonds with Ser330 and Asn376. The addition of another feature atom-Q at the sulfur atom of the sulfonamide moiety enhanced the specificity of the model, resulting in a more restricted model. Further, both occluded and exclusion volumes were added to improve the quality of the model.

Figure 1

Finally selected structure-based pharmacophore models 2, 5, and 10 based on the crystal structures of lmNMT; PDB IDs 4A2Z, 4A31, and 5AG6, respectively. Panel (a) depicts the features of pharmacophore model based on 4A2Z (model 2). The features of 4A31 and 5AG6 (model 4 and 10) are depicted in subfigures (c,e). In the right panels, subfigures (b,d,f) show the flattened structure of cognate inhibitor and limited definition of the surrounding residues. The khaki spheres show hydrophobic or aromatic, blue, and magenta spheres depicts hydrogen bond donor and acceptor, respectively. The gray sphere exhibit exclusion volumes, while translucent red dotted spheres in figure (a) exhibit the occupied volumes. The magenta dashed lines in (a,c,e) depict hydrogen bonds. The water-mediated bridges are presented as blue dashes.

Model 4 consisted of six features with partial match enabled for one feature. The selected acceptor and donor features mapped to hydrogen bonding interaction with Ser330, Asn376, and a water-mediated bridge with Tyr80 (Figure c,d). The two aromatic features depict π-stacking interaction with the aromatic amino acids in serine and peptide binding pockets. The pharmacophore model retrieved 32 active compounds and 2 inactive compounds.

The pharmacophore model 10 consisted of four features complementary to the architecture of the binding site. An aromatic feature depicts the interaction of inhibitor with Phe232 and Tyr217. The model also contains two hydrogen bond features depicting interactions with Asn376 and Asp396 (Figure e,f).

A comparative overview of interaction profiles of cognate ligands in case of pharmacophore models 2, 4, and 10 is tabulated as Table . The nitrogen atom of the pyrazole ring in case of models 2 and 4 and oxygen substitution at the thiazole ring of thiazolidinone in model 10 mediate hydrogen bonding interaction with Ser330. Similarly, all of the three models mediate hydrogen bond with Asn376. In model 2, the ligand mediated multiple interactions with Tyr217, including hydrogen bond, while the other models emphasize on the hydrophobic interactions with Tyr217. In model 10, multiple interactions were observed between the ligand and Phe88, whereas the other two ligands do not interact with Phe88. Thus, the three models provide a slightly different definition of the active site, which might incorporate diversity in the obtained virtual hits (VHs). To further investigate the binding mode of the obtained VHs, molecular docking simulations were carried out. Figure provides a sequential summary of the multistep VS carried out in the present study.

Table 2

The Comparative Interaction Pattern of the Models; 2, 4, and 10a

residue	4A2Z (viq)	4A31 (2cb)	5AG6 (5pe)
Tyr80		Hw
Asn167		Hw
Asn376	Hb	Hb	Hb
Asp396			Hb
Asp83	Hd
Gly205		Hb
Gly397			Hb
His219			Hw
Leu341	Hd	Hd	Hd
Phe232	Hd	Hd	Hd
Phe88			Hm
Ser330	Hb	Hb	Hb
Tyr217	Hd, Hb, π-stacking	Hb, π-stacking	Hd, π-stacking
Tyr345		Hw

Hd = hydrophobic, Hb = hydrogen bond, Hw = water-mediated bridge, Hm = multiple interactions.

Figure 2

Overview of the multistage VS employed in the present study.

Molecular Docking

The binding modes of compounds screened by pharmacophore were developed using molecular docking studies. Molecular docking is a fast and reliable method for large-scale drug screening. The accuracy of a docking program to be used in a particular VS campaign can be evaluated in two ways. First, the software should be able to reproduce the near-native; that is, the experimentally determined pose and secondly the scoring function should be able to differentiate between the active and inactive compounds. To evaluate the accuracy of docking software and to establish the binding mode of novel VHs, benchmarking studies were carried out. The experiments revealed that the top-ranked poses generated by Vina were similar to the native binding mode (Figure S2). Further experiments with multiple receptors indicate that Vina in combination with at least three different structures, 2WSA, 4A2Z, and 4A31, was able to score active and inactive compounds differently. In the case of PDB: 2WSA, the ROC lies in the left upper quadrant of the graph pronouncing the high specificity of the protocol (Figure ). The observation is further supported by the AUC value closer to 1 (0.98). In the light of the above experiments, the binding modes of the shortlisted hits were established using Vina and PDB: 2WSA as a receptor.

Figure 3

Results of score-based ranking experiments carried out using Vina and crystal structures of NMT; mentioned below as four-letter PDB codes. The data consisted of the members of active and inactive datasets. The AUC values and ROC curves were obtained using ROCR function in R (Sing et al. 2005; R Development Core Team 2008). The higher AUC values in case of 2WSA, 4A2Z, and 4A31 depict a good binary classifier.

The VHs screened during the pharmacophore screening were subjected to molecular docking studies. The top-ranked pose of the shortlisted VHs was extracted using an in-house script, and the hits with scores higher than −7.7 kcal/mol were included in the further calculations. The selection cut off was based on the results of the score-based ranking experiments carried out earlier.

The compounds were further clustered according to the nature of the core structure, using the MACCS structural keys fingerprinting approach implemented in SAR report utility in MOE. As a result, 13 large clusters were obtained (Figure ). A cursory look at Figure shows that our pharmacophore models were able to identify diverse chemotypes, pronouncing the scaffold hopping ability of selected pharmacophores. The selected hits were subjected to protein–ligand interaction fingerprinting (PLIF) implemented in MOE. The interaction pattern of each VH was analyzed. Based on the interaction with hotspot residues, 12 hits were shortlisted. The chemical structures of the selected 12 compounds are presented in Figure . The IUPAC names, three-dimensional coordinates, and docking scores are present in the Supporting Information; Table S2 and S2.zip.

Figure 4

Substructures of the most prominent clusters identified in selected VS hit.

Figure 5

Chemical Structures of the finally selected hits identified during multistage VS against NMT along with their Maybridge product codes.

Analysis of the Binding Mode

The binding mode of selected hits was established using Vina. Figure a depicts the simulated binding pose of the compound HTS05346. The methyl groups at the terminal benzene ring exhibit a “methyl effect”. However, Tyr217 is stacked against benzene ring, exhibiting face-to-face interaction. The oxygen of isoxazole ring demonstrates a hydrogen bond with a backbone amide group of His398. The protein–ligand contacts are further stabilized by apolar interaction between the ligand and Glu82.

Figure 6

Binding poses of VHs: HTS05346, SPB03259, CD09798, HTS11181, SPB08283, and HTS05445 depicted in panels (a–f) in alphabetical order. The protein residues are depicted as thick gray sticks (carbon atoms colored in light grey). For the clarity of the picture, only residues interacting with ligand were included in the graphic. The dashed green line depicts hydrogen bonding interaction. Images were rendered using UCSF Chimera.

In case of compound SPB03259, the chlorobenzene ring exhibits a parallel π-stacking interaction with the aromatic ring of Tyr 217, as well as a hydrogen bond (distance 2.97 Å). However, the terminal isoxazole ring captures the hydroxyl group of Tyr80 (Figure b). The chlorobenzene moiety of CD09798 interacts with Tyr217, with additional contact mediated by hydrogen bond interaction via a side chain hydroxyl group. The nitrogen atoms in the pyrazole ring of the ligand demonstrate a hydrogen bond with Thr203 and a water bridge with Asn167 (Figure c).

The compound HTS11181 exhibit hydrogen bond with the backbone amide of His398 (Figure d). Another hydrogen bond is observed between the ligand and Glu82. Further anchorage is provided by face-to-face interaction between the ligand and Tyr217. In case of SPB08284, the heteroatoms (oxygen and nitrogen) of 5-methylisoxazole demonstrate hydrogen bonding interaction with backbone amide groups of the Asn396 and Gly397 (Figure e).

Figure f presents the top-ranked simulated pose of compound HTS05445. As evident from figure, the edges of the benzene ring of the ligand exhibits weak interactions with Tyr217. The chlorobenzene ring in ligand delves deeper into the ligand binding site, producing a “hydrophobic effect” which in the absence of polar interactions act as a compensatory stabilizing factor. The compound also exhibits hydrogen bonds with Asp83 and Asn376.

Figure a presents the binding mode of compound AW00765. The compound exhibits hydrophobic contacts with Val81 and Phe90. The methyl group at the nicotinic acid establishes a hydrophobic contact with edge of Tyr217. There is an apparent dearth of hydrogen bonding with the ligand mediating hydrogen bond interaction with only Asp83.

Figure 7

Binding poses of VHs: AW00765, SPB03712, HTS11453, RH06193, BTB05308, and KM04091 depicted in panels (a–f) in alphabetical order. The protein residues are depicted as thick gray sticks (carbon atoms coloured in light grey). For the clarity of picture, only residues interacting with ligand were included in the graphic. The dashed green line depicts hydrogen bonding interaction. Images were rendered using UCSF Chimera.

The compound SPB03712 mediates hydrophobic contacts with Phe90, Leu341, Ala343, and Leu399. The thiazole ring of compound is sandwiched between Tyr217 and Phe232. Apart from hydrogen bonding interaction with Asn376, the chloride moiety is involved in halogen bonding interaction with Asn167 (Figure b). The plausible binding mode of compound HTS11453 is depicted in Figure c. The dichlorobenzene moiety of the ligand exhibits parallel π-stacking interaction with Tyr217. One of the chloride substitutions mediates a halogen bond with Tyr345 side chain hydroxyl group. The backbone amide of Asp396 interacts with oxygen substitution at the purine ring.

The compound RH06193 resides comfortably in the binding site extending hydrophobic interactions with Val81, Phe90, and Tyr345 and hydrogen bond with Asn376 (3.4 Å). The aromatic ring of the Tyr217 mediates face-to-face stacking interaction against the benzene ring of ligand (Figure d). The analysis of the top-ranked binding mode of the compound BTB05308 highlights hydrophobic contacts as the most dominant interaction between the ligand and protein. As evident from the Figure e, the ligand presents hydrophobic contacts with Val81, Phe88, Phe90, Leu341, and Leu399. The aromatic rings of the Phe90 and Tyr217 are positioned perpendicular to the ligand, reinforcing a T-shaped π-stacking interaction. The polar atoms of the ligand mediated hydrogen bonds with the side chains of the Tyr80 and Asn376. Moreover, a water-mediated bridge connects the ligand with Asn167.

Figure f depicts the binding mode of compound KM04091. The compound also exhibits hydrophobic contacts with the Val81, Phe88, Phe90, and Phe232. An essential feature of the ligand is water-mediated interaction with Gly397. The ligand 646 mediates similar interaction with Asn396 (PDB: 2WSA). Nonetheless, the compound 646 exhibits water-mediated interaction with Gly205 via a different moiety. The simulated binding pose of KM04091 demonstrates a direct hydrogen bond between Gly205 and ligand. The observation highlights the scaffold hopping ability of our protocol.

ADME-ToX Prediction

The selected hits were subjected to ADME-ToX prediction using the QikProp module in Schrödinger. The molecular weights and volumes of all of the 12 compounds were found to comply with the ranges established using the 95% of the reported drugs (Tables and S3). All of the compounds were found to observe the Lipinski rule of five of drug likeliness[29] with zero violations, except for SPB08284 which was found to violate one of the filters. As presented in Table most of the compounds were also in accord with the Jorgensen rule of three. All of the compounds were predicted to have moderate to low blood brain barrier permeability with high oral absorption rate except SPB03712. The observation is further supported by the good Caco-2 permeability of the selected compounds. Caco-2 is often employed as an in vitro model of human intestinal mucosa for prediction of the oral bioavailability of a drug.[30]

Table 3

Predicted ADME Properties of the Selected Hits

molecule	CNSa	mol_MWb	volumec	QPlogPo/wd	QPlogHERGe	QPPCacof	QPlogBBg	percent human oral absorptionh	rule of fivei	rule of threej
aw00765	–1	364.341	1007.017	4.435	–3.718	293.05	–0.077	100	0	1
btb05308	–2	361.313	1053.058	1.925	–6.377	86.882	–1.783	72.922	0	0
cd09798	0	260.682	801.816	2.423	–4.898	791.38	–0.311	93.01	0	0
hts05346	0	321.378	1076.456	3.454	–5.927	1259.41	–0.569	100	0	0
hts05445	0	323.741	956.885	3.077	–5.688	497.446	–0.539	93.23	0	0
hts11181	–1	388.428	1208.834	3.407	–6.096	957.396	–0.85	100	0	0
hts011453	0	433.296	1216.372	4.264	–4.892	1278.487	–0.173	100	0	1
km04091	–1	450.386	1140.141	2.819	–5.01	695.371	–0.545	94.322	0	0
spb03219	0	318.716	970.456	2.683	–5.385	831.922	–0.495	94.917	0	0
spb03712	–2	459.924	1311.448	3.923	–6.572	266.051	–1.106	93.319	0	1
spb08284	1	447.335	1208.13	5.543	–5.098	3446.525	0.261	100	1	1

Predicted CNS activity (−2/+2).

Molecular weight (130.0/725.0).

Molecular volume (A3 500.0/2000.0).

log octanol/water partition (−2.0/6.5).

HERG K+ channel blockage: log IC50 (concern below −5).

Apparent Caco-2 permeability (nm/s <25 poor, >500 great).

log brain/blood permeability (−3.0/1.2).

Percentage human oral absorption in GI (<25% is poor).

Lipinski rule of 5 violations (maximum is 4).

Jorgensen rule of 3 violations (maximum is 3).

MD Simulation

In order to understand the dynamic behavior of ligands, a short production run of selected hits was carried out using AMBER. The protein–ligand complex was solvated and equilibrated, followed by the production run. The MD simulation provided insights into the overall stability of system. The quality of the simulation was evaluated by analyzing the root mean square deviation (RMSD), radius of gyration (RoG), and root mean square fluctuation (RMSF).

RMSD serves as a key factor to describe the stability of a system undergoing MD simulation. The internal motion of protein backbone atoms is analyzed to investigate the stability of the system throughout simulation. The RMSD of the main chain heavy atoms of system of protein was obtained using the “rms” command in CPPTRAJ. Figure depicts the RMSD deviation of heavy atoms of protein backbone. The assembly was not fully equilibrated until 10 ns into the production run (as determined by the protein RMSD); therefore, only the remaining 70 ns was considered for following calculations.

Figure 8

RMSD of the systems calculated as a function of time.

As evident from Figure , all of the systems displayed variable deviations with an average RMSD value of 2.5 Å. All of the systems presented stable internal motion. The observation was further supported by analysis of RoG, which suggests that all of the systems were well-converged (Figure S3). The residue contact information was retrieved to observe the immediate environment of each ligand. To assess the dynamics of side chains of residues, RMSF of the protein residues was calculated as a function of time. The data indicate that the amino acid residues in all of the protein–ligand complexes were fairly stable when compared to the apo protein (Figure ).

Figure 9

RMSF of the protein residues calculated as a function of time.

The visual analysis of trajectories suggests that HTS11453 mediates hydrophobic and π-stacking interaction with Tyr217 and His219. Further, methyl substitutions at the pyrazole ring are near Phe88, Phe90, and Phe232 which validate the previously hypothesized methyl effect. Moreover, the other methyl group extends weak interactions with aliphatic side chains of Leu341, Ala343, and Val374 (Figure S4). Taken together, these pieces of evidence suggest the critical role of methyl groups in lowering the solvation energy. During simulation, the benzene ring of ligand was found to be sandwiched between the Tyr217, aromatic group and His219 (imidazole) consistent with the docking pose.

The analysis of time-dependent trajectories of SPB08284 demonstrated some interesting results. The polar atoms of the isoxazole moiety remained in contact with the receptor via polar interaction with Gly397 and Asp396. The minimum distance between isoxazole polar atoms and amide groups of interacting residues is depicted in Figure S4. Another interesting observation was found in case of the 1,2-dimethoxybenzene moiety. After 16 ns into the production (8 ns of the reported production run used for the analysis), ring rotates (approx. 10°) to expose the other methoxy group to Ser330. However, the minimum distance between the Ser330 hydroxyl group and ligand remained the same.

During simulation, the benzene ring of compound AW00765 penetrated more in-depth into the binding side extending interaction with Val81, Tyr217, Tyr345, and Leu399. It has been reported that these aromatic–aliphatic interactions are amongst the most commonly occurring interactions at the nonpolar interfaces.[31] After the removal of the first restraint, fluorine atoms were reoriented such that one of the fluorine atoms remained in contact with the amide backbone of Gly205 (Figure S4). Of particular importance is the role of Tyr217 in the stability of compound. The side chain hydroxyl group of Tyr217 maintained polar contact with ligand supporting the hydrophobic contact between ligand and its aromatic ring.

Binding Free Energy

The MM(GB/PB)SA free-energy contributions were calculated to get insights into the forces involved in protein–ligand binding. To evaluate the binding ability of all three complexes, 100 conformations were extracted from last 10 ns of MD trajectories with 0.1 ns time interval. Table enlists the contribution of individual energetic components to the binding of the inhibitor. The calculated binding free energies (ΔGb) of HTS11181, SPB08284, and AW00765 complexes with NMT are −51.16, −38.99, and −63.82 kcal mol–1, respectively. The data indicate that the intermolecular van der Waals interactions (ΔEint vdW) contribute significantly to the binding free energies. However, the electrostatic interaction energies (ΔGele) were positive in case of complex AW00765 and the data indicate favorable contributions to the binding free energies, which are in good agreement with MD result (electrostatic interaction of ligand with Gly205 perturbs after 40 ns; showed in md graphs). Further, the nonpolar salvation terms, corresponding to the burial SASA upon inhibitor binding, contribute slightly favorably. However, the total binding affinity arises from a more complex interplay between all these components.

Table 4

Results of Binding free Energy Calculations of Selected Hits

	HTS11181	SPB08284	AW00765
van der Waals	–51.1676	–38.9934	–63.8209
EEL	–233.6963	238.1473	–7.7965
EPB/EGB	252.866	–219.888	23.2525
ESURF/ECAVITY/ENPOLAR	–5.0245	–3.8207	–6.0056
ΔG binding (GB)	–37.0223	–24.5548	–54.3706
ΔG binding (PB)	–5.148	6.4545	–15.2767

Conclusions

Myristoylation, a post-translational modification process, is essential for the growth and survival of protozoa. In the present study, we have carried out multistage VS to identify novel NMT inhibitors as antiprotozoal agents. In this connection, multimodel pharmacophore screening was carried out using the available crystal structures. The subsequent optimization and statistical evaluation highlight the ability of generated hypotheses to identify chemical entities with reliable sensitivity and specificity (Table and Figure ). The compounds screened from the pharmacophore-based VS were subjected to docking simulation for binding mode analysis. The benchmarking results suggest that the global search and scoring functions of AutoDock Vina were able to reproduce the experimental pose and differentiate between the active and inactive compounds. The analysis of molecular docking and subsequent analysis identified 12 significant hits (Figure ). Finally, selected compounds from the obtained hits were subjected to a short production run of 80 ns. The analysis of the quality metrics suggests that all of the systems remain stable throughout simulation (Figures and 9).

Methodology

Protein Preparation

The crystal structures of the lmNMT bound to different classes of inhibitors with high resolution (<3 Å) were retrieved from the RCSB Protein Data Bank (Table S1).[32] The structures were prepared using the protein preparation wizard in MOE v. 2016.0801 [Molecular Operating Environment (MOE), 2016.08; Chemical Computing Group Inc., 1010 Sherbrooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7, 2016].

During preparation, the structures were subjected to addition of missing atoms and residues, bond order, formal charge correction, and tautomer adjustment. Hydrogens were added using the protonate 3D algorithm. For protonation, generalized Born volume integral[33] (GB/VI) was used as the electrostatics function, with a dielectric value of 80 (for solvent). The van der Waals and electrostatics cut off were set as 15 and 10 Å, respectively. Partial charges were applied using AMBER10:EHT force field[34,35] implemented in a MOE software suite. The water molecules involved in bridging ligand and protein were retained, which were kept rigid during energy minimization.

Pharmacophore Modeling

Pharmacophore modeling is a sophisticated computational chemistry method used for the elucidation of intermolecular interactions stabilizing the protein–inhibitor complexes.[36] Pharmacophore-based screening allows rapid screening of millions of compounds, without compromising the overall accuracy. We used the unified annotation scheme for pharmacophore modeling implemented in MOE, as it is the most comprehensive annotation scheme. A series of different pharmacophore models were developed using structures mentioned in Table S1. Each pharmacophore model was refined to increase its specificity by manipulating the exclusion volumes. The partial matches were also allowed in some instances to improve the quality of model.

Data Set Preparation

In order to evaluate the selectivity of developed protocol, a data set of structurally diverse compounds with potent NMT activity was manually curated. Similarly, a data set of inactive compounds was also created using structurally similar inactive compounds. To access the efficiency of generated pharmacophore hypotheses, we developed a decoy library using DUD-E web server.[37]

For each active ligand, 50 decoys with similar physical properties but different topological characteristics were created. The inclusion of the decoy data set further challenges the specificity of the hypothesis. The prime objective behind using the decoy data set is to avoid the artificial enrichment.[37] The decoys were prepared according to the protocol established earlier.[38]

For VS, Maybridge Screening Collection was downloaded from the Maybridge website (https://www.maybridge.com). As of December, 2018, Maybridge screening collection consisted of 52 066 compounds. All of the compounds were charged using the MMFF94x force field implemented in MOE, allowing the adjustment of hydrogens and lone pairs, when required. The database was subjected to energy minimization using default RMS gradient of 0.1 kcal/mol/Å2. Following protonation and minimization, the compounds were saved in the MOE database format.

Theoretical Validation

For theoretical validation, the developed pharmacophore models were subjected to three different subsets; active, inactive, and decoys. A modest pharmacophore should be able to differentiate between the active and inactive data set, identify the members of active data set, and reject the subset of inactive compounds. In the present study, we have calculated different statistical metrics to validate our pharmacophore models. The details of different statistical metrics used for the validation are mentioned below.where TP = true positive, TN = true negative, FP = false positive, and FN = false negative.

Molecular docking simulation is commonly employed in structure-based VS campaigns to explore the behavior of small molecules in binding site of target protein.[39] Molecular docking studies were carried out to identify the binding mode and binding scores of compounds shortlisted by the pharmacophore hypotheses. The selection of an optimal docking software with suitable sampling and scoring algorithm is instrumental in successful VS campaigns. The accuracy of a docking software depends on the ability of the software to quantitatively reproduce the crystal pose and to differentiate between the active and inactive compounds. In the present study, we have evaluated the performance of three different docking softwares; MOE [Molecular Operating Environment (MOE), 2018.01; Chemical Computing Group ULC, 1010 Sherbrooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7, 2018.], FRED[40] and AutoDock Vina[41] (hereafter referred to as Vina). The ability of the software to produce the crystal pose was assessed using the redocking experiments, whereas the ability of the software to differentiate between the active and inactive compounds were assessed by score-based ranking.

Redocking

To access the ability of the software to reproduce the crystal pose, we subjected the cognate ligands to docking simulation. For Vina, the receptor and ligand were prepared using the Python-based prepare_receptor and prepare_ligand scripts. The grid was prepared centering at the ligand. The search spaces and grid coordinates of each complex are tabulated in Table S4. The global search exhaustiveness was set to the default value (8), and a total of 10 poses were generated for each complex.

In case of MOE, the docking site was defined using the coordinates of cognate ligands. We used “Triangle Matcher” and London dG as placement and scoring methods, respectively. The rigid receptor scheme was used for refinement of top 30 poses generated by placement algorithm, followed by scoring using the GBVi/WSA dG scoring method. In every instance, 10 poses were recorded.

Fast Rigid Exhaustive Docking (FRED) is implemented in OpenEye software suite.[40] FRED performs ultrafast docking by symmetrically searching all rotations and translations in ligand using an exhaustive search space in receptor. The receptors were converted into the OpenEye binary format (oeb) using the prepare_receptor application in OpenEye. A conformation database of the cognate ligand was generated using Omega2 (OpenEye Software, SantaFe, NM, USA). We used ChemGauss4 scoring function and the hit list size was set to 150. For each complex, a total of 20 poses were recorded. The results were analyzed using VIDA.[42]

The quantitative measurement of difference in the cartesian coordinates of simulated pose and experimental pose (RMSD) serve as gold standard to evaluate the reliability of a docking experiment.[43] For a satisfactory redocking experiment, the RMSD value between the cartesian coordinates of simulated and crystal pose should be less than 3 Å.

Bench-Marking and Scoring

Apart from the regeneration of crystal pose, a docking protocol should be able to differentiate between the active and inactive compounds. Herein, we challenged the ability of Vina to differentiate between the active and inactive compounds. In this connection, the available lmNMT crystal structure was utilized (Table S1). The data sets of both active and inactive compounds were subjected to molecular docking studies using Vina. The predicted binding free energy of the top-ranked pose of every compound was retrieved in all instances. For statistical analysis of results, a receiver operator curve was plotted using ROCR package in R.[44,45]

Screening of Maybridge Database

The hits obtained from the pharmacophore-based screening of the Maybridge database were subjected to molecular docking using Vina using the protocol mentioned in the redocking section. The prepare_ligand script was used for the format conversion (from mol2 to PDBQT). For each compound, 10 poses were generated. However, only the top-ranked pose for further calculations.

Protein–Ligand Interaction Fingerprinting

The PLIF approach converts the molecular interaction data into the bit array conforming to the residue ID and interaction type. In the present study, we calculated the interaction fingerprints of the top-ranked pose of the VHs obtained from Vina. For format conversion from PDBQT to mol2, open babel v. 2.4.0 was used.[46] The molecules were then imported into MOE as MDB, and PLIF was generated using only a weak cut off value of 0.5 kcal/mol for each type of interaction. The compounds were shortlisted on the basis of interaction with the crucial residues; Tyr80, Phe90, Thr203, Tyr217, His219, Phe232, Ser330, Asn376, Asp396, and Leu421. The selected hits were further subjected to MD simulation.

The QikProp module implemented in Schrödinger was used to calculate the ADME-ToX properties of the selected compounds. The compounds were imported in Maestro workspace in mol2 format and were subjected to QikProp calculations in the fast mood. For each molecule, top five similar drugs were also evaluated. The results were exported in CSV format for further processing.

MD simulations were carried out using the pmemd.CUDA module in Amber18 using explicit water model under periodic boundary conditions.[47] The force field ff14SB force field was used for the parametrization of protein.[48] The nonpolar hydrogen atoms of protein were trimmed using the reduce module in AMBER followed by subsequent protonation via tleap.[49] The cofactor Myristoyl-CoA was parameterized using the second generation GAFF2 force field.[50] However, the atomic partial charges of the ligand were derived by AM1-BCC charge calculation formalism using Antechamber module in Amber14.[51−53] The parmchk2 utility in Antechamber was employed to generate missing dihedral parameters. The system was then solvated with the explicit TIP3P water molecules into a periodic box that extended at a distance of 10 Å from the protein atoms.[54]

To ensure electroneutrality of the system, counter ions were introduced into the solvated system, when required. The topology, coordinates, and the AMBER parameters were generated using the tleap module in AMBERTOOLS-14. The box dimensions, number of ions required for neutralization, number of water molecules added, and the total number of atoms in the systems are tabulated as Table S5.

Energy Minimization, Equilibration, and MD Production Run

All of the complexes were relaxed before MD simulations to remove the bad contacts and possible steric clashes. By performing a short energy minimization, the system was relaxed before MD simulations, to remove initial bad contacts and steric clashes by performing energy minimization. Briefly, energy minimization was carried out to relax the system by decreasing the harmonic force restraints to the solute atoms (protein–ligand-co-factor complex) involving the steepest descent and conjugate gradient methods.[55] The system was equilibrated in a stepwise manner, by gradually heating the restrained system, from 0 to 300 K, in an NVT ensemble. The system was then relaxed gradually in an NPT ensemble (pressure = 1 atm). During relaxation, a harmonic restraint of 25 kcal/mol/Å2 was applied on the solute atoms for 25 000 steps, which was released gradually.

The production run was carried out for 100 ns using the isothermal–isobaric (NPT) ensemble at 300 K temperature and 1 atm pressure with periodic boundary conditions. The temperature was controlled using the Berendsen weak-coupling algorithm.[56]

The particle mesh Ewald method[57] integrated into AMBER was utilized to simulate the long-range electrostatic interactions. The cutoff values for the nonbonded (electrostatics and van der Waals) interactions were set to 8 Å, and the dielectric constant was taken as 1. The integration step (dt) was set to 2 fs to integrate the Newtonian equation of motion. The bond involving the hydrogen atoms were constrained using the SHAKE algorithm.[58]

The mass-weighted RMSD, RMSF, and RoG calculations were performed using CPPTRAJ[59] module in AMBERTOOLS18. The distances between the heavy atoms were calculated using the CPPTRAJ. All of the graphs were generated using gnuplot.[60]

Molecular Mechanics/Generalized Born Surface Area Calculations

The molecular mechanics/generalized Born surface area (MM/GBSA) calculations were performed to assess the electrostatic contribution of the solvation free energy using generalized Born (GB) model implemented in AMBER18.[61] For solute protein (εin), the dielectric constant value was set to 1 and the exterior was set to 80. The GB model developed by Onufriev et al. (GBOBC1) with igb = 2 and mbondi2 parameterized radii implemented in Amber software package were used.[62] The nonpolar contributions (ΔG nonpolar) were calculated with the linear combinations of the pairwise overlaps method.[63] The surface tension (γ) and offset (β) values were set to 0.00542 kcal/(mol Å2) and 0.92 kcal/mol, respectively.