Identification of novel inhibitors of Neisseria gonorrhoeae MurI using homology modeling, structure-based pharmacophore, molecular docking, and molecular dynamics simulation-based approach

Abstract MurI is one of the most significant role players in the biosynthesis of the peptidoglycan layer in Neisseria gonorrhoeae (Ng). We attempted to highlight the structural and functional relationship between Ng-MurI and D-glutamate to design novel molecules targeting this interaction. The three-dimensional (3D) model of the protein was constructed by homology modeling and the quality and consistency of generated model were assessed. The binding site of the protein was identified by molecular docking studies and a pharmacophore was identified using the interactions of the control ligand. The structure-based pharmacophore model was validated and employed for high-throughput virtual screening and molecular docking to identify novel Ng-MurI inhibitors. Finally, the model was optimized by molecular dynamics (MD) simulations and the optimized model complex with the substrate glutamate and novel molecules facilitated us to confirm the stability of the protein-ligand docked complexes. The 100 ns MD simulations of the potential lead compounds with protein confirmed that the modeled complexes were stable. This study identifies novel potential compounds with good fitness and docking scores, which made the interactions of biological significance within the protein active site. Hence, the identified compounds may act as new leads to design and develop Ng-MurI inhibitors. Communicated by Ramaswamy H. Sarma


