Identification of novel natural inhibitor for NorM – a multidrug and toxic compound extrusion transporter – an insilico molecular modeling and simulation studies

The emergence of bacterial multidrug resistance is an increasing problem in treatment of infectious diseases. An important cause for the multidrug resistance of bacteria is the expression of multidrug efflux transporters. The multidrug and toxic compound extrusion (MATE) transporters are most recently recognized as unique efflux system for extrusion of antimicrobials and therapeutic drugs due to energy stored in either Na+ or H+ electrochemical gradient. In the present study, high throughput virtual screening of natural compound collections against NorM – a MATE transporter from Neisseria gonorrhea (NorM-NG) has been carried out followed by flexible docking. The molecular simulation in membrane environment has been performed for understanding the stability and binding energetic of top lead compounds. Results identified a compound from the Indian medicinal plant “Terminalia chebula” which has good binding free energy compared to substrates (rhodamine 6 g, ethidium) and more favorable interactions with the central cavity forming active site residues. The compound has restricted movement in TM7, TM8, and TM1, thus blocking the disruption of Na+ – coordination along with equilibrium state bias towards occlude state of NorM transporter. Thus, this compound blocks the effluxing pathway of antimicrobial drugs and provides as a natural bioactive lead inhibitor against NorM transporter in drug-resistant gonorrhea.


