Helicobacter pylori, a Gram negative, micro-aerophilic bacteria that colonize in human stomach leads to many diseases. The statistical analysis also revealed that more than 50% of the world’s population was infected with H. pylori [1,2]. Colonization of H. pylori in gastric mucosa results in the development of chronic gastritis with intolerable gnawing pain. Patient with chronic gastritis progresses to complications like gastric ulcers, gastric neoplasia and some distinct extra gastric disorders [1,2]. Vacuolating cytotoxin autotransporter is a 140 kDa protoxin with 1287 amino acids released by H. pylori in lumen of gut. It undergoes N- and C-terminal cleavage during the secretion process to yield a mature 88 kDa toxin consisting of 1-821 residues length (p88) [3-5]. The N-terminal p33 domain (residues 1–311) contains a hydrophobic sequence residues from 6–27 involved in pore formation, whereas the p55 domain residues from 312–821 contains one or more cell-binding domains [6].

VacA makes anionic channels on the plasma membrane, these pores are very important for entry of pathogens into the cells and forming vacuoles in H. pylori epithelial cells [7,8]. VacA can cause membrane depolarization, alteration of mitochondrial membrane permeability, detachment of cells from the basement membrane, activation of mitogen-activated protein kinases, inhibition of antigen presentation, inhibition of T cell activation and proliferation and VacA-induced cell death [8,9]. This process confers development of pathogenesis in the human gut such as gastritis, peptic ulcerative syndromes and other pleiotropic effects linked to adenocarcinoma [7,10]. A receptor-like protein-tyrosine phosphatase RPTPα and A receptor-like protein-tyrosine phosphatase β (RPTPβ) were reported as important VacA-binding proteins located in epithelial membrane but the more information about binding residues of VacA was unclear [11,12]. However, amino acid sequences within the p-55 domain β-helices mediate the binding of VacA to host cells [7,13-16]. Many research groups have indicated that binding region on p55 domain considered as an essential site in making attachment to epithelial cells results in oligomeric anionic channels and vacuoles, specifically 111 amino acids from the N terminus of p-55 domain are required for vacuolating toxin activity [13,17-19]. In the present study this pore forming protein factor p55 domain was therefore considered as an important drug target to propose inhibitors which would interfere p55 domain attachment with epithelial cells as well as with p33 [8,20,21]. About the cleavage of VacA into p55 domain and p33 domain was unclear yet some research groups came up with in vitro studies regarding p55/p33 binding [13,14]. However, their studies reported that p55 domain is more crucial in adherence than p33 domain of VacA to epithelial cells which lead to further maladies [7,13,14]. In this scenario, the p55 domain of functional VacA became an attractive target of interest and attained great attention towards finding of inhibitors through in silico approaches. Despite the unavailability of complete VacA crystal structure, the p55 domain structure made it possible to study and design novel inhibitors against VacA toxin as it is significant in adherence to epithelial cells as well as to bind p33 domain if it cleaved at any circumstances [13,14].

VacA employs cellular pleiotropic effects which worsen the gastrointestinal tract towards severe uncontrollable maladies. To evade such conditions there is a great need to design novel inhibitors against VacA. So, the goal of the study is to investigate the novel inhibitors against p55 domain of VacA toxin produced by H. pylori.

Therefore, VacA p55 domain was explored with docking and molecular dynamic simulations to propose potential lead molecules and also to check the stability of target-lead complex. Moreover, proposing inhibitors could act against gastritis, peptic ulcers and aid in obstruction of pathogen entry, colonization and survival of H. pylori.

Hardware and software

The work was carried out in HPZ800 workstation (12 processors; 12 GB RAM; NVIDIA graphics card; Cent OS 6.0) and state of art HP personal computers. The systems were installed with vendor licensed Schrodinger software suite (Schrödinger LLC, 2014.2). The computers were connected to 4 Mbps dedicated broadband connection to access major Bioinformatics resources.

Binding site selection and protein preparation

