Analyzing a Potential Drug Target N-Myristoyltransferase of Plasmodium falciparum Through In Silico Approaches.

BACKGROUND: Despite concerted global efforts to combat malaria, malaria elimination is still a remote dream. Fast evolution rate of malarial parasite along with its ability to respond quickly to any drug resulting in partial or complete resistance has been a cause of concern among researcher communities.
MATERIALS AND METHODS: Molecular modeling approach was adopted to gain insight about the structure and various analyses were performed. Modeller 9v3, Protparam, Protscale, MEME, NAMD and other tools were employed for this study. PROCHECK and other tools were used for stereo-chemical quality evaluation. RESULTS AND
CONCLUSION: It was observed during the course of study that this protein contains 32.2% of aliphatic amino acids among which Leucine (9.5%) is predominant. Theoretical pI of 8.39 identified the protein as basic in nature and most of the amino acids present in N-Myristoyltransferase are hydrophobic (46.1%). Secondary structure analysis shows predominance of alpha helices and random coils. Motif analyses revealed that this target protein contains 2 signature motifs, i.e., EVNFLCVHK and KFGEGDG. Apart from motif search, three-dimensional model was generated and validated and the stereo-chemical quality check confirmed that 97.7% amino acid residues fall in the core region of Ramachandran plot. Molecular dynamics simulation resulted in maximum 1.3 Å Root Mean Square Deviation (RMSD) between the initial structure and the trajectories obtained later on. The template and the target molecule has shown 1.5 Å RMSD for the C alpha trace. A docking study was also conducted with various ligand molecules among which specific benzofuran compounds turned out to be effective. This derived information will help in designing new inhibitor molecules for this target protein as well in better understanding the parasite protein.

INTRODUCTION

Every year, millions of people pay the cost of residing in malaria-endemic regions of the world with their life. It is speculated that 40% of population face the risk of malaria in 109 countries across the globe.[1] Most vulnerable population figures among the most downtrodden and is further weakened by poverty, malnutrition, unhygienic conditions, and above all with limited or no access to healthcare and life-saving drugs. Malaria now claims more life than ever before and is tightening its grip in newer regions. In 2006, 247 million cases were reported and 3.3 million live under the shadow of malaria.[1] A glance at malaria figures shows clearly that it is the biggest killer of children below 5 years of age and the number is bound to increase in next few years according to a recent estimate.[1]

Malaria is caused by four species of protozoan parasite Plasmodium vivax, Plasmodium falciparum, Plasmodium ovale, and Plasmodium malariae transmitted by vector Anopheles mosquitoes. Plasmodium falciparum has emerged as one of the most successful and dreaded parasite of our times. Hopes raised by initial success of insecticide Dichloro Diphenyl Trichloroethane (DDT) have been marred by increasing resistance of vectors to insecticides and drug resistance in parasite. Emergence of multi–drug resistant parasite has been the cause of failure of our efforts to control malaria. Efficacy of many conventional anti-malarial drugs has been compromised.[2-3] Malaria control still remains a far fetched dream and big claims of malaria elimination have been a major letdown despite the enormous funding for research,[4] efforts by national as well as international agencies and advances in basic and applied sciences. Widespread resistance of Plasmodium to chloroquine and sulfadoxine/pyrimethamine poses a major threat to malaria control programs.[5] Failure of conventional drugs for treatment of malaria has warranted the need for exploring new drug targets. Malaria drug pipeline is nearly empty due to lack of interest from big pharmaceutical giants owing to low profit margins. Few anti-malarial drugs that were registered during past several decades can be counted on fingers.[6] The woefully small number of available anti-malarial drugs in arsenal for battling malaria is a cause of concern. In order to replenish the drug pipeline, efforts for identification of novel chemotherapeutic targets must be intensified.