Introduction
Bacterial multidrug resistance to antibiotics is a significant obstacle in clinical treatment of many infectious diseases caused by Mycobacterium Tuberculosis, Salmonella typhi, Vibrio cholerae, Staphylococcus aureus, and many other micro-organisms. Different virulent target proteins of these pathogens have been identified, which showed resistance to antibiotics such as β-lactams, kanamycin, trimethoprim, methicillin, vancomycin, etc. Many drug targets were used to understand inhibition of drug-resistant micro-organisms using screening of herbal leads by structure-based virtual screening (VS) approach, and explore the structural insight in binding (Skariyachan et al., 2013). There are several mechanism of resistance of these genes and their corresponding gene product, such as acquired or transmitted through horizontal gene transfer or transposable elements and modification in enzymatic activity by mutations. The Klebsiella pneumoniae carbapenemases (KPCs) is one of the β-lactamases which is present in plasmid and transferred through horizontal gene transfer. Most of the second-generation antibiotics such as cefotaxime and ceftazidime showed resistance against KPC β-lactamases by hydrolysis. Inhibition of this enzyme by structure-based VS had been carried out and identified potential drug candidates (Danishuddin et al., 2014).
Drug-resistant tuberculosis is a one of the mostly studied diseases. Many computational approaches have been used in identifying chemotherapeutic agents against many potential drug targets. One of the known drug therapy agent pyrazinamide (PZA) converts pyrazinoic acid by pyrazinamidase enzyme in Mycobacterium Tuberculosis and loss of this conversion by mutation in this enzyme leads to resistance into this antibiotic. Based on interaction analysis of PZA in wild and mutated pyrazinamidase using computational docking followed by simulation approaches, drug-resistance mechanism has been analyzed and discussed (Rajendran & Sethumadhavan, 2014). One of other key mechanism of resistance in bacteria is the extrusion of cytotoxic as well as antibiotics molecule across cell membrane, being executed by functionally diverse family of multidrug resistant (MDR) transporters possessing overlapping substrate specificity. MDR transporters can be categorized into five main groups such as resistance-nodulation-cell division, ATP binding cassette (ABC) family, major facilitator family, small multidrug resistance family (SMR), and multidrug and toxic compound extrusion (MATE) family. ABC transporters use the energy obtained from ATP hydrolysis to transport drugs, whereas the other groups use H+/Na+ electrochemical gradient for substrate export (He et al., 2010;Song, Ji, & Zhang, 2014). The mycobacterium multidrug resistant protein belonging to SMR efflux pump had been studied for identifying novel lead molecule by computational VS and ADME approaches (Malkhed, Mustyala, Potlapally, & Vuruputuri, 2014). The NorMa MATE transporter has been recently identified in 1999 which extrudes dissimilar lipophilic and cationic drugs such as fluroquinolones, ethidium bromide, rhodamine 6G, crystal violet, berberine, etc. across cell membrane by means of preexisting H+/Na+ gradient. This extrusion of antimicrobial and therapeutic drugs leads to multidrug resistance, a serious and growing public health threat (Brown, Paulsen, & Skurray, 1999). The crystal structure of the NorM -MATE transporter from Vibrio cholera (NorM-VC) has been determined to 3.7 Å resolution in 2010 and revealed an outward facing conformation with two portals opened to the outer leaflet of the membrane. It also showed a unique topology of the predicted 12 transmembrane (TM) helices which are distinct from any other known MDR transporters (He et al., 2010). Molecular dynamics (MD) simulations of MATE transporter had revealed a loss in the V-shape of extracellular vestibule due to the presence of hydrophobic residues in outer parts of TM1, TM7, and TM8 and due to the extrusion of water molecule in the lumen, which provide the access of water and Na+ ions to the cation binding site within the protein (Vanni, Campomanes, Marcia, & Rothlisberger, 2012). MATE proteins of NorM-NG and NorM-VC are Na+ dependent and both having specificity towards high diverse poly aromatic and poly cationic substrates. The crystal structure of NorM transporter from Neisseria gonorrheae (NorM-NG) bound with substrates such as rhodamine 6G (R6G), Tetraphenyl phosphonium (TPP), and ethidium showed that all these substrates bind at the same location on the periplasmic site of the cleft between the two halves of NorM-NG in an apparently outward state. It also showed that crystal structure of NorM-NG contains 12 TM helices which are connected by 11 intracellular and extracellular loops denoted by L1-L2 to L11-L12 ( Figure 1). Among them, L3-L4, L9-L10, and L6-L7 loops are longer. The amino and carboxyl halves of the transporter are distinguished by intracellular loops L6-L7, the extracellular L3-L4 and L9-L10 loops which extend into the drug-binding central cavity formed between amino and carboxyl terminal halves. The central cavity is defined by Asp41, Ser61, Phe265, Ile292, and surrounded by Thr42, Ala57, Leu58, Val269, Gln284, Val285, Ile287, and Ser288. The cavity is also capped with Ser129, Asp130, Asp355, Asp356, and Pro357 from outer leaflet of bilayer (Lu et al., 2013). It has also been found that helices TM 7 and TM 8 moved 6 Å away from the central cavity in analysis from superimposition of the X-ray crystal structures of NorM-VC representing cation bound/ drug-free state and NorM-NG representing cation free/ drug-bound state (Lu et al., 2013). Further, TM7 and TM8 move towards TM10 to engage Asp377 to form Na+ coordination during drug extrusion thus making Na+-driven TM rearrangement which is important in efflux of drugs. The TM helices rearrange such as pulling of TM7 (Phe265) and TM8 (Gln284 and Ser288) away from central cavity to disrupt drug binding. NorM-NG is Na+ coupled and hence conformational changes associated with drug release are likely triggered by Na+ binding. The simulation study revealed that NorM-NG is capable of binding Na+ in cation binding site even when it is in drug-bound state and stabilization of ion in cation binding site is likely to release the drug (Leung, Holdbrook, Piggot, & Khalid, 2014). Thus, conformational rearrangement of TM helices formed by residues Glu261, Tyr294, and Asp377, combined from the NorM-VC and NorM-NG structures enlightened the antiport transport mechanism (Lu et al., 2013;Tanaka et al., 2013). Recent research report on the independence of substrate and Na+ binding at their respective binding site conflicts the classical model of antiport, where binding of the substrate and coupling ion is mutually exclusive (Steed, Stein, Mishra, Goodman, & Mchaourab, 2013). Evidences suggest that transporter adopts multiple states in equilibrium rather than complete transformation between distinct states (Steed et al., 2013). In the current era of drug discovery, computer-aided drug design provides in silico methods for identification of drug based on complementarily favorable interaction in binding pocket of protein. Due to enormity of chemical space, finding new drugs against particular drug target using virtual high-throughput screening approach of large chemical database such as Traditional Chinese Medicine (TCM), Asinex, Pubchem, etc. is of utmost importance and provides top hit list for lead optimization, and scaffold hopping. This will finally lead to the identification of bioactive compounds and new chemical starting points through computational means and thus becomes an important and integral part of modern drug discovery (Reddy, Pati, Kumar, Pradeep, & Sastry, 2007, Ripphausen, Nisius, Peltason, & Bajorath, 2010Shoichet, 2004).
The extrusion of antimicrobial and therapeutic drug by MATE transporter has resulted in problem of drug resistance in bacteria. In our study, high throughput virtual screening (HTVS) of natural compound collections has been carried out against Na+ bound NorM transporter from Neisseria gonorrhea (NorM-NG) for identifying lead compounds for inhibiting efflux pathway of antimicrobial drugs. Flexible-induced fit docking was used in sortlisting top hit compounds from hierarchical filtered hit lists (Farid, Day, Friesner, & Pearlstein, 2006). Further, five best identified lead molecules and two substrates have been used for molecular dynamic simulation and their conformational stability, binding interaction, and conformational changes associated with active tunnel forming TM helices were compared. The identified molecule showed significant binding affinity in terms of theoretically calculated binding free energy from simulation approach which are more favored compared to substrates and provide lead inhibitors against MATE transporter for further in vitro and in vivo experimental studies.

