New molecular insights into the tyrosyl-tRNA synthase inhibitors: CoMFA, CoMSIA analyses and molecular docking studies

Drug resistance caused by excessive and indiscriminate antibiotic usage has become a serious public health problem. The need of finding new antibacterial drugs is more urgent than ever before. Tyrosyl-tRNA synthase was proved to be a potent target in combating drug-resistant bacteria. In silico methodologies including molecular docking and 3D-QSAR were employed to investigate a series of newly reported tyrosyl-tRNA synthase inhibitors of furanone derivatives. Both internal and external cross-validation were conducted to obtain high predictive and satisfactory CoMFA model (q 2 = 0.611, r 2 pred = 0.933, r 2 m = 0.954) and CoMSIA model (q 2 = 0.546, r 2 pred = 0.959, r 2 m = 0.923). Docking results, which correspond with CoMFA/CoMSIA contour maps, gave the information for interactive mode exploration. Ten new molecules designed on the basis of QSAR and docking models have been predicted more potent than the most active compound 3-(4-hydroxyphenyl)-4-(2-morpholinoethoxy)furan-2(5H)-one (15) in the literatures. The results expand our understanding of furanones as inhibitors of tyrosyl-tRNA synthase and could be helpful in rationally designing of new analogs with more potent inhibitory activities.

Infectious diseases caused by bacteria are known as one of the most life-threatening health problem all over the world, whose chemotherapy using antimicrobial agents and antibiotics has been a critical public health tool for nearly a century, saving millions of lives around the world 1 . However, due to the indiscriminate usage of antibiotics in particular, surviving bacteria have evolved resistance against several antibiotics in recent decades 2 . Especially with the growth of multidrug resistance in bacteria, finding new antibacterial drugs becomes increasingly crucial in global researches 3 . Interruptions of protein synthesis have long been recognized as an attractive target of anti-bacterium, and its crucial enzyme aminoacyl-tRNA (aaRS) involved in protein synthesis catalyzed the bond between specific amino acid and its cognate tRNAs 4 . Hence, an increasing number of researches focused on aaRS inhibitors as a potent antibacterial agent. The inhibitors of leucyl-tRNA synthase of icofungipen and AN-2690 are both in clinical development for the treatment of onychomycosis 5 . Mupirocin, an inhibitor of isoleucyl-tRNA, shows good effect on infectious diseases as antibacterial agents 6 . Recently, a series of furan-2(5H)-one derivatives have been found remarkable inhibitory activities against tyrosyl-tRNA synthase 5 .
The three-dimensional quantitative structure-activity relationship (3D-QSAR) models, calculated via the most widely used comparative molecular field analysis (CoMFA) 7 and comparative molecular similarity indices analysis (CoMSIA) 8 incorporating the information of conformation or spatial orientation of molecules, are used for the rational design of most potent novel inhibitors. In this study, they were used to extract the structural features favored for tyrosyl-tRNA synthase inhibitors based on the skeleton of furan-2(5H)-one. The best model, which was developed from a dataset consisted of a 44 molecule training set and a 8 molecule test set, have been validated appropriately, and 10 novel compounds designed on the basis of the model have been predicted better activity than compound 15, the most active molecule in the literatures 5,9 . Therefore, the established 3D-QSAR models of fifty-two molecules by CoMFA and CoMSIA could not only give the key structure requirements for the antimicrobial activity but also serve as a helpful guidance in design of novel antibiotics, especially for the control of drug-resistant superbugs.