The crystallized apo structure of p55 domain (2QV3) of VacA was downloaded from the protein data bank [8]. Understanding the structure and exploiting the function of active site and residues is a cornerstone of virtual screening and drug design. SiteMap v3.1 module of Schrodinger software was used to explore binding site on p55 domain of VacA. The protein was further analyzed with FINDSITE to verify binding site residues on p55 domain. Torres et al. [13,21]; Ye et al. [18,17] reported that 313-478 interacting residues on p55 domain were important for p33 domain attachment in case of cleavage occur, Wang et al. [16] reported that total p55 domain is important for VacA binding to cell studied from HeLa cells, so it was also taken into consideration in this study. VacA p55 domain structure was imported to Maestro v.9.8 and protein preparation was performed using protein preparation wizard. Hydrogen atoms were added to the receptor subsequently minimized by optimized potentials for liquid simulations (OPLS) 2005 force field. Minimization was performed by restraining heavy atoms with hydrogens, to allow free rotation of hydrogens setting the maximum root mean square deviation (RMSD) of 0.30Å [22].

Ligand preparation

Ligands from small molecule databases such as ChemBank, ChemPDB, Drug likeliness NCI, KEGG Ligand, Anti-HIV NCI, AKos GmbH, Unannotated NCI, Asinex Ltd. and TimTec subsets consisting of 2344, 4009, 192,323, 10,005, 42,689, 544,391, 15,237, 348,276 and 7,500 entries respectively were prepared using Schrodinger software [23,24]. LigPrep v3.0 [25] and Epik v2.8 modules of Schrödinger were utilized for ligand preparation with attaining chemical correctness for all ligands by elevating the stereo chemical nature, protonation and tautomeric states. Energy of the ligands was minimized at pH 7.0 ± 2.0 and utilized for docking studies.

Docking and scoring

The prepared protein was imported to Glide v.6.3 to generate a grid box of 10 × 10 × 10 Å around the predicted interacting site on p55 domain. Three different docking strategies such as rigid receptor docking (RRD), quantum polarized ligand docking (QPLD) and induced fit docking (IFD) were performed to calculate the scores and binding affinity of the interactions between p55 domain and leads [27]. The in-house library of prepared ligands and protein were imported to dock and to screen the best ligands by high throughput virtual screening (HTVS), standard precision (SP) and extra precision (XP) docking methods [28]. Prime/MM–GBSA (molecular mechanics/generalized born surface area) method was used to calculate the binding-free energy (ΔG) of the docking complex [29]. The docked poses were minimized using the local optimization feature in Prime v3.6, and energy of the complexes were calculated using the OPLS-AA 2005 force field and generalized-born/surface area (GB/SA) continuum solvent model [26]. The best ligands based on ΔG scores obtained from RRD were subjected to QPLD using Q-Site v6.3 module of Schrödinger. Quantum mechanics and molecular mechanics (QM/MM) were calculated and polarizable ligand charges were induced by protein [26].

The best lead resulted from QPLD with lowest binding free energy (ΔG) was further assigned to induced fit docking. IFD protocol was adopted to consider the flexibility of both lead and receptor [27]. In this docking stage, Glide docking parameters were set to the hardpotential function was used with van der Waals radii scale limiting of 0.5 Å for receptor lead docking calculations [27,30]. Binding affinity of the complex was reported in terms of binding free energy (ΔG). The lead molecules with reactive functional groups and high energy ionization states and compounds disobeying the Lipinski’s rule of five were removed from the study. Physicochemical descriptors and pharmaceutically relevant properties (ADME/T) for the leads were predicted by using QikProp v.4.0 [28].

Molecular dynamics simulations