Preparation of ligand
The lead ligands for identifying inhibitors against NorM transporter, few phytochemicals of nine Indian medicinal plants (Albizia lebbeck, Terminalia chebula, Syzgium cumin(L.), Solanum nigrum, Picrorhiza kurrooa Royle ex Benth., Butea monosperma, Saraca indica, Aegle marmelos and Withania somnifera) which have antimicrobial properties against MDR bacteria were downloaded from from DUKE database (http://www.ars-grin.gov/duke/) (Acharyya, Patra, & Bag, 2009). Natural herbs collections from TCMnatural compound database (Chen, 2011) were also used as lead ligands. The ligands were prepared for docking using LigPrep module of Schrodinger suite of tools (LigPrep 2.3, Schrödinger Suite 2009;Sastry, Adzhigirey, Day, Annabhimoju, & Sherman, 2013). First, all the hydrogen atoms were added to ligand molecules as they had implicit hydrogen atoms. The bond orders of these ligands were fixed. The ionization states of the ligands were generated in the pH range of 5.0-9.0 using EpiK, and most probable tautomers and all possible stereo isomers were generated. In the final stage of Ligprep, ligands were minimized with OPLS-2005 force field (LigPrep 2.3, Schrödinger Suite 2009; Sastry et al., 2013). The prepared leads were filtered based on Lipinski's rule of five for their drug-like properties.

Preparation of protein
Three-dimensional protein structure of Na+ bound Norm-NG transporter (PDB ID: 4HUL) was downloaded from RCSB Protein databank (Bernstein et al., 1977). For in silico drug designing, high resolution of atomic structure of protein with correct bond order in minimum energy state is very important to get the correct interactions with the ligands in any successful drug discovery. Thus, protein preparation should be performed (Sastry et al., 2013). The protein structure was prepared by multi-step process through the "Protein Preparation Wizard" of the Schrödinger Suite (2009) (Protein Preparation Wizard, Epik 2.2). Firstly, the bond orders in the protein structure were assigned, hydrogen atoms were added and optimized. Further, the protein-Na+ complex was subjected to energy minimization using the Schrodinger implementation of steepest descent algorithm with OPLS-2005 force field in implicit solvation. The entire complex was minimized and the minimization terminated when the root mean square deviation (RMSD) of the heavy atoms in the minimized structure relative to X-ray structure exceeded 0.3 Å.

VS based on semi flexible docking
Semi-flexible docking studies on prepared ligands were carried out in the substrate recognition site of NorM transporter using docking program, Glide 5.5 Program Halgren et al., 2004). The shape and properties of the receptor were represented on a grid by several set of fields that help progressively in more accurate scoring of ligand poses. The Norm-Na+ complex prepared as described above were employed to build energy grids using default values of protein atom scaling (1.0 Å) within a cubic box dimensions 20 Å centering the residues of the central cavity such as Asp41, Glu261, Ser288, Asp355, Asp377, etc. In this docking, semi-flexible docking protocols were used. The ligands being docked were kept flexible, in order to explore an arbitrary number of torsional degrees of freedom in addition to the six spatial degree of freedom spanned by the translational and rotational parameters. The process of VS was followed in three phases, using three different protocols, i.e. HTVS, Standard Precision, and Extra Precision docking protocols described in Flow Chart 1.

Flexible docking studies of best leads from VS
Induced fit docking studies of the top-50 ligands were carried out by semi-flexible docking studies, wherein induced fit models have been obtained to fit ligands in non-cognate structures. In other words, the protein structure was induced to fit the ligands (Sherman, Day, Jacobson, Friesner, & Farid, 2006). Twenty poses were chosen to be saved after initial Glide docking, which was carried out with the van der Waals scaling of 0.4 Å for both protein and ligand non-polar atoms. After obtaining initial docking poses, Prime side chain and backbone refinement together with the minimization of the docked pose was carried out within a sphere of 5 Å from each pose saved. Glide redocking was performed on Prime refined structures having Prime energy values within 20 kcal/mol of the lowest energy value. Best ligand binding mode was visualized using Maestro 10.1.013 (2014).