Methodology
Methods and Data sets. Fifty-two furan-2(5H)-one tyrosyl-tRNA synthase inhibitors were collected from papers published by a certain research group 5,9 . The biological data of the compounds and their structures were showed in Table 1. The pIC 50 (−logIC 50 ) values converted from the origin IC 50 was required as dependent variables in 3D-QSAR analysis. The pIC 50 values covering 3 log units were considered as a homogenous and wide range dataset for 3D-QSAR studies 10 . Eight molecules of structural variety and broad range of activity in the data were randomly chosen as test set to assess predictive ability of the resulting models, and therefore the remaining forty-four molecules were selected as a training set to generate the 3D-QSAR models.
Molecular modeling and alignment. The 3D structures of the furan-2(5H)-ones were build in SYBYL 8.1 (Tripos, Inc, St. Louis, MO, USA) molecular modeling package. The force field of standard Tripos molecular mechanics along with Gasteiger-Hűckel charge were employed to perform structure energy minimization 11 . The quality of molecular alignment was considered as a key factor for the robustness and predictive power of CoMFA and CoMSIA models 4, 12 . Here we applied molecular alignment to align all the molecules by using furan ring as the common skeleton and the most active compound 15 as the template molecule. To ensure the energy level of the conformations was reasonable, a further conjugate gradient method minimization of compound 15 was conducted and the convergence was reached. The aligned molecules were shown in Fig. 1.
CoMFA and CoMSIA analysis. The descriptor fields of both methods were calculated in a three-dimensional cubic with one angstrom grid spacing. The frontier of the box extended extra 4 angstrom units from the border of aligned structures in each direction.
For CoMFA method, incorporating steric and electrostatic fields 13 , the the probe atom of a charged sp 3 hybridized carbon atom was applied to compute electrostatic and steric fields. The cutoff value was 30 kcal•mol −1 14 . As to CoMSIA approach, three extra fields including hydrogen bond acceptor, hydrogen bond donor and hydrophobic were considered 15 . A Gaussian function was also applied in calculating the similarity indices, making it accounts for all grid points 16 . The equation (1) for the similarity indices calculation is as below: where A q means the similarity index of point q, and k represent the physicochemical properties of electrostatic and steric descriptors; W probe,k means the probe atom and the attenuation factor is a default value of 0.3; i means summation index of the molecule j, and W ik is the observed value k of a specific property of the atom i 17 .
Partial least squares (PLS) analysis. In order to build a statistically significant model, PLS method 18 was introduced to correlate the both fields to the pIC 50 values linearly. Leave-one-out (LOO) method was performed first in the cross-validation in which each compound is deleted from the dataset and the activity of the "leave-one" molecule is predicted using the model build based on the remaining molecules in the dataset 19 . With the default value of column filtering at 2.0 kcal•mol −1 , the optimum number of components (ONC) was obtained on the lowest standard error of prediction (SEP), which is usually corresponds to the highest cross-validated squared coefficient (q 2 ). To avoid the model over-fit, higher numbers of component will not be accepted unless the q 2 could have risen by 10% or more 20 . Non-cross-validation was then conduct to build the final 3D-QSAR models. The correlation coefficient of LOO method (r 2 cv ) was defined by equation (2) as follows: where Y mean , Y pre and Y obs , means average, calculated and actual pIC 50 , respectively 21,22 .
Sensitivity of a PLS model. Most molecules of the data set may have "twins", which make a near twin of each left-out molecule likely remain in the training data and lead to a good prediction, therefore, the q 2 statistic may give you a false sense of confidence. Progressive scrambling was often used to determine the sensitivity of a QSAR model to small systemic perturbations of the response variable for the model's stability 23,24 .
Predictive correlation coefficient (r 2 pred ). Normally q 2 is considered as a productive but not sufficient parameter in validating the model. Sometimes models with high r 2 cv and r 2 values may be improper in many cases. The external predictive correlation coefficient (r 2 pred ) 25 was calculated to estimate the predict ability. The predictive correlation coefficient (r 2 pred ) was calculated by the following equation (3):

