Identification of potential inhibitors for LLM of Staphylococcus aureus: structure-based pharmacophore modeling, molecular dynamics, and binding free energy studies

Abstract Staphylococcus aureus causes various life-threatening diseases in humans and developed resistance to several antibiotics. Lipophilic membrane (LLM) protein regulates bacterial lysis rate and methicillin resistance level in S. aureus. To identify potential lead molecules, we performed a structure-based pharmacophore modeling by consideration of pharmacophore properties from LLM-tunicamycin complex. Further, virtual screening of ZINC database against the LLM was conducted and compounds were assessed for Lipinski and ADMET properties. Based on pharmacokinetic, and molecular docking, five potential inhibitors (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) were identified. Molecular dynamics simulation (MDS) of these five molecules was performed to evaluate the dynamics and stability of protein after binding of the ligands. Several MDS analysis like RMSD, RMSF, Rg, SASA, and PCA confirm that identified compounds exhibit higher binding affinity as compared to tunicamycin for LLM. The binding free energy analysis reveals that five compounds exhibit higher binding energy in the range of −218.76 to −159.52 kJ/mol, which is higher as compared to tunicamycin (–116.13 kJ/mol). Individual residue decomposition analysis concludes that Asn148, Asp151, Asp208, His271, and His272 of LLM play a significant role in the formation of lower energy LLM-inhibitor(s) complexes. These predicted molecules displayed pharmacological and structural properties and may be further used to develop novel antimicrobial compounds against S. aureus. Communicated by Ramaswamy H. Sarma


Introduction
Staphylococcus aureus is a Gram-positive coccal bacterium, which is a major cause of diseases in the hospital and community (Lowy, 2003). It is frequently discovered on human skin and in the respiratory tract. It causes several skin infections to life-threatening diseases such as endocarditis, meningitis, pneumonia, and toxic shock syndrome, etc. (Lowy, 2003). S. aureus resistant to multiple-drug resistant such as methicillin-resistant S. aureus (MRSA) causes it difficult to treat and leads to a mortality of 20 to 40% (Lowy, 2003). The need for new antibiotics to treat S. aureus infection is urgent as it was in 1940 when it was treated with penicillin (Bassetti et al., 2013;Walsh, 2003).
S. aureus acquired methicillin resistance to b-lactam antibiotics due to the presence of penicillin-binding protein 2 0 (PBP2A), a penicillin-binding protein (PBP) encoded by mecA gene present within the mec locus (Fontana, 1985;Pinho et al., 2001;Reynolds & Fuller, 1986;Ubukata et al., 1989). PBP2A exhibits a very low affinity for b-lactam antibiotics and mainly responsible for synthesis of peptidoglycan in the presence of antibiotics. The expression of methicillin resistance can be homogenous and heterogeneous (Tomasz et al., 1991). The mecI-mecR1 gene, present upstream of mecA regulates the expression of PBP2'A (Suzuki et al., 1993;Tesch et al., 1990). Moreover, blaI gene also regulates the transcription of mecA gene in S. aureus (Hiramatsu et al., 1990). The llm and fem genes, autonomous of mec locus, are involved in resistance to methicillin for S. aureus (Maki et al., 1994;Matthews & Tomasz, 1994). It has been found that inactivation of llm results in a lower level of methicillin resistance in S. aureus (Maki et al., 1994). The llm gene encodes a lipophilic membrane (LLM) protein, which affects bacterial lysis rate and methicillin resistance level while its function is not determined yet. Maki et al. observations show that this protein is involved in the metabolism of peptidoglycan (Maki et al., 1994).
In the current study, the 3-dimensional structure of LLM from S. aureus was predicted using homology modeling. The protein sequence of LLM was taken from NCBI database and BLAST tool was used to find the homologous structures. Molecular docking of Tunicamycin, drug blocks N-linked glycosylation, was done along with LLM. Structure-based pharmacophore modeling was done by considering the LLMtunicamycin complex. Drugs like molecules from ZINC database were retrieved by using pharmacophore features of LLM-tunicamycin complex. Virtual screening of downloaded molecules was performed to check their binding affinities with LLM. The molecules showed higher binding affinity than tunicamycin were considered and checked on the basis of pharmacokinetic properties. Further, molecular docking of compounds with LLM was done to cross-check the interactions of ligands with protein. Molecular dynamics was employed to understand the stable behavior of the proteinligand complex in a dynamics environment. Binding free energy calculations and per residue decomposition method was used to confirm the binding stability of ligand in the protein-ligand complex. The summary of the workflow done for the identification of novel potent molecules for LLM is shown in Figure 1. Our study suggests that five identified molecules have the potential to interact with LLM.