MD simulation of best leads from VS
MD simulations were carried out to explore the conformational space. The aim was to uncover differences between nine systems such as X-ray crystal structure of Na+ bound NorM transporter in complex with two substrates (Lu et al., 2013) [by superimposition of Na+ bound NorM (PDB ID: 4HUL) with Ethidum (PDB ID: 4HUM) and Rhodamine 6G (PDB ID: 4HUN) NorM substrate complexes], best five leads from VS and free NorM (by removing Na+ in PDB ID: 4HUL) (apo), and Na+ free rhodamine 6 g complex. A very long simulation (~75 ns) was run in all-atom (AA) molecular simulation in membrane environment for allowing the system to rearrange. Before insertion of all nine protein systems in membrane, TM helices region of NorM was estimated from OPM server (Lomize, Pogozheva, Joo, Mosberg, & Lomize, 2012) and was oriented along the z-axis of lipid bilayer. Based on cross-sectional area of protein, 128 POPC (palmitoyl-oleyl-phosphatidyl-choline) lipid was inserted in the top and bottom leaflet of bilayer and water thickness of 20 Å on both sides of bilayer was maintained. Further, solvation with TIP3P water molecule and neutralization with 0.15 M NaCl ion concentration have been performed. All simulation system was constructed using CHARMM-GUI Membrane Builder web tool (Jo, Kim, & Im, 2007;Jo, Kim, Iyer, & Im, 2008). The simulation was performed with GPU enabled NAMD2.9 program (Kalé et al., 1999;Phillips et al., 2005;Stone et al., 2007) using CHARMM36 force field parameter for protein, lipid molecule (POPC) (Brooks et al., 2009;Klauda et al., 2010) and TIP3P parameter for water molecules (Jorgensen & Jenson, 1998). Ligands were also parameterized with CGENFF force field (Vanommeslaeghe et al., 2010). All simulation system was minimized in solvated environment using conjugate gradient minimization of 5000 cycles, while keeping the protein and ligand fixed by harmonic restraint of 100 kcal/mol. This was followed by an AA conjugate gradient minimization of the entire system for 10,000 steps. Further, 25 ns of MD simulations were carried out for gradual equilibration of all nine assembled systems in six consecutive steps with various restraints applied on the protein, water, ions, and lipid molecules such as (i) harmonic restraints to ions and heavy atoms of the protein, (ii) repulsive planar restraints to prevent water from entering into the membrane hydrophobic region, and (iii) planar restraints to hold the position of head groups of membranes along the Z-axis (Jo et al., 2007). These restraint forces were slowly reduced as the equilibration progressed. Finally, 50 ns of MD simulations with the 2 fs time step size were carried out for all equilibrated system. Langevin dynamics was used for temperature control with the thermostat set at 303 K. The Nose-Hoover Langevin piston pressure control was used to control fluctuations in the barostat, which was set at a pressure of 1 bar (1 bar = 100 kPa) (Berendsen, Postma, van Gunsteren, DiNola, & Haak, 1984). Here, the periodic cell was constrained to remain tetragonal, but the cell parameters were allowed to vary. The SHAKE algorithm was used to constrain all bonds between hydrogen and heavy atoms. A dielectric constant of 1 was used for the electrostatic interactions, which was calculated using the particle mesh Ewald (pme) method (Essmann et al., 1995). The grid in the x, y, and z directions used for the particle mesh Ewald method was calculated based on rectangular simulation box size. The van der Waals interactions were described using a Lennard-Jones function multiplied by a cubic spline switching function starting at 10 Å and stopping at 12 Å. The cut-off radius for including atoms in the nearest-neighbor list was 16 Å. All 1-2 and 1-3 interactions were excluded, and 1-4 interactions were scaled by multiplication with a predefined factor. The bonded interactions were calculated every time step, the non-bonded interactions were calculated every other time step, and the electrostatic interactions were calculated every fourth time step. The nearest-neighbor list was updated every 20 time steps. For every 2 ps interval, snapshots were written to the trajectory file for subsequent analysis. All MD simulations were performed on 8-GPU cluster built in CAS in Crystallography & Biophysics, University of Madras, Chennai, India.

Analysis of MD trajectories
Simulated trajectories of NorM systems were analyzed for membrane bilayer properties such as area per lipid using VTMC program for last 10 ns (Mori, Ogushi, & Sugita, 2012) and bilayer thickeness for overall 50 ns using MEMBPLUGIN of VMD program (Guixà-González et al., 2014;Humphrey, Dalke, & Schulten, 1996). Further, structural deviation with respect to equilibrated conformation in terms of overall Cα atoms, RMSD, RGYR, and RMS fluctuations in residues using EUCB tool (Tsoulos & Stavrakoudis, 2011) were analyzed. Chimera (Pettersen et al., 2004) and VMD (Humphrey et al., 1996) were used for visualization of MD trajectories. For understanding of energetic and dynamic behavior of binding of ligands, last 10 ns harmonic simulated trajectories have been used for analyzing hydrogen and hydrophobic interaction, Dynamic cross-correlation residual matrix and clustering of conformation based on RMSD and binding free energy calculation for protein-ligand complexes.
For inter molecular hydrogen bond and hydrophobic contact analysis, last 10 ns simulated trajectory was used. Hydrogen bond analysis was performed using 3.6 Å (D-H…A) distance cutoff and 30°cutoff for the (D-H… A) angle. Hydrophobic contact was analyzed using 6.5 Å cut-off around ligand binding region of NorM transporter. Principal component analysis (PCA) of simulated trajectory was performed to understand collective structural motion in NorM transporter using Prody software package (Bakan, Meireles, & Bahar, 2011) and projected using NMWIZ in VMD (Bahar, Lezon, Bakan, & Shrivastava, 2010). Due to higher flexibility in loop region, only Cα atoms of TM helices were used for PCA calculation. Calculation of Dynamic cross-correlation matrix (DCCM) of TM helices was performed using tcl script (Cunningham, 2012), helical movement of TMs was analyzed using Trajelix analysis module (Mezei & Filizola, 2006) available in Simulaid program (Mezei, 2010) and clustering of conformation has been performed based on RMSD in Chimera mdtool packages (Pettersen et al., 2004).
The binding free energy calculations have been carried out for converged trajectories of last 10 ns using Molecular Mechanics Generalised Born Implicit Membrane solvation (MMGBIM) approach implemented in CHARMM (Brooks et al., 2009;Spassov, Yan, & Szalma, 2002) without entropy. The binding free energy referred as enthalpies of ligand and protein was computed in terms of membrane solvation energy, calculated using Generalised Born Implicit Membrane salvation model and desolvation (gas) phase energy using molecular mechanics force field. A dielectric constant of 2 was assigned to protein, and protein-solvent surface was defined as a set of atomic radii (Nina, Beglov, & Roux, 1997) and membrane slab thickness of 30 Å. Further, decomposition of binding free energy was carried out in residual basis. Nonpolar contribution in residual free energy was not considered due to negligible changes.