Introduction
There is a growing concern worldwide regarding the problem of infectious diseases, especially sexually transmitted diseases (STDs). STDs are a group of communicable diseases that are transmitted predominantly through sexual contact. Amongst STDs, gonorrhea continues to be a major health problem with serious health, social and economic consequences. Gonorrhea is a commonly and universally encountered sexually transmitted bacterial disease, caused by gram-negative diplococcus, and Neisseria gonorrhoeae. It involves the mucosae of the genital, anorectal, pharyngeal, and ocular regions. The clinical spectrum of gonococcal infection extends from minimally asymptomatic mucosal colonization to frank inflammatory mucosal disease (urethritis, cervicitis, proctitis, pharyngitis, conjunctivitis) that may lead to local complications such as epididymitis, salpingitis, and pelvic inflammatory disease as well as serious systemic complications such as arthritis, carditis, and dermatitis. Apart from causing disastrous sequelae such as infertility and ectopic pregnancy, it constitutes a major risk factor for HIV transmission through increased viral shedding from the inflamed mucosae (Donnell et al., 2010). In 2019, a total of 616,392 cases of gonorrhea were reported to the CDC, making it the second most common notable condition in the United States for that year. During 2018-2019, the overall rate of reported gonorrhea increased by 5.7%. Rates increased among both males and females, in all regions of the United States, and among all racial/Hispanic ethnicity groups (National Overview -Sexually Transmitted Disease Surveillance, 2019, 2021). Globally, approximately 40% of all reported cases were from the Western Pacific WHO reporting region. This region comprises 37 countries including developed countries such as Australia, Japan, New Zealand, the Republic of Korea, China, and Singapore, and developing nations such as Vietnam, Malaysia, and the Philippines. Therefore, there are cleaner and more cogent reasons than ever before to ensure that gonococcal disease is timely diagnosed and treated effectively.
Monitoring of antimicrobial susceptibilities of Neisseria gonorrhoeae is important to investigate treatment failure and to evaluate the efficacy of currently recommended therapies. The continuing emergence of antibiotic-resistant strains of Neisseria gonorrhoeae and the preponderance of disease in patients in resource-poor settings where many recommended antibiotics are either unavailable or expensive are the major limiting factors. Antibiotics are recommended at the initial clinical visit, before any knowledge of the organism's susceptibility pattern or local sensitivity pattern. This has led to an alarming increase in the number of isolates of Neisseria gonorrhoea resistant to commonly used drugs (Bala & Sood, 2010;Bhalla et al., 2002;Edwards & Apicella, 2004;National Overview -Sexually Transmitted Disease Surveillance, 2019, 2021Nizamuddin et al., 2011;Olivieri, 2002). Another major reason for increasing resistance across the world is due to injudicious use of antibiotics and increased global travel  which increases the selective pressure on bacterium and transfer of resistance genes via horizontal gene transfer (Hamilton & Dillard, 2006;Spratt et al. 1992). The greater concern is the recent emergence of resistant strains to extended-spectrum cephalosporins (ESCs) namely, cefixime and ceftriaxone in several countries (Alirol et al., 2017) Therefore, Neisseria gonorrhoeae has been classified by the CDC as a Superbug and the prospect of untreatable gonorrhea was voiced in 2012 by both the CDC and WHO. Currently, WHO has recommended the use of multiple drugs for the treatment of gonorrhea (Fenton et al., 2003). Consequently, it has become imperative to initiate sustained national and international efforts to reduce infection and the misuse of antibiotics to prevent the further emergence and spread of antimicrobial resistance. It is necessary not only to monitor drug resistance and optimize treatment regimens but also to find new drug targets in Neisseria gonorrhoeae and the possibility of identification of new drugs for treating gonorrhea. It is important as antibiotic therapy remains the primary line of treatment since no effective vaccine is currently available for Neisseria gonorrhoeae.
The above-mentioned global concern as well as the advancement in computational research prompted us to understand the biology and drug discovery potential of one such enzyme, glutamate racemase (MurI), which is involved in the early phases of peptidoglycan biosynthesis, an essential component of the cell wall. The bacterial cell wall, in general, is a highly cross-linked polymeric structure consisting of repeating peptidoglycan units of disaccharides N-acetyl glucosamine and N-acetyl muramic acid (joined in b-1-4 linkage) which contain a novel pentapeptide substitution (Heijenoort, 2001). Cross-linking occurs through the transpeptidation of the peptide linkages between adjacent glycan strands. Such a structural mesh imparts a cellular skeleton and protects the bacterial cell from osmotic lysis. Inhibitors that target proteins/enzymes of cell wall synthesis will therefore be deleterious to the bacteria and hence bactericidal.
Based on the location of the peptidoglycan biosynthesis machinery, three distinct phases have been categorized. Phase III involves the extracellular cross-linking and final maturation of the cellular envelope and is most frequently targeted. Unfortunately, resistance continues to rise and limits the therapeutic utility of both the existing and future compounds within these classes. The first committed step to peptidoglycan biosynthesis is the conversion of glucosamine-1-phosphate to UDP-N-acetylglucosamine (UDP-GlcNAc) and eventually production of UDP-N-acetylmuramic acid-pentapeptide (UDPMurNAc-pentapeptide during the Phase I of cell wall synthesis. An important feature of peptidoglycan is the incorporation of D-amino acids, D-glutamate by the ligase Mud, and D-alanine as a dipeptide formed by D-ala-Dala ligase, which is a substrate of the ligase MurF. It is believed that the presence of these unusual residues helps the bacteria to evade the proteolytic cleavage by host proteases. While some variation exists in the peptide composition of the pentapeptide intermediate, the incorporation of D-glutamate as the second amino acid is highly conserved (Heijenoort, 2001). Hence, it is important to maintain levels of D-glutamate in the cell and this is achieved by the racemization of L-glutamate to its D isomer by glutamate racemase (MurI).
There are various similarities in the primary amino acid sequence of several glutamate racemase enzymes from different organisms. More importantly, glutamate racemase enzymes share similar active site architecture and catalytic mechanisms. Recent high-resolution structural data from a wide range of species including B. subtilis (Ruzheinikov et al., 2005), B. anthracis (May et al., 2007), E. coli, S. aureus, H. pylori (Lundqvist et al., 2007), and S. pyogenes (Kim et al., 2007) further support the highly conserved structural topology of glutamate racemase. The enzyme is comprised of two symmetric domains that are joined by a two-stranded hinge region. Since a number of these structures contain substrate bound at the active site, it was predicted that the active site is formed at the interface of these two domains. The structure also assigned the participating amino acid residues contributing to the catalytic machinery required for enantiomeric proton abstraction/donation. It was suggested that the conserved residues of the N-terminal domain are involved in the deprotonation of D-glutamate while the C-terminal domain encodes the residues involved in L-glutamate deprotonation. This hypothesis is supported by sitedirected mutagenesis of the conserved amino acids (Cys70 and Cys181) in the active site (Glavas & Tanner, 2001). Since D-glutamic acid is an essential building block for peptidoglycan in bacterial cell walls, and D-glutamate is synthesized from L-glutamate by glutamate racemase, it is considered an attractive target for the novel drug discovery process of the potential antibacterial agents.
Despite the importance of this enzyme in the cell wall biosynthesis pathway, the 3D structure for glutamate racemase of Neisseria gonorrhoeae is still not available. Since there is considerable homology between Neisseria gonorrhoeae glutamate racemase with enzymes found in other bacteria, the knowledge of similarity between the template and the query sequences can be translated into a threedimensional model of the query protein based on the identity between their sequences. In this study, we constructed the homology-based model of Ng-MurI, validated it, and then used it as a starting point for the structure-based drug design to identify the lead compounds via high throughput molecular docking. The pharmacophore-based approach was utilized to perform the high throughput virtual library screening of the natural compounds to identify the most potential hits and to identify the lead compounds.

Protein homology model
For homology modeling, the protein sequence of Neisseria gonorrhoeae Glutamate racemase (Ng-MurI) with 270 amino acids was retrieved from the NCBI with GenBank ID: BBM75780. Using PDB structural blast, multiple sequence alignment, proteins with high degrees of query coverage, and homology, we selected 1ZUW (Crystal structure of B. subtilis glutamate racemase (RacE) with D-Glutamate) (Ruzheinikov et al., 2005) as the structural template. Homology models were generated using the 'Build homology model' module of Biovia DS 2020. The generated models were ranked based upon their discrete optimized protein energy (DOPE) scores (Shen & Sali, 2006) which is a statistical potential used to assess homology models in protein structure prediction. The target model having the least DOPE score was selected. Missing loops and side chains were filled in with Protein Preparation Wizard and missing hydrogens were also added (Meslamani et al., 2012). To assess model accuracy in the predicted experimentally derived binding pose, the co-crystallized ligand 'D-glutamate' was manually removed from the selected model and then re-docked into the binding site.

Model validation
The Ng-MurI model was superimposed with the template structure using the Align and Superimpose Protein protocol to calculate the root mean square deviation (RMSD) of coordinates between the homology model and the template structure. Different web servers like SWISS-MODEL, Phyre2, Bhageerath-H, and Modeller (Jayaram et al., 2014;Kelley et al., 2015) were also used to validate and compare our generated model. The model was then optimized using an energy minimization protocol of DS 2020 with a gradient of 0.01 using the CHARMM force field (Vanommeslaeghe et al., 2010) to remove any steric clashes within the amino acid side chains. Using the PROCHECK server, the accuracy of the predicted model and its stereochemical properties were assessed using the Ramachandran plot and the overall goodness factor (G-factor). In addition, the model was analyzed by SAVES and ERRAT online servers (Messaoudi et al., 2013).

Structure-based pharmacophore generation
Pharmacophore is an abstract description of molecular features that are necessary for the molecular recognition of a ligand by a biological macromolecule. The pharmacophore model explains how structurally diverse ligands can bind to a common receptor site (Qing et al., 2014). For the generation of structure-based pharmacophore, the protein of interest MurI-Ng was selected as the receptor and D-Glutamate as the ligand. The substrate 'D-glutamate' was docked into the conserved binding site of modeled Ng-MurI receptor. The best docking pose of D-Glutamate was selected based on highest CDOCKER energy as well as the lowest RMSD as compared to D-Glutamate in the crystal structure of B.subtilis glutamate racemase. Finally the structure-based pharmacophore model was generated using the 'receptor-ligand pharmacophore generation module' of Biovia DS 2020 (Meslamani et al., 2012). This software interprets ligandreceptor interactions based on pharmacophore features like hydrogen bond donor, hydrogen bond acceptor, hydrophobic, and hydrophilic regions. 'Excluded volume' was added to the active site which is necessary to maintain the steric circumference of the macromolecule to identify a better and optimum pharmacophore (Seidel et al., 2010).

Database preparation, high-throughput virtual screening, and molecular docking studies
A library of drug-like compounds containing approximately 1,95,399 compounds from the ZINC database was employed and mapped on the selected pharmacophore model using the screen library module of Biovia DS 2020. After the initial round of screening based on pharmacophore mapping, and fit value cutoff, compounds were further filtered using ADMET, TOPKAT, VEBER, and Lipinski's filter (Moroy et al., 2012).Virtual screening is of cardinal importance for in silico drug discovery process to speed up drug development (Ekins et al., 2007), as it helps in the selection of the best candidate drugs. In the process of virtual screening, a pharmacophorebased approach was required for the most energetically favorable site that finds matching pharmacophores in the ligands. For filtering the database molecules, a minimum of 4-5 sites are required to match hypotheses with 8-10 sites. This criterion was followed in the present work to screen the ZINC database. Based on their fitness score, database hits were first ranked to measure how well the aligned ligand conformer matched the hypothesis based on RMSD, and site matching (Jha, Saluja, et al., 2022). Molecules filtered from HTVS (High-throughput virtual screening) were subjected to molecular docking. Before molecular docking, all the hit compounds were energy minimized using the steepest descent (3000 steps) followed by the conjugate gradient (5000 steps) using the energy minimization tool of Biovia DS 2020 (Helal et al., 2011). The CDOCKER docking module of Biovia DS 2020 was used and compounds with the best docking scores in terms of CDOCKER energy and fit value were shortlisted (Jha, Singh, et al., 2022).
To further improve the accuracy of the molecular docking study, we performed flexible docking by keeping active site residues (Ser15, Thr79, Cys189, Thr190, and His191) flexible by using the Flexible docking module of Biovia DS 2020 (Koska et al., 2008).

Molecular dynamics simulations
MD simulation studies were carried out to gain a better understanding of the stability and dynamical properties of protein-bound to final 5 lead molecules. The Ng-MurI model complexed with the shortlisted five lead molecules from the ZINC library was subjected to MD simulation for 20 ns using GROMACS 2020.1 (Abraham et al., 2015;Rashidieh et al., 2015) package. Next, after trajectory analyses, we performed a 100 ns MD run for the control compound, and the best Lead 3 and Lead 5 complexed with Ng-MurI. Gromos96 43a1force field (Pol-Fachin et al., 2009) was used to generate the protein force fields, while Prodrg server (Sch€ uttelkopf & van Aalten, 2004) was used to generate the force fields for the reference and lead molecules.
Briefly, these complexes were immersed in a cubic box of extended simple point-charge (SPC) water molecules. The complexes were neutralized by adding Na þ ions and Cl À ions. To relieve the short-range bad contacts, energy minimization was performed using the steepest descent method for 50,000 steps. Initially, the position-restrained simulations were carried out at 298 K for 100 ps. The top five complexes were initially subjected to a 20 ns MD simulations production run followed by 100 ns for selected compounds (Lead 3 and Lead 5) and the control at 298 K temperature and 1 bar pressure, using a 0.002 ps time step. The Parrinello-Rahman method (Marton ak et al., 2003) was used to control the pressure and a V-rescale thermostat was used to maintain the temperature. The long-range electrostatics were handled using the particle mesh Ewald (PME) method (Essmann et al., 1995) with a real-space cut-off of 10 Å, PME order of six, and a relative tolerance between long and short-range energies of 10 À6 kcal/mol. Short-range interactions were evaluated using a neighbor list of 10 Å updated every 10 steps, while Lennard-Jones (LJ) interactions and the real-space electrostatic interactions were truncated at 9 Å. Hydrogen bonds were constrained using the LINCS algorithm (Hess et al., 1998). The final models in all five complexes were obtained by averaging the snapshots from the trajectory generated by MD simulations after the structure stabilization was achieved.
To study the effect of inhibitor binding on the stability of MurI, we also performed the MD simulation studies of protein-bound with the substrate (i.e. glutamate racemase) under similar conditions. To study the conformational variations in the structures of MurI-substrate and MurI-inhibitor complexes, the root-mean-square fluctuations (RMSF) and RMSD of the Ca atoms were calculated. The convergence of MD simulations was analyzed in terms of the potential energy, RMSD, RMSF, H-bond analysis, and Radius of gyration (RoG) (Jha, Singh, et al., 2022).

Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA)
To gain more insight into the structural dynamics of the protein-ligand complex and to further validate the MD simulation findings, free energies of docked complexes from the top two lead compounds (Lead 3 and Lead 5) along with the control were calculated using molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) using g_mmpbsa tool. From the last 20 ns of each MD run (i.e. from 80 ns to 100 ns) MM-PBSA calculations were done with all the default parameters. Also, the binding energies of these three protein-ligand complexes were enumerated using the following equation as reported earlier (Kollman et al., 2000).
The receptor's binding free energy with its ligand in a solvent can be expressed as where DG complex is receptor-ligand total free energy and DG Rec and DG lig are isolated total free energies of receptor and ligand, respectively in a solvent. Furthermore, each individual entity's free energy can be provided by where G is the free energy for ligand, receptor or receptorligand complex, E gas is molecular mechanics (MM) potential energy in gaseous state and G sol is the solvation free energy of the respective entity. TS represents entropic contributions, where T and S refer to temperature and entropy, respectively. Equation (3) given below calculates various interactions.
where the E gas includes both bonded interaction (E int ) as well as non-bonded interaction (E ele þ E vdw ) terms based on MM force field parameters. The solvation free energy was calculated by Eqs. (4) and (5) given below.
In Eq. (4), G sol is solvation free energy calculated in terms of electrostatic (G PB(GB )) and non-electrostatic (G sol-np ) contributions, whereas in Eq. (5), SAS represents solvent accessible surface area of the model.