Molecular docking and MD simulation. Molecular docking technique is an important method in discov-
ering novel small-molecule drugs [28][29][30][31][32] . In our study, Surflex-Dock was used to perform the molecular docking. Tyrosyl-tRNA synthase's crystal structure 3P0H which contains a ligand was obtained from the PDB database 33 . The hypothetical protomol was used to probe steric and electrostatic interactions of the active pocket. H-bond acceptor and donor substituent as well as hydrophobic fields were also investigated 34,35 . Before docking, hydrogen atoms of the receptor were added first, and the Kollman-All force field was applied in the prepared structures [36][37][38] .
The definition of active site definition was performed based on the original ligand in the crystal. We chose compound 15 as the subject to dock into the active pocket under the conditions previously optimized. As to the newly designed molecules, the docking study along with MD simulation was also conducted to validate the calculated activities and the interaction mechanisms. The docking method was similar as the docking process above and the binding conformations with highest score were subjected to MD simulations in AMBER molecular dynamics package 39 . Each MD time was 10 ns in explicit solvent. Antechamber tool was applied to generate the partial atomic charges of each molecule and the force field ff12SB was loaded to the receptor. Truncated octahedral water box by the TIP3P water model was chosen to encompass the complexes before minimization and MD. After the MD, analysis was performed with the ptraj analysis tool based on the last 2 ns trajectories of the simulation. The binding free energies of complexes were calculated using the MM-PBSA 40 .

Results and Discussion
CoMFA and CoMSIA results. Both models were obtained on the basis of a 44 molecules training set of whose pIC 50 values ranging from 4 to 7. Table 2 listed the statistical parameters of CoMFA and CoMSIA models.
And the results of the progressive scrambling demonstrated that the value for the slope in the five component model is admissible (Fig. 2), and with a minimum cSDEP and maximumQ 2 the optimum statistics are also seen for six components in Table 3    CoMFA contour map analysis. CoMFA contour maps were drawn vividly to explore the areas in three-dimensional space around the compounds where modifications would alter activity. The contour maps are revealed in Fig. 5, with compound 15 as the template molecule. In Fig. 5a, the green blocks mean a bulky group favored area, while the yellow blocks indicate that minor substituent are preferable to enhance the activity. In Fig. 5b, the electron-donating group and electron-withdrawing group favored region are represented by blue and red contours, respectively. As shown in Fig. 5a, a giant green block located near the hydroxyl group substituted to phenyl ring along with the large green contours around the non-aromatic hexatomic ring of compound 15 indicate that bulky groups here can increase the activity. These factors may explain why compounds 11, 12, 41, 42 and 47 with halogen substituents in this area are more potent than molecules without any substituents at this particular position. The yellow contour on the top or bottom of the non-aromatic hexatomic ring of compound 15 suggested the unfavorable influence of bulky groups. This might be the reason why compound 40 whose methoxyl substituent is not on the plane of the molecule showed significantly decreased activities.
As shown in Fig. 5b, red contours around substituent group R 1 or R 4 , and para-and meta-positions of hexatomic ring of group R 3 or R 5 revealed that electron-withdrawing substituents were considered to be beneficial to the activity. The compounds 29, 41, 42 and 47, with electron-withdrawing groups at R 1 and the meta-position of hexatomic ring of group R 3 or R 5 , displayed the good bioactivity. On the contrary, the blue contour around meta-position of hexatomic ring of group R 3 or R 5 told us that electron-donating groups here would be expected. Therefore, compound 28 with the methyl group at the meta-position of benzene ring displayed good IC 50 values.
CoMSIA contour map analysis. CoMSIA contour maps of steric, electrostatic, H-bond donor and acceptor field are revealed in Fig. 6. Generally the CoMSIA contour maps of electrostatic and steric field are similar to CoMFA models. However, the electrostatic field in CoMSIA model appended a blue contour around group R 2 , which means that electron-donating substituents were also preferred in this region.
The hydrogen-bond donor field was presented in Fig. 6c. Cyan contour around hydroxyl substituent of compound 15 revealed that hydrogen bond donor was preferred in this region. In hydrogen bond acceptor field (Fig. 6d), the purple contour was around group R 3 , which revealed that hydrogen acceptor was preferred in this region. Therefore, compounds 6 and 11 with hydrogen acceptor at group R 3 showed good IC 50 values. Furthermore, the hydrogen acceptor was also preferred in the side chain of furan ring. Molecular Docking Analysis. The Figs 7 and 8 demostrated the binding modes of tyrosyl-tRNA synthase with compound 15. The hydroxyl groups at R 1 position formed H-bonds with Tyr46 and Asp170 residue by acting as H-bond acceptor and donor simultaneously. The morpholine nitrogen in R 3 , which served as a H-bond acceptor, also formed a hydrogen bond with Glu40. These interactions summarized from Fig. 7 coincided satisfactorily with the contour maps derived from CoMSIA method above. Moreover, as can be seen in Fig. 8, a cavity located beneath the morpholine ring where steric contour map hints us a bulky group favored verified our results.
Summary of the structure-activity relationships. The SAR derived from the present work were illustrated in Fig. 9. To be specific, the bulky, electron-withdrawing and h-bond donor groups such as hydroxyl group at R 1 or R 4 position are essential to improve the antibacterial activity. The electron-withdrawing, H-bond acceptor or bulky groups like ether bond at para-position or meta-position of substituent group R 3 or R 5 may increase their activities. Introducing electron-donating groups at meta-position of substituent group R 3 or R 5 may also enhance the activity. Furthermore, the hydrogen acceptor was also preferred in the side chain at furan ring.