Results and discussion
Docking studies of diverse substrate of NorM transporter From the crystallographic studies carried out in literature it is confirmed that the NorM drug resistant substrates like rhodamine-6G [R6G] (MW: 442), ethidium (MW: 314), TPP (MW: 339) belonging to neoflavonids, benzenoids, and aminoglycosides, respectively, bind to Asp41, Ser61, Phe265, Gln284, and Ser288 substrate binding site of central cavity. Before understanding of binding interaction of any molecules to NorM transporter, validation of docking protocol has been carried out by cocrystallised substrates such as rhodamine 6 g and TPP. Molecular docking studies revealed that orientation of binding mode of cocrystallised substrate was almost similar in central cavity with little bit structural drift in terms of RMSD (heavy atoms) which has been given in Supplementary Figure S1. Docked complex of rhodamine 6 g shows rmsd of 2.89 Å due to more rotational flexibility in ethoxycarbonyl phenyl substitutent in main core of rhodamine 6 g which showed only 0.56 Å only. Similarly in docked TPP complex also showed rmsd of 1.40 Å, we believe due to different in position of arsenic atom (cocrystal) and phosphate atom ( (Table 1). But high molecular weight substrates such as Streptomycin (MW: 581), Tigecyline (MW: 585) belonging to aminoglycosides and glycylcyclines, respectively, confirm the binding of these to the bottom of central cavity located towards cation binding site and having interaction with Glu261 and Note: Bold colored residues are known substrate binding residues from co-crystallised complex study.
Tyr294 residues (Supplementary Figure S2(b)). In addition, similar movement of TPP towards cation binding site was observed in the dynamic simulation studies (Leung et al., 2014). Summary of docking scores, glide energy, and binding interaction residues of NorM substrates are tabulated (Table 1).

VS of NorM transporter
VS against NorM transporter has been carried out in two phases, first in semiflexible docking and second in flexible docking with natural compound collections. In semiflexible docking, three different docking methods with increasing accuracy for predicting binding affinity have been used for filtering ligands from VS as described in Flow Chart 1. Furthermore, flexible docking has been performed for top-50 lead compounds and docking results have been analyzed in terms of docking score, glide energy, and binding interaction. The top-five lead compounds which have the most negative docking score and glide energy and interactions with central cavity residues are presented ( Table 2) and their interacting residues were compared with substrates of NorM (Table 1). All lead compounds show better docking score and similar hydrogen and hydrophobic interactions with central cavity residues compared to the substrates bound transporter. The spatial binding of identified leads in central cavity tunnel is represented in Supplementary  Figure S3(a). Further chemical similarity of these top hit lists were analyzed and given below. Leads 3, 4, and 5 contain common naringenin flavanone chemical moiety and only difference in substituent of rhamnose sugar moiety or conjugated with gallate moiety via glycosidic bond (Figure 2). Leads 4 and 5 binds more favorable in substrate binding region in central cavity via common hydrogen bonding pattern with Asp41, Ser60, Gln284, Asp355, and Asp356. They also have common hydrophobic interaction with Ser61, Ala64, Ile68, and Phe265 residues in naringenin chemical moiety (Supplementary Figure S3(b)). In particular, lead-5 have more favorable docking score (−11.94 kcal/mol) and glide energy (−73.15 kcal/mol) compared to lead 4 as well as substrates (Tables 1 and 2).
Leads 1 and 3 have common gallate moiety substituted with rhamnose sugar (Figure 2) which have similar binding towards cation binding site in bottom of central cavity and have common hydrophobic interaction with Ile68, Tyr258, Glu261, Tyr294, Met295, Asp377, and Tyr385 (Supplementary Figure S3(b)). Naringenin moiety of lead 3 binds with Gln34, Ile37, Val40, Met44, Tyr67, Met71, Leu156, and Tyr159 via hydrophobic interaction (Supplementary Figure S3(b)). In particular, lead 3 has good docking score (−11.22 kcal/mol), glide energy (−72.81 kcal/mol) compared to other leads and substrates, and it binds more favorably to bottom of central cavity towards cation binding site with Tyr159, Tyr258, and Glu261 via hydrogen bonding and has several hydrophobic interactions with central cavity covering both N-and C-domain halves residues. Thus, lead 3 showed to be the best compared to others due to interaction with cation binding residues which are important for coordinating Na+ in coupling of Na+/drug efflux. Above docked substrates mostly interacted with substrate binding site residues such as Asp41, Ser61, Phe265, Gln284, and Ser288 in central cavity. But few of them such as norfloxacin, ciprofloxacin, streptomycin, and tigecycline also interacting with cation binding residues Glu261, Tyr294. Due to lackness of any reported inhibitor against NorM transporter, in our computational study for identifying lead compounds, comparison of binding interaction with substrate only had been performed and it reveals Asp41, Ala45, Ala57, Ser61, Ala64, Ile68, Glu261, Ala262, Ala264, Phe265, Gln284, Gly291, Ile292, Met295 Note: Bold colored residues are known substrate binding residues from co-crystallised complex study.
the importance of N-term cavity lining residues such as Glu34, Ser60, Ile68, and C-term cavity lining residues Tyr258, Glu261, Tyr294, Met295, Asp377, and Tyr385 in inhibitor binding which were not reported in known substrates.

MD simulation for validation of binding in top-five identified leads
The docked conformation of lead molecule and crystallographic substrate bound NorM complex gives only static binding event in which ligand and protein molecule are fixed. But molecules adopt accordingly to optimize maximum binding affinity and fluctuate in open and close forms in physiological dynamic conditions. Thus, for studying the dynamic properties, in terms of binding energetic and interactions of the five lead compounds with Na+ bound NorM transporter, MD simulations were carried out in membrane phase. For constructing the membrane phase, 128 POPC lipid were assembled in top and bottom leaflet of bilayer and further protein, ions, and water were inserted. For equilibration of simulation system, long 25 ns simulation was performed in NPT ensemble with various force restraints due to transgauche isomerization of hydrocarbon tail of lipid as well as fluids environment. Finally, 50 ns simulation in NPT ensemble was carried out. For understanding of comparison in dynamics, simulations were also carried out with Na+ free and Na+ bound NorM substrate complexes (rhodamine6G and ethidium) and free (apo protein) NorM for 50 ns. All the nine simulations showed thermodynamic stability in terms of potential energy during 50 ns simulation time. Before understanding of any structural as well as functional properties of membrane protein in simulation, physical properties of lipid bilayer are most important. Bilayer membrane characteristics such as area per lipid and thickness of lipid bilayer were analyzed at 30°C temperature and summarized in Table 3. It suggests that area per lipid in POPC bilayer during equilibration and production phase of simulation are almost equal to 64 Å 2 observable in experimental area per lipid of POPC lipid (Kučerka, Nieh, & Katsaras, 2011) along with little bit changes (±1) due to adaptation of lipids in presence of protein. Bilayer thickness of POPC lipid bilayer also matches to experimental and calculated from simulation trajectories of all nine systems.

Exploring the stability of NorM transporter and its complexes using MD simulation
For understanding the conformational stability of complexes of NorM transporter, all the nine simulation systems were analyzed in terms of structural deviation of overall fold (RMSD), radius of gyration (RGYR), and  Figure S4). It suggests that lead 1, 3, and 4 have more conformational restriction. RMSD of ligands displayed small conformational variation with respect to simulation time and showed adaption in structure for optimizing maximum affinity in central cavity. Further, superimposition of initial equilibrated and average structures (Cαatoms) from clustering of last 10 ns simulated trajectories of substrate complexes and best leads of NorM transporter, was carried out and their RMSD is summarized in Table 4. Lead 1 and 3 complexes showed  less deviation in ligand as well as overall protein compared to substrate bound complexes. Ligand surrounding residues (in tunnel lining residues) were also more conformation restricted in RMSD observation. Small perturbation in conformational stability of complexes causes smaller variation in the interaction partner residues of leads binding due to conformational adaptation in simulation ( Supplementary Figures S5 (a) and (b)). All simulation systems of NorM transporter have very least deviation in RGYR~0.3 Å suggesting that binding of substrate as well as lead compounds with protein do not disrupt the overall folding, thus confirming that the complex is stable during simulation time.

