Binding mode exploration of LuxR-thiazolidinedione analogues, e-pharmacophore-based virtual screening in the designing of LuxR inhibitors and its biological evaluation

Master quorum sensing (QS) regulator LuxR of Vibrio harveyi is a unique member of the TetR protein superfamily. Recent studies have demonstrated the contribution of thiazolidinedione analogues in blocking QS by decreasing the DNA-binding ability of LuxR. However, the precise mechanism of thiazolidinedione analogues binding to LuxR is still unclear. In the present study, molecular docking combined with molecular dynamics (MD) simulations was performed to understand the mechanism of ligand binding to the protein. The binding pattern of thiazolidinedione analogues showed strong hydrogen bonding interactions with the amine group (NH) of polar amino acid residue Asn133 and carbonyl (C=O) interaction with negatively charged amino acid residue Gln137 in the binding site of LuxR. The stability of the protein–ligand complexes was confirmed by running 50 ns of MD simulations. Further, the four-featured pharmacophore hypothesis (AHHD) consists of one acceptor (A), two hydrophobic regions (HH) and one donor (D) group was used to screen compounds from ChemBridge database. The identified hit molecules were shown to have excellent pharmacokinetic properties under the acceptable range. Based on the computational studies, ChemBridge_5343641 was selected for in vitro assays. The 1-(4-chlorophenoxy)-3-[(4,6-dimethyl-2-pyrimidinyl)thio]-2-propanol (ChemBridge_5343641) showed significant reduction in bioluminescence in a dose-dependent manner. In addition, ChemBridge_5343641 inhibits biofilm formation and motility in V. harveyi. The result from the study suggests that ChemBridge_5343641 could serve as an anti-QS molecule.


Introduction
Quorum sensing (QS) is a bacterial cell-cell communication process that controls the expression of target genes in response to population density (Han, Li, Qi, Zhang, & Bossier, 2010;Koch et al., 2005;Vannini et al., 2002). The phenomenon of QS in bacteria was originally described in the luminescent marine bacterium Vibrio fisheri (Lupp & Ruby, 2004;Nelson & Hastings, 1979;Yildiz & Visick, 2009). The QS process is mediated by N-acyl-l-homoserine lactone (AHL) signalling molecules (also called as autoinducers) (Marketon, Gronquist, Eberhard, & Gonzalez, 2002). Accumulation of AHL molecules at an optimal concentration activates a transcriptional factor which results in the regulation of gene expression profiles. The defined chemical signals have a crucial role in many cellular processes that include biofilm formation, virulence factor expression, antibiotic resistance and bioluminescence (Ahumedo, Díaz, & Vivas Reyes, 2010;Miller & Bassler, 2001).
In Vibrio harveyi QS system, the information encoded in the autoinducers (AIs) control the production of LuxR (master regulator of QS gene expression). Three types of AIs are used by V. harveyi QS system, that include HAI (3-hydroxybutanoly homoserine lactone), AI-2 [2s, 4s]-2-methyl, 2, 2, 3, 4-tetrathydroxytetrachydrofuran borate and CAI-2 [(3)-3-hydroxytridecan-4-one] (Kim et al., 2010). LuxR is a unique member of TetR superfamily of transcriptional regulators. Recently, the structure of many LuxR type proteins has been defined, which provides an evidence that the family of proteins contains a helix-turn-helix DNA-binding motif in its N-terminal domain and small molecule or signalling molecule binding region in the C-terminal domain (Ramos et al., 2005). The binding of small molecule in the C-terminal domain unwinds DNA from N-terminal domain and allows transcription. This type of proteins can act as both an activator and a repressor, which control the expression profile of 625 genes. However, the majority of the TetR family of proteins typically acts as transcriptional repressor that controls the expression of their own gene and as well as other genes. Additionally, the LuxR protein can bind to the multiple binding sites in its target promoters (van Kessel, Rutherford, Shao, Utria, & Bassler, 2013;Waters & Bassler, 2006). Many LuxR homologues have been found in Vibrio species that includes HapR of V. cholera (De Silva et al., 2007), SmcR of V. vulnificus (Kim et al., 2010), LitR of V. fisheri (Fidopiastis, Miyamoto, Jobling, Meighen, & Ruby, 2002) and OpaR of V. parahaemolyticus (Yiquan et al., 2012). Among these groups, HapR possess 71% identity with LuxR of V. harveyi. Functional analysis suggests that these proteins possess similar functions. The crystal structure of HapR and SmcR proteins provide information's about the ligand binding pocket (De Silva et al., 2007).
Recently, researchers have been focusing on QS systems as a therapeutic target for drug development (Galloway, Hodgkinson, Bowden, Welch, & Spring, 2011). So far, no potent ligand molecules were identified for LuxR type of proteins. Thiazolidinedione analogues were reported to be most active against AI-2 QS system of V. harveyi (Brackman et al., 2013). But the mechanism of binding is still poorly understood. In the present study, LuxR-thiazolidinedione docked complexes were used to generate the energy-optimized pharmacophore (e-pharmacophore) hypothesis. The best pharmacophore model was used as a query in virtual screening to identify potent hits from ChemBridge database. Further, density functional theory (DFT) and ADMET calculations were performed to determine the reactivity and drug like properties of the compounds. Finally, based on the above analysis the bestpotent compound (ChemBridge_5343641) was selected for in vitro assays. Overall, the results of the present study expected to be useful in designing of novel anti-QS drug of V. harveyi.