Molecular modeling
LLM protein sequence from Staphylococcus aureus was retrieved from NCBI database (BBA23401.1) and BLAST tool was utilized to determine the homologous structures (Altschul et al., 1990). The multiple sequence alignment (MSA) of LLM from S. aureus with other homologous structures was done by using Clustal Omega webserver (http://www.ebi.ac.uk/Tools/ msa/clustalo/) and visualized using Espript 3.0, as done by Dalal et al. (http://espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi) Gouet et al., 1999;Sievers et al., 2011). The three-dimensional structure of LLM was generated using RaptorX on windows based system (Wang et al., 2016). RaptorX is a webserver used a powerful in-house deep learning model DeepCNF (Deep Convolutional Neural Fields) to predict secondary structure content, solvent accessibility (ACC) and disordered regions (DISO) of a protein. ModLoop was used for the refinement of the disordered loops of the predicted model, as done by Singh et al. (Fiser & Sali, 2003;Singh et al., 2017). The structure was energy minimized using Swiss PDB viewer, as done by Saini et al. (Guex & Peitsch, 1997;Saini et al., 2019). The structure was finally optimized by 100 steps of steepest descent and conjugate gradient minimization in chimera (Pettersen et al., 2004). The model was validated by PROCHECK and ProSA integrated webserver (https://prosa. services.came.sbg.ac.at/prosa.php) to estimate the deviation of Z-score from the high-resolution structures, as done by Pandit et al. (Laskowski et al., 1993;Pandit et al., 2018;Wiederstein & Sippl, 2007). Additionally, molecular dynamics using GROMACS was run to check the stability of the model, as done by Singh et al. (Gunsteren et al., 1996;Singh et al., 2019;Van Der Spoel et al., 2005). Furthermore, root means square deviation (RMSD) between predicted structure and template structure was calculated using PyMOL (DeLano, 2002).

Molecular docking of tunicamycin
Tunicamycin was retrieved from MraY tunicamycin complex (PDB ID: 5JNQ) (Hakulinen et al., 2017). Molecular docking of tunicamycin (TUM) with LLM was done using AutoDock Vina and AutoDock Tools, as done by Singh et al. (Morris et al., 2009;Singh et al., 2018;Trott & Olson, 2010). AutoDock Tools was used for the addition of Kollman charges (12.055) and the receptor was saved in .pdbqt format. Polar hydrogen atoms and gasteiger charges (-0.0002) were applied to TUM. Potential grid around active site residues (Lys97, Asn148, Asp151, Gly152, Leu153, Asp154, Asn198, Phe205, Asp208, His271, and His272) of LLM was generated using AutoGrid4. The grid dimension and centre point coordinates are as 66 Å, 70 Å, 60 Å, and X ¼ 154.87, Y ¼ 144.16, and Z ¼ 264.74. The lamarckian genetic algorithm was utilized to generate the 50 poses of ligand. In AutoDock Vina, configuration file consists of grid box properties along with receptor and ligand files were used to generate the 20 conformations of ligand. Additionally, HADDOCK webserver was also used for molecular docking of TUM with LLM and active site residues of receptor were mentioned as active residues (De Vries et al., 2010;Kurkcuoglu et al., 2018). In HADDOCK, 198 structures in 2 clusters were generated with 99.0% of water refined models. All the generated docking confirmations were visualized and analyzed in PyMOL2.0.5 and Maestro11.2, as done by Saini et al. (DeLano, 2002;Saini et al., 2021).

Pharmacophore modeling
Pblast results showed that several homologous structures of LLM are in complex with TUM: MraY tunicamycin complex (PDB ID: 5JNQ), Human GPT (DPAGT1) in complex with tunicamycin (PDB ID: 6BW5), and Crystal structure of human UDP-N-acetylglucosamine-dolichyl-phosphate N-acetylglucosaminephosphotransferase (DPAGT1) (V264G mutant) in complex with tunicamycin (PDB ID: 5O5E) (Dong et al., 2018;Hakulinen et al., 2017;Yoo et al., 2018). The structure of tunicamycin is similar to the substrate of the protein. Therefore, structure-based pharmacophore modeling of TUM with LLM was performed using the Pharmit webserver, as done by Sunseri & Koes, 2016). The pharmacophore model shows the common chemical environmental features that are required for optimal interactions of a ligand with a receptor. In pharmit, LLM protein and TUM were considered as receptor and ligand features, respectively. Pharmit generally displays the aromatic, hydrogen donor, hydrogen acceptor, hydrophobic, negative, or positive ion of a ligand. 123399574 conformers of 13,190317 molecules from ZINC database were considered as ligands for pharmacophore modeling (Sterling & Irwin, 2015). Pharmacophore mapped features and receptor was set as inclusive and exclusive shapes in pharmit. Molecules were filtered by specifying molecular weight (200 to 500 Da), hydrogen bond acceptor (<10), and hydrogen bond donors (<5). Further, ligands were screened to 687 molecules having RMSD less than 1.5 Å with queried features.