Analysis of binding interactions in NorM transporter and best leads
Hydrogen bond occupancy was calculated from last 10 ns simulation time between protein and ligand with donoracceptor distance (D…A) less than 3.6 Å and donor and acceptor angle (D-H…A) more than 30 and presented in Table 5. Due to adaptation in simulation, lead 5 complex which binds in substrate binding region showed hydrogen bond interaction with cation coordinating residue such as Glu261 (41%) as well as other central cavity residues Ser60 (94%), Ser61 (84%), Gln72 (62%), and Ile287 (88%). Lead 3 complex also has stable hydrogen bond interaction with cation binding residues (Glu261 with 99%, Tyr294 with 70%, Asp377 with 77%) as well as other residues in tunnel such as Glu34 (60%), Ser60 (99%) which makes tight binding to naringnin chemical moiety. Further, based on analysis of stable interaction with cation coordinating residues in central cavity tunnel, lead 3 compound shows to be the best drug molecule among others.
Hydrophobic interaction statistics calculated from ligand around 6.5Å cut-off distance from last 10 ns simulated trajectories ( Figure 5) revealed that substrates (Ethidium and Rhodamine 6G) mainly bind with Val269 in TM7, Glu284, Ile287, and Ser288 in TM8, Asp41, Thr42 in TM1 and Ala57 in TM2. The binding is more favorable for lead 5 with the residues Glu261, Phe265 in  TM7, Ile287 in TM8, Thr42 in TM1 and Ser61 in TM2 compared to lead 2 and 4 and substrate (Ethidium) in terms of hydrophobic interaction ( Figure 5). Similarly, lead 3 binds more favorably with Glu261, Phe265 in TM7, Ile287, Tyr294 in TM8 and Asp377 in TM10 along with Asp41 in TM1 compared to substrates and lead 1 which binds mainly with TM7 and TM10 only. Lead 3 compound formed favorable hydrophobic interactions with substrate binding region such as Asp41, Phe265, and Ile287 along with cation binding residues (Glu261, Tyr294, and Asp377) ( Figure 5). Thus, hydrogen and hydrophobic interaction studies also suggest lead 3 to be the best drug candidate.