Molecular dynamics (MD) simulations were performed to VacA p55 domain-lead1 dock complex for 50 ns using Desmond v3.8 [31]. Simple point charge (SPC) water molecules were embedded with OPLS_ AA-2005 force field parameters to depict a solvated model system [31]. The system was neutralized by adding minimum number of sodium and chloride ions in order to balance the net charge in the solvated system.0.15 mol/L sodium and chloride were then added to mimic the osmotic effect of water [32]. Electrostatic interactions were applied using particle mesh Ewald (PME) [31,33]. System was minimized before simulating and was carried out with the periodic boundary conditions in the NPT ensemble. The temperature and pressure were kept at 300 K and 1.01325 bar pressure using Nose-Hoover temperature coupling and isotropic scaling. The minimized system was passed through NVT ensemble for short 12 ps (picoseconds) simulations at 10 K temperature, additionally non-hydrogen solute atoms were restrained for 24 ps at 300 K temperature. Further, the system was simulated for 24 ps in NPT ensemble at 300 K temperature without restrains so as to attain equilibrium state [34]. The minimized system without restrains was subjected to 50 ns NPT simulation production [35,36]. Trajectories were recorded for every 4.8 ps time interval. The stability of VacA p55 domain-lead1 docked complex was examined through potential energy and root mean square deviation (RMSD), root mean square fluctuations (RMSF) and hydrogen bond interactions.

Binding site selection and ligand preparation

As a result of the sitemap, the best binding site with site score 0.928 was selected, in which the predicted sitemap binding residues are A431, S434, N436, F437, E438, K440, F461, N463, K465, F483, N484, T485, D487, S504, T505, N506, A508, G527, E528, Y529, H531, S533, E534 and K562. As per software module, a site score with 0.80 has been found to accurately distinguish between drug binding and non-drug binding sites (Schrodinger LLC, 2014-2). The interacting residues from the literature study [13,16,17,36] were also correlated and present within the predicted binding region on p55 domain of VacA considered for docking study. So the interruption of this binding site makes VacA unable to form oligomerization as well as adherence to the RPTPα and RPTPβ receptors on epithelial cells. FINDSITE analysis also revealed to be made sure that the same residues whether important for VacA p55 domain binding or not. The prepared in-house library consisting of 1,166,774 compounds were used for docking study.

Molecular docking analysis

Docking is a method which predicts the preferred orientation of ligand to receptor that forms a stable complex and used to predict the binding affinity between two molecules using scoring functions [37]. A wide range of different docking programs are available, most of which use semi-rigid docking, where the ligands are treated as flexible and the receptors as rigid. It offers full spectrum of speed and accuracy from VHTS of millions of compounds to extremely accurate binding mode predictions, providing consistently high enrichment at every level. Docking of 1,166,774 prepared ligands with p55 domain through highthroughput virtual screening had given only 10% good conformers of total screened outcomes with setting option and allowed to standard precision. With the refinement, SP resulted 8439 the best poses as the 10% outcomes from the screening. The dataset was further condensed to 4309 based on best scoring ligands through XPG score along with Prime-MM/GBSA. XPG Score is an empirical scoring function that approximates the ligand binding energy and the parameters such as force fields, penalties for the interactions that have the influence of ligand binding with the receptor. The XPG Score is given by:

GScore = a*vdW+ b*Coul + Lipo + Hbond + Metal + BuryP + RotB + Site

Whereas vdW denotes van der Waals energy, Coul denotes Coulomb energy, Lipo denotes lipophilic contacts, HBond indicates hydrogen-bonding, Metal indicates metal-binding, BuryP indicates penalty for buried polar groups, RotB indicates penalty for freezing the rotatable bonds, Site denotes polar interactions with the residues in the active site and the a = 0.065 and b = 0.130 are coefficient constants of van der Waals energy and Coulomb energy respectively [38,39].

Finally, 10 ligands were identified as RRD best ligand candidates through careful inspection of the docking poses, free energy scores and possible interactions with the binding site for all of the active compounds. Here, HTVS, SP and XP docking strategies has given ligands selectively from lower stringency to higher stringency at every level of glide docking. Among which top ten ligands were selected based on XP Glide score for QPLD to calculate ΔG scores. The QPLD protocol gives more accurate treatment of electrostatic interactions by using quantum mechanics, which results in an improvement of the docking accuracy. Top ten leads from QPLD with lowest binding free energies were ranked based on ΔG scores. Natarajan et al. [27,39], reported that quantum polarized ligand docking interactions are reliable for ranking the leads than the regular Glide XP scoring Quantum polarized ligand docking (QPLD) was employed to compute the atomic partial charges (APC) of the leads through quantum mechanical and molecular mechanical (QM/MM) calculations. Natarajan et al. [27] reported that absolute closer experimental affinities can be derived by the including the entropic effect on the binding free energy (ΔG) calculated using MM-GBSA. The binding free energy (ΔG) derived through MM-GBSA calculations had much more consistency than the docking scores and improves the ranking of potential leads [39].