Virtual screening
687 molecules downloaded from ZINC database were considered as ligands for virtual screening. All the molecules were energy minimized using Universal Force Field (UFF) in PyRx, as done by Malik et al. (Dallakyan & Olson, 2015;Malik et al., 2019). Open Bable was utilized to convert the minimized ligands into pdbqt format (O'Boyle et al., 2011). Virtual screening was done using AutoDock Vina in PyRx0.8 and nine distinct poses were generated for each molecule, as done by Trott & Olson, 2010). The molecular grid for virtual screening was considered around Lys97, Asn148, Asp151, Gly152, Leu153, Asp154, Asn198, Phe205, Asp208, His271, and His272 of receptor. The dimension and grid points for screening are as: 62 Å, 65 Å, 56 Å, and X ¼ 154.8, Y ¼ 144.2, and Z ¼ 264.7. Further, conformations showing the highest binding affinity and corresponding interactions with protein were saved and visualized in PyMOL (DeLano, 2002). Five molecules (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) showing higher binding affinity as compared to TUM were considered for further study.

In-silico pharmacokinetic and ADMET studies
The pkCSM webserver was utilized to calculate the physiochemical and pharmacokinetic parameters like molecules, polar surface area, number of hydrogen bond donors, and acceptors, as done by Pires et al., 2015). The in-silico absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties have been determined by using pkCSM server. The ADMET profile of molecules was predicted by various factors like water solubility, cytochrome P450, human intestinal absorption, skin permeability, volume of distribution (VDss), blood-brain barrier (BBB) permeability, AMES toxicity, total clearance, and skin sensitization.

Molecular dynamics
Molecular dynamics of LLM-Native, LLM-TUM, and LLM-inhibitor(s) complexes were done to check the stability at an atomic level. Molecular dynamics simulation was run using GROMOS96 43a1 force field in GROningen Machine for Chemical Simulations (GROMACS), hosted on NMRbox computing platform (Gunsteren et al., 1996;Maciejewski et al., 2017;Van Der Spoel et al., 2005). Topology of TUM and inhibitors (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) were generated using PRODRG webserver, as done by Kesari et al. (Kesari et al., 2020;Sch€ uttelkopf & Van Aalten, 2004). PRODRG produces the parameters for bond, angles, and charges on GROMOS96 43a1 force field. The protein complexes were solvated with simple point charges (SPC) in a triclinic box of volume 472.67 nm 3 with 1.5 nm distance from the nearest atom of protein. Systems were neutralized by the addition of counterions (Cl-ions) by using 'genion' tool and minimized by using the steepest descent algorithm to remove the steric clashes. Minimization was performed for 50,000 iteration steps up to an energy cut of 10.0 kJmol À1 . Particle Mesh Ewald (PME) method was used for the calculations of long-range electrostatic interactions, as done by Singh et al. (Abraham & Gready, 2011;Singh et al., 2020). Systems were equilibrated for two steps equilibration: a constant number of particles, volume, and temperature (NVT) and a constant number of particles, pressure and temperature (NPT) at 300 K for 1 ns with Linear constraints Solver (LINCS) algorithm to constrain the covalent bonds, as done by Dalal et al. (Dalal et al., 2021;Hess et al., 1997). NVT was performed using Parrinello-Rahman barostat pressure coupling method while NPT was done using Berendsen thermostat, as done Kumari et al. (Berendsen et al., 1984;Parrinello & Rahman, 1981;R. Kumari et al., 2021). Final, molecular dynamics production run for 100 ns and the coordinates were updated at every 10 ps, as done by Gupta et al. (2021). The trajectory retrieved from the molecular dynamics was utilized for the analysis of different results like Root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), hydrogen bond numbers, and principal component analysis (PCA).