Homology modeling and validation
In the absence of crystal structures, homology modeling has become the primary method for obtaining a 3D representation of the target in recent years. Since no crystal structure of Ng-MurI glutamate racemase (RacE) was available, it was determined using homology modeling in Biovia DS 2020. From the structural PDB blast results, it was observed that the Ng-MurI sequence showed the best alignment with the crystal structure of B.subtilis glutamate racemase (RacE) with D-Glu (PDB Id: 1ZUW) (Ruzheinikov et al., 2005) with a query coverage (99%), sequence similarity (40%), e-value (5e À 30) and bit score of 169.483. Gaps between the template and query proteins were found to be 1%. The 3D structure of Ng-MurI was modeled based on the sequence alignment with (RacE) using PDB Id: 1ZUW (Figure 1).
The modeled structure with the lowest DOPE score of À31.936 was superimposed on the template structure i.e. 1ZUW (Figure 2a), and the RMSD value of 0.2475 A was calculated, which indicated significant backbone similarity between the modeled protein and the template structure. The modeled structure was then validated, with the PROCHECK server producing an overall G-factor value of 0.4, and the Ramachandran plot (Oberholser, 2010) (Figure 2b) revealed that 93.1% of residues were in the allowed region, and 6.5% in the fairly favored region.; no residue was present in disallowed regions. These values strongly suggested good overall stereochemical quality and stability of the model. Furthermore, according to ERRAT analysis, the model's overall quality factor was 70.23 ( Figure S1) and the value of verify-3D was 96.24 percent ( Figure S2), indicating that the generated model is of extremely high-quality 3D structure and reliable quality. All the evidence suggested that the model's backbone conformation, non-bonded interaction, and energy score were well within the range of the highquality model.