All lead molecules satisfied pharmacological properties of 95% drugs from the Lipinski’s filter. All the ten leads were well interacted with p33 domain binding site in the same orientations and with better binding affinities as well as with the lowest binding free energies (ΔG). Among them lead1 has the lowest binding free energy (ΔG) of –48.286 kcal/mol (Table 1).

Table 1: Gscores and ΔG scores of the ten proposed leads.

Lead1-receptor docking complex showed one hydrogen bond with T485, one π-π stacking interaction with F483 and good van der Waals interactions with S434, N463, N484, S504, T505, N506, G527, E528, Y529 and H-531 residues (Figures 1 and 2). Lead1 has also better XPG docking score of –3.860 kcal/mol as well as good interactions with p55 domain among the ten leads was shown in Table 1. Pharmacological features of the leads were correlated favorably with more than 95% approved drugs showed in Supplementary Tables S1-S3. In IFD, the molecular interactions between p55 domain and lead1 docking complex showed that one hydrogen bond with polar residue T485 and one π-π stacking interaction with hydrophobic Phe483. Some other polar residues S434, N463, N506, N484, S504, T505, H511; nonpolar and hydrophobic residues F461, F483, Y529 as well as negative charged E528 were involved in intermolecular van der Waals contacts were considered as significant (Figure 3). Side chain orientations were automatized in IFD within flexible mode. P55-lead interaction energies and total energy of the system was calculated as IFD score. Based on the IFD scores the poses generated were ranked. The side chains of binding site residues were solely refined and reoriented without disturbing the lead of docking complex in IFD, so the binding mode of lead will not undergo much conformational change. In RRD and QPLD, the receptor and ligand do not undergo structural changes. In IFD, due to applying side chain refinement and rearrangement of the receptor, the docking complex showed additional interaction with F461. The same interactions found in IFD indicated that binding residues were well correlated with QPLD. Therefore three different docking strategies (RRD, QPLD and IFD) and MM-GBSA re-scoring further affirmed the ability of ten leads as potent inhibitors of VacA p55 domain.

Figure 1: VacA p55 domain-lead1 complex.

Figure 2: VacA p55 domain-lead1 from QPLD.

Figure 3: Induced fit docking of VacA p55 domain-lead1.

MD simulations analysis