mmpbsa binding free energy calculation
The g_mmpbsa tool was used to calculate the binding free energy of LLM-TUM and LLM-inhibitor(s) complexes (R. Kumari et al., 2014). The binding free energy of the proteinligand complex in a solvent was calculated as:- Where DG complex is the total free energy of the protein-ligand complex, DG protein and DG ligand are the total energy of protein and ligand in a solvent, respectively. Molecular Mechanics/Position-Boltzmann Surface Area (MMPBSA) method to determine the binding free energy of the proteinligand complex, as done by . MMPBSA performs in three steps: a) calculation of potential energy in a vacuum, b) polar solvation energy, and c) non-polar solvation energy. MMPBSA do the molecular mechanics (vdw and electrostatic) vacuum energy with energy decomposition and further compute the polar and non-polar solvation energy. In the current work, the trajectory was retrieved from the last 20 ns (80 to 100 ns) of molecular dynamics. Total 2000 snapshots produced at every 10 ps during an interval of 80 to 100 ns were used to do MMPBSA calculations of protein-ligand complexes. Per-residue decomposition was done to determine the quantitative analysis of the energetic contribution of each amino acid. The binding energy decomposition study is crucial for elucidating the contribution of each residue to the ligand.

Molecular modeling
An amino acid sequence of LLM from S. aureus was retrieved from NCBI database. Multiple sequence alignment of LLM with its homologous structures: MraY tunicamycin complex (PDB ID:5JNQ), human GPT (DPAGT1) in complex with tunicamycin and crystal structure of polyprenyl-phosphate N-acetyl hexosamine 1-phosphate transferase (PDB ID:4J72), was generated using Clustal Omega and ESPript3 as shown in Figure  S1. RaptorX server was utilized for the prediction of the three-dimensional model of LLM. The model was energy minimized by using Swiss PDB viewer and disordered loops were refined by using ModLoop. Ramachandran plot generated by PROCHECK shows that 88.0% residues are in the most favored region, 10.5% and 1.1% in additionally allowed and generously allowed and 0.4% residues in disallowed regions (Table S1 and Figure S2A). The generated graphical plot having scored for all experimentally determine NMR and X-ray structures were used to compare the quality of the model. ProSA results depict that the refined model and its template Z score are À3.28 and À2.58, respectively, as shown in Figure S2B and 2 C. The Z score shows that the predicted model is in good agreement with the homologous structure. Additionally, the refined model was subjected to molecular dynamics using GROMACS. The RMSD plot shows that model got convergence at 0.45 nm at 15 ns and is stable throughout the simulation Figure S2D. Moreover, the structural superposition of the predicted model with homologous structure showed an RMSD of 0.25 Å for 221 Ca atoms. All these parameters affirm that the predicted model is good and can be considered for further study. The validated model consists of 304 residues along with 11 a helices and 2 b sheets, as shown in Figure S3.

Molecular docking of tunicamycin
Molecular docking was done to analyze the binding affinity of the tunicamycin with LLM. Docking of tunicamycin with LLM was done using AutoDock Tools, AutoDock Vina, and HADDOCK. Tunicamycin shows hydrogen bonding with Asp32, Lys33, Asn35, Lys97, Asn148, Asp151, Asp208, and His272, as shown in Figure 2. Molecular docking performed using AutoDock Vina and AutoDock Tools shows the binding affinity of À8.1 and À7.49 kcal/mol, respectively, as shown in Table 1. HADDOCK docking results show that tunicamycin show a haddock score of À67.3 ± 3.5 with LLM as shown in Table S2. Further, the LLM-TUM complex was considered for further study for pharmacophore modeling.

Pharmacophore modeling
LLM-TUM complex predicted by molecular docking was fed to the pharmit for pharmacophore modeling, as shown in Figure 3. The pharmacophore features are two hydrogen bond donors (white, radius: 1.0 Å), two hydrogen bond acceptors (orange, radius: 1.0 Å), hydrophobic feature (green, radius: 1.0 Å), and aromatic feature (purple, radius: 1.0 Å) as shown in Figure 3. Shape filters for the protein and ligand were applied after performing the pharmacophore search. Shape of ligand and receptor is kept as inclusive and exclusive features for the pharmacophore modeling. The constructed pharmacophore models were utilized to screen the 123399574 conformers of 13,190317 molecules from ZINC database. The results were further filtered and refined with molecular properties, RMSD, and energy minimization. All total of 687 filtered hits were considered for virtual screening against LLM.

Virtual screening
687 filtered molecules from pharmacophore modeling were retrieved from the ZINC database in the sdf format. Virtual screening was done using AutoDock Vina in Pyrx 0.8 using Universal Force Field (UFF) to minimize the ligands. Openbable was utilized to convert the sdf into pdbqt and nine generated poses were analyzed in PyMOL. 5 molecules (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) exhibit higher binding affinity in range of À8.1 to À8.8 kcal/mol than TUM (-7.6 kcal/mol) for LLM as shown in Table 1. All these selected molecules were bound at the active site of LLM.