Active site analysis and validation
The identification of the binding site is essential for structure-based virtual screening of compound libraries. The binding of the substrate D-glutamate in the template protein (1ZUW) requires a twisted conformation at the binding site pocket, to access the binding site. For the validation of the active site cavity, substrate D-glutamate bound to the protein was removed and prepared using the 'Ligand Preparation' module of Biovia DS 2020 software and redocked with the active site residues of the Ng-MurI. The re-docked ligand D-glutamate exhibited an acceptable docking score of À30.5649 and was found in the vicinity of amino acids in the substrate-binding pocket that is Ser15, Tyr46, Gly47, Asn78, Thr79, and Thr190.
The modeled structure of the Ng-MurI protein (RacE) is shown in Figure 3. Ng-MurI complex shows that the D-glutamate moiety mostly made contact with amino acid residues in the deep cleft between the enzyme's two domains. Similar to the template, B.subtilis, the substrate-binding pocket of Ng-MurI protein (RacE) was also found to be formed by residues from the b1-a1, b2-a2, and b3-a3 loops from domain I as well as important counterparts from domain II comprising b5-a5, b6-a6 and b7-a8. Substrate Dglutamate present in the vicinity of the binding site residues Ser15, Tyr46, Gly47, Asn78, Thr79, Thr190 amino acid residues in Ng-MurI is indicated in Figure 4. Conserved residues Ser15, Tyr46, Gly47, Asn78, Thr79, and Thr190 characteristic of cell wall biosynthesis proteins form a part of the active site of Ng-MurI. The RMSD between the original conformation of substrate D-glutamate and docked conformation of substrate D-glutamate in Ng-MurI was found to be 0.2475 Å. Based on the lowest CDOCKER score (À30.5649) and RMSD value (0.2475 Å), the best docking pose was obtained for the generation of the receptor-ligand based pharmacophore model (Figure 4a).