Materials and methods
Molecular docking and e-pharmacophore analysis were carried out in a single machine running on Intel Core TM 2 Duo processor with 2 GB RAM and 160 GB hard disk with Red Hat Linux Enterprise v5 as operating system. Molecular dynamics (MD) simulations were carried out with an Intel Xenon processor in a 2-node cluster system with 64 GB RAM.

Homology modelling and structure validation
The amino acid sequence of V. harveyi LuxR protein was retrieved from UniProt database (Accession Number Q8GBU2). The 205 amino acid residue of LuxR was subjected to Blast (Altschul, Gish, Miller, Myers, & Lipman, 1990) search in the Protein Data Bank to identify suitable templates. Sequence alignment was performed using ClustalW program (Larkin et al., 2007). MODELLER 9v7 (Eswar, Eramian, Webb, Shen, & Sali, 2008) was used to build the three-dimensional structure of LuxR using the crystal structure of SmcR, a transcriptional regulator from V. vulnificus (PDB: 3KZ9). The protein models containing all non-hydrogen atoms were obtained automatically using methods implemented in MODEL-LER. The generated models are ranked on the basis of Discrete Optimized Protein Energy (DOPE), MOLPDF and GA341 scores. Model with least DOPE and MOLPDF scores was selected for further studies. The quality and stereo-chemical property of the model were evaluated using PROCHECK (Laskowski, MacArthur, Moss, & Thornton, 1993). Reliability of the model was assessed by calculating Root Mean Square Deviation (RMSD) between target and template by structural superimposition using PyMol.

MD simulation
MD simulation was performed using GROMACS v4.5.5 package (Groningen Machine for Chemical Simulation) with GROMOS96 43A1 force field (Scott et al., 1999). The calculated charges and topology files for the ligand were generated using PRODRG web server. Initially, the system was relaxed to eliminate unwanted atomic contacts that may lead to unstable MD simulation. The protein-ligand complex was solvated in a cubic simulation box under Simple Point Charge 216 (SPC216) water environment (Berendsen, Grigera, & Straatsma, 1987). Five appropriate Na + counter ions were added by replacing the water molecules to obtain a neutral system. Under equilibrium period, a free MD simulation of 10 ps was simulated with a force constant of 1000 kJ mol −1 mol −2 . In order to remove the unwanted van der Waals contacts, the system was minimized with steepest descent minimization followed by position restrained (PR) of the ligands and convergence criterion of 100 kcal mol −1 Å −1 , followed by steepest descent without PR, conjugate gradients and finally, quasi Newton Raphson until an energy of 1.0 kcal mol -1 .Å -1 . After energy minimization, the system was equilibrated with isothermal-isochoric (NVT) by applying the reference temperature of 300 K using Berendsen thermostat with a coupling time of .1 ps. Then in the isothermal-isobaric (NPT), the pressure was maintained with the reference pressure of 1 bar and time constant of .5 ps. The minimized structures were subjected to MD simulations in two steps. Initially for 500 ps with PR was performed at 300 K for the entire system, except the water molecules, in order to enable a balance of the solvent molecules around the residues of the protein. After equilibration, the final MD simulation was carried out for 50 ns without restriction. All covalent bond lengths involving hydrogen atoms were constrained with LINCS algorithm (Hess, Bekker, Berendsen, & Fraaije, 1996). The Particle Mesh Ewald (PME) method (Darden, York, & Pedersen, 1993) was employed to enlighten the electrostatic interactions. SETTLE algorithm was used to constrain the geometry of all covalent bonds containing water molecules. Analyses such as, RMSD, Root Mean Square Fluctuation (RMSF), Solvent Accessible Surface Area (SASA), Radius of gyration (Rg) and number of inter-molecular hydrogen bonding (H-bond) was calculated using scripts included with the GROMACS distribution. The trajectories were graphically visualized using OriginPro for graphical representation of the results.

Protein preparation
The modelled protein was prepared using Protein Preparation Wizard utility from Schrödinger, LLC, New York, NY, USA (Protein Preparation Wizard, 2014). Polar hydrogen atoms were added to the parent carbon atoms. Side chains that do not participate in the formation of salt bridges were neutralized. Correct bond orders and formal charges were assigned. The protein was subjected to energy minimization using the Optimized Potentials for Liquid Simulations (OPLS)-2005 force field with an implicit solvation model (Jorgensen, Maxwell, & Tirado Rives, 1996). Positions of hydroxyl, thiol hydrogens, protonation states, tautomers of His residues, Chi "flip" assignments for Asn, Gln and His residues were selected by the protein assignment script. Restrained minimization was terminated until the average RMSD of non-hydrogen atoms reaches .3 Å.

Ligand preparation
Six thiazolidinedione analogues which showed anti-QS activity against V. harveyi with EC 50 value in the low micro molar range was used in this study (Brackman et al., 2013). All the selected ligands and hits retrieved from ChemBridge database were prepared using LigPrep module of Schrödinger (LigPrep, 2013). Hydrogen atoms were added and bond orders of the ligands were fixed. The most relevant ionization and tautomeric states were generated between pH 6.8 and 7.2 using Epik module. Appropriate charges were assigned and the ligands are then neutralized. Energy minimization was performed using OPLS-2005 force field with a default setting. For each ligand, low energy ring conformations were generated and optimized ligands were further used for docking studies.