In-silico pharmacokinetic and ADMET studies
Pharmacokinetic properties are considered to be the most significant challenging topics in structure-based drug designing.
These properties depend on Lipinski rule, absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties for consideration of a molecule as a drug candidate. Lipinski rule states that a drug candidate exhibits molecular weight 500 Da, logP 5, hydrogen bond donor 5, and hydrogen bond acceptors 10 (Lipinski, 2004). LogP is the ratio of a molecule in octanol to water, which plays a vital role in the distribution and excretion of compound. The pharmacokinetic properties of the selected molecules (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) were predicted by pkCSM webserver. All of the identified compounds fulfill the Lipinski rule of five, as shown in Table S3. pkCSM was used to predict the ADMET properties of the molecules, as shown in Table S4. All of the compounds shown good water solubility in the range of À3.18 to À4.59 log mol/L. All the molecules exhibit good human intestinal absorption of more than 65%. Likewise, skin permeability less than À2.7 log Kp indicates that these molecules are skin permeable. The steady-state volume of distribution (VDss) and unbound plasma fraction of all the molecules are more than 0.55 and 0.1, respectively. Blood-brain barrier (BBB) permeability and central nervous system (CNS) permeability lies in the range of À1.01 to À1.53 and À2.24 to À3.39, respectively. The total clearance of all the compounds is below 0.9. Negative AMES toxicity indicates that the molecules are non-mutagenic and non-carcinogenic.

Molecular docking of compounds
AutoDock Tools and HADDOCK were used for the molecular docking of potent molecules fulfilling the pharmacokinetic properties. Molecular docking results shows that 5 molecules (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) shows higher binding affinity in range of À7.49 to À8.83 kcal/mol as compared to TUM (-7.2 kcal/mol) as shown in Table 1. ZINC000072380005 forms a low energy complex with LLM as compared to other predicted molecules, as indicated by AutoDock Tools, AutoDock Vina, and HADDOCK. It forms seven hydrogen bonds with Asn35, Asn148, Gly207, Asp208, Ser268, His269, and His272 with an energy complex of À8.83 kcal/mol as shown in Figure 2 and Table 1. ZINC000072380005-LLM complex is also stabilized by four hydrophobic interactions with Leu149, Phe205, Ala211, and Ala251, as shown in Figure S4. HADDOCK docking results show score (-102.0 ± 4.6) with electrostatic energy of À236.4 ± 18.1 kcal/mol. Molecular docking performed with AutoDock Tools and HADDOCK shows that binding of ZINC000257219974 to LLM form a low energy complex. ZINC000257219974 stabilized by seven hydrogen bonds via Asn32, Asn148, Asp151, Asp208, Arg254, Ser268, and His269, as shown in Figure 2. It forms hydrophobic interactions with residues Val46 and Phe205 of LLM, as shown in Figure S4. Additionally, polar and charged interactions are contributed by Asn35, Arg255, Lys267, and His272. AutoDock shows a binding affinity of À7.82 kcal/mol, as shown in Table 1. Molecular docking was done using HADDOCK show score À94.5 ± 0.9 along with electrostatic energy of À206.1 ± 20.4 kcal/mol, as shown in Table S2. ZINC000176045471 also formed a low energy complex as compared to TUM with LLM. Molecular docking performed using AutoDock Tools exhibit affinity of À8.3 kcal/ mol, as shown in Table 1. It forms six hydrogen bonds with amino acids Asn35, Asn148, Asp151, Gly207, Asp208, and His272 of LLM, as shown in Figure 2. Further, the molecule was stabilized at the active site of LLM by hydrophobic interactions contributed by Phe36, Val46, Leu149, and Phe 205, as shown in Figure S4. Moreover, docking done using HADDOCK shows a score of À76.5 ± 2.3, as shown in Table S2. ZINC000035296288 exhibits a binding affinity of À8.12 kcal/mol as per AutoDock Tools, as shown in Table 1. The ZINC000035296288-LLM complex was stabilized by three hydrogen bonds (Asn35, Asp208, and His272) and two hydrophobic interactions by Val46 and Phe205, as shown in Figures 2 and S4. Other polar and charged interactions to ZINC000035296288 from LLM are contributed by Asp32, Asp88, Asn148, Asp151, Arg254, Lys267, Ser268, and His269.

Molecular dynamics
Molecular dynamics (MD) simulation was done to study the effective conformational stability of protein and ligand. In several studies, MD simulation has been carried out to determine the variations at the atomistic level in the protein-ligand complex in the dynamic environment Dhankhar, Dalal, Kotra, et al., 2020;Dhankhar, Dalal, Mahto, et al., 2020;Dhankhar, Dalal, Singh, et al., 2020;N. Kumari et al., 2020;R. Kumari et al., 2021). In the present study, we investigated the different parameters such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), hydrogen bond formation, and principal component analysis (PCA) for LLM-TUM and LLM-inhibitor(s) complexes during the molecular simulation of 100 ns.

