Comprehensive screening of marine metabolites against class B1 metallo-β-lactamases of Klebsiella pneumoniae using two-pronged in silico approach

Abstract The emergence of antibiotic resistance is one of the major global threats in healthcare. Metallo-β-Lactamases (MBL) are a class of enzymes in bacteria that cleave β-lactam antibiotics and confer resistance. MBLs are further divided into subclasses B1, B2 and B3. Of these, subclasses B1-MBLs (including NDM-1, VIM-2 and IMP-1) constitute the clinically prevalent lactamases conferring resistance. To date, no effective drugs are available clinically against MBLs. In this work, we aim to identify potent inhibitors for the B1 subclass of MBL from available marine metabolites in Comprehensive Marine Natural Product database through integrated in silico approaches. We have used two methods, namely, the high-throughput strategy and the pharmacophore-based strategy to identify potential inhibitors from marine metabolites. High-throughput virtual screening identified N-methyl mycosporine-Ser, which had the highest binding affinity to NDM-1. The pharmacophore-based approach based on co-crystallized ligands identified makaluvic acid and didymellamide with higher binding affinity across B1-MBLs. Taking into account of the advantage of a pharmacophore model-based approach with higher binding affinity, we conclude that both makaluvic acid and didymellamide show potential broad-spectrum effects by binding to all three B1-MBL receptors. The study also indicates the need to take multiple in silico approaches to screen and identify novel inhibitors. Together, our study reveals promising inhibitors that can be identified from marine systems. Communicated by Ramaswamy H. Sarma