The dynamics process and interaction mechanisms of receptorligand system were obtained by simulating their internal motions. A number of marvelous biological functions in proteins with their profound dynamics mechanisms, such as switching between active and inactive states, cooperative effects, allosteric transition and intercalation of drugs can be revealed by studying their internal motions [32]. An analysis of the interaction patterns in the production simulations allowed to identify the key residues of P55 involved in the binding with potential lead [38]. Molecular dynamics simulations were performed in the system embedded with water molecules, temperature, and pressure in closer to the physiological environmental condition [35,36]. The dynamical properties of system were analyzed from trajectory data obtained from 50 ns MD simulations [35,38]. The final system contained 69,565 atoms and passed through a six step relaxation protocol prior to MD simulations. The energy plot indicated that the system has lower energy in the complex was fairly stable throughout the 50 ns simulations period; it was shown in Figure 4. The energy plot showed that the complex was stable with constant lowest energy level during the entire run of MD simulations. RMSD plot for backbone and-heavy atoms showed stability of the complex after a small rearrangement from the initial conformation. The RMSD of lead1 was observed to stay upto 1Å. The backbone atoms mostly remained within 2.4Å in most of the trajectories and the graph was represented in Figure 5. Therefore, the lower RMSD of the system indicates lesser conformational changes and stable nature of the complex during entire MD simulation period [23]. RMSF determines fluctuations of each residue in VacA p55 domain during 50 ns simulation period; it was shown in Figure 6. The average RMSF for VacA p55 domain residues backbone and side chain was 2.0 Å and 2.5 Å respectively. One intermolecular hydrogen bond was observed in lead1-receptor complex during MD simulations. More in detail, the carbonyl oxygen of the lead molecule interacted with three water molecules. None of the water bridge or water mediated interaction was found between the receptor and ligand. Transient variations of set of interactions were also observed in binding cleft of p55 domain along with ligand individually. Upto 200 pico seconds simulation a few number of water molecules were involved in interaction with more rigidity. After 5000 ps to 10 ns the water molecular contacts were gradually decreased in the system. But three water molecules interacted with oxygen (O) of lead1 were constantly involved in the system upto 10416 trajectories. The hydrogen bond between T485 and lead1 (atomic distance 2.5 Å in ~95% trajectories was conserved in ~95% of MD simulation time. So, the consistent water molecules binding to lead1 was found in the trajectory analysis throughout 50 ns simulations which concerns that may increase complex stability [35]. The interactions of lead1 and residues of receptor were shown in Figure 7. The intermolecular interactions of lead1-receptor complex obtained from MD simulations were well corroborated the QPLD and IFD docking results. The residues include polar residue T485 with H-bond, S504, T505, N506, G527, E528, Y529 and H531with van der Waals interactions were attained importance in the study due to significantly involved in correlated docking and MD simulations. In addition, the identified residues were all within the boundary of binding site and identical. The energy plot, RMSD and RMSF analysis revealed that VacA p55 domain–lead1 docking complex was conformationally stable during 50 ns MD simulations.

Figure 4: Energy plot of VacA p55 domain-lead1 complex.

Figure 5: RMSD plot of p55 domain and lead1.

Figure 6: RMSF plot of VacA p55 domain-lead1 docking complex.

Figure 7: Interactions of lead1 with VacA protein p55 domain during 50 ns MD simulations.

Trajectory analysis of docking complex in entire MD simulations showed that high conformational stability and less structural flexibility in physiological environment. In the context of inappropriate drugs against H. pylori and its related uncontrollable diseases rose the demand to find novel leads against H. pylori and its virulent factors. Here, VacA was preferred as a significant target due to its role in causing intolerable disease like chronic gastritis and having no direct inhibitors or drugs against. Due to unavailability of inhibitors against VacA the study needed to design prominent novel leads against VacA or its domains like P55/p33. Therefore, proposed leads with further experimental validation would be a starting point towards designing potent inhibitors molecules for treatment of gastritis and peptic ulcers.

The goal of this study was to find novel and more powerful lead molecules against painful chronic gastritis caused by H. pylori. In the study virtual screening (HTVS, SP, XP), docking (RRD, QPLD, IFD), MD simulation studies were applied successfully to propose ten novel lead molecules against p55 domain of VacA toxin of H. pylori. In comparison, the lead1 has the highest binding affinity than other leads. The interactions in docking complex of lead1-p55 domain were well correlated with 50 ns MD simulations and showed that consistently lower fluctuations as well as good binding affinity at physiological environment. As the proposed leads are obeying good (ADMET) pharmacological properties with greater binding affinity as well as the lowest binding free energy would impede the VacA binding onto epithelial plasma membrane if synthesized and tested in animal models. Therefore, the proposed leads could be potential inhibitors against H. pylori causing chronic gastritis and other ailments by deliberately obstructing the anionic channels formation, vacuoles and entry of pathogenic H. pylori into epithelial cells.

PC is thankful to ICMR, New Delhi for the award of Junior Research Fellowship (3/1/3/JRF-2014/HRD-8) and authors are grateful to DBT, Ministry of Science and Technology, Government of India, New Delhi for supporting the work through BTISnet BIF program (BT/BI/04/055/2001).