RMSD
The dynamics variation in Ca backbone of the LLM-TUM and LLM-inhibitor(s) complexes were determined by monitoring the RMSD during the molecular simulation. The RMSD plot shows an initial sharp increase in RMSD during the first 5 ns. RMSD plot signifies that most of the protein-ligand complexes achieved equilibration at 22 ns and stable up to 100 ns as, shown in Figure 4(A). The comparison of average RMSD of LLM-TUM and LLM-inhibitor(s) complexes are shown in Table S5. The RMSD fluctuations of LLM-ligand(s) complexes are in the range of 0.3 to 0.41 nm. RMSD plot indicates that there is no large conformational changes in Ca backbone of LLM after binding of inhibitor(s). Ligand RMSD graph shows that identified compounds exhibit lesser RMSD as compared to TUM during the molecular simulation, as shown in Figure 4(B). All identified compounds showed  average RMSD in a range of 0.08 to 0.11 nm, which is lesser as compared to TUM (0.12 nm). Overall, RMSD results suggest that the binding of inhibitor(s) at the active site of LLM is stable and does not affect the stability of the protein.

RMSF
RMSF was determined to calculate the flexibility of the backbone of protein after fitting to the frame of the reference. The higher RMSF values suggest the loosely organized loops and turns, while lower flexibility indicates the secondary structures like helix and sheets. In this study, we determined the Ca atoms flexibility of each residue of LLM-TUM and LLM-inhibitor(s) complexes. The RMSF plot of LLMinhibitor(s) complexes is almost comparable to LLM-TUM, as shown in Figure S5. Results show a higher RMSF in residues range of 240 to 280, which contains a partial helical region between loops. The average RMSF of LLM-inhibitor(s) complexes are in the range of 0.13 to 0.2 nm, which is lesser than LLM-TUM (0.22 nm), as shown in Table S5. So, overall, RMSF results imply that the inhibitors are well fitted at the active site of LLM and form the stable LLM-inhibitor(s) complexes.

Radius of gyration (Rg)
The radius of gyration (Rg) depicts the overall compactness of protein during the molecular simulation. The radius of gyration demonstrates the moment of inertia of atoms of protein from their centre of mass during a particular time interval. The radius of gyration plot of LLM-TUM and LLMinhibitor(s) complexes are shown in Figure 5. Rg results show that the compactness of LLM-inhibitor(s) complexes are higher as compared to the LLM-TUM complex. Overall, the findings of Rg indicate that the binding of inhibitors to LLM results in the creations of stable complexes of LLM-inhibitor(s) compared to LLM-TUM complex.

Solvent accessible surface area (SASA)
Solvent Accessible Surface Area (SASA) of a protein is the surface area determined by its polar and non-polar interactions. SASA of a protein decrease with an increase in compactness of a protein so change in SASA value can be used to determine the stability of a protein. SASA values of LLM-inhibitor(s) complexes were lesser as compared to LLM-TUM complex, as shown in Figure S6. Average SASA values of LLM-TUM, LLM-ZINC000072380005, LLM-ZINC000257219974, LLM-ZINC000176045471, LLM-ZINC000035296288, and LLM-ZINC000008789934 were 141. 3, 131.1, 133.3, 137.9, 135.4, and 135.1 nm 2 , respectively. SASA results affirm that LLM-inhibitors complexes are more compacted and stable than LLM-TUM complex.

Hydrogen bond analysis
Hydrogen bond numbers and distribution were studied to determine the stability of LLM-TUM and LLM-inhibitor(s) complexes during the molecular simulation of 100 ns. The g_hbond utility of GROMACS was used to generate the number and distribution of hydrogen bonds in the LLM-TUM and LLM-inhibitor(s) complexes. Intra-protein hydrogen bond plot of LLM-TUM and LLM-inhibitor(s) complexes are shown in Figure 6(A,B). The average number of intra protein hydrogen bond numbers of LLM-inhibitor(s)  complexes are in the range of 209.7 to 218.5, which is more than LLM-TUM complex (208.9), as shown in Table S5. Inter-molecular hydrogen bonding results show that LLM-inhibitor(s) complexes maintained at least four numbers of hydrogen bonds during the simulation, as shown in Figure 6(C). The hydrogen bonding length suggested that the affinity of hydrogen bonds of LLM-inhibitor(s) complexes are similar to LLM-TUM complex as shown in Figure  6(B,D). Overall, hydrogen bonding results suggested that binding of inhibitors to LLM results in stable LLM-inhibitor(s) complexes as compared to LLM-TUM complex.