Binding free energy analysis of NorM complexes
Binding free energy from converged simulated trajectory offers accurate estimation using free energy perturbation, thermodynamic integration, and MMPBSA/MMGBSA approaches and provides comprehensive understanding of molecular recognition in protein-ligand. Here, binding free energy calculations have been carried out for converged trajectories of last 10 ns using MMGBIM approach without entropy and tabulated ( Table 6). The total binding free energy is decomposed into components such as electrostatic, van der Waals, polar, and non-polar contribution. The substrates such as rhodamine 6 g and ethidium in presence of Na+ ion have binding free energy of −16.74 (±5.09) and −38.35 (±6.23) kcal/mol, respectively. Binding free energy of substrates was more favorable due to van der Waals contribution. The lead 5 complex has the most favorable binding free energy −50.52 (±5.39) kcal/mol compared to lead 2 having −35.71 (±5.02) kcal/mol and lead 4 having −31.93 (±3.72) kcal/mol and substrates (Table 6) due to electrostatic, van der Waals and polar contribution in binding energy. NorM complex with lead 3 has the most negative binding free energy −81.67 (±4.17) kcal/mol compared to lead 1 which has −66.79 (±6.19) kcal/mol due to electrostatic and van der Waals contribution. Thus, based on binding free energy, lead 3 compound is more favorable compared to other leads as well as substrates. Further, decomposition of binding free energy apart from non-polar contribution for central cavity lining residues has been carried out ( Figure 6). In lead 3, the important substrate binding residues such as Phe265, Ligands-NorM complexes Hydrogen bonding residues (Occupancy %)