Pharmacophore generation
The structure-based pharmacophore generation approach was explored for Ng-MurI protein to compare the important pharmacophore features. The results of docking studies were used to determine the pharmacophore features of the protein. A total of fourteen features were found in the docked glutamate racemase (RacE) where five features match the receptor-ligand interactions which were used to generate six pharmacophore models. Based on the selectivity score and number of features, the best pharmacophore was finalized with five features (4 hydrogen bond acceptors and 1 positive ionizable, AAAAP) and a selectivity score of 8.9323 ( Figure 5).

Virtual screening
The virtual screening protocol reported in this study is based on the application of sequential filters to identify a limited number of compounds that could be potential leads for the target protein. The best pharmacophore model (Pharm1) was subjected to virtual screening of the ZINC database. Build 3D database module and Screen library module of Biovia DS 2020 were used to create and screen the data set of the libraries. The cutoff value of 2.5 was decided concerning the fit value of our positive control ligand i.e. glutamate. A total of 1,95,399 compounds from a library of drug-like compounds from the ZINC database were employed. After the initial round of screening based on pharmacophore mapping, 2214 compounds were shortlisted and mapped successfully against Pharm1 (AAAAP) and were sorted based on their fit values. The pharmacophore fit value is typically displayed as the geometric fit of features to the 3D-structure-based pharmacophore model. The molecules with the highest fit score in the validated pharmacophore model should be active against our target MurI-Ng (RacE) protein. Therefore, 594 hit compounds were chosen that showed a fit value of more than 1.5 ( Figure 6).
Next, 378 hit compounds were filtered using ADMET, TOPKAT, VEBER, and Lipinski's filter and were able to pass all the pharmacokinetic parameters and were subjected to docking studies.