Myristoyl-CoA protein N-myristoyltransferase (NMT: EC 2.1.3.97) is a ubiquitous cytosolic[7] enzyme, which follows a bi–bi catalytic reaction mechanism to bring about co-translational transfer of the rare cellular fatty acid myristate (C14:0) from myristoyl CoA to the N-terminal glycine residue in a variety of eukaryotic cellular protein substrates.[8-10] N-terminal N-myristoylation results in significant change in important properties like lipophilicity of the protein and thus facilitates interactions of protein with hydrophobic domains and membranes.[11-15] NMT of P. falciparum is known to express in its asexual blood stage and has been cloned.[16] NMT of P. falciparum is an attractive drug target as established by comparative biochemical studies involving human and parasite enzyme.[17] Difference in properties can be exploited for development of specific P. falciparum NMT inhibitors. So, the logical step for characterization of protein is making use of available bioinformatics tools and techniques for determining important aspects of the enzyme. Crystal structure of P. falciparum NMT has not been solved yet. Comparative modeling provides a way for obtaining structural information in absence of experimentally derived structure. We have undertaken this study to characterize P. falciparum. NMT in silico as knowledge of key properties and structural features of NMT will aid in development of new chemotherapeutic agents for treatment of malaria.

MATERIALS AND METHODS