Principal component analysis (PCA)
Essential dynamics was utilized to check the significant motion after ligand binding to the protein. Principal component analysis (PCA) was used to generate the dynamical motion along the 911 eigenvectors of all the proteinligand(s) complexes. The first eigenvector (PC1) and second eigenvector (PC2) contributes the major global motion to protein. So, the dynamical difference between LLM-TUM, LLM-ZINC000072380005, LLM-ZINC000257219974, LLM-ZINC000176045471, LLM-ZINC000035296288, and LLM-ZINC000008789934 can be easily understood from Ca motion along the first two eigenvectors. According to 2 D projection map, it can be easily observed that the binding of inhibitor(s) results in the formation of stable complexes, as shown in Figure S7. From the 2 D projection, we can confirm that LLM-ZINC000072380005 complex is most stable as it exhibits a very stable cluster and occupied the least phase space as compared to other compounds. While LLM-ZINC000008789934 complex showed the least stability among identified compounds.

Secondary structure content analysis and
The active site of LLM is generally lined by loops regions. These loop regions were shown to be maintained during the molecular simulation, as shown in Figure 7. Dictionary of Secondary Structure of Protein (DSSP) was used to predict the secondary structure of LLM-TUM and LLM-inhibitor(s) complexes during the molecular simulation. The secondary structure content of the protein remains constant throughout the molecular simulation, as shown in Figure S8.  TUM showed interactions with Asn148, Leu149, Asp151, Phe205, Asp208, His271, and His272, as shown in Figure S9. The identified inhibitor(s) exhibited non-polar interactions with the active site residues of LLM (Figures S10-S14). Overall, snapshots generated at different time periods of MD simulation confirmed that binding of TUM and inhibitors were stable to LLM.

mmpbsa binding free energy calculation
The quantitative estimation of the binding free energy of ligand to a protein was determined using the Molecular Mechanics/Position-Boltzmann Surface Area (MMPBSA). The molecular dynamics trajectories of the last 20 ns were used to predict the binding free energy of LLM-TUM and LLM-  MMPSA results indicate that all the LLM-inhibitor(s) complexes are highly stable. It affirms that inhibitors bind efficiently at the active site of LLM and form the stable LLMinhibitor(s) complexes. The contribution of individual residues actively participated in the binding of ligands was quantified using MMPBSA. The per residue interaction profile of Asn148, Asp151, Asp208, His271, and His272 were calculated. Binding energy contributed by important residues to inhibitors is higher as compared to TUM, as shown in Table  3. Overall, MMPBSA results conclude that Asn148, Asp151, Asp208, His271, and His272 residues play an important role in the binding of inhibitor(s) to LLM and results in the formation of stable complexes.