Molecular docking
Molecular docking (Glide, 2014) was performed to determine the probable binding mode of thiazolidinedione analogues in the binding site of LuxR. It examines the complementary of protein-ligand interactions using a grid-based method, after the empirical ChemScore function. Amino acid residues such as Phe-75, Phe-78,  were reported to play a crucial role in ligand recognition (Kim et al., 2010). In the present study, Sitemap protocol of Schrödinger was used to assess the druggability site for the LuxR protein. From the result, best site was selected based on site and druggability scores. The amino acid residues which lie in the selected active site perfectly coincide well with the predicted active site residues of SmcR. Thus, the grid was generated at the centriod point of the active site using Receptor Grid Generation panel with default parameters. Two docking protocols were used in this study. First, in receptor rigid docking the internal geometry of the receptor is fixed while the ligand is flexible to explore an arbitrary number of torsional degrees of freedom in addition to the six degrees of freedom spanned by the translation and rotational parameters. Second, a mixed molecular docking/dynamic protocol called Induced Fit Docking (IFD) was used (Wang, Aslanian, & Madison, 2008). In this method, both receptor and ligand are considered flexible during the docking process. The protocol generates a diverse ensemble of ligand poses; the procedure uses reduced van der Waals radii and an increased Coulomb-vdW cutoff, and can temporarily remove highly flexible side chain during docking. For each pose, a prime structure prediction was used to accommodate the ligand by reorienting nearby side chains. The residues and the ligand are then minimized. Finally, each ligand was redocked into its corresponding low energy protein structure and the resulting complexes were ranked according to the glide score. Accuracy is ensured by Glide's superior scoring function and prime's advanced conformational refinement. Scoring functions such as glide score and glide energy were employed to estimate the proteinligand complexes. Protein-ligand interactions were analysed using Glide pose viewer tool. The pose with better glide score and maximum number of hydrogen bond interactions was considered for MD simulation.

Density functional theory
DFT computations were performed using Jaguar v8.3 module, implemented in Schrödinger (Jaguar, 2014). Thiazolidinedione analogues and hits retrieved from Chem-Bridge database were used as inputs for DFT calculation. Complete geometry optimization was carried out using Hybrid DFT with Becke's three-parameter exchange potential and Lee-Yang-Parr correlation functional (B3LYP), using basis set 3-21*G level (Becke, 1993;Lee, Yang, & Parr, 1988). Energy calculations were performed using the Poisson-Boltzmann Finite solvation method. Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) and its corresponding energy gaps were calculated.

e-pharmacophore analysis
The e-pharmacophore hypothesis was generated using the protein-ligand complexes via Pharmacophore Alignment and Scoring Engine (PHASE) module of Schrödinger (Dixon, Smondyrev, & Rao, 2006). The common pharmacophoric features such as hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic region (H), positive inoizable (P), negative inoizable (N) and aromatic ring (R) were utilized to rank the pharmacophoric hypothesis. The final hypothesis was selected based on ranking and quantification process. Sites less than half of the heavy atoms contributing to the pharmacophoric features were excluded from the final hypothesis. Thus, if only two heavy atoms in a six member ring exhibit energetic interactions, the ring is not considered as a pharmacophoric feature. In order to identify potent hits against LuxR, the best hypothesis was selected to screen compounds from ChemBridge database. The final hits were ranked based on the fitness score, ligand conformer that matches the hypothesis based on RMSD site matching, vector alignments and volume terms. The pharmacokinetic properties of the hits were analysed (Qikprop, 2011). Finally, two successive hits were subjected for toxicity prediction using PROTOX web server (Drwal, Banerjee, Dunkel, Wettig, & Preissner, 2014).

Prime MM-GBSA analysis
The prime MM-GBSA (Molecular Mechanics, The Generalized Born Model and Solvent Accessibility) was performed to calculate the ligand binding energies, ligand strain energies for a set of ligand molecules and a receptor (Prime, 2014). It works with the combination of OPLS-AA, Molecular Mechanics Energies (EMM), SGB solvation model for polar solvation (GSGB) and a nonpolar solvation (GNP), that composed of different non polar SASA and van der Waals interactions. The free energy changes upon binding were calculated using the following equations.
The G complex represents complex energy, G protein is the receptor energy and G ligand is the unbound ligand energy. EMM represents molecular mechanics energies, GSGB is an SGB solvation model for polar solvation and GNP is a nonpolar solvation term. The Spearman's rank correlation co-efficient is used to quantify the correlation that exists between the predicted binding free energy (ΔG bind ) against the docking score of six thiozolidinedione analogues. The Spearman's rank correlation co-efficient is calculated as: 2.9. In vitro assays 2.9.1. Bacterial strain The aquatic bacterial pathogen V. harveyi (MTCC 3438) was kindly gifted by Dr. A. Veera Ravi, Department of Biotechnology, Alagappa University, Karaikudi, India. The V. harveyi was cultured in Luria Bertani (LB) media (pH 7.5 ± .2) and incubated overnight at 30°C. For experimental purpose, the overnight culture was inoculated into fresh LB liquid medium to an initial OD 600 of .5 and the culture was aerated on an orbital shaker at 300 rpm at 30°C.