Amino acid sequence of P. falciparum NMT (Accession Number: AAF18461.1) was obtained from the protein database of NCBI. In silico characterization of physiochemical properties of the protein was performed using Protparam[18] and Protscale available at EXPASY(http://www.expasy.ch/). Secondary structure predictions were made using multi-prediction server at Network Protein Sequence Analysis at PBIL. CONSEQ was used for determining important residues conserved during the course of evolution.[19] Multiple Em for Motif Elicitation (MEME) available at (http://meme.nbcr.net/meme4_1_1/cgi-bin/meme.cgi) was used to predict motifs in the PfNMT keeping default options (minimum width: 6, Maximum width: 50, Motifs to find: 3 and minimum sites >=2). For identifying potential domains, PROSCAN was used. Quasi motifinder was employed to predict signatures and motif like patterns keeping default cut-off P value for motif similarity of 0.05 and number of possible pseudomotifs predicted was kept 3.

CYSREC (http://linux1.softberry.com/) was used to predict the cysteine pairing pattern. As some proteins are known to contain several unstructured or disordered regions, metaProtein disorder prediction system (metaPrDOS)[20] was used for prediction of disorder regions applying PrDOS, DISOPRED2, DisEMBL, DISPROT (VSL2P), DISpro and IUPred predictors. Transmembrane regions were predicted by DAS,[21] TMPRED.[22] SOSUI,[23] HMMTOP,[24] TMHMM,[25] and SPLIT.[26] A search was made using target protein sequence as query against PDB database by Position Specific Iterative (PSI) BLAST[27] keeping default parameters so as to find homologous structure, which can serve as template for theoretical model construction. Sequence alignment between target and template sequence was created using CLUSTALX,[28] a multiple sequence alignment program.

Homology model was constructed using a computer program MODELLER 9v3,[29] which relies on satisfaction of spatial restraints. MODELLER works by accepting target–template alignment as an input and employs a series of model building steps for construction of the model. Obtained model was further refined and subjected to molecular dynamics simulations using strategy followed in our earlier work.[3031]

Structural diagram of the model was prepared and viewed using VMD.[32] The model was subjected to various tests for assessing the quality. For checking the consistency as well as validity of the model, stereochemical evaluation was performed using PROCHECK,[33] VERIFY3D,[34] and ERRAT[35] server. As low folding energy indicates stability of the model, energetic properties were also determined using PROSA.[36] POLYVIEW[37] was used to predict secondary structure profile and solvent accessible area.

Molecular dynamics simulation was performed on the molecule showing lowest MODELLER objective function. These studies employed NAMD 2.5 (Nano-scale Molecular Dynamics)[38] by applying CHARMM27 force field[39] for lipids and proteins[4041] along with the TIP3P model for water. Energy minimization of the structure was done using 10,000 steps. Multiple time-stepping algorithm was used[4243] with an integration time step of 2 fs. Various interactions were computed in 1, 2, and 4 time steps for covalent bonds, short-range non- bonded interactions and long range electrostatic forces, respectively. For every ten time steps, non-bonded interactions were with a pair list distance of 13.5 Å. Van-der-Waals and electrostatics interactions were defined as interactions between short-range non-bonded interactions between particles within 12 Å. For Van-der-Waals interactions at a distance of 10 Å, a smoothing function was employed. Simulations were performed on the equilibrated system for 1 ps while maintaining a restraint of 500 kcal/mol/A° 2 on the protein backbone under constant pressure and temperature of 1 atm and 310 K, respectively, with Langevin damping coefficient set to 5 ps. Structure showing least energy with converged root mean square deviation was used for subsequent exercise. Final structure was visualized using VMD.[32]

Binding pockets were determined and explored using CASTp (Computed Atlas of Surface topology of Proteins)[44] and residues lining the functional site were also identified. The results were compared with highest scoring amino acids (score=9) predicted by Conseq results so as to identify the amino acids conserved in other species. Based on such observation, atom 3739 of Phenylalanine 226 was selected as target atom for docking. A total of 41 ligand molecules were selected based on literature[45-47] for docking studies to determine binding affinities of the ligand molecules toward the modeled NMT.

AMPAC software (http://www.semichem.com/ampac/default.php) was used for molecular mechanical calculation where AM1 calculations with SCF were performed applying restricted Hartree-Fork method. GOLD 2.0 package, (Cambridge Crystallographic Data Centre, Cambridge, UK)[48] which performs exhaustive and exclusive search for different conformations efficiently maintaining the flexibility of the ligand molecule provided by variation in dihedral angles was employed for docking studies. Default annealing parameters were considered for Van der Waals force and hydrogen bond calculation. Parameters considered for genetic algorithm (GA) are as follows: Population size 100, Selection pressure 1.1, No. of operations: 100000, No. of Islands: 5, Niche Size: 2, Migrate: 10, Mutate: 95, and Crossover: 95. Active site radius was set to 15 Å for docking calculations. Molecular interactions between the ligands and the protein were analyzed using SILVER.

RESULTS AND DISCUSSION

Sequence and structural analysis of the P. falciparum NMT was carried out using bioinformatics tools. Amino acid composition of a protein reveals a lot about its nature. Amino acid composition was calculated using Protparam[18] from Expasy (http://us.expasy.org/tools/protpar/) and the obtained results are summarized in Table 1.

Table 1

Amino acid composition of P. falciparum N-myristoyltransferase

It was found that majority of amino acids present in NMT are hydrophobic (46.1%). Polar amino acids and charged amino acids constitute 23.2% and 27.5%, respectively, while glycine constitutes only 3.2% of all amino acids. Figure 1a shows the detailed representation of amino acid composition of NMT.

Figure 1

(a) Pie chart diagram representation of composition of P. falciparum N-Myristoyltransferase; (b) Protein disorder predicted using metaProtein disorder prediction system. 2 peaks clearly visible above threshold value of 0.5 denote the disordered regions; (c) Transmembrane region predicted using DAS server. The peak above 1.7 indicates the transmembrane region

Other properties calculated using Protparam are summarized in Table 2.

Table 2

Properties determined using Protparam

Instability index of NMT indicates about the instability of the protein as a value above 40 is hallmark of unstable proteins. The aliphatic index (AI), which denotes the relative volume of a protein occupied by aliphatic side chains was found to be 94.37. GRAVY value of 0.327 tells about its hydrated state. Theoretical pI of 8.39 classifies the protein as basic. Other important physico-chemical properties were calculated using Protscale.[18] Protscale assigns value to each amino acid using a scale and the results are presented below [Table 3].

Table 3

Important physicochemical properties calculated using protscale

Numerous NMT sequences were retrieved by the CONSEQ server[19] using P. falciparum NMT protein sequence as the query keeping the default options of BLAST E-value threshold: 0.001, maximum number of homologs: 50, iteration: 1. 32 out of 33 PSI-BLAST hits were found to be unique and the calculation was performed on the unique hits. The conservation scores versus residue number were determined and are shown in Figure 2a. An unrooted phylogenetic tree was constructed using the tree building facility of CLUSTAL-W employing the multiple sequence alignment obtained from MUSCLE [Figure 2b]. It was found that among the secondary structure elements, alpha helices were found to be predominant followed by random coils, extended strands and beta turns in descending order as predicted by DPM, DSC, GOR3, HNNC, SOPMA, Sec. Consensus. Results of three servers viz. MLRC, PHD, Predator indicated that random coils outnumbered alpha helices and extended strands in the protein [Table 4].

Figure 2

(a) Conservation scores of amino acids on a scale varying from 0–9 indicating variable to conserved amino acids where e-An exposed residue according to the neural-network algorithm, b-A buried residue according to the neural-network algorithm, f-A predicted functional residue (highly conserved and exposed), s-A predicted structural residue (highly conserved and buried), X-Insufficient data-the calculation for this site was performed on <10% of the sequences (b) Phylogenetic tree obtained using CONSEQ

Table 4

Secondary structure prediction using NPS server

MEME was used for the elucidation of motifs in P. falciparum NMT with the parameters set to their default values. Three motifs predicted using MEME server along with their positions are shown in Table 5.

Table 5

Detailed information of motifs predicted using Meme

5 motifs were predicted using PROSCAN and 2 signature motifs viz Myristoyl-CoA: Protein N-myristoyltransferase signature 1(EVNFLCVHK) and Myristoyl-CoA: Protein N-myristoyltransferase signature 2(KFGEGDG) were found. Pattern, probability and description of the motifs are presented [Table 6]. 3 pseudomotifs were also predicted using Quasimotifinder apart from those predicted using PROSCAN [Table 7].

Table 6

Output of PROSCAN along with the probability and patterns of motifs predicted

Table 7

Pseudomotifs predicted using QuasiMotif

Disuphide bridges known as “switches for protein function”,[49] result from covalent bonding of suphur from cysteine residues. Disulphide bridges play an important part in folding of protein and are also responsible for stabilization of protein structure. Keeping this in mind, disulphide bridges were calculated using CYS-REC. The search resulted in eight cysteine residues and the most probable pattern of pairs are 74–400 and 76–397 [Table 8 in supplementary].

Table 8

Details of patterns of Cystine–Cystine binding

Some regions in protein occur as dynamic and unstructured ensembles and are called disordered regions. Identification of protein disorder is important for understanding protein function.[50] The disorder in protein facilitates its molecular interaction with multiple partners and is implicated in provision for various modification sites.[51] Disorder prediction by metaProtein disorder prediction system showed the presence of two such regions viz. MNDDKKDFVGRD (1-12) and DEFDENVNEPFISDN (41-55) with a value higher than the threshold of 0.5 [Figure 1b].

DAS[21] predicted a transmembrane region spanning from 139–146 at a cutoff value of 1.7 [Figure 1c], while TMPRED predicted three helices (inside to outside 102(102)-122(120); 186(188)-209 (206); 329(329)-349(349) and two outside to inside helices (191(193)-209(209; 329(329)-345(345)). Results from SOSUI[23] server classified the protein as a soluble protein. HMMTOP, TMHMM, and SPLIT server did not predict any region spanning the membrane in the protein.

Protein structure often reflects its function. This fact makes the protein structure prediction a lucrative exercise in wake of unavailability of experimentally derived protein structures. Comparative modeling involves assigning the structure to a protein for which the structure has not been determined based on its sequence similarity to already known protein structures. This is based on the assumption that similarity at primary structure level is often indicative of structural similarity. This method assures reasonably good protein structure prediction which can provide deep insight in the mechanism of protein function in absence of its crystal structure.

Chain A of myristoyl-CoA: Protein N-myristoyltransferase (PDB ID: 1RXT_A) determined by X-ray diffraction at a resolution of 3.0 Angstrom was selected as a template as the protein shared 44.3% identity with the target protein sequence. Sequence alignment between target and template was generated using CLUSTALX and used as input for model construction [Figure 3a].

Figure 3

(a) Sequence alignment of template and target protein. Blocked regions represent the conservation between the target and template; (b) DOPE score of models generated using Modeller; (c) Modeller objective function of models generated using Modeller; (d) RMSD variations in dynamics calculation

A three-dimensional model was constructed using comparative modeling employing MODELLER9v3 based on probability density functions (PDFs). Total 25 models were generated using Modeller and evaluated based on DOPE score and Modeller objective function [Figure 3b–c]. Model having least objective function, i.e., 4486.498 was selected for further refinement. Model was refined by NAMD simulations using strategy discussed elsewhere.[3031]

Obtained model was further refined using molecular dynamics. A graph was plotted between the trajectories generated as function of time (ps) and RMSD of C alpha trace of the protein molecule. During the MD simulation run, a total of 109 frames were generated and it was observed clearly that the obtained RMSD data showed a variation from 0.47 to 1.30 and thereafter, attained a plateau state with minor variation far away from the decimal points. Model was found to be stable above 1.3 ps of molecular dynamics simulation [Figure 3d]. A schematic presentation of the modeled protein is shown in Figure 4a.

Figure 4

(a) Three-dimensional structure of P. falciparum NMT where alpha helix are shown in purple, 3_10_helix in blue, Pi_helix in red, extended beta sheets in yellow, bridge beta in tan, turn in cyan and coil in white color (b) Superimposition of target and template structure (Target in cyan color, template in mauve color) (c) Ramachandran plot analysis (d) PROSA curves representing the residue interaction energies

Superimposition of target and template indicates the structural similarities and differences and a close homology between template and target is expected as revealed by RMSD of 1.5 [Figure 4b]. This further reinforced the reliability of our model and this model was used for subsequent analysis. The protein belongs to alpha class of protein.

Stereochemical properties of the model were evaluated using PROCHECK. Ramachandran plot analysis of Psi and Phi dihedral angles showed that 97.7% fall in the allowed regions of the plot while 2.3% of the residues fall in disallowed regions [Table 9 and Figure 4c]. These results are comparable to the template structure. Several residues Asn 19, Gln 94, Lys25, Asp22, Tyr65, Tyr28, Ser33, Phe345, Phe51, Leu16, Lys23, Thr134, Tyr393, Phe226, and Ala198 fall outside energetically favorable regions. The ERRAT score of 76.294 (which is greater than 60.542 score of template) confirms about the reliability of the structure as a score of greater than 50 indicates the high quality of the model (data not shown). About 77.41% of residues in the model have a score >0.2 as revealed by VERIFY3D. This also underlines the reliability of structure (data not shown). PROSA score is a measure of model quality and calculates the deviation of the total energy of the structure with respect to an energy distribution derived from random conformations. PROSA score of model protein is -6.55, which is negative and is comparable to -9.36 [Figure 4d]. As negative value of residue interaction energy is the criterion for the correctness of the model, this further indicates good quality of model.

Table 9

Stereochemical evaluation of obtained protein model in Ramachandran Plot analysis using Procheck

Secondary structure of a protein is considered as the local spatial arrangement of main chain atoms. Obtained model was analyzed for determining the secondary structure and physiochemical profiles, solvent accessible surface area using POLYVIEW.[37] This provided the information on presence of secondary structure elements at specific locations [Figure 5a]. It was found that the protein contains 11 α helices, 46 turns, and 20 β strands altogether. A saddle shaped mixed β sheet was found in the core, which was further surrounded by several other α helices.

Figure 5

(a)Secondary structure, physiochemical profile and solvent accessible surface area as predicted by POLYVIEW H-α and other helices (view 1), H-α and other helices (view 2), E-β-strand or bridge, C–coil, Relative solvent accessibility (RSA) 0-completely buried (0-9% RSA), 9-fully exposed (90-100% RSA), where H-hydrophobic: A,C,F,G,I,L,M,P,V; A-amphipathic: H,W,Y; P-polar: N,Q,S,T and N/C-charged: D,E–negative, R,K–positive; (b)Top 10 possible binding sites as predicted by CASTP where Pocket 1=green, Pocket 2=blue, Pocket 3=cyan, Pocket 4=yellow, Pocket 5=magenta, Pocket 6=pink, Pocket 7=orange, Pocket 8=purple, Pocket 9=brown, Pocket 10=gold. (c) Active site of modeled protein. (d) Surface occupied by Pf NMT in spacefill model

The local environment inside binding pockets and functional groups of amino acid lining the cavities on protein determines the function of protein and influence substrate binding. Active site of the modeled structure was determined using CASTP. Total 94 binding sites were predicted. Top 10 binding sites based on area are shown [Figure 5b]. Binding site having area of 754.9 and volume of 979 was further explored [Figure 5c]. Active site was found to consist of L16, I 17, N19, A20, K21, W73, D83, R84, Y95, V96, E97, D98, D99, N101, F103, F105, Y107, K167, T197, A198, G199, V200, Y211, U213, F226, L316, S318, L330, A332, F334,V363, N365, L367, D385, G386, S387, L388, Y390 and L410 [Figure 5c]. Volume and surface of the protein model are represented in Figure 5d.

Based on GOLDSCORE, which is dependent on the protein-ligand hydrogen bond energy (external H-bond), protein-ligand van der Waals (VDW) energy (external VDW), ligand internal VDW energy (internal vdw), and ligand torsional strain energy (internal torsion), the ligand protein-binding efficiency was computed. Molecule 7 in our study, which is a benzofuran compound with CH2CH2Ph and CH3 side chains[45], showed highest fitness score (65.93). It showed the presence of two hydrogen bonds, i.e., N-37 of ligand molecule with Serine 387(Bond length=2.422) and O12 of the ligand molecule with Tyrosine 211 (Bond length=2.154). Besides, it showed close contacts with C19 of ligand molecule with Glutamic acid 97 (Bond length=(2.689) and H52 of ligand molecule with Leucine 410 (Bond length=1.846)) [Figure 6].

Figure 6

Top scoring benzofuran compound obtained in the docking study

Molecule 9 in the study,[45] belongs to the similar compound group with CH2SPh and CH3 side chains showed the second best fitness score of 65.74 in best rank file and displayed several close contacts which did not come up to the range of hydrogen bond. Following close contacts were observed: S53 of the ligand molecule with ALA198 of protein (Bond length=2.529), C43 of the ligand molecule with TYR95 of protein (Bond length=2.477), C41 of the ligand molecule with TYR95 of protein (Bond length=2.589), C5 of the ligand molecule with GLU97 of protein (Bond length=2.588), C34 of the ligand molecule with SER319 of protein (Bond length=2.586), N37 of the ligand molecule with SER319O of protein (Bond length=2.451) [Figure 7a]

Figure 7

Second (a) and third (b) best benzofuran molecule found in the docking calculations

Molecule 14 in the study (A benzofuran compund with a Phenyl R group)[45] displayed score of 63.67 and showed the formation of 2 hydrogen bonds; one between O24 of the ligand molecule with Tyrosine 95 (Bond length=2.644) and another between O13 of the ligand molecule with TYR211 (Bond length=2.192) and 4 close contacts [Figure 7b]. Ranks of the rest of the molecules based on the GOLDSCORE fitness values are represented in Table 10.

Table 10

Fitness scores of ligand molecules considered for the docking studies

CONCLUSION

Despite advent of high throughput methods, there exists a huge gap in number of available protein sequences and experimentally derived protein structures. Comparative modeling is often envisioned to play a role in providing details about protein structure in absence of crystal structures. We have applied bioinformatics tools for determining important features and properties of the NMT of Plasmodium falciparum.

The extent of reliability of structure prediction is dependent totally on the degree of sequence similarity between the target and template sequences. The methodology adopted by us resulted in a good quality model as the model was constructed based on a sequence homology of 44.3%. Model obtained was assessed using various structure validation servers and was found to be reasonably good. We have performed a thorough in silico characterization of N-myristoyltransferase of malarial parasite. The outcome of this study will enhance our knowledge about this enzyme and will boost drug designing process as various aspects of structural characteristics of the enzyme explored in this study can provide a basis for effective inhibitor design. The model so obtained can be used as a potential target and results of this study can be exploited for the development of effective treatment plans for combating this dreaded disease.