Molecular docking
Molecular docking is a vital step in the drug design process and is used to assess the binding ability of the hit compounds to the target (RacE) protein. 378 hit compounds obtained above were docked with modeled (RacE) protein using the CDOCKER tool of Biovia DS 2020 to evaluate their binding capacity, which satisfied the pharmacophore model's characteristics. Out of 378 successfully docked compounds, 86 hit compounds exhibited favorable binding interactions with MurI-Ng (RacE). Based on structural and biological properties, 20 hit compounds (Table ST1) were selected and redocked using the Flexible docking protocol of Biovia DS 2020 as the comparatively rigid protocol of C-DOCKER allows only limited discovery of novel conformational space. Six active site residues were selected as a group for flexible docking: Ser15, Gly47, Cys77, Thr79, Cys189, and Thr190. Among the selected top hits, best five lead compounds are shown in Figure 7.
These final five lead compounds were shortlisted and further analyzed based on their earlier known biological activity, Figure 2. (a) Structural superimposition of Ng-MurI and 1ZUW. The predicted structure of Ng-MurI (Yellow color) was superimposed on 1ZUW (Orange color) using Align and Superimpose Protein protocol of Biovia DS 2020 software. The binding site was predicted, and D-Glutamate was docked into the binding site (shown in red color). (b) The Ramachandran plot shows the percentage of residues in allowed and disallowed regions. (The plot shows 93.1% of residues were in the allowed region, 6.5% in the additionally allowed region, and 0.4% in the generously allowed region). interaction with active site residues as well as other docking parameters. The chemical structures of these shortlisted compounds are illustrated in Table 1. These top lead compounds showed good interactions with some of the substrate-binding amino acids, such as Ser15, Tyr 46, Gly47, Asn78, Thr79, and Thr190, and fit well in the active site cavity of the protein. The docking scores (CDOCKER Energy), Fit value, H-bond, and important interactions of these lead compounds are provided in Table 1.