Determination of minimum inhibitory concentration of compound
The selected compound (1-(4-chlorophenoxy)-3-[(4, 6dimethyl-2-pyrimidinyl) thio]-2-propanol) (Chem-Bridge_5343641) was purchased from ChemBridge online chemical store (http://www.hit2lead.com). The stock was prepared by dissolving 5 mg of compound in 155 μl of 100% sterile dimethyl sulfoxide (DMSO), whereas the working standard was freshly prepared on the day of the experiment by diluting (20 μl) the stock solution with 80 μl of DMSO. Determination of Minimal Inhibitory Concentration (MIC) of the compound was performed using the method of Nithya & Pandian, 2010. LB broth medium was added to the 24-well microtitre plate to which the bacterial suspension (adjusted OD 600 of .5) was inoculated and then supplemented with twofold serially diluted solution of the compound. The assay plate was incubated at 30°C for 24 h. The complete suppression of visible growth achieved through the lowest concentration of the compound was noted as MIC. Further, investigation was performed with sub-MIC concentration of the compound.

Bioluminescence inhibition assay
The V. harveyi cells enriched in Alkaline Peptone Water (APW) was treated in the presence (.5-8 μg/ml) or absence of test compound for 12 h at 30°C. The bioluminescence intensity was measured in Relative Light Units (RLU) by a luminometer. The untreated control culture of V. harveyi was used as positive control, whereas the sterile APW was considered as negative control.

Effect of compound on biofilm development
The biofilm formation of V. harveyi was initially quantified using crystal violet staining as described (Musthafa, Saroja, Pandian, & Ravi, 2011). Briefly, the overnight culture of V. harveyi (OD adjusted to .4 at 600 nm) was treated with the presence (.5-8 μg/ml) or absence of test compound for 16 h. After incubation, the plate was sequentially washed thrice with sterile deionized water to remove loosely associated bacteria. Subsequent to washes, the cells were stained with 200 μl of .4% (v/v) crystal violet solution for 15 min at room temperature. After staining, the plates were rinsed thrice with deionized water to remove the excess stain and then 1 ml of absolute ethanol was added and incubated to extract the crystal violet solution. The absorbance was measured at 600 nm using UV-Visible spectrophotometer.

Percentage inhibition of biofilm biomass
¼ ½Control OD -Test OD=Control OD Â 100 2.9.5. Disintegration of mature biofilm The biofilm was incubated in the presence or absence of the test compound to evaluate the biofilm disintegration. The glass slides were incubated for 5 h at 30°C. After incubation, the biofilm formed on the glass slides was stained with 200 μl of .4% (v/v) crystal violet solution for 15 min at room temperature. Then, the glass slides were washed thrice with deionized water and air dried at room temperature. The disintegration of biofilm was confirmed through light microscope.
2.9.6. Swimming and swarming assays The motility assays such as swimming and swarming were performed according to the procedure described . For swimming assay, the plates were prepared with .3% agar, supplemented with 8 μg/ml of the test compound. For swarming assay, the plates were prepared with .5% agar, supplemented with 8 μg/ml of the test compound. About, 3 and 5 μl overnight culture of V. harveyi (OD 600 nm of .5) were incubated at the centre of the swimming and swarming agar plates and incubated for 16 h at 30°C to observe the reduction in swimming and swarming migration zones.
2.9.7. Growth curve and bioluminescence kinetic assay The overnight culture of V. harveyi was diluted with LB broth to an OD 600 nm of .5. About 250 μl of the diluted culture was dispensed in 25 ml of LB broth and then various concentrations (.5-8 μg/ml) of the test compound were added and incubated on a rotary shaker under 180 rpm at 30°C. The cell density was measured using a UV-Visible spectrophotometer over 1-h time intervals up to 16 h. LB broth containing no bacteria was used as a negative control and bacteria without compound treatment serve as a positive control. In bioluminescence kinetic assay, 250 μl of bacterial cell culture was inoculated in 25 ml of APW, supplemented with various concentrations (.5-8 μg/ml) of the test compound and incubated by following the procedure as stated above. The bioluminescence intensity was calculated every 1 h up to 16 h using luminometer. APW containing no bacteria was used as a negative control and bacteria without compound treatment serve as a positive control.

Structure determination and validation
The three-dimensional structure of LuxR was successfully predicted using homology modelling (Srinivasan, Sudha, Umamaheswari, Vanajothi, & Ramesthangam, 2012). Blast search of the LuxR protein showed maximum homology with the crystal structure of SmcR, a transcriptional regulator from V. vulnificus (PDB: 3KZ9). The sequence alignment results from ClustalW is shown in Supplementary Figure 1S. Sequence identities and similarity between LuxR and template were found to be 46 and 90%, respectively. Sequence alignment results showed that some divergent evolution exist between the sequences. Such differences exist among the transcriptional regulator since each system has been optimized to provide survival in the specialized niche in which a particular species of bacteria reside (Winans, 2006). The target and template belongs to the TetR protein superfamily. It has been proposed that TetR family of proteins consist of two domain architecture. The N-terminal domain of LuxR is required for DNA binding and a large C-terminal region for ligand binding (Kim et al., 2010). Based on the template structure (PDB: 3KZ9), a suitable model for LuxR was predicted using MODEL-LER program. A total of five models were generated and ranked based on their internal scoring functions. Structural analysis indicated that the constructed model entirely composed of nine alpha-helices. The model 5 with the least DOPE score of −26,614.88 was selected and the quality of the model was further evaluated. Structural superimposition of LuxR and template showed good RMSD of .47 Å, confirmed the quality of the predicted model. The predicted structure of LuxR is shown in Figure 1. The figure reveals that, the first three alpha-helices (14-31, 39-46, 50-56) forms a three helix bundle that contains a canonical helix-turn-helix motifs. The successive 4-9 helices (60-82, 89-106, 109-119, 125-150, 159-179, 183-196) contribute to the dimerization domain with the helix 8 and 9 in each monomer forming an anti-parallel four helix dimerization motif. Nine helices are connected through loops and overall protein folds into a compact conformation in a Ω-shaped structure. As reported previously, the LuxR protein family forms a putative ligand binding site within the dimerization domain (Ramos et al., 2005). Structural parameters of the constructed model were found to fold in a similar way as that of template. Geometric evaluations of the model performed using PROCHECK represent the distribution of phi and psi angles for the amino acid residues. Ramachandran plot showed that a total of 94.2% of residues were distributed within the most favourable region and 5.8% of residues in the allowed region. None of the residue presents in the generous and disallowed region indicates that bond length, bond angle, main chain, side chain, torsion angles and clashes between non-bonded pairs of atoms are located in the reasonable position in the constructed model.

Molecular docking
Six thiazolidinedione analogues ((Z)-5-octylidenethiazolidine-2,4-dione (Compound 1), (Z)-5-decylidenethiazolidine-2,4-dione (Compound 2), (Z)-5-undecylidenethiazolidine-2,4-dione (Compound 3), (Z)-5-dodecylidenethiazo-lidine-2,4-dione (Compound 4), 5-decylthiazolidine-2,4dione (Compound 5) and (E)-3-decylidenepyrrolidine-2,5dione (Compound 6)) exhibiting potent activity against AI-2 QS system in V. harveyi with EC 50 values in the lowest micro molar range. Hence, these six compounds were considered for molecular docking studies. The 2D chemical structures of thiazolidinedione analogues are illustrated in Supplementary Figure 2S. Molecular docking was carried out to understand the probable binding mode of thiazolidinedione analogues in the binding site of LuxR. Docking results of thiazolidinedione analogues are shown in Table 1. All the six analogues were successfully docked into the binding site of LuxR and exhibited glide score in the range from −5.08 to −5.62 kcal/mol, and glide energy in the range from −27.54 to −35.82 kcal/mol, respectively. The formations of hydrogen bonds indicate the strong interactions between protein and ligands. From the binding mode analysis, it was found that six analogues produced almost similar binding mode. The carbonyl group present in the six analogues were established hydrogen bond interaction with the help of polar amino acid residue Asn133. Flexibility of protein is considered as the most challenging issue in molecular docking due to the large size and degree of freedom of protein. In order to validate the receptor rigid docking results, an additional sophisticated and computationally intensive method referred as IFD protocol was used in this study. The IFD protocol mainly considers the flexibility of protein to predict the binding modes of analogues. The method keeps ligand and receptor more flexible during docking. Compared with receptor rigid docking results, IFD protocol showed highest glide score ranging from −10.77 to −14.04 kcal/mol for the six analogues. The glide energy and IFD score ranges from −38.43 to −49.12 kcal/mol, −468.31 to −472.40 kcal/mol, respectively. From the IFD results, it was observed that the carbonyl group present in the six analogues was established hydrogen bond interaction with the polar amino acid residue Asn133. Most interestingly, IFD protocol exposed additional interaction with the negatively charged amino acid residue Gln137. Both the residues were co-ordinated within the dimerization domain. From the results, it was observed that the carbonyl and amine groups present in the analogues plays a major role in forming hydrogen bond interactions with the binding site amino acid residues of LuxR. It was also observed that two amino acid residues, namely Asn133 and Gln137 play a key role in the binding of six analogues. Figure 2 simulations (Karplus & McCammon, 2002). In the present study, MD simulation of 50 ns was performed for the protein-ligand complexes. The exclusive use of RMSD, RMSF and number of protein-ligand contacts represents the realistic and stable state of the proteinligand complexes in the timescale studied. Analysis of  backbone RMSD from the starting structure over the course of the MD simulations was used to check the structure and dynamic properties of protein-ligand complexes. Figure 3 illustrates the RMSD profile of proteinligand complexes. In case of complex 1, the RMSD showed initial drift in the first 25 ns of the simulation period, this is due to the absence of restraining the production phase of MD simulation. After 25 ns, the protein-ligand complex attains stable conformation with an RMSD of .3 nm with the SD of .04 nm. It was noted that the complex 2 does not show much deviation during the simulation time and the average RMSD is .4 nm with the SD of .04 nm. Complex 3 shows initial fluctuation in the first 15 ns simulation with the RMSD value of .38 nm and then the RMSD slightly falls back to .22 nm, which may be due to unpacking and repacking. However, after 35 ns of MD simulation, the complex attained stable conformation with the RMSD around .4 nm with the SD of .06 nm. In complex 4, constant RMSD was observed throughout the whole simulation. The complex showed stable RMSD around .3 nm with the SD of .03 nm. Complex 5 and complex 6 showed initial fluctuation in the first 20 ns of MD simulation. In complex 5, the RMSD curve shows stability near 25-40 ns with RMSD deviation around .3 nm with the SD of .05 nm, whereas in the case of complex 6, the system relaxes stable after 25 ns with the constant RMSD of .3 nm with an SD of .03 nm. The RMSD result confirmed the stability of the protein-ligand complexes.

RMSF
In order to analyse the fluctuation rate with respect to each residue, RMSF of C-α atoms are analysed during the simulation time is shown in Figure 4. The study mainly focus on the fluctuation rate of amino acid residues located at the binding site region of LuxR, especially 7 th helix which is mainly involved in ligand binding process. The catalytic residues present in the active site of the protein were unchanged upon binding of six analogues. In complex 1, the catalytic amino acid residues such as Asn133, Gln137 has shown .19 and .18 nm of fluctuation. However, a few other amino acid residues including Met1 (.60), Arg11 (.51 nm), Arg202 (.61 nm), Glu203 (.57 nm), His204 (.77 nm) and His205 (.96 nm) were showed higher fluctuations. The residual fluctuation was found pronounced in loop regions that correspond to the N-terminal and C-terminal domains of LuxR. In complex 2, the catalytic amino acid residues have shown reduced fluctuation rate of .18 nm for Asn133 and .12 nm for Gln137 and the major RMSF can be seen in the amino acid residues such as Met1 (.6 nm), Arg36 (.45 nm), Glu203 (.43 nm), His204 (.48 nm) and His205 (.54 nm). The result indicates that the binding of inhibitors with the protein plays a crucial role in the flexibility of the protein. In complex 3, major conformational changes were observed in the amino acid residues such as Met1 (.71 nm), Gly2 (.51 nm), Ser3 (.43 nm), Arg9 (.42 nm), Arg11 (.61 nm), Arg202 (.48 nm), Glu203 (.57 nm) and His204 (.76 nm). These residues are insignificant for the study since the major dynamic changes with the ligand were absent in the catalytic site of the protein. The results indicated that these amino acid residues present in the protein are more unstable during MD simulation. The interacting amino acid residues including Asn133 and Gln137 were shown the fluctuation rate of .15 and .14 nm, respectively. In complex 4 and complex 6, it was noted that the catalytic The overall fluctuation rates of 7 th α-helix for the complex 1-complex 6 are .04, .03, .03, .02, .03, and .03 nm, respectively. The result clearly indicates that all the major fluctuations were located far away from the active site. The fluctuations observed between protein-ligand complexes demonstrate the overall stability of protein was not affected upon ligand binding in the binding site of LuxR.

SASA
The SASA was calculated to predict the magnitude of ligand-induced conformational changes in the proteinligand complexes. The overall free energy solvation (ΔG) of the protein-ligand complexes is shown in Figure 5(a). It clearly defined that the binding of ligand molecules into the protein binding site does not produce any significant changes in the SASA value. Complex 1complex 6 shows the average solvation value of 306. 45,306.33,308.19,306.34,313.37 and 305.96 nm with the SD of 15.27,15.25,11.46,11.80,11.20 and 13.63 nm, respectively. Comparing the SASA value of the proteinligand complexes, complex 5 has shown higher value of SASA with time, while other complexes (complex 1-4 & complex 6) have shown lower SASA value. Low variation in the SASA plot indicates higher thermodynamic stability. The results also indicate that the complex 5 undergoes large conformational changes upon ligand binding. The SASA plot is more comparable with the RMSF plot.

Hydrophobic interaction
The nature of hydrophobic interactions between proteinligand complexes is investigated and shown in Figure 5(b). The figure describes that the protein-ligand complexes have shown similar hydrophobic contacts during the simulation. The average hydrophobic contacts of the protein-ligand complexes are found to be 58. 06, 58.21, 57.89, 58.06, 58.70 and 57.84 nm with the SD of 2.3, 2.1, 1.5, 1.6, 1.5 and 2.2 nm, respectively. From the result, it was noted that hydrophobic contacts observed between protein-ligand complexes are similar in nature and these contacts make protein-ligand complexes more stable.

Rg
In order to determine the compactness of the proteinligand complexes during the simulation time, the Rg was calculated. Figure 5(

H-bond interaction
Since, hydrogen bond plays an important role in structural conformation and enzyme inhibition in the macro molecules. It also plays a major role on controlling the steady conformation of protein. Time-dependent hydrogen bond formations during simulation were observed in order to comprehend the relationship between hydrogen bond formation and flexibility. The constant hydrogen bond interactions observed throughout the simulation period of 50 ns are displayed in Figure 6. H-bond analysis revealed that a maximum of four hydrogen bond interactions are observed in the protein-ligand complexes during the simulation time. The H-bond interactions observed between the amino acid residues (Asn133 & Gln137) and analogues in docking studies were found stable. Overall, these bonds makes protein-ligand complex more stable.

Distance analysis
The distance of thiazolidinedione analogues from the key amino acid residues (Asn133 & Gln137) in the binding pocket were measured throughout the simulation. The hydrogen bond distance between protein-ligand complexes was plotted and is shown in Supplementary Figure 3S. Amino acid residue Asn133 binds with the carbonyl group present in the six thiazolidinedione analogues with an average distance of 2.53, .66, .90, .99, .94 & 1.54 nm and SD of 1.37, .12, .12, .17, .15 and .24 nm with minor periodic deviations. Whereas, the amino acid residue Gln137 binds with the carbonyl group of thiazolidinedione analogues with an average distance of 2.70, .56, .92, .56, .89 & 1.32 nm and SD of 1.51, .07, .10, .12, .13 and .19 nm, respectively. The hydrogen bond interactions observed in the docking studies were found stable throughout the whole simulation. The results indicated that the observed hydrogen bond interactions are essential for protein-ligand binding and stabilization of the complex.

Prime MM-GBSA
Prime MM-GBSA was calculated to evaluate the affinities between the selected hits with LuxR protein (protein-ligand complexes) generated from molecular docking. The results are tabulated in Tables 1 and 5. The Spearman's correlation co-efficient of r s = .71 demonstrates good correlation exists between the two variables. Thiazolidinedione analogues produced a favourable binding free energy in the range from −93.79 to −119.06 kcal/mol and the hit molecules obtained from ChemBridge database showed the free energy of binding in the range from −35.64 to −78.53 kcal/mol. The result provides the additional evidence indicating the potent application of hits towards LuxR of V. harveyi.

Density functional theory
Frontier orbital energies including HOMO and LUMO are used to characterize the electronic properties of the compounds. LUMO is directly associated with the electron affinity and system's tendency to accept electron density, whereas HOMO is directly associated with the ionization potential (AI-Sabagh et al., 2013). There are many reported evidence on using 3-21*G basis set and the results obtained are more reliable (Reddy, Singh, & Figure 6. H-bond variations in protein-ligand complexes during the simulation time. Singh, 2014; Sindhu & Srinivasan, 2015;Tripathi & Singh, 2014). The HOMO and LUMO distribution, energies and energy gaps for the hit molecules were computed. Energy gap between the HOMO and LUMO is used to determine the chemical stability of a molecule. The HOMO, LUMO and energy gap values for the thiazolidinedione analogues and hits molecules are given in Tables 2 and 3. The energy gap observed between HOMO and LUMO becomes smaller, ranging from .26 to .20 eV and .06 to .00 eV respectively, indicating the fragile nature of the bound electron. A small gap between the two frontier orbitals implies high reactivity, high polarizability and low stability of the compound. However, the calculated HOMO_LUMO gap values are small for thiazolidinedione analogues and hit molecules. Thus, rapid electron transfer and exchange are equally possible by making these molecules extremely more reactive. The molecules with small gap are more polarizable than the molecules with the large gap. This helps in the reaction of these compounds with the active site residues of the receptor.

e-pharmacophore mapping and virtual screening
The six protein-ligand complexes were selected to develop the e-pharmacophore hypothesis. Results revealed four pharmacophoric features, namely two hydrophobic regions (H), one acceptor (A) and one donor (D) group. Based on the docking results, the common features such as two hydrophobic, one acceptor and one donor regions were selected and their distance is shown in Figure 7. The pharmacophoric distances are shown in Angstrom. The ranking of hypothesis and pharmacophoric features are shown in Table 4. All these pharmacophoric features are considered to discover potent LuxR hits using ChemBridge database. After screening, 14 compounds were selected based on the fitness value of ligand to the pharmacophore hypothesis.
The fitness values for the identified hits range from 1.6 to 1.5. The docking result suggest that strong proteinligand complex was observed with the glide score and glide energy of −7.0 to −7.9 kcal/mol and −29.4 to −48.6 kcal/mol, respectively. The ligands were further analysed for their interactions with key amino acid residues in the LuxR binding site. The glide score, glide energy and their interacting residues of the hits are shown in Table 5. The binding mode of the top two hits with their hydrogen bonding interactions are shown in Supplementary Figure 4S   The interaction analysis revealed that the two amino acid residues such as Asn133 and Gln137 interact with almost all hits molecules identified from e-pharmacophore-based virtual screening. It was concluded that the two amino acid residues play a significant role in the protein-ligand recognition process. The ADME properties of the hits are presented in Table 6. Through computational prediction, all the hits are in the acceptable range of the pharmacokinetic parameters indicating their potential as drug like molecules. The identified hits provide an excellent idea to select the compounds for the development of drug against LuxR of V. harveyi. PROTOX server revealed no toxicity-related fragments present in the hit molecule. On the basis of pharmacokinetic and toxicity properties, the hit (ChemBridge_5343641) was taken for experimental studies. The 2D structure of Chem-Bridge_5343641 is shown in Supplementary Figure 5S.

Minimal inhibitory concentration
In drug designing, the first relatively unexplored is the identification of molecular drug target that is relevant to the disease (Greer, 2000). While the second relatively unexplored area is to regulate the QS among pathogenic microorganisms using small molecules (Choi et al., 2012). In this study, the ability of the small molecule on QS was tested using V. harveyi. The test compound inhibited the growth of V. harveyi with the least concentration of 16 μg/ml. Therefore, all further experiments were performed at mean sub-MIC concentration (.5, 1, 2, 4, 8 μg/ml) of the compound. DMSO which was used as negative control in all the bioassays does not show any inhibition effect in the test organism.

Bioluminescence inhibition assay
In V. harveyi, QS circuits control the production of several virulence factors and the expression of bioluminescence (Tala et al., 2014). Therefore, the anti-QS properties of the test compound were assessed using V. harveyi. The inhibitory effect of the test compound on the bioluminescence of V. harveyi is shown in Figure 8. The result demonstrated that the compound exhibited a concentration-dependent inhibition of bioluminescence in V. harveyi. From the results, it was also observed that the test compound found to inhibit the enzyme involved in the QS-based biofilm formation in V. harveyi. Our result fall in line with report (Defoirdt, Pande, Baruah, & Bossier, 2013), in which pyrogallol showed significant   (Federle & Bassler, 2003). As discussed in multiple reports, the interspecies auto inducer AI-2 has been involved in the regulation of bacterial physiology, production of virulence factors and formation of biofilm in various bacteria (Guo, Gamby, Zheng, & Sintim, 2013). Also, AI-2 plays a promotive role in the formation and maturation of biofilm (Huang, Li, & Gregory, 2011).    Therefore, it is hypothesized that the compound disrupts bioluminescence due to the inhibitory effect on biofilm formation. The effect of test compound on biofilm formation of V. harveyi was investigated and shown in Supplementary Figure 6S. Inhibition of biofilm formation by test compound was found to be in a dose-dependent manner. The compound produced an inhibition effect on biofilm (64%) at a concentration of 8 μg/ml and has no apparent effect on bacterial growth rate. The results of the current analysis are consistent with the previous report (Vikram, Jayaprakasha, Jesudhasan, Pillai, & Patil, 2010), wherein the flavonoid compounds such as quercetin and naringenin have been shown to possess potent anti-biofilm efficacy. Our results also correlate with the report (Uckoo, Jayaprakasha, Vikram, & Patil, 2015); the polymethoxyflavones isolated from the peel of miaray mandarin was shown to have biofilm inhibitory activity in V. harveyi without affecting the bacterial growth rate.

Disintegration of mature biofilm
The biofilm growing bacteria can be up to 1000 fold more resistant to antibiotics, host immune defence system than that of planktonic bacteria (Sun et al., 2013). Therefore, biofilm disintegration plays a crucial role to overcome infections that are highly resistant to antibiotics. After treatment, the compound showed a thin and disintegrated matrix of biofilm on glass slides. In addition, the control without compound treatment represents a continuous matrix with membranous structures spanning the biofilm. The light microscopic analysis has explained the ability of test compound to disintegrate the mature biofilms of V. harveyi and is shown in Figure 9.

Swimming and swarming assays
It has been demonstrated that there is link between swimming, swarming motility and biofilm formation in V. harveyi (Brackman et al., 2008). Vibrio species switches from swimming cell to swarmer cell to move in high viscous media or in surfaces. Also, it uses polar flagellum to promote biofilm formation in Vibrio species (Yildiz & Visick, 2009). Therefore, inhibiting the bacterial biofilm formation and motility by a small molecule represent an important strategy to control bacterial colonization. The effect of the test compound on swimming and swarming motility was assessed using V. harveyi. The result clearly indicated that the reduction in the swimming and swarming motility of V. harveyi was observed at the higher concentration of 8 μg/ml (Figure 10(a)-(d)). The result of current finding is consistent with the report (Paciavathy, Agilandeswari, Figure 11. The effect of ChemBridge_5343641 on the growth curve of V. harveyi. Notes: Cell density was quantified by turbidity measurement at an absorbance at 600 nm. Musthafa, Pandian, & Ravi, 2012), wherein the curcumin was found to inhibit the mobility of the reporter V. harveyi strain MTCC 3438.

Growth curve and bioluminescence kinetic assay
To determine the growth curve in the presence of the test compound, V. harveyi bacterium was used. No significant influence on the growth of V. harveyi was observed in both treated (.5-8 μg/ml) and untreated control after 16 h incubation. Figure 11 shows the antibacterial activity of the test compound at various concentrations (.5-8 μg/ml). The result clearly demonstrates that the test compound does not produce significant antibacterial activity at the test concentration. In bioluminescence kinetic assay, it was observed that the test compound reduced the bioluminescence of V. harveyi in a dose-dependent manner. The Figure 12 represents the reduction in the bioluminescence of V. harveyi by the treatment of the test compound.

Conclusion
The current study investigated the exact binding mode of thiazolidinedione analogues with LuxR by combining protein homology modelling, MD simulations and molecular docking. Based on the high-resolution structure of V. vulnificus SmcR, we generated the homology model for LuxR. Molecular docking was performed to study the binding mode of thiazolidinedione analogues in the binding site of LuxR. These analogues bind in the binding site residues and utilize hydrogen bonding interactions with two amino acids, namely Asn133 and Gln137 resulting in clue about the inhibitory mechanism about the analogues. MD simulation analysis, such as RMSD, RMSF, SASA, Rg and H-bond analysis showed small fluctuations, indicating that the docked complexes are stable over time. Furthermore, e-pharmacophore mapping was performed for the protein-ligand complexes. Based on the best hypothesis (AHHD), screening was performed against ChemBridge database. The final hits were selected based on the fitness score, hydrogen bond interactions and ADME properties. The identified potent hit molecule was validated using in vitro assays. The findings provide new insights into ligand recognition by LuxR and it can be used as a reference to elucidate the structure of protein-ligand complexes. The results can also offer useful references for designing potent inhibitors targeting QS system of V. harveyi. In this study, we report for the first time the anti-quorum sensing properties of ChemBridge_5343641.