PLS Statistics
CoMFA CoMSIA Design of Novel Derivatives. After gaining the SAR revealed by the foregoing study, ten new furan-2(5H)-ones (D1-D10) were designed and predicted. Molecular alignment was conduct on these new molecules and their activities were calculated to be better than compound 15 (pIC 50 = 6.941). Such results indicated that these 3D-QSAR models with considerable predictive ability could be prospectively used in structure modification and optimization. The structures and calculated pIC 50 of these newly designed molecules were listed in Table 4.    Validation of Newly-designed Derivatives. These ten new furan-2(5H)-ones (D1-D10) were re-docked to the receptor 3P0H and based on which MD simulations were also performed to check the potencies by assess the binding affinity and free energy. The docking results was listed in the following Table 5. The ΔG calculated from the MD simulation results were also shown below. Compound D7 showinged significant affinity and the interactions between ligand and receptor were depicted in Fig. 10, where compound 15 was also exhibited as the reference. We can find that more hydrogen bonds with stronger bond energy formed based on the modified group, for instance the acylamino in R 1 connecting Asp170 with a new hydrogen bond and interacting Tyr36 with   a hydrogen bond of −3.3 K cal/mol is more favored than the original hydroxyl in this position whose hydrogen bond was −2.9 K cal/mol. However, comparing with the predicted pIC 50 values, compounds with trichloromethyl or bromine demonstrated less potency, especially the compound D6 are of poor affinity. By analyzing the conformations optimized from the MD simulation, we found that although the substituent group of R 3 are bulky favored, the group size is not unlimited, substituent group larger than furan ring with trifluoromethyl may cause crash of protein and ligand thus changed the equalized conformation after the MD simulation. Compound D1 suffered the same situation that showed dissatisfactory inhibitory activity in the docking and optimization. In Fig. 11, the conformation of compound D6 altered hugely to extend the trifluoromethyl out of the surface to reduce the tension. Additionally the amide group in R 1 served as a H-bond donor by forming a stable hydrogen bond with His210 and verified the SAR obtained above. According to the docking and MD simulation study the compound D7 was also proven to be most promising in tyrosyl-tRNA synthase inhibition.

Conclusion
The tyrosyl-tRNA synthase inhibitors provide a promising approach in fighting against drug-resistant bacteria. In our present studies, we have established CoMFA model (q 2 = 0.611, r 2 = 0.940) and CoMSIA model (q 2 = 0.546, r 2 = 0.905) with satisfactory correlation and predictive abilities from fifty-two tyrosyl-tRNA synthase inhibitors. The CoMFA and CoMSIA contour maps provided information to summarize the SAR and characterized key features affecting the antibacterial activity. Moreover, the prediction ability of the model validated by the test set turned out to be satisfactory which means these models could be applied in predicting the activities of new compounds. Thus 10 new molecules were designed and predicted to be more potent than compound 15, these results were further validated by MD simulation study.All the work above confirmed that our models can provide a ponderable clue in designing novel antimicrobial agents.