Molecular dynamics simulations of top potential lead compounds
To gain insights into the stability and dynamic properties of the protein-ligand complexes, explicit solvent MD simulations were performed. MD simulations provided detailed insight into protein-ligand interactions in motion contributing to their stable bound conformation and visualizing the effect of ligand binding on protein conformational changes. Initially, MD simulation for 20 ns was conducted for all the top potential lead compounds complexed with Ng-MurI. Out of the total five shortlisted lead compounds, the MD simulation results for Lead 3 and Lead 5 were taken further due to the stability and consistency shown in the plots. The RMSD, RMSF, and Radius of Gyration plots for the rest of the lead compounds i.e. Lead 1, Lead 2, and Lead 4 are depicted in supplementary figure ( Figure S3). The Ng-MurI complexed with final lead compounds (Lead 3 and Lead 5) along with the control (glutamate bound with Ng-MurI) was studied further using a 100 ns MD simulations trajectory. The time evolution of the root means square deviation (RMSD) during MD simulation is used to monitor protein stability. The distributional probability of RMSD up to 100 ns trajectories is shown in Figure 8(A). The mean RMSD values of Ng-MurI-control   complex, Ng-MurI-Lead 3, and Ng-MurI-Lead 5 complex were around 0.2 nm which shows the structural stability throughout the MD run, where the RMSD values of Ng-MurI-Lead 3 were found to be more stable ( Figure 8A). Root-mean-square fluctuation (RMSF) was assessed to analyze the impact of lead-compounds binding on the flexible portion of the targeted protein. From the RMSF plots, we observed that there are slightly higher fluctuations in the amino acid regions (66)(67)(68)(69)(70)(71)(72) in the case of the control compound (glutamate), whereas Lead 3 and Lead 5 show a minimal amount of fluctuation. These regions of fluctuations for Ng-MurI-control, Ng-MurI-Lead3 complex, and Ng-MurI-Lead5 complex are shown in Figure 8(B). We also calculated the radius of gyration (Rg) value, a parameter directly associated with the overall conformational changes in the structure of the enzyme upon ligand binding, to further validate our results. It also reveals the stability, compactness, and folding behavior of the structure. We calculated the Rg values of all the selected compounds and the reference complex to determine their compactness. The average Rg values for the Ng-MurI-control complex, Ng-MurI-Lead 3 complex, and Ng-MurI-Lead 5 complex was below 1.775 nm and were stable after 30 ns implying increased compactness and improved binding ( Figure 8C). In addition to RMSD, RMSF, and Rg, the number of hydrogen bonds generated between protein and ligand throughout the simulation duration was determined. Figure 8(D) displays the graphs of these H-bonds formed between the protein and the corresponding ligand throughout the 100 ns simulation run. The average number of H-bonds formed between MurI-Ng-Lead 3 is highest followed by MurI-Ng-Lead 5 and MurI-Ng-control.

Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA)
Conformations from the last 20 ns MD simulation run were used to gain more insight into the structural dynamics of protein-ligand complexes. To validate the MD simulation findings free energies of docked complexes from two top lead compounds (Lead 3 and Lead 5) along with control were calculated using molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) of g_mmpbsa tool ( Figure 9). All the top three complexes exhibited negative binding energy in the MM-PBSA threshold. As shown, out of the two lead compounds, Lead 3 (ZINC000014657962) exhibits a better net binding energy score which may be comparable with the flexible CDOCKER energy score (À7.76499 kcal/ mol) ( Table 1). The binding energy score for Lead 3 depicted by MM-PBSA exhibits better protein-ligand binding as compared to Lead 5.

Discussion
Ng-MurI is an indispensable enzyme that acts in bacterial peptidoglycan biosynthetic pathway and is also a nonhuman analog. Therefore, it is an important target for developing new antibiotics and novel drug molecules. In this    current study, we demonstrated that Ng-MurI, a potential drug target, may be employed to develop anti-gonococcal drugs. This study also signifies the promising candidature of the lead compounds identified through virtual screening of the libraries from ZINC and the stability of the same was confirmed with the help of Molecular Dynamics simulation studies performed on the lead compounds along with the control molecule for 100 ns. The structural dynamics of the protein-ligand complex were further elucidated with the help of MM-PBSA calculations.
Since the crystal structure of glutamate racemase for Neisseria gonorrhoeae was not available, therefore, we took the path of generating the homology model for the same using various authenticated web servers. The template to generate the model for the Ng-MurI was selected based on the query coverage and percentage similarity score and 1ZUW (Crystal structure of B. subtilis glutamate racemase with D-Glutamate) was selected as the template. The purpose behind developing the model of this target protein was to identify the ligand binding active sites in the modeled protein structure. In this study, we also characterized the model of the Neisseria gonorrhoeae glutamate racemase (Ng-MurI) using a conventional drug discovery pipeline. Virtual screening followed by docking studies ascertained that ZINC000014657962 (Lead 3) is one of the best ligands which exhibited good docking energy scores, and showed maximum stability and binding energy during simulation studies.
A series of molecular docking steps retrieved top leads with good fitness, docking scores, fit values, and interaction patterns among residues. The overall analysis of the results suggests that the structure-based pharmacophore modeling provides useful information required for a proper understanding of the important structural and binding features for designing novel Ng-MurI inhibitors.

Conclusion
The emergence of resistant strains of Neisseria, in the absence of any vaccine, is a threatening public health crisis as it could significantly increase reproductive health complications (pelvic inflammatory disease, ectopic pregnancy, epididymitis in men, increased infertility) neonatal blindness, leading to increased socioeconomic cost. Hence, there is an urgent need to identify new drug targets and develop new drugs/antibiotics to effectively treat the infection.
This fact has prompted us to identify a set of inhibitors for Ng-MurI using the standardized computational approaches which include homology modeling, pharmacophore-based virtual screening, molecular docking, and MD simulations. We report the identification of two potential compounds which may be taken further for Ng-MurI in-vitro inhibition studies. These identified lead compounds can provide a platform for further detailed efficacy analysis, leading to discovering and designing novel drug molecules against Ng-MurI specifically for further studies in multidrug-resistant gonococcal infection. This study provides momentum for future investigations in finding novel compounds that can act as antibacterial agents against resistant microbes.