Analysis of conformational dynamics of best identified leads and substrates in NorM transporter
The aim of conformational dynamics in NorM substrate complex is to understand the substrate triggering mechanism in terms of conformational rearrangement of TM helices in the presence of both Na+ and substrate and to compare the movement of TM helices in identified lead compounds. Normal mode analysis/essential dynamics analysis (EDA) is the foremost simulation techniques used to explore the substrate/lead recognition, conformational transitions and adaptability, etc. by identifying essential degrees of freedom in protein.
To understand the conformational dynamics in NorM transporter complexes, EDA or PCA was performed on only Cα atom of TMhelices and the complex motion and elasticity derived from the MD simulations were analyzed using ProDy program (Bakan et al., 2011). PCA analysis was omitted for the loops due to higher flexibility. The principal modes represent the direction (eigenvector) and magnitude ( (Table 7). The projection of top-PCA modes on trajectory of lead molecules were evaluated for understanding relative movement in all TM helices and compared with substrates and apo NorM (Figure 7). PCA projection depicts that TM3 forces TM4 to move towards the central cavity in case of substrate with Na+ and to move away from the central cavity in substrate only. Similarly TM1 and TM5 also moved away from the central cavity in substrate with Na+ and moved towards cavity in substrate only. Further, TM7 has less movement towards TM10 in case Figure 6. Decomposition of binding free energy of major active tunnel residues in substrates and best leads with active tunnel residues.  site of TM7 and TM8 towards TM10 tries to form/disrupt the coordination of Na+ (Figure 7). This results in opening of two halves by TM1 and TM5 in N-terminus and TM11 and TM12 on C-terminus which move away from central cavity so that drug release event occurred in the presence of Na+. The TM7 and TM8 moved in lateral axis and did not move towards TM10 in lead 3 complex. TM1, TM11, and TM12 do not move away from central cavity as the movement occurred only in helical axis. The correlated movement of TMs like 2, 3, and 4 occurred towards central cavity. In lead 5, movement of TM7 was fully arrested and movement of TM7, TM8, and TM10 (in lateral axis) occurred in independent direction (no correlated direction). TM1 andTM5 have movement in their helical axis. Further, analysis of displacement, tilt and rotation in TM helices was calculated (Table 8) and it reveals that TM1 and TM5 of lead 3 have less helix rotation compared to the abovementioned leads and substrates. Helix rotation in TM1 of lead 3 was found to be in anticlockwise while it was found clockwise in substrate and apo form. The average global helixtilt angle with respect to z-axis in TM1 and with respect to y-axis in TM2 of lead 3 was found to be less compared to substrate and free NorM-NG. Lesser helix rotation of TM7, TM8, and TM10 in lead 3 also occurred and rotation was in opposite direction compared to substrates and other leads, i.e. TM7 and TM10 showed clockwise rotation compared to substrate where it was found to be anticlockwise. The displacement with respect to y-axis (lateral axis) of TM7 and TM8 was lesser compared to all substrates and leads. Lead 5 complex also has lesser helical tilt in TM2 and TM7 compared to substrates. The DCCM which serves as dynamical fingerprint, identifies persistent interactions between residues in protein. In DCCM, diagonal elements are united, and are depicted in red diagonal line ( Figure 8). Significant cross-correlated movement was observed in apo NorM and its complex with substrates and compared with leads. The DCCM was visualized in Figure 8 where negative (anti-correlated motion) matrix element was plotted in the lower half of matrix and positive matrix elements (correlated motion) in upper half.
As above PCA projection suggest that TM3 forced TM4 for enforcing TM2 to move, correlated movement was also observed between these TM helices in DCCM matrix of substrate with Na+ or without Na+. Similarly, correlated movement was also observed in TM helices such as TM1 vs. TM5, TM11 vs. TM12, TM7 vs. TM8, TM8 vs. TM10 in case of substrate/Na+ coupled NorM complexes ( Figure 8). Matrix plot also revealed that lead 3 complex has less cross-correlated movement in TM helix pair such as TM9 vs. TM12 and TM10 and TM12 compared to Na+ bound substrate complex (R6G_NA). Similarly correlated movement was observed for TM7 vs. TM8 and TM10, TM8 vs. TM10 as suggested in the above PCA projection. It also gives a strong correlated movement of TM8 vs. TM9 and TM10 vs. TM11, which may provide the conformational restriction in TM 7, 8 and 10 and TM11 and TM12 in lead 3 binding for Norm transporter. In lead 5 complex, correlated movement of TM7 vs. TM8 was found to be negligible and it has lesser cross-correlated movement of TM10 with respect to TM7 and TM8 compared to substrates and apo form. But we observe little bit of movement of lead 5 ligand towards cation binding site during simulation whereas lead 3 doesn't move at all in 50 ns simulation time. Thus, it can be suggested that lead 3 will be better drug candidate in binding.

Analysis of Na+ coordination with cation binding residues in NorM
Previous studies suggest that during drug release, flipping of Glu261 brings Na+ ions in cation binding site from central cavity and makes coordination with Glu261, Tyr294, and Asp377 which provide stabilizing interactions that anchor the cation within the binding site (Lu et al., 2013). Comparison of coordination distance distribution function of these three residues with Na+ has been performed for Na+ bound complex (substrates (rhodamine 6 g, ethidium) and best three compounds lead 1, 3, and 5) from 50 ns simulation trajectory. Cation binding residues like Glu261,Tyr294, and Asp377 in lead 3 is more distributed in coordination distance (2-3 Å) compared to substrates (Figure 9). Disrupting the Na+ coordination with above mentioned residues causes influx of ions resulting in the disruption of substrate interactions with binding site residues for its efflux. In the above simulations, Na+ coordination was maintained due to strong electrostatic interaction of lead 3, thus this inhibitor arrests the transport of Na+ and forms stable coordination in cation binding site.

Conclusion
The extrusion of antimicrobial and therapeutic drugs by NorM -MATE transporter has resulted in increased drug resistance in bacteria. Design and evaluation of chemotherapeutics still remain a major challenge. In this study, we identified five lead molecules for Na+-dependent NorM-transporter based on the structure-based approaches with the combination of HTVS, flexible molecular docking, and molecular simulation approach. Binding free energy calculation from simulation approach showed lead 1, 3, and 5 as best compounds. Further, based on binding free energy decomposition and hydrogen and hydrophobic statistics, Lead 3 was summarized as better compared to other leads. For understanding of false positive namely, whether this compound will be effluxed out, conformational analysis of TMs helices was performed using PCA, DCCM, and Na+ coordination distance distribution approach and results were compared with known substrates. Lead 3 was shown to be stable during the simulation and the conformational restriction in binding interaction comprising TM helices was more compared to substrate/Na+-coupled NorM complexes. Based on conformational analysis of TM helices in substrates, it can be hypothesized that the major conformational rearrangement of TM1 and TM5, TM11, and TM 12 upon Na+-induced movement of TM7, TM8 towards TM10 and in the case of lead 3 movement of these TM helices was observed in opposite direction. This lead 3 compound is a Prunin 7″-O-gallate (CID: 42607908) and has been identified from the Indian medicinal plant "Terminalia chebula" which has been studied for MDR inhibitory properties in crude extract (Acharyya et al., 2009). Finally, overall studies suggest that the compound which interacts towards cation binding residue in central cavity will have better binding for arresting of substrate efflux pathway compared to substrate binding region (as observed in lead 2, 4, and 5 which bind in substrate binding region but binding is not stable during the simulation). Good inhibitory activity of any NorM inhibitor should have interaction with Gln34, Val35, Ile37, Gly38, Phe63, Tyr67, and Ile68 as well as with cation binding residues such as Glu261, Tyr294, and Asp377 paving the way for maximum stability and affinity with NorM transporter.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2015.1132391.