Introduction
One of the main concerns in the healthcare sector is the development of bacterial resistance against antibiotics (Mancuso et al., 2021).The main clinical pathogens that exhibit antibiotic resistance include Acinetobacter spp., Pseudomonas spp., Klebsiella pneumoniae, Escherichia coli and Enterobacter spp.These pathogens cause severe and often fatal infectious diseases such as bloodstream infections and sepsis (Mancuso et al., 2021).Cephalosporin and carbapenems are the main class of antibiotics with antimicrobial effects against Enterobacter class of bacteria.Among these, the carbapenem-resistance strains of K. pneumoniae are the most clinically prominent resistant strains (De Oliveira et al., 2020).Carbapenems possess a broad-spectrum antibiotic effect against b-lactamase-producing Gram-positive and Gram-negative bacteria (Meletis, 2016).Carbapenems bind to and inhibit penicillin binding proteins inhibiting the formation of peptidoglycans in the cell wall (Papp-Wallace et al., 2011).However, over-prescription of carbapenems as the antimicrobial agent has resulted in bacteria gaining resistance through multiple mechanisms including production of b-lactamases, efflux pumps etc., conferring resistance to carbapenems (Papp-Wallace et al., 2011).Overall, the development of antibacterial resistance is estimated to lead to about 10 million deaths per year by 2050 (Mancuso et al., 2021).A multi-prong approach is essential to tackle this alarming development of bacterial resistance against antibiotics.Among this, the development of novel inhibitors for specific enzymes that confer antibacterial resistance is one approach.b-lactamases are one class of enzymes that confer resistance to b-lactam antibiotics such as penicillin and cephalosporins.These enzymes cleave the b-lactam rings and thus render the antibiotic ineffective.Broadly, b-lactamases are classified into Metallo-b-lactamases (MBL) and Serine b-lactamases (SBL) based on their hydrolysis activity mechanism.Among the groups of MBLs, diverse range of sequences with 25% is observed.However, they are structurally similar with characteristic ab/ba sandwich fold with the active site located at the interface between domains including either one or two zinc ions in the catalytic domain (Palzkill, 2013).MBLs are classified into three major subclasses, B1, B2 and B3, based on their sequence identity and zinc ion dependencies.B1 subgroup is further classified into three subclasses based on their amino acid homology viz., New Delhi metallo b-lactamases (NDM), VIM and IMP (Palzkill, 2013).Strains harboring MBLs, particularly NDM, have been reported to be a serious threat because there is no standard detection test, and their plasmids are highly transferrable causing widespread transmission among bacterial populations, and a higher probability of unrecognized asymptomatic carriers (Khan et al., 2017;Yang et al., 2020).Kang et al. identified novel scaffolds for the development of both enzyme-specific and broad-spectrum MBL inhibitors in Serratia marcescens and Pseudomonas aeruginosa (Kang et al., 2018).Consequently, many experimental and in silico works have attempted to discover inhibitors against NDMs, and in general against MBLs.Chandar et al. showed that specific plant extracts can act against NDM-1 producing Escherichia coli (Chandar et al., 2017).Khan et al. identified four compounds as potential inhibitors that exhibit high affinity for NDM-1 producing Klebsiella pneumoniae through virtual screening and microbiological analysis (Khan et al., 2017).Wang et al. studied the thermokinetic profile of NDM-1 and its inhibition by carboxylic acids (Wang, 2018).They also identified a specific NDM-1 inhibitor, isoliquiritin from Glycyrrhiza uralensis that enhances the activity of meropenem against NDM-1-positive Enterobacteriaceae through in vitro experiments.Their results showed that isoliquiritin in combination with meropenem could be used as a potential 'lead compound' for the development of NDM-1 inhibitors (Wang et al., 2020).Spyrakis et al., identified three compounds exhibiting broad-spectrum b-lactamase inhibitor activity from a series of non-covalent scaffolds targeting SBLs and concomitantly MBLs including NDM-1 and VIM-2 producing P. aeruginosa, through a large virtual screening study and in vitro validation (Spyrakis et al., 2020).B. Liu et al. identified taniborbactam, a boronic-acidcontaining compound as a broad spectrum Serine-and Metallob-lactamase inhibitor for carbapenem-resistant infections.The combination of taniborbactam with fourth-generation cephalosporin, cefepime, against Pseudomonas aeruginosa and carbapenem-resistant Enterobacteriaceae is currently in phase 3 clinical studies (Liu et al., 2020).Hinchliffe et al. identified novel and efficient cross-class metallo-b-lactamase inhibition by bisthiazolidines revealing multiple binding modes in Escherichia coli-producing MBLs (Hinchliffe et al., 2016).Palacios et al. described a number of MBL substrate-mimicking inhibitors that contain combinations of hydroxyl, ketone, ester, amide, or sulfonyl groups.The study also reported that bicyclic boronates as a promising new scaffold for MBL inhibitors (Palacios et al., 2020).Pedroso et al. showed that the anti-Helicobacter pylori drug called colloidal bismuth subcitrate (CBS) and its related Bi (III) compounds irreversibly inhibited different types of MBLs under bacterial phyla: Proteobacteria, Actinobacteria, Bacteroidetes and Firmicutes (Pedroso et al., 2020).Abbas et al., showed that inhibition of carbapenemase producing Enterobacteriaceae activity by the combined therapies with carbapenems and tigecycline/polymyxins (Abbas et al., 2019).Sofia Maraki et al., reported efficient combinatorial drug treatment against multidrug resistant MBLproducing Klebsiella pneumonia (Maraki et al., 2021).Taken together, these studies indicate that there have been several efforts to identify novel inhibitors against MBLs in various pathogenic organisms including K. pneumoniae.However, in the case of K. pneumoniae, so far, efforts have been made to use combination therapies against MBL-producing K. pneumoniae, which have been shown to have varying efficacy.Of the subclasses of MBLs, the B1 classes of MBLs are significantly important enzymes as they are resistant to almost all the known antibiotics.Such resistance is an alarming threat because the genetic element of such MBLs spreads over a wide spectrum compared to other b-lactamases.Therefore, there is an urgent need to identify inhibitors with high efficacy.In this work, we focus on identifying inhibitors against the B1 subclass of MBLs in K. pneumoniae.
Marine natural products are a high-volume source of bioactive secondary metabolites with the characteristic chemical novelty to treat several diseases.The enormous novelty in the chemical structure and function makes marine natural compounds an important drug component for consideration (Montaser & Luesch, 2011).Many studies showed the potential of marine natural compounds as drug molecules as it possesses antioxidant, antibacterial, antitumor, antivirus, anticoagulant, anti-inflammatory properties (Barzkar et al., 2019).Kong et al. showed that marine natural products are superior to terrestrial natural products with chemical novelty through comparative analysis (Kong et al., 2010).Ziconotide, the first marine-derived peptide drug on the market, has a potent analgesic effect, inhibiting the activity of a subset of neurons including pain-sensing primary nociceptors.Taking lead from these findings, in this work, our objective is to identify the effective inhibitors against B1-MBLs through molecular docking analysis using both High Throughput Virtual Screening (HTVS) and pharmacophore-based-screening from the marine metabolites made available in the Comprehensive Marine Natural Product Database (CMNPD) (Lyu et al., 2021).In this work, we identified potent inhibitors from CMNPD using two approaches, viz., High Throughput Virtual Screening (HTVS) and pharmacophore-based-screening.We validated the identified inhibitors using Güner-Henry's Goodness of Hits approach and molecular dynamics simulation.

Ligands for high throughput virtual screening & pharmacophore generation
Workflow of the methodology used for identification of inhibitors for B1-MBLs from marine sources is given in Figure 1.The marine metabolites in the CMNP database were downloaded in sdf format (Lyu et al., 2021).There were approximately �32000 molecules.These ligands were used for high throughput virtual screening and pharmacophore-based screening.For generation of the pharmacophore model, we used two different sets of ligands to generate two different pharmacophore models-(i) the FDA-approved drugs and (ii) the co-crystallized ligands that are reported to inhibit MBLs.A total of 32 drugs (Table 1S, Supplementary material) including higher order carbapenems and cephalosporins that act against b-lactamases and are approved by the FDA for clinical applications against MBLs were collected from the literature.These were used to generate the FDA-based pharmacophore model.The co-crystallized ligands of the Metallo-b-lactamase class of protein structures in all available organisms (Table 2S, Supplementary material) were collected from the b-lactamases database (BLDB) (Naas et al., 2017).For this, we first searched for MBLs that were solved using the X-ray diffraction method and had no mutations.The co-crystallized ligands from the filtered protein structures with reported inhibitory properties were collected.This resulted in a total of 139 co-crystallized ligands for pharmacophorebased screening.

Protein & ligand preparation
Three representative structures of the B1 class of MBLs of K. pneumoniae were obtained from the Protein Data Bank -NDM-1 (PDB ID: 4EXS), VIM-2 (PDB ID: 5O7N) and IMP-1 (PDB ID: 6C6I).These proteins structures have no mutations and were solved using X-ray diffraction method.Furthermore, Chain A from each of the structures with the zinc metal in the ligand binding site excluding other heteroatoms was considered for docking and simulations.The protein preparation wizard in Maestro of the Schr€ odinger suite was used to prepare the input protein structure (Sastry et al., 2013).This procedure includes fixing the missing hydrogen atoms, ensuring a proper charge/ionization state, protonation state for histidine residues, capping the termini, loop refinement through Prime, optimizing the hydrogen bonds and energy minimization using the OPLSe force field.Water molecules in the protein beyond 4 Å were removed.
LigPrep module of Schr€ odinger suite was used to generate the three-dimensional coordinates and conformers for all the compounds from the CMNPD.Tautomeric states were generated with Epik which is based on the more accurate Hammett and Taft methodology (Sastry et al., 2013;Shelley et al., 2007).

HTVS-based molecular docking
HTVS was performed using the Glide module of Schr€ odinger Suite (v2021-3) against all the three receptors for the compounds from CMNPD.For this, we first generated the grid using the Glide module of the Schr€ odinger suite.The grid center was constructed based on the position of the ligand in the x-ray crystal structure of the protein.In case of IMP-1, the grid was generated using the sitemap of the Schr€ odinger suite (Halgren, 2009).Docking was then performed to this grid using the Glide module-Schr€ odinger (v2021-3).The Glide ligand docking procedure is known to reliably docks ligand compounds in the correct binding site with correct pose (Friesner et al., 2004;Halgren et al., 2004;Friesner et al., 2006).

Pharmacophore-based screening and molecular docking
Pharmacophore-based feature generation was performed using Phase module of Schr€ odinger Suite (v2021-3).Pharmacophore modeling involves determining the spatial arrangement of chemical features which interact with the binding site of the receptors (Dixon et al., 2006;Dixon et al., 2006).Pharmacophore model is generated for the set of co-crystallized ligands and FDA-approved drugs with hypothesis threshold of 30% and 50%, respectively.The default pharmacophore features in the module include -hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic (H), negative (N), positive (P) and aromatic ring (R).Default settings were used to generate the pharmacophore hypothesis and the maximum number of sites was set to 5 and minimum to 3. The pharmacophore box size was adjusted to 2 Å and the molecules were scored for a given pharmacophore using the default weights of the scoring parameters.The generated model was used for screening against CMNPD (Lyu et al., 2021).Screening for potential hits against generated pharmacophore model of co-crystallized ligands and FDA-approved drugs was performed using the Schr€ odinger Suite (v2021-3) Shape screening (Sastry et al., 2011).This method works by identifying the similarity from superimposition between two structures.The similarity score is calculated as the ratio of the occupied volume to the total volume in the range between 0 and 1.
where A and B are the two structures; V A\B is the shared or jointly occupied volume of the structure and V A[B is the total volume of the structures.Molecular docking was performed using Glide module of Schr€ odinger Suite (v2021-3) after pharmacophore-based shape screening for both sets of co-crystallized ligands and FDA-approved drugs (Friesner et al., 2006).

Validation of the docking based on Güner-Henry's goodness of hits approach
The Güner-Henry's Goodness of Hits approach is useful to validate the docking and the pharmacophore model.In this approach, a set of decoys are constructed based on the active molecules (that are known to have the desired property).These decoys are then subjected to virtual screening/ pharmacophore-based screening.In our case, we used the 32 FDA-approved antibiotics (Table 1S, Supplementary material) as the actives for the virtual screening procedure.
The decoy sets for these molecules were generated were generated using the dUD-E Webserver (Mysinger et al., 2012).This resulted in 1320 decoy molecules.The docking score for these decoy molecules were calculated using the same approach as that of ligands from CMNPD.To validate the pharmacophore-based model, we generated the decoy set assuming the co-crystallized ligands as the actives (Table 5S, Supplementary material).The dUD-E webserver was used to generate the decoy set.A total of 1672 decoy molecules were generated.This was then followed by pharmacophore screening of the decoy set.The Güner-Henry score (goodness of hit -GH) in each case was then calculated using Equation (2) (Sirous et al., 2020): where D is the total number of molecules including both actives and decoys, A is the number of actives considered for decoy generation, Ht represents the total number of topranked hits identified from screening considered, Ha is the total number of active molecules (A) with docking score better than that of the top-ranked hits.

Molecular dynamics simulation
Molecular dynamics simulations of the protein-ligand complexes of the highest-binding ligands were performed using Desmond (Schr€ odinger, 2020; Bowers et al., 2006).The details of the simulated systems are given in Table 3S (Supplementary material).Each protein-ligand complex was placed in a cubic box.The size of the cubic box was set such that the distance between the outermost atom of the protein and the box was 10 Å.The TIP3P water model was used for solvation.The system was neutralized by adding appropriate number of counter ions (Naþ/ClÀ ) positioned randomly in the box.Minimization of the protein-ligand complexes were performed using minimization of the system using a hybrid method of the steepest decent and the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithms.Energy minimization was carried out for 500ps with maximum steps of 5000 and convergence threshold of 1.0 kcal/ mol/Å.Minimization was followed by NPT simulations for 10 ns and final production run for 100 ns were performed using OPLSe force field parameter.The temperature of the system was maintained at 300 K and pressure of 1.01325 bar was maintained using Nose Hoover coupling and isotropic scaling, respectively.
Trajectories were written every 10 ps.

High throughput virtual screening (HTVS)
HTVS of the marine metabolites was done against the three target proteins in the class B1-MBLs.The ligand interaction diagram of the top five molecules based on the Glide score is given in Figures 2-4 for each of the proteins respectively.The results showed a novel scaffold of inhibitors against MBLs.In case of the NDM-1, top docking score of À 11.31 kcal/mol was observed for N-Methyl mycosporine-Ser.Mycosporine-like amino acids (MAA) are typically present in red & green algae and corals.MAAs are the secondary metabolites found in different organisms inhabiting ecosystems with a high concentration of sunlight, such as marine and freshwater environments, to protect against solar radiation (Chrapusta et al., 2017).MAAs are also considered to be potential anti-cancer agents due to their anti-proliferative activities on neoplastic cells (Wada et al., 2015).Other than N-Methyl mycosporine-Ser, Isonamidine D (À 10.831 kcal/mol), 3-(4-hydroxyphenyl)-4phenyl-1H-pyrrole-2,5-dicarboxylice acid (À 10.726 kcal/mol) and Leucyl 4-hydroxyproline (À 10.632 kcal/mol) displayed higher docking score with novel chemical scaffold.In NDM-1, Gln123 showed water-mediated interactions with top-scoring ligands and zinc ions exhibited interactions with carbonyl oxygen of top-scoring ligands.Zinc ions at the binding site of the MBLs were known to activate both hydrolytic water and the b-lactam carbonyl.In addition, these ions also participate in substrate binding via the carboxylate present in most b-lactam antibiotics (Pettinati et al., 2016).
In case of IMP-1, HTVS results showed Sydowiol B with the highest docking score (À 8.66 kcal/mol).However, it displayed no interaction with zinc ions in the binding pocket.Similarly, other top scoring marine metabolites such as spartinol (À 8.343 kcal/mol), penic acid (À 7.961 kcal/mol), nitrosporusine (À 7.519 kcal/mol) and sysowin A (À 7.332 kcal/mol) did not show zinc-mediated interactions at the binding site.However, we found that all these ligands showed interaction with Lys216 in the binding pocket.Reports on the Sydowiols A and C extracted from the marine-derived fungus, Aspergillus sydowii in the East China Sea showed inhibitory activity in the crude fungal extract against protein tyrosine phosphatase A of Mycobacterium tuberculosis (Ji et al., 2021).

Validation using Güner-Henry (GH) score
For the HTVS validation, the dataset with the FDA-approved drugs as actives and the decoys generated based on these FDA drugs was used.HTVS screening was carried out on this dataset along with the top twenty hits from the CMNP database for each protein.Similarly, for the pharmacophore-based screening, the dataset consisting of co-crystallized ligands was used as the actives and the decoys generated based on these ligands were used along with the top twenty hits from the CMNP database for each protein.The number of actives that scored more than the hits from the CMNPD were identified and was used to calculate the GH score.The percentage of decoy molecules that had a docking score greater than that of the top hits was also calculated.The results are provided in  Table 1.We observe that the GH score ranges between 0.63 and 0.75.A GH score of 0 is considered to be a null model, while a GH score 1 denotes the generation of an ideal model.Since our GH scores are above 0.5, we can consider the models to be valid.Furthermore, the percentage of decoy molecules that scored higher than the top hits is less than 2% in all but one case and has a maximum of 8.73%.This again indicates that the screening models are reliable.

Molecular dynamics simulations for top-scoring ligands from HTVS
In addition to the docking scores, we evaluated the conformational stability and flexibility of the ligands in the binding pockets using molecular dynamics simulations.The proteinligand complexes with the top-scoring ligand for each target were simulated.Figure 5 gives the Root Mean Squared Deviation (RMSD) of the protein C-a atoms, Root Mean Square fluctuation (RMSF) of residues and RMSD of ligands as a function of time for each complex simulated for 100 ns each.In case of RMSF, residue numbers 25-32 and 173-186 showed maximum fluctuations.
These residues form loops around the binding pocket of the target proteins and hence show higher fluctuation.However, the overall complex is stable (does not undergo large conformational change) as evident from the RMSD graph.
Visual inspection (Figure 6) showed that Anguibactin changed its conformation at about 40 ns, explaining the change in the RMSD.However, the complex remained stable with ligand in the binding pocket throughout the simulation.
In addition to analyses of conformational stability and flexibility of protein and ligands, the percentage time of the ligand   interactions at the binding site was calculated and given in Figure 7. Interactions that occur more than 30.0% of the simulation time throughout the trajectory (0 through 100.0 ns) were shown.The interaction profile showed a water-mediated interaction of the binding site residues with the ligands.In addition to water-mediated interactions, the top-scoring ligand of NDM-1 and VIM-2 displayed interactions with zinc that existed more than 99% of the simulation period.In both NDM-1 and VIM-2, the binding of ligand with HIS and ASP residues in the binding interface was observed.
Figure 8 shows the interaction fractions of the ligand with the amino acid residues of the proteins such as H-bond interactions, hydrophobic interactions, ionic interactions and water bridge interactions.The stacked bar charts are normalized over the course of the simulation trajectory (100 ns).Hydrogen bond interactions and water mediated interactions were found in all the three receptors-ligand complexes.In cases of NDM-1 (HIS189, CYS208 & HIS250) and VIM-2 (HIS179, CYS198 & HIS240), histidine and cysteine contributed for interactions between ligands and the metal ions in the binding pocket.However, in case of IMP-1, hydrophobic interactions along with hydrogen bond interactions and water mediated interactions were observed.Since hydrogen bonded interactions were observed in all the three receptors and its complexes, we further analyzed the number of hydrogen bonds between protein and the ligands (Figure S10, Supplementary material).Hydrogen bond interaction of identified compound with its receptor shows their strong influence on drug specificity and adsorption.Thus, analyses from the molecular dynamics simulation of the top-scoring compounds show that the identified compounds were stable and can be further used considered for lead optimization of the drug discovery process.

Pharmacophore-based screening and docking
Pharmacophore-based feature generation following screening CMNPD against the three protein targets was performed.Pharmacophore features were generated with two sets of  molecules-(1) ligands in the co-crystallized structures of MBLs; (2) FDA-approved known antibiotics.Figure 9 shows the pharmacophore model generated for co-crystallized ligands and FDA-approved drugs with a hypothesis threshold of 30% and 50%, respectively.Once the feature generation was performed, the generated features were used to screen the CMNPD to identify potential inhibitors.The screened molecules were docked against three protein targets using Glide-Schr€ odinger.The list of PDB IDs of the structures with co-crystallized ligands and FDA-approved drugs were given in Tables S1 and S2 (Supplementary material).

Pharmacophore-based screening with co-crystallized ligands
The results for the top-scoring ligands (both co-crystallized ligands and FDA-approved drugs) docked with protein targets are summarized in Tables 2-4.From the features generated with co-solved ligands, Makaluvic acid was observed to occur as one among the top five ligands binding to each of the three protein targets of NDM-1, VIM-2 and IMP-1.In addition to makaluvic acid, Didymellamide was observed to show a higher docking score with NDM-1 and IMP-1, respectively.Studies show that derivatives of makaluvic acid isolated from Strongylodesma aliwaliensis, exhibiting in vitro cytotoxicity against esophageal cancer cells (Schr€ odinger, 2020).Similarly, pyridine alkaloid didymellamide isolated from Stagonosporopsis cucurbitacearum showed significant activities against azole-resistant Candida albicans (Haga et al., 2013).
The ligand interaction diagram of pharmacophore-based screening and docking of three protein targets were given in Supplementary information.Pharmacophore-based screening against co-crystallized ligands showed interactions with HIS, ASP, TRP and VAL along with Zn metal ions at the binding site of NDM-1 and VIM-2.For IMP-1, HIS, ASP, TRP and VAL in the binding site were observed to have interactions with the ligands (Figure S1-S3, Supplementary material).Validation of the docking function using Güner-Henry Goodness of hit approach was detailed in Section 3.2 for reference.

Pharmacophore-based screening with FDA-approved drugs
In the case of pharmacophore-based screening for FDAapproved drugs, the interaction diagram showed that the interactions were primarily with the Zn ion and HIS residue  at the binding site of NDM-1 and VIM-2.From the features generated with FDA-approved drugs, Damirone showed higher binding affinity with NDM-1 and VIM-2 respectively.Similar to the case of IMP-1 in the pharmacophore screening based on co-crystallized ligands, no Zn based interactions were observed.However, HIS and SER interactions with the ligands at the binding sites were observed (Figure S4-S6, Supplementary material).

Molecular dynamics simulations for top-scoring ligands from pharmacophore-based approach
We performed molecular dynamics simulation for all three protein structure with makaluvic acid to check the conformational stability and flexibility of the protein-ligand complexes.Makaluvic acid B and Makaluvic acid C were visualized and observed to have similar structure.It was considered as same molecule and used as makaluvic acid in the following sections.Figure 10 shows the conformational stability and flexibility of three protein structures with makaluvic acid.In all the cases of RMSD of protein C-a atoms (Figure 10A), RMSD of ligand (Figure 10C) and RMSF of the protein residues (Figure 10B) remained stable throughout the simulations of 100 ns.As discussed in the Section 3.3, RMSF residue numbers between 25 and 32(GLU, VAL, ASN, GLY, TRP, GLY, VAL, VAL) & between 173 and 186 (LYS, ALA, LYS, SER, LEU, GLY, ASN, LEU, GLY, ASP, ALA, ASP, THR, GLU) showed maximum fluctuations along with the fluctuations in the N-terminal regions.These regions were identified as the loops around the binding pocket of the receptors which generally fluctuates more compared to the secondary structures of the protein.
Thus, the complexes are stable throughout the simulations.
In addition to analyses for conformational stability and flexibility, interaction profile of makaluvic acid with NDM-1, VIM-2 and IMP-1 with their percentage of contact was given in Figure 11.The result showed an HIS interaction with ligands in all three protein structures.Along with HIS, Zn metal ion interaction was also observed in the case of NDM-1 and VIM-2, respectively.Similarly to the interaction profile of IMP-1 based on the HTVS approach, IMP-1 in the pharmacophore-based approach also did not show Zn ion interaction at the protein binding site.However, HIS in the binding pocket was consistent in making interactions with makaluvic acid in all three protein receptors.
Figure 12 shows the interaction fractions of the ligand with the amino acid residues of the proteins such as H-bond interactions, hydrophobic interactions, ionic interactions and water bridge interactions.The stacked bar charts are normalized throughout the simulation trajectory (100 ns).Here, ionic interactions between ligand and metal ion in the binding pocket are dominant in cases of NDM-1 and VIM-2, whereas the hydrogen bond interactions and water mediated interactions are dominant for IMP-1 receptors (Table S10, Supplementary material).

Molecular docking of FDA-approved drugs
In addition to both approaches, we also performed docking of FDA-approved drugs with the three protein structures.The results were summarized in Table S3 (Supplementary material).
Identified top-scoring molecules obtained from HTVS approach and pharmacophore-based screening approach drugs showed better ligand interactions and docking score than with FDA-approved.This shows that the identified compounds can be used as inhibitors against antimicrobial resistance.In case of NDM-1 and VIM-2, the carbapenem class of antibiotics showed higher binding affinity.Top-scoring drugs with NDM-1 were ertapenem with docking score of À 9.407 kcal/mol.In case of VIM-2 and IMP-1, top-scoring drugs showed docking score of À 12.737 kcal/mol (Tomopenem) and À 10.075 kcal/mol (Lenopenem), respectively.The ligand interaction diagram of the binding was given in Figure S7-S9 (Supplementary material) for each protein target.The profile of NDM-1 and VIM-2 binding showed a zinc interaction with drugs that was observed with both approaches -HTVS and pharmacophore-based screening.This shows that even though sequence similarity is observed among NDM-1, VIM-2 and IMP-1, the patterns of binding and interactions are unrelated.The prediction accuracy depends on the similarity of the input compound to compounds with known LD50 values as well as the hit rates obtained in a cross-validation study.
3 logP is the Octanol/water partition coefficient ranges between -3 (very hydrophilic) and þ10 (extremely lipophilic/hydrophobic). QPlogS -Predicted aqueous solubility, log S where S is the concentration of the solute in a saturated solution that is in equilibrium with the crystalline solid with accepted values in the range -6.5 -0.5.
3 QPPCaco -Predicted apparent Caco-2 cell permeability in nm/sec.Caco-2 cells are a model for the gut-blood barrier.QikProp predictions are for non-active transport in which value of 500 is great. 4 QPlogBB -Predicted brain/blood partition coefficient for orally delivered drugs and accepted range is between -3.0 and 1.2.
5 HumanOralAbsorption (HOA) -Predicted qualitative human oral absorption: 1, 2, or 3 for low, medium, or high.The assessment uses a knowledge-based set of rules, including number of metabolites, number of rotatable bonds, logP, solubility and cell permeability. 6 Percent Human Oral Absorption (% HOA) -Predicted human oral absorption on 0 to 100% scale.The prediction is based on a quantitative multiple linear regression model.This property usually correlates well with Human Oral Absorption, as both measure the same property with accepted value at >80%.Compounds with fewer (and preferably no) violations of these rules are more likely to be orally available with maximum of 3.

Drug-like property of the potential ligands
The identified molecules from CMNPD with B1-MBL class of Klebsiella pneumoniae were further checked for their drug-like properties using the QikProp of Schr€ odinger Suite (2022) (Schr€ odinger, 2021).The drug-like properties from QikProp for the molecules were given in Table 5. QikProp predicts the Absorption, Distribution, Metabolism, Elimination (ADME) properties of all natural compounds.Drug-like properties for all the top ligands were found to be in the accepted range for all the QikProp properties.Along with the QikProp properties, toxicity of the identified molecules was predicted using the PROTOX web server (Banerjee et al., 2018).Table 6 shows that the lead molecules either fall into class 4 or class 5. Class 4 denotes that the LD50 value (LD50 is the dose at which 50% of test subjects die upon exposure to a compound) lies between 50 and 2000, whereas class 5 denotes that the LD50 value lies between 2000 and 5000.Both these classes are said to be nontoxic, however, class 4 is harmful if swallowed and class 5 may be harmful if swallowed.Hence the identified molecules can be taken as lead molecules for further lead optimization.

Conclusions
In summary, molecular docking of B1-MBLs via two approaches, namely HTVS and pharmacophore-based screening with marine metabolites, reported promising ligands to consider in lead optimization.N-Methyl-mycosporine-Ser, Anguibactin and Sydowiol B from the HTVS approach based molecular docking are novel molecular scaffolds.However, molecules with broad range effect (binding to all three proteins) from the HTVS-based approach is not observed.From pharmacophore-based screening using co-crystallized ligands, Makaluvic acid B/Makaluvic acid C showed binding with all three protein structures indicating its potential to be a broad-range inhibitor.In the case of pharmacophore-based screening using FDA-approved drugs, Damirone is observed as a potential inhibitor for the NDM-1 and VIM-2 proteins.Molecular Dynamics Simulations of molecules obtained from both approaches showed stability throughout the simulation.Even though the identified top-scoring ligands showed lower docking score compared to the known antibiotics, this study has its advantages: (1) identification of marine metabolites in the context of antibacterial property is less explored, (2) Pharmacophore-based approach identified broad spectrum effect against B1-MBLs, (3) HTVS approach identified novel scaffold of molecule against B1-MBL inhibition expanding the library for antimicrobial drugs.From our study, Makaluvic acid could be used as the lead molecule to test the broad-spectrum effect against MBL.So far, this is the first study to identify inhibitors from marine metabolites with broad-spectrum effect for B1-MBL against antibiotic resistance.Thus, two approaches resulted in two sets of ligand molecules that could be explored further to be used against inhibition of class B1-MBL.

Figure 1 .
Figure1.The two-pronged approach for the identification of inhibitors for B1-MBLs from marine sources.

Figure 5 .
Figure 5. (A) Conformational deviation of Protein C a atoms as a function of time, (B) Fluctuations of amino acid residues, and (C) conformational deviations of ligand molecules as function of time in the protein-ligand complexes simulated for 100 ns.

Figure 6 .
Figure 6.VIM-2 (PDB ID: 5O7N) complexed with anguibactin showing changes in the conformation during molecular dynamics simulation.Green color represents conformer of anguibactin at 35 ns of simulation and red color represents conformer of anguibactin at 45 ns of simulation.

Figure 7 .
Figure 7. Ligand Interaction Profile for top ligands from three target protein structures as a function of time (100 ns) using HTVS approach.

Figure 8 .
Figure 8. Interaction fractions for top ligands with amino acid residues throughout the simulation of 100 ns for HTVS approach.

Figure 9 .
Figure 9. Pharmacophore model generated from co-crystallized ligands and FDA-approved drugs.

Figure 10 .
Figure 10.(A) Conformational deviations of Protein C a atoms as a function of time, (B) Fluctuations of amino acid residues, and (C) conformational deviations of ligand molecules as function of time in the protein-ligand complexes simulated for 100 ns.

Figure 11 .
Figure 11.Simulation interaction diagram of inhibitors obtained from Pharmacophore-based screening and docking against B1-MBL throughout 100 ns.

Figure 12 .
Figure 12.Interaction fractions for top ligands with amino acid residues throughout the simulation of 100 ns for Pharmacophore-based approach.

78
Rule of five -Number of violations of Lipinski's rule of five.The rules are: molecular weight < 500, QPlogPo/w < 5, donor H-Bond � 5, and acceptor H-Bond � 10.Compounds that satisfy these rules are considered drug-like with maximum is 4.Rule of three -Number of violations of Jorgensen's rule of three.The three rules are: QPlogS > -5.7, QP PCaco > 22 nm/s, # Primary Metabolites < 7.

Table 1 .
Validation of the identified hits using Güner and Henry's 'goodness of hits' approach.

Table 2 .
Glide score for top ligands obtained from pharmacophore model generation and screening for NDM-1 (PDB ID: 4EXS) against CMNPD.

Table 3 .
Glide score for top ligands obtained from pharmacophore model generation and screening for VIM-2 (PDB ID: 5O7N) against CMNPD.

Table 4 .
Glide score for top ligands obtained from pharmacophore model generation and screening for IMP-1 (PDB ID: 6C6I) against CMNPD.

Table 6 .
Toxicity prediction for the identified hits using PROTOX server.
1 QPlogPw -Predicted water/gas partition coefficient with accepted values in the range 4.0 -45.0.2