Discussion
S. aureus is a Gram-positive pathogen which causes several skin infections and diseases like pneumonia, meningitis, and toxic shock syndrome, etc. penicillin-binding protein 2 0 (PBP2A) is one of the main factors of methicillin resistance in S. aureus and has a very low affinity for b-lactam antibiotics (Fontana, 1985;Pinho et al., 2001;Reynolds & Fuller, 1986;Ubukata et al., 1989). Several studies have reported that various factors such as llm, fmta, fema, femC, femD, femE, or sigma factor involved in expression of the high level of methicillin resistance in S. aureus (Fan et al., 2007;Gustafson et al., 1994;Komatsuzawa et al., 1997;Kornblum et al., 1986;Maki et al., 1994;Ornelas-Soares et al., 1994;Wu et al., 1996). Lipophilic membrane (LLM) protein encoded by llm gene affects the bacterial lysis rate and inactivation of LLM causes the reduction of methicillin resistance in S. aureus (Maki et al., 1994). It has been found that LLM is responsible for the biosynthesis of peptidoglycan (Maki et al., 1994).
In the present study, the LLM structure from S. aureus was predicted and energy minimized using RaptorX and chimera. The quality of the refined model was assessed by utilizing PROCHECK and ProSA. Molecular dynamics confirmed that the predicted model is reliable and can be considered for further study. Molecular docking results showed that tunicamycin bind with a good binding affinity at the active site of LLM. Structure-based pharmacophore modeling was done using pharmacophore properties from LLM-TUM complex. Total of 687 molecules were filtered from 13,190317 molecules from ZINC database by six pharmacophore features (two hydrogen bond donors, two hydrogen bond acceptors, hydrophobic feature, and aromatic feature). Filtered molecules were screened against a validated model of LLM using AutDock Vina in PyRx. 5 compounds (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) bind at the active site of LLM and showed a higher binding affinity than tunicamycin.
The pharmacokinetic properties of the identified molecules were determined to screen the drug-like molecules from the selected compounds. All identified compounds exhibit molecular weight less than 500 Da, H-bond donor less than 5, H-bond acceptor less than 10 and cLogP less than 5, fulfilled the Lipinski rule of five. ADMET profile of identified compounds predicted by several factors such as water solubility, cytochrome P450, Human intestinal absorption, skin permeability, volume of distribution (VDss), bloodbrain barrier (BBB) permeability, and AMES toxicity, etc. confirmed that these compounds satisfied all the parameters for a potent inhibitor.
Further, molecular docking was done using AutoDock 4.2.6 and HADDOCK to study the interactions of identified molecules (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) with LLM. All the compounds showed a higher binding affinity than tunicamycin for LLM. All the above-mentioned compounds possess interactions with Asn148, Asp151, Asp208, His271, and His272 of LLM, as shown in Figures 2 and S4. ZINC000072380005 exhibits the highest binding energy via the formation of seven hydrogen bonds with Asn35, Asn148, Gly207, Asp208, Ser268, His269, and His272 of LLM as shown in Figure 2 and Table 1. The molecular docking results inferred that five compounds are stabilized by hydrogen bonding and hydrophobic interactions in LLM-inhibitor(s) complexes. The best binding confirmations of compounds generated by AutoDock Vina, AutoDock Tools, and HADDOCK were almost identical.
Molecular dynamics simulation study of the LLM-TUM and LLM-inhibitor(s) complex(s) were done to explore the atomistic variations in the protein-ligand complex. Lesser RMSD values of LLM-inhibitor(s) complex(s) as compared to LLM-TUM indicate that binding of inhibitor(s) to LLM results in the formation of stable protein-ligand complexes. Identified compounds exhibit RMSD in a range of 0.05 to 0.11 nm, which is lesser than TUM (0.12 nm). RMSF results suggested that molecules are well fitted at the active site of LLM without causing structural and conformational changes in protein-ligand complex. The smaller Rg values of LLM-inhibitor(s) complex in comparison to LLM-TUM complex infer that protein-inhibitor(s) complexes are compactly packed and stable. SASA results revealed that the binding of predicted molecules to LLM tends to make a more stable complex as compared to LLM-TUM complex. The intermolecular hydrogen plot confirms that the minimum four hydrogen bonds stable proteinligand complex during the molecular simulation. Further, PCA results depict that LLM-inhibitor(s) complexes showed less flexibility as compared to LLM-TUM complex. Secondary structure analysis and snapshots generated at different time intervals of molecular simulation confirm that identified molecules remain bind at the active site of LLM and sustain throughout the molecular simulation.
The LLM-TUM and LLM-inhibitor(s) complexes were utilized for calculations of binding free energy by using MMPBSA. All the identified compounds show higher binding affinity (-159.52 to À218.76 kJ/mol) compared to TUM (-116.13 kJ/mol). Further, the per residue decomposition method was used to determine the contribution of individual residues participated in the binding of the ligand in proteinligand complex. Overall MMPBSA results conclude identified compounds interact with Asn148, Asp151, Asp208, His271, and His272 of LLM and results in stable LLM-inhibitor(s) complex as compared to LLM-TUM complex.

Conclusion
LLM protein affects the methicillin resistance level and is involved in the metabolism of peptidoglycan in S. aureus. In the current study, structure-based pharmacophore modeling, molecular docking, molecular dynamics simulation, and MMPBSA were used to identify the potent molecules for inhibition of LLM. 687 molecules were filtered from 13,190317 molecules by structure-based pharmacophore modeling from the pharmacophore features of LLM-tunicamycin complex. Molecular docking results show that 5 molecules (ZINC000072380005, ZINC000257219974, ZINC000176045471, ZINC000035296288, and ZINC000008789934) bind at the active site of LLM with higher binding affinity in comparison to tunicamycin. Molecular dynamics simulations were done to explore the interactions of inhibitor(s) with LLM. Molecular dynamics results reveal that binding of selected molecules to LLM results in the formation of stable LLM-inhibitor(s) complexes. Moreover, MMPBSA results confirm that inhibitor(s) molecules interact with active site residues (Asn148, Asp151, Asp208, His271, and His272) of LLM and tend to form the lower energy complexes are compared to LLM-tunicamycin complex. Further, these identified molecules can be tested in vitro and provide a path for the development of antimicrobial compounds against S. aureus.