Mutations in passive residues modulate 3D-structure of NDM (New Delhi metallo-β-lactamase) protein that endue in drug resistance: a MD simulation approach

Abstract The ability of antimicrobial resistance developed by bacteria enhanced the complexity of bacterial treatment leading a serious threat to human health. Production of β-lactamase by bacteria that inactivates β-lactam is a generic cause of resistance. One such β-lactamase enzyme is New Delhi Metallo-β-lactamase (NDM) which is recently reported to have clinically more importance and recognized as an antibiotic resistance marker. Mutations in active and passive residues of NDM protein play a fateful role in the substrate and inhibitor specificity. In this study, in silico point mutations of residues near the active site and flexible regions of protein were investigated. Hybrid modelling and molecular dynamics (MD) simulations were carried to build up the mutant models and monitored structural stability. Molecular docking results articulated that mutant proteins had lesser binding affinities with methicillin, oxacillin and doripenem drugs. Further, to scrutinize the structural alterations and rescore the binding energies per-residue basis, MD simulations of wildtype (WT) and mutant (MT) NDM proteins with methicillin, oxacillin and doripenem were performed. Our results demonstrated that mutations in N193A, S217A, G219A and T262A residues led to protein destabilization and amend their binding affinities with methicillin, oxacillin and doripenem. The present study exploited computational approaches which displayed differential binding of drugs with WT and MT NDM proteins that confer resistance to oxacillin and doripenem. The study features the significance of passive residues, thus provides a clue to accelerate the process of designing an ergastic antibiotic against NDM protein. Communicated by Ramaswamy H. Sarma


Introduction
Antibiotics resistance is the ability of microorganisms or bacteria to resist the effects of antibiotics. Higher concentration of same drug is required when microbes become resistant or less sensitive (Singer et al., 2016) and emergence of antibiotic resistance appeared shortly after the administration of new drugs (Richardson, 2017). Antibiotics cease the growth of microbes by interfering various processes such as b-lactams and glycopeptide agents interfere cell wall synthesis, macrolides and tetracycline inhibit protein synthesis, fluoroquinolones and rifampin interfere nucleic acid synthesis (Drlica & Zhao, 1997;Neu, 1992). Genetic exchange with mutation and selection enables bacteria to adapt the environment quickly in response to antimicrobial agents (de Sousa et al., 2017). Methicillin-resistance Staphylococcus aureus (MRSA) is the most common antibiotic resistance bacteria that affects human. SCCmec is a mobile genetic element containing mecA resistance gene which is responsible for methicillin and oxacillin resistance (Okuma et al., 2002). Similarly, Escherichia coli is a causative agent of urinary tract infections in human and frequently becomes resistant to narrow-spectrum of cephalosporins and ampicillin (Ili c et al., 2011;Yang et al., 2020). b-lactam antibiotics such as penicillins, cephalosporins, carbapenems and methicillin are the most broadly used chemotherapeutic agents due to its low toxicity and high efficacy against bacterial infections (Drawz & Bonomo, 2010). However, the efficacy of these antibiotics is destroyed by b-lactamase enzymes and the genes code for these enzymes significantly contribute in multidrug resistance (Bush, 2010).
Based on the protein sequence homology, there are four classes of b-lactamases viz; A, B, C and D. Among these A, C and D classes are serine hydrolases that recruit serine in the active site and catalyse the b-lactam hydrolysis (Wang et al., 1999). Class B contains zinc ions for their activity and known as metallo-b-lactamases (MBLs). On the basis of sequence conservation, MBLs are further classified into three subclasses viz; B1, B2 and B3. Imipenemase (IMP), New Delhi metallob-lactamases (NDMs) and Verona integrin-encoded metallob-lactamase (VIM) belong to subclass B1 enzyme and are clinically more significant (Lauretti et al., 1999;Ooi et al., 2021;Yong et al., 2009). New Delhi metallo-b-lactamase 1 (NDM-1), firstly reported in Klebsiella pneumoniae from India, consists of single chain protein in which N-terminal domain has a signal peptide and a core of 270 amino acids with enzyme activity (Rehman et al., 2019;Yong et al., 2009). The core of protein harbours active site or binding pocket surrounded by two flexible loops that help in increasing its catalytic activity and also responsible for drug resistance. Various studies have been done regarding the participation of active residues in drug binding and resistance (Liang et al., 2011;Ma et al., 2021). The current study was carried out to understand the molecular mechanism of resistance caused by highly conserved passive residues. Passive residues which are not the part of active site but are essential during drug binding and present in flexible region of the protein that also provide stable conformational regime during drug binding process (Khan & Rehman, 2016;Teague, 2003).
To explore the molecular basis of NDM resistance, in silico substitution mutations of highly conserved passive residues (Aspharagine193, Serine217, Glycine219 and Threonine262) present near to active site and loop region of NDM were created (Ali et al., 2019;Zhang & Hao, 2011). Experimental and computational studies of above mutations was previously performed to elucidate the carbapenem (imipenem and meropenem) resistance towards NDM-1 protein (Ali et al., 2019;Jackson et al., 2020). NDM also hydrolyse other antibiotics such as cephalosporins, methicillin and oxacillin. Here, we have selected methicillin and oxacillin drugs complexed with NDM protein which are the most common and highly resistant antibiotics in human and demonstrated the mechanism of their resistance (King et al., 2012;Okuma et al., 2002). In addition, a carbapenem compound such as doripenem was also included by the fact that NDM also known to hydrolyse the carbapenem (Chahine et al., 2010). Doripenem is a member of b-lactam class of antibiotic which is having bactericidal activity and noninferior to both imipenem and meropenem in treatment of pneumonia and other urinary tract infections (Chahine et al., 2010). Wildtype (WT) NDM protein with methicillin (PDB id: 4EY2) and oxacillin (PDB id: 4EYB) complexes were taken from PDB (Protein Data Bank) and tertiary model of mutant (MT) proteins N193A, S217A, G219A and T262A were constructed using hybrid modelling approach (King et al., 2012). Models were refined by molecular dynamics (MD) simulation technique. A pile of validation tools was applied to monitor the local and global qualities of tertiary structures and ligand binding modes and affinities were obtained through molecular docking. The stability and conformation of protein-ligand complexes were further examined using MD simulation and rescoring of binding energy was elucidated by calculating different energy components using MM-PBSA (Molecular Mechanics-Poisson Boltzmann Surface Area) method. Based on the modelling and MD simulation results, 4 mutants (MTs) N193A, S217A, G219A and T262A showed similar spatial arrangement of a-helices, b-sheets and loops while minor structural differences were observed during simulation. The significant decreased in binding energies were observed in mutant-drug complexes as compared to WT during docking and MD simulation. We observed mutations in the passive residues of NDM protein contribute in drug resistance by modulating the tertiary structure of protein. The current study provides an opportunity for structure-based inhibitor design.

Tertiary structure modelling and molecular dynamics simulation
Experimentally solved structure of WT NDM (PDB id: 4EY2) was obtained from PDB (https://www.rcsb.org/) (King et al., 2012;Rose et al., 2011) and different MT proteins N193A, S217A, G219A and T262A were constructed using hybrid modelling (homology and de novo) approach through I-TASSER server service (https://zhanglab.ccmb.med.umich.edu/ I-TASSER/) (Zhang, 2008). Selection of best models was chosen on the basis of appropriate C (confidence)and TM (template modelling)scores in I-TASSER (Iterative-Threading ASSEmbly Refinement) tool as described previously (Kumar et al., 2018). Constructed MT proteins along with WT were followed for optimization through molecular dynamics (MD) simulation. Molecular dynamics (MD) simulations for all proteins (WT, N193A, S217A, G219A and T262A) were performed by using GROMACS (Groningen Machine for Chemical Simulation) 5.0 suite in conjunction with GROMOS53A6 force field (Lemkul et al., 2010;Van Der Spoel et al., 2005). Briefly, each protein system was solvated with the extended simple point charge (SPC/E) water model in the cubic box and neutralized by adding counter ions which were followed by energy minimization. Each system was equilibrated with a constant temperature of 300 K for 100ps and pressure of 1 bar for 500ps, were maintained by Berendsen and Parinello-Rahman methods, respectively. Lennard-Jones potential and Particle Mesh Ewald (PME) calculations were used to handle the van der Waals and electrostatic interactions, respectively. All bonds were constrained by using LINCS algorithm (Hess et al., 1997). Finally, 100 ns production run was carried in which 2fs time steps were followed and coordinates were saved every10ps.

Structure validation and molecular docking
Structure validation of protein models was assessed through different quality check tools. Initially, the correctness of generated models was validated by superimposing or aligning them with the reference template and structural identity (global and local) were measured by RMSD using SuperPose online tool (http://wishart.biology.ualberta.ca/SuperPose/) (Maiti et al., 2004). Geometry and stereochemistry of models were monitored by the inspection of phi/psi distribution using Ramachandran plot obtained from PROCHECK analysis in SAVES (Structure Analysis and Verification Server) server (http://servicesn.mbi.ucla.edu/SAVES/) (Laskowski et al., 1993). ERRAT and Verify 3D programs of SAVES and ProSA (Protein Structural Analysis) web servers were further used to assess the overall quality of protein models (Luthy et al., 1992;Wiederstein & Sippl, 2007). The degree of nativeness of different models was examined through QMEAN (Qualitative Model Energy ANalysis) server (https://swissmodel.expasy. org/qmean/) (Benkert et al., 2008). Refinement of protein models was also done by ModRefiner tool that refines models in 2 steps (https://zhanglab.ccmb.med.umich.edu/ ModRefiner/) (Xu & Zhang, 2011). Initially it constructs the main chain model from hydrogen bonding network of main chain and tracing the alpha carbon with acceptable backbone topology. Finally, additions of side chain atoms to the backbone conformation were optimized by combining physical and knowledge-based force field. The refinement of global topology in ModRefiner was measured by mean of RMSD (root mean square deviation) and TM-score to native structures. Refined protein models were further used for molecular docking.
Molecular docking was carried by employing 3 different offline and online tools viz; AutoDock, AutoDock Vina and SwissDock (http://www.swissdock.ch/) (Grosdidier et al., 2011;Norgan et al., 2011;Santos-Martins et al., 2014;Trott & Olson, 2010). Chemical structure of methicillin (DB01603), oxacillin (DB00713) and doripenem (DB06211) drugs were obtained from Drugbank (https://www.drugbank.ca/) (Wishart et al., 2018). Open Babel tool was used for different file format conversion of receptor and drug molecules (O'Boyle et al., 2011). Prior performing docking in AutoDock and AutoDock Vina, both receptors and ligands were prepared in Auto Dock tools (Morris et al., 2009). Polar hydrogens were added to receptor molecules and non polar hydrogens were merged whereas ligand molecules were prepared by adding Gasteiger charges. Since NDM is metal protein and has two divalent zinc cations (Zn 2þ ), therefore, two Zn 2þ ions were also included during docking. However, the standard AutoDock force field unable to handle the metal proteins, therefore, a specialized force field such as AutoDock4 Zn that have improved parameters to dock the zinc protein was utilized. The force field (AutoDock4 Zn ) has a potential to describe both geometric and energetic components for the interaction of zinc coordinating ligands which is already implemented in AutoDock tools. Docking parameters involve 100 genetic algorithm (GA) runs, having 150 population size with 25 Â 10 5 evaluations and 27 Â 10 5 generations. The binding site of ligand was obtained from the available protein-ligand complex in PDB which was used as a reference to set the dimensions of grid box in all receptor molecules (King et al., 2012). Size of grid box contains 28, 28, 28 points of dimensions with x, y, z coordinates (À4. 903, 13.284, 21.541) to cover the binding site of protein and kept constant for all WT and MT proteins. Binding poses were generated by Lamarckian genetic algorithm and root mean square deviation (RMSD) were used to cluster the docked conformations (Morris et al., 1998). The best conformation of each receptor-ligand complex was obtained on the basis of binding energy of clusters. Docking procedure used in AutoDock Vina had been previously described (Kumar et al., 2017). Online docking was accomplished by uploading the receptor and ligand molecules on SwissDock server. The binding energies were expressed in kilocalorie per mol and best binding poses for each receptor-ligand complex obtained from all docking tools were selected for downstream processes.

MD Simulations and binding free energy calculation
The high scored docking complexes of WT, N193A, S217A, G219A and T262A MTs with methicillin, oxacillin and doripenem were selected as initial conformations for MD simulation. MD simulations were carried out using GROMACS 5.0 suite. Receptor and ligand topologies were derived from GROMOS53A3 force field in GROMACS suite and PRODRG server, respectively (Schuttelkopf & Van Aalten, 2004). The procedure of MD simulation used for protein-ligand docking complexes was same as used in case of receptor molecules. Binding energies of protein-ligand complexes were calculated using MM-PBSA (Molecular Mechanics -Poisson Boltzmann Surface Area) method as described previously (Kumar et al., 2017). Initially, protein-ligand interaction energy was calculated in GROMACS suite using gmx energy module as given in Equation (1) E Int is the interaction energy calculated by the short range Lennard-Jones E LJ and Coulomb energy E Coul : Later, binding energy was calculated using g_mmpbsa module of GROMCAS (Kumari et al., 2014) as given in Equation (2) where G complex is the total free energy of the protein-ligand complex and G protein and G ligand are the total free energies of proteins and ligands, respectively. All energies were expressed in kilojoule per mol. The binding free energy of protein-ligand complexes were measured from last stabled 10 ns period of MD simulation run. Finally, energy contributions from each residue of the receptor molecule in ligand binding were calculated by energy decomposition method (Kumari et al., 2014).
bond) and interaction energy were calculated using gmx rms, gmx rmsf, gmx gyrate, gmx hbond and gmx energy modules of GROMACS utility. Secondary structures plots were obtained from MD simulation trajectories using DSSP (Dictionary of secondary structure of protein) method (Kabsch & Sander, 1983). All 2D and 3D figures of proteins and ligands were illustrated in LigPlot (Wallace et al., 1995), UCSF Chimera v1.7 (Computer Graphics Laboratory, University of California, San Francisco) and PyMOL (The PyMOL Molecular Graphics System, Version 1.3 Schrodinger, LLC), respectively. 2D graphs were plotted in Prism 6 (GraphPad Software, CA, USA, www.graphpad.com).

Tertiary structure modelling
Knowledge of protein structure is advantageous for structure-based drug screening, provide information regarding disease-related mutations and ultimately determine the function of protein (Gao et al., 2015;Lionta et al., 2014). Experimental methods such as X-ray diffraction and NMR are used to determine the structure of proteins, but these methods are expensive, time-consuming and also limited. While the computational methods are faster, cheaper and provide an accurate way of protein structure prediction (Schames et al., 2004). Here, we harnessed latest computational tools and techniques to construct the protein models. Initially, in silico substitution mutations were accomplished by replacing the targeted residue with alanine and model were constructed through I-TASSER server. Alanine is widely used as substitution residue because it eliminates side chain interactions without altering the main chain conformation and does not introduce electrostatic and steric effects which is preferred for examining the contribution of specific side-chains while preserving native protein structure (Kasturi et al., 1992). I-TASSER predicted 5 models for each MT and final models were selected on the basis of appropriate values of C-and TM-scores. C-score is a confidence score that should be in the range of À5 to 2, with higher scores indicating a model of better quality and TM-score is template modelling score that should be >0.5, for a model indicated correct global topology. The correctness of models was also validated by superimposing the model with reference template from where it was constructed as discussed below. C-and TMscores of N193A and S217A MTs were À0.83, 0.61 and À0.99, 0.59, respectively (Supporting Figure S1(A,B); Table S1). MT proteins G219A and T262A exhibited À1.25, 0.56 and À1.17, 0.57 C-and TM-scores, respectively (Supporting Figure  S1(C,D); Table S1). Above results showed that all generated models had C-score in the range between À5 to 2 and TMscore greater than 0.5, indicated that all models were accurate with correct global topologies. Generated MT models showed similar spatial arrangement of secondary structures (a-helices, b-sheets and loop and turns) and conformations to WT protein.

Tertiary structure optimization by MD simulations
Molecular dynamics (MD) simulation has been commonly applied to refine and check the stability of protein models (Raval et al., 2012). All generated models along with WT protein were optimized through MD simulation in a solvent to imitate the real physiological environment. Each system was subjected to 100 ns MD simulations and data were collected in the form of trajectories. Prior to actual analysis, temperature of 300 K, pressure of 1 bar, total energy (potential and kinetic) of all the systems and minimum distance between periodic images (>2 nm) were examined. In MD simulation, the stability of all proteins was measured by its deviation from the initial structures in terms of root mean square deviation (RMSD) and observed that a stabled conformation was maintained during entire simulation time of 100 ns ( Figure  1 Mobility and flexibility of residues and structures were monitored by the assessment of root mean square fluctuation (RMSF) (Supporting Figure S2). Fluctuation of Phy70 (phenylalanine) with RMSF value of $0.59 nm was observed in N193A (Supporting Figure S2(A)). Inspection of N193A tertiary structure suggested that the above residue located at the boundary of the sheet that occupied turn region (Supporting Figure S2(C)). In S217A mutant, fluctuation of Gly153 (glycine), Met154 (methionine), Ala156 (alanine) and Gln157 (glutamine) with RMSF values of $0.41, $0.39, $0.31 and $0.34 nm were observed (Supporting Figure S2(A)). The analysis of S217A tertiary structure suggested that the corresponding residues located at loop or random regions of protein (Supporting Figure S2(D)). While all other MT proteins (G219A and T262A) exhibited no significant fluctuations compared with WT (Supporting Figure S2(A,E,F)). In conclusion, the RMSF plot confirmed the fluctuations were more pronounced to loop regions and turns that connected secondary structures. Results of MD simulation indicated that all MT protein models along with WT showed stable and consistent behaviours of RMSD, Rg and RMSF with no significant conformational changes existed during entire simulation period.

Estimation and evaluation of stereochemical properties
Validation and evaluation of geometrical and stereochemical properties of tertiary models are necessary prior to prepare the models for further experiments. Initially, validation of generated models was performed by superimposing the generated models with the reference templates and both global and local structural similarities were measured using SuperPose tool. We found that all models exhibited similar and low values of global and local RMSD thus models showed less deviation and higher similarities to its templates (Table S1). The stereochemical properties of WT and MT proteins were deduced from SAVES server. In WT protein, 92.6% and 0.5% residues were placed in the most favoured and disallowed regions of the Ramachandran plot, respectively (Table S1). While in N193A, S217A, G219A and T262A MT proteins showed 94.1%, 94.1%, 94.6%, 94.1%, and 0.5%, 1.0%, 0.5%, 0.5% residues were placed in the most favoured and disallowed regions of Ramachandran plots, respectively (Table S1). A larger number of protein residues in most favoured regions suggested that all MT and WT protein structures were stereochemically stabled. The overall quality of WT and MT proteins were further monitored by ERRAT program and observed that all proteins exhibited fair qualities located with the range of $90.5 to 94.4% (Table S1). Total quality index of protein models was monitored by measurement of Z-score in ProSA web server. The results demonstrated that the total quality index (Z-score) was found in the range of À8.0 to À8.6 for WT as well as MT models (Table S1). Low values of QMEAN score (À1.27 to À5.87) indicating that all models displayed native-like structures (Table  S1). Global and local qualities of WT and MTs were also examined through calculating the RMSD and TM-score of ModRefiner tool and found that all protein models exhibited low RMSD and better TM-scores implying protein models had good local and global topologies (Table S1). Verify3D results displayed about 99.5 to 100% of the residues had 3D-1D compatibility for all protein models. Above results indicated that all protein (WT and MTs) models had good structure qualities with optimized geometry and also all models have stereochemically stabled structures.

Binding modes generated by molecular docking
Binding modes of protein and ligands (methicillin, oxacillin and doripenem) were studied by measuring the binding energy obtained from 3 docking tools (AutoDock, AutoDock Vina and SwissDock) and inspection of physical and chemical bonding between protein and drug molecules. 3D structures for all WT and MTs obtained from last stabled MD simulation     (Ali et al., 2019). Furthermore, docking scores of methicillin and oxacillin were more or less similar to meropenem but higher than imipenem and doripenem. Binding modes of protein and ligands were also monitored through assessing the different types of physical and chemical interactions by plotting them on 2D plot (Supporting Figures S3-S5). Hydrophobic and hydrophilic interactions were monitored during protein-ligand complexes (Table S4). WT-methicillin complex exhibited 7 hydrogen bonds with His122, Gln123, Asp124, His189, Cys208, Asn220 and Asp223 and 6 hydrophobic bonds with Leu65, Met67, Val73, Trp93, Gly219 and His250 (Supporting Figure S3 and Table S4). MT proteins such as N193A-methicillin and S217Amethicillin exhibited 6, 8 hydrogen bonds and 6 hydrophobic interactions, respectively (Supporting Figure S3 and Table  S4). G219A-methicillin and T262A-methicillin complexes exhibited 5, 7 hydrogen bonds and 6 hydrophobic interactions (Supporting Figure S3 and Table S4). On the other hand, the WT-oxacillin complex exhibited 4 hydrogen bonds with Gln123, Asp124, Lys211 and Asn220 and 8 hydrophobic interactions (Supporting Figure S4 and Table S4). N193A and S217A-oxacillin complexes formed only 3, 4 hydrogen bonds and 7, 8 hydrophobic interactions, respectively, whereas G219A-and T262A-oxacillin complexes exhibited 2, 4 hydrogen bonds and 11, 8 hydrophobic interactions (Supporting Figure S4 and Table S4). Binding strength of WT and MTs with doripenem were also elucidated by 2 D plots and found that WT-doripenem exhibited 3 and 13 hydrogen bonds and hydrophobic interactions, respectively (Supporting Figure S5 and Table S4). N193A and S217A-doripenem complexes formed 5 hydrogen bonds and 5, 10 hydrophobic interactions, respectively, whereas G219A-and T262A-doripenem complexes exhibited 7, 3 hydrogen bonds and 8, 12 hydrophobic interactions (Supporting Figure S5 and Table S4). Above results indicated that methicillin exhibited stronger binding affinities with WT as compared to all MT proteins while oxacillin exhibited stronger binding affinity with WT and T262A MT as compared to N193, S217A and G219A MT proteins. Also, doripenem strongly bind with WT, G219A and T262A as compared to N193A and S217A MTs. Unlike, other carbapenem such as imipenem and meropenem which displayed weaker binding with N193A, S217A and G219A MTs compared with WT and T262A MT as shown previously (Ali et al., 2019). Overall docking results suggested that oxacillin, imipenem and meropenem showed similar binding affinities as compared to methicillin and doripenem.

Binding modes refined by MD simulation
Prediction of binding pose by molecular docking are generally reliable, but it also fails to estimate the drug binding affinity because it treats receptor as a rigid molecule in order to save computational time that limits the conformational movement of receptor during docking (Sousa et al., 2006). The assumption of rigidity works well in case of rigid protein or proteins which are less flexible, but NDM protein is more flexible and the flexibility of this protein is assisted in ligand binding as also reported previously (Zhang & Hao, 2011). MD simulation is extensively applied to monitor the stability of protein-drug/protein/DNA complexes, to rescore the binding energy of ligands, to study the conformational changes occur during ligand binding and unbinding and to study the protein folding and protein dynamics (Kumar & Saran, 2018;Lindorff-Larsen et al., 2011;Piana et al., 2013). In account of NDM protein flexibility, MD simulation is necessary to monitor the conformational changes and stability existed in protein during ligand binding.
Here, MD simulation was executed for all protein-ligand complexes to refine and rescore the binding modes which were obtained through molecular docking. In addition, the effect of drugs on protein secondary structures formation during simulation were also examined. RMSDs of WT and MTs complexed with methicillin were measured and found that no significant drift was originated from WT and mutant structures and steady and similar behaviours were found (Figure 2(A)). WT and N193A, G219A, and S217A MTs were aligned with average RMSD values of $0.21, $0.22, $0.22 and $0.20 nm, respectively, except T262A MT which was having $0.25 nm of an average RMSD (Figure 2(G)). RMSDs of methicillin complexed with WT and MT proteins showed stable behaviours except WT and S217A MT which showed deviations throughout the simulation period (Figure 2(B)). RMSDs of methicillin complexed with WT and S217A MT were stabilized after $92 and $85 ns, respectively ( Figure  2(B)). Average RMSD values of WT, N193A, S217A, G219A and T262A were $0.19, $0.27, $0.17, $0.21 and $0.22 nm, respectively (Figure 2(G)). On the other hand, RMSDs of WT and MT proteins complexed with oxacillin were monitored   Higher RMSD values in MT-methicillin and MT-doripenem complexes compared to WT-methicillin and -doripenem were observed suggesting that all MTs led to structural destabilization. Similar results to doripenem were also observed in case of imipenem and meropenem where higher RMSDs were observed in MTs as compared to WT (Ali et al., 2019). In contrast, RMSDs of methicillin and oxacillin complexed with N193A, S217A and T262A MTs were higher as compared to rest of MTs and WT. Above results indicating that MTs greatly affected carbapenem and other antibiotics such as methicillin and oxacillin-NDM complexes stabilities.
Shape and compactness of protein-drug complexes were monitored by measuring the Rg (Figure 3). Throughout the simulation period, WT and MT proteins complexed with methicillin ( Figure 3(A)), oxacillin (Figure 3(B)) and doripenem (Figure 3(C)) displayed similar pattern of Rg with almost constant values in the range of $1.67 (minimum) to $1.72 nm (maximum) (Figure 3(D)). However, no such changes in Rg of imipenem and meropenem drugs were existed, suggesting that these drugs had no effect on overall compactness and globularity of WT and MT proteins as previously observed (Ali et al., 2019). Contrary to that, methicillin, oxacillin and doripenem drugs complexed with MTs had altered the compactness of protein, suggesting that MTs led to destabilization of NDM-methicillin, -oxacillin and -doripenem complexes.
The mobility of different residues of a protein in presence of drug was examined by measuring the RMSF values (Supporting Figures S6-S8). RMSF results indicated that higher fluctuations were observed in apo protein as compared to complex form. Fluctuations in some residues were observed which were located at loop region in WT-methicillin complex with $0.42 nm RMSF value (Supporting Figure  S6(A,B)), whereas no fluctuations were found in N193Amethicillin (Supporting Figure S6(C,D)). S217A-methicillin and G219A-methicillin exhibited $0.34 and $0.32 nm RMSF values, respectively (Supporting Figure S6(E,F,G,H)), whereas no fluctuations were found in T262A-methicillin complex (Supporting Figure S6(I,J)). Results of RMSF demonstrated that WT, S217A and G219A MT proteins underwent some conformational changes in order to accommodate the drug for making a stable complex with methicillin. Similar results of RMSF were also observed in case of protein-oxacillin complex (Supporting Figure S7). Leu221 in S217A-(Supporting Figure S7(E,F)) and T262A-oxacillin complexes (Supporting Figure S7(I,J)) showed high fluctuations with RMSF values $0.57 and $0.59 nm indicating the structural alterations in those MT proteins, while rest of protein-drug complexes showed minor or negligible fluctuations (Supporting Figure S7(A-D,G,H)). RMSF analyses of WT and MT-doripenem complexes showed no major fluctuations and it is restricted only the flanking or extreme end of the proteins (Supporting Figure S8) except T262A MT where RMSF of >0.5 nm were observed particularly, in Gly219, Asn220 and Asn223 (Supporting Figure S8(I,J)). Hence, these results were quite similar to results of imipenem and meropenem as previously noted but the structure levels fluctuations of WT and MTs have not been investigated in presence of these drugs (Ali et al., 2019). Results of proteindrug simulations demonstrated that all WT and MT complexes exhibited stable and similar dynamics behaviours with some fluctuations were notable in S217A-and T262Aoxacillin and -doripenem complexes.

Secondary structures and hydrogen bond analyses during MD simulation
Structural stability of WT and MT proteins were carried out to trace the alterations in secondary structures formation (Coils, b-sheets, Bend, Turn and a-helices) in apo and complex form of proteins with various drugs during the course of simulation period. The secondary structures of proteins were deduced from MD simulation trajectories and analysed by do_dssp module which is an external tool, but GROMACS provide an interface to do_dssp utility. Variations in bend and turn regions of proteins were observed in presence of drug which led to changes in the percentage of helices and sheets ( Figure 4 and Table S5). The total coil contents of WT and N193A, S217A, G219A and T262A MTs in complexed with methicillin were 22.5%, 23%, 23%, 23.5% and 24%, whereas coil contents of WT and MTs complexed with oxacillin and doripenem were 23%, 23.5%, 23.5%, 23.5%, 23.5% and 23%, 23.5%, 23.5%, 23.5%, 24%, respectively ( Figure  4(A)). Secondary structure particularly, sheet contents of WT and N193A, S217A, G219A and T262A MTs complexed with methicillin and oxacillin were 26.5%, 26.5%, 26.5%, 26%, 26.5% and 25.5%, 26%, 25.5%, 25.5% and 26%, respectively (Figure 4(B)), whereas sheet contents of WT and MTs complexed with doripenem were 25.5%, 26.5%, 25.5%, 25% and 25.5%, respectively (Figure 4(B)). Bend contents of WT and N193A, S217A, G219A and T262A MTs complexed with methicillin and oxacillin were 15%, 13.5%, 15%, 14%, 15% and 15%, 13.5%, 15%, 14.5% and 14.5% respectively ( Figure  4(C)), whereas bend contents of WT and MTs complexed with doripenem were 14.5%, 13.5%, 15.5%, 14.5% and 14%, respectively (Figure 4(C)). Turn contents of WT and N193A, S217A, G219A and T262A MTs complexed with methicillin and oxacillin were 12.5%, 12%, 12.5%, 12%, 11.5% and 11.5%, 12.5%, 12%, 12.5% and 12.5%, respectively ( Figure  4(D)), whereas bend contents of WT and MTs complexed with doripenem were 12.5%, 12%, 12%, 11% and 12%, respectively (Figure 4(D)). Secondary structure such as helix contents of WT and N193A, S217A, G219A and T262A MTs in complexed with methicillin were 22.5%, 23.5%, 21%, 23.5% and 22.5%, whereas helix contents of WT and MTs complexed with oxacillin and doripenem were 22.5%, 22%, 21.5%, 21.5%, 20% and 21.5%, 22%, 21%, 22%, 22.5%, respectively (Figure 4(E)). Above results demonstrated that the coils remained almost similar in all MT and WT proteins while rest of the secondary structure entities such as b-sheets, bend, turn and a-helices were altered in case of MT proteins and marginal alterations in WT protein. However, bend and turns appeared to slightly increased in N193A, S217A and T262A MT proteins when complexed with drugs as compared to WT, resulting more flexible nature of proteins ( Figure 4 and Table S5). In contrast to other carbapenems such as imipenem and meropenem, increased content of sheets and helices were observed in WT and MTs indicating proteins acquired rigid conformations in presence of these drugs (Ali et al., 2019). This has been not observed in case of protein-methicillin and -oxacillin complexes as protein acquired flexibility which led to unstable protein secondary structures. The results of secondary structure analysis were in agreement with results obtained through RMSF investigations and hence we conclude that drug molecules modulated the structure of MT proteins by shifting the formation of secondary structure components.
Further stabilities of WT-and MT-ligand complexes were assessed by examining the protein-ligand interaction and hydrogen bond (H-bond) formation during simulation. Protein-ligand binding interactions and hydrogen bonding plots of all methicillin, oxacillin and doripenem were shown in Figures 5-7. Protein-ligand binding complexes were extracted from last trajectory of MD simulation and H-bond plots were calculated by utilizing entire MD simulation trajectories based on the geometric criteria that H-bond formation should have <0.30 nm distance between H-bond acceptor and donor and angle should larger than 120˚. The study revealed that methicillin drug was well accommodated in the binding pocket of all MT and WT proteins, suggesting that drug remained stable in the binding surface of protein during entire course of simulation ( Figure 5(A,C,E,G,I)). An average of 4-5 H-bonds was observed in WT-methicillin complex ( Figure 5(B)) whereas 3-4 H-bonds were found in S217Amethicillin and G219A-methicillin complexes ( Figure 5(F,H)). An average of 3-4 H-bonds was observed in N193A-methicillin complex during the end of simulation period ( Figure 5(D)) while 5 H-bonds were formed during entire period of simulation in T262A-methicillin complex ( Figure 5(J)). H-bonding capacity of methicillin remained similar in both WT and MTs.
Similar results were also observed in case of protein-oxacillin complexes as the perfect orientation of drug was found at binding surface of all MT and WT proteins suggested that ligand remained stable during entire simulation time ( Figure  6(A,C,E,G,I)). An average of 4-5 H-bonds was observed in WTand G219A-oxacillin complexes (Figure 6(B,H)) while 3-4 Hbonds were found in N193A-oxacillin, S217A-oxacillin ( Figure  6(D,F)) and T262A-oxacillin complexes (Figure 6(J)).
On the other hand, doripenem was well oriented in the binding cleft of WT and all MT proteins again indicating the drug remained stable during simulation (Figure 7(A,C,E,G,I)). An average of 4-6 H-bonds was observed in WT-, S217A and T262A-doripenem complexes (Figure 7(B,F,J)) while 3-4 Hbonds were found in N193A-doripenem and G219A-doripenem complexes (Figure 7(D,F)). H-bonding results demonstrated that methicillin constituted lesser hydrogen bonds with MT proteins as compared to WT, while oxacillin showed more hydrogen bonds with S217A and T262A MT proteins compared with WT, G219A and N193A MT proteins. Doripenem formed maximum H-bonds with WT, S217A and T262A while it formed minimum H-bonds with N193A and G219A MTs. Unlike other carbapenems such as imipenem and meropenem, slightly lesser number of H-bonds were found in MTs as compared to doripenem, indicating H-bonding between protein and drugs were reduced in both MTimipenem and -meropenem complexes (Ali et al., 2019). Hbonding were found to be highly altered in MT-methicillin, -oxacillin and -doripenem complexes.

Binding free energy calculation
Rescoring of binding energies were accomplished by using MD simulation trajectories of last 10 ns simulation time through g_mmpbsa utility where all free energy binding components (van der Waal's interaction, electrostatic interaction, polar solvation energy, SASA energy, and binding energy) were computed (Tables 1-3).
Above data suggested that oxacillin has less binding capacity with S217A as compared to WT and other MT proteins whereas binding affinities of methicillin towards MT were higher than WT. Moreover, binding affinities of doripenem with all MT proteins were less than the WT. On the other hand, binding energies of others carbapenems like imipenem and meropenem were higher in N193A, G219A, S217A and T262A MTs as compared to other MTs and WT, respectively, indicating that different carbapenems expressed differential binding affinities to NDM protein (Ali et al., 2019). In other words, oxacillin showed resistance to S217A MT and doripenem showed resistance to all MT proteins, respectively. Hence, we conclude that mutant proteins displayed resistance to available drugs, particularly oxacillin and doripenem by modulating the spatial arrangement of secondary structures and conformation of the tertiary structure of proteins.

Conclusions
The objective of the current study was to elucidate the role of passive residues in drug resistance towards NDM protein.
To fulfil above objective, we manifested different substitution mutations such as N193A, S217A, G219A and T262A, and 3D models were constructed through a hybrid modelling approach. Stability of models was examined by MD simulations and stereochemical properties were estimated via various structure validation tools. We observed that all MTs showed low and consistence behaviours of RMSD, Rg and RMSF except minor structural alterations in loop and turn regions. Assessment of stereochemical properties of all models showed stable local and global protein geometry. The contribution of protein residues in methicillin, oxacillin and doripenem binding were monitored through molecular docking that displayed low binding affinities towards MT as compared to WT. Binding affinity and mechanism were further confirmed by MD simulations suggested that oxacillin had less binding capacity to S217A and doripenem showed low binding affinities to all MTs as compared to WT. Moreover, all drugs displayed lesser hydrogen bonds with almost all MTs as compared to WT. Additionally, drug molecules also enhance the flexibility of secondary structure of MT proteins that provide the conformational instability of protein during ligand binding. Taken all together, this study sheds light on the molecular mechanism of resistance exerted by oxacillin and doripenem on NDM protein. Hence, current study highlights the importance of passive residues that facilitate the significant insights into structure-based development of potent inhibitor against NDM protein.

Disclosure statement
No potential conflict of interest is reported by the author.

Funding
Author thanks Indian Council of Medical Research (ICMR-India) for financial support.