Identification of natural inhibitors against Acinetobacter baumannii d -alanine-d -alanine ligase enzyme: A multi-spectrum in silico approach

Accepted Manuscript Identification of natural inhibitors against Acinetobacter baumannii d-alanine-d-alanine ligase enzyme: A multi-spectrum in silico approach

Sajjad Ahmad, Saad Raza, Sumra Wajid Abbasi, Syed Sikander Azam PII: DOI: Reference:

S0167-7322(17)36043-9 doi:10.1016/j.molliq.2018.04.124 MOLLIQ 9027

To appear in:

Journal of Molecular Liquids

Received date: Revised date: Accepted date:

19 December 2017 22 March 2018 24 April 2018

Please cite this article as: Sajjad Ahmad, Saad Raza, Sumra Wajid Abbasi, Syed Sikander Azam , Identification of natural inhibitors against Acinetobacter baumannii d-alanined-alanine ligase enzyme: A multi-spectrum in silico approach. The address for the corresponding author was captured as affiliation for all authors. Please check if appropriate. Molliq(2017), doi:10.1016/j.molliq.2018.04.124

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Title page Manuscript Title: Identification of natural inhibitors against Acinetobacter baumannii D-alanineD-alanine ligase enzyme: A multi-spectrum in silico approach Authors: Sajjad Ahmad, Saad Raza, Sumra Wajid Abbasi, Syed Sikander Azam* Computational Biology Lab, National Center for Bioinformatics, Quaid-i-Azam University,

T

Islamabad-45320, Pakistan.

IP

*Corresponding Author:

E-mail: [email protected]; [email protected] Phone: 0092-51-906 44130

CR

Syed Sikander Azam

US

Computational Biology Lab, National Center for Bioinformatics, Quaid-i-Azam University, Islamabad-45320, Pakistan.

AN

Emails:

Saad Raza ([email protected])

M

Sajjad Ahmad ([email protected])

Sumra Wajid Abbasi ([email protected])

AC

CE

PT

ED

Syed Sikander Azam ([email protected]; [email protected])

1

ACCEPTED MANUSCRIPT Abstract

CE

PT

ED

M

AN

US

CR

IP

T

D-alanine-D-alanine ligase (Ddl), an enzyme that catalyzes the D-ala-D-ala dipeptide formation in UDPMurNAc pentapeptide, is a part of peptidoglycan biosynthesis machinery. Inhibition of enzyme leads to bacterial growth arrest making it a viable and attractive target for screening of potent antibacterial drugs. Combination of virtual screening, molecular docking, in silico pharmacokinetics, molecular dynamics (MD) simulation, and binding free energy calculations based on Molecular Mechanics Generalized Born and Surface Area (MMGBSA) and WaterSwap were applied in the current framework for the detailed analysis of potent natural inhibitors against Ddl enzyme. Comparative molecular docking supported with computational druglikeness revealed compound-331 (6-(4-((3-methoxyphenylsulfonamido) methyl) phenyl)-2-methylnicotinamide) as the best docked inhibitor. The inhibitor has Genetic Optimization for Ligand Docking (GOLD) fitness score of 84.2 and AutoDock Vina binding energy of -7.2 kcal/mol. The inhibitor exhibited to have a better druglikeness by adhering to Lipinski rule of five, Ghose rule, Veber filter, Egan filter, and Muegge filter in contrast to D-cycloserine, which violates Ghose and Muegge rule. MD simulation unravels that in complex over the course of 100-ns, the enzyme remained highly stable with mean Root Mean Square Deviation (RMSD) of 1.4 Å when compared to an undocked structure having RMSD of 1.6 Å. Root Mean Square Fluctuation (RMSF) predicted stable behaviour of the active site residues in both undocked and docked system with average RMSD values of 1.2 Å and 0.7 Å, respectively. Radial Distribution Function (RDF) and Axial Frequency Distribution (AFD) demonstrated Lys176 and Trp177 as critical residues of enzyme for binding, anchoring and bridging strong hydrogen and hydrophobic contacts between enzyme and the inhibitor. The estimated MMGBSA based binding free energy for the complex is -47.36 kcal/mol, signifying its stable nature. WaterSwap further indicated a worthy agreement on the affinity of inhibitor towards enzyme active pocket. This newly discovered natural inhibitor might serve a parent structure for the development of more potent derivatives with encouraging biological activity.

AC

Keywords: D-alanine-D-alanine ligase; Molecular docking; Molecular dynamics simulation; Radial distribution function; Axial frequency distribution; MMGBSA; WaterSwap.

2

ACCEPTED MANUSCRIPT 1. Introduction

CR

IP

T

Acinetobacter baumannii is a Gram-negative bacilli that has emerged as a multidrug-resistant (MDR) pathogen, tough to control and treat [1]. The organism produced several types of hospitalacquired infections including hospital-acquired pneumonia, bacteremia, endocarditis, meningitis, soft and skin tissue infections, and urinary tract infections [2]. According to Center for Disease Control and Prevention (CDC), an estimated 12,000 A. baumannii associated infections occurred in the United States each year. Approximately, 7,000 (63%) of these are MDR, which leads to 500 deaths [3]. The clinical isolates of A. baumannii have been shown to confer resistance to majority of the currently used antibiotics, including aminoglycosides, β-lactams, diaminopyrimidines, macrolides, phenicols, quaternary amines, sulfonamides, and streptothricins [4]. Moreover, A. baumannii is recently placed in “Critical” category based on the urgency of need for the development of novel antibiotics [5].

CE

PT

ED

M

AN

US

The pathways that govern cell wall biosynthesis are generally a target of great interest against bacterial pathogens for designing of novel antibiotics since they are essential and specific for bacterial survival [6, 7]. One of the integral structures of the cell wall is peptidoglycan that preserves cell integrity and ensures defined cell shape [8]. Peptidoglycan is a common feature of both Gram-positive and Gram-negative bacteria and is an attractive target for antibacterial [7]. Peptidoglycan biosynthesis is broadly divided into three sequential steps; (i) synthesis of nucleotide precursors in the cytoplasm, (ii) synthesis of lipid-linked intermediated on the inner side of cytoplasm membrane, and (iii) polymerization reactions that occurred at outside of cytoplasmic membrane [8]. Each of these steps encompasses a series of enzymatically catalyzed reactions, which occurred at different locations [7, 9]. D-alanine-D-alanine ligase (Ddl) enzyme acts in the first stage and utilized D-alanine as a substrate [10, 11]. Ddl functions by ligating two D-alanine units to produce dipeptide D-alanyl-D-alanine (D-ala-D-ala) as shown in S-Fig.1. Dala-D-ala is a fundamental component of bacterial cell that is essential for maintaining cell wall stability by cross-linking peptide chains of peptidoglycan [12].

AC

So far, four main kinds of Ddl inhibitors are reported: (i) D-Ala analogues, (ii) D-ala-D-ala analogues, (iii) transition state analogues, and (iv) more recently those discovered by screening and modeling methods [9, 12]. Among these, D-cycloserine (D-4-amino-3-isoxazolidine) is the first discovered and most important [13]. D-cycloserine is an analogue of D-alanine having a Ki value of 27 µM and is the only clinically used inhibitor mostly employed in a combination therapy for treatment of Mycobacterium tuberculosis associated infections [13]. The use of D-cycloserine is now almost completely abounded due to its high Minimum Inhibitory Concentration (MIC) value and neurological side effects [14]. Over the last decade, virtual screening [15-19] and de novo structural-based drug design [20-21] have reported several new scaffolds of Ddl inhibitors. Most of these inhibitors show no similarity with the substrate, reaction intermediates, and product [9]. Due to lack of potent Ddl inhibitors, the search for novel inhibitor scaffolds against this target enzyme is urgently required for the routine use in antibacterial therapy. In this current investigation, we documented the potential of natural inhibitors against Ddl enzyme based on the 3

ACCEPTED MANUSCRIPT

CR

IP

T

applications of in silico drug discovery. The library of natural inhibitors was subjected first to Lipinski’s rule of five to screen only drug-like compounds. Structure based virtual screening was then followed to identify compounds with the best orientation in enzyme active site and docking score. The functionality and physical behavior of the enzyme were further deciphered to understand its dynamics in docked and undocked form. This was especially important for the behavior of active pocket of the enzyme in the presence and absence of inhibitor. Further, in complex the enzyme was analyzed for the first time through Radial Distribution Function (RDF) [22] and Axial Frequency Distribution (AFD) [23-24] to examine the hydrogen-bonding patterns of vital residues from enzyme active site. Lastly, binding free energy calculations were performed for enzyme complex to understand the inhibitor affinity for active site. Molecular Mechanics Generalized Born and Surface Area (MMGBSA) [25] and a more recent and sophisticated approach of WaterSwap [26-27] were applied.

US

2. Materials and methods 2.1. Protein and inhibitors preparation

CE

PT

ED

M

AN

Since the crystal structure of A. baumannii Ddl is resolved at a resolution of 2.2 Å and is available in protein data bank (PDB id: 5D8D) [28], it was retrieved and subjected to initial preparation phase. The enzyme is present in homodimer stereochemistry; each monomer with the independent identical active site, so only one monomer was retained [29]. The protein structure was relaxed by assigning Gasteiger charges and steric clashes were removed by minimizing the protein through UCSF Chimera for 1500 rounds of minimization (750 steepest descent steps and 750 conjugate gradient steps). The steepest descent and conjugate gradient step size were set to 0.02 Å under tripose force field (TFF) [24]. For inhibitors, Asinex antibacterial library (http://www.asinex.com/antibacterial_compound_library-html/) of natural compounds was used. Virtual screening of the library based on Lipinski rule of five filters [30] was first accomplished through Ligand scout 4.1[31]. The screened compounds were minimized by the MMFF94 force field [32]. 2.2. Binding cavity prediction

AC

Binding cavity prediction for docking of inhibitors was done through multiple approaches. First, the catalytic cavity was predicted through online tools (Metapocket [33], DoGSiteScorer [34]). The predicted cavities were reconfirmed thorough literature search [28]. After confirmation, coordinates of a highly conserved catalytic residue, as point of inhibitors docking, was identified via multiple sequence alignment (MSA) approach [35]. 2.3. Molecular docking Molecular docking was performed to unveil the preferred binding mode of drug-like inhibitors into protein active cavity and to understand their affinity toward the active pocket. Docking was carried out using Genetic Optimization for Ligand Docking (GOLD) [36] and AutoDock Vina [37]. In both cases, the receptor protein was kept rigid while the binding site for compounds was defined 4

ACCEPTED MANUSCRIPT to include all the residues for binding present within 10 Å of Glu247. GOLD docking was accomplished through the automated GOLD wizard on Intel Xeon QuadTM Core processor of 3.0 GHz with Linux workstation. Default parameters were employed for genetic algorithm while number of iterations was set to 10. GOLD fitness score was assigned to each compound, which is based on the following equation. GOLD fitness score = Shb-ext + SvdW-ext +Shb-int + Sint

IP

T

In the above equation, Shb-ext and SvdW-ext are hydrogen and van der Waals interaction score, respectively for the protein-inhibitor complex. Shb-int is the intramolecular hydrogen bonding contributions to overall GOLD fitness score while Sint is an intramolecular strain of the inhibitor.

US

CR

AutoDock Vina calculations were done on Intel Core (TM) i5 CPU M 540 @ 2.53 GHz processor. The binding site coordinates were the same used in GOLD. The inhibitors binding affinity was calculated in terms of binding energy in kcal/mol. The selection of a compound for molecular dynamics (MD) simulation was based on GOLD fitness score, AutoDock Vina binding energy, and druglikenss predicted by SwisADME [38] and preADMET [39].

AN

2.4. Molecular dynamics simulation

AC

CE

PT

ED

M

To further gain dynamic and mechanistic insights of enzyme and its docked complex, a 100-ns MD simulation was carried out [23]. Simulation protocol was accomplished using Assisted Model Building with Energy Refinement version 14 (AMBER14) [40]. Initial libraries and parameters for the inhibitor were generated through antechamber program of the AMBER. The integration of complex in TIP3P water box (size of 12 Å) was assisted with Leap program using ff14SB force field to explain molecular interactions of the system [41]. The energy of systems was minimized in a gradual manner [42] started with relaxing the entire system hydrogens for 500 rounds of minimization followed by water box minimization for 1000 cycles with a restraint of 200 kcal/mol. The systems were then subjected to restraint of 5 kcal/mol-Å2 on alpha carbon atoms for 1000 steps. In the last phase of minimization, all non-heavy atoms were minimized for 300 cycles with restraint of 100 kcal/mol- Å2. Followed by minimization, each system was heated gradually to 300 K for a time scale of 20-ps with restraint on carbon alpha atoms of 5 kcal/mol- Å2 and a time step set to 2-fs. To maintain temperature of the systems, Langevin dynamics [43] was used. Running Langevin dynamics at constant pressure (NPT) is quite useful to maintain a constant density and let the box reduce to the "proper" size for the number of solvent molecules in system. The gamma value was set to 1.0, while NVT ensemble and SHAKE algorithm [44] were used for heating and applying constraints on hydrogen bonds of the system. Systems equilibration was achieved for 100-ps and time step of 2-ns [45]. Systems pressure was maintained using NPT ensemble along with isotropic position scaling and restrained on carbon atoms of 5 kcal/mol- Å2. The pressure phase was extended for another 50-ps with restraint on carbon atom, which is reduced by 1 kcal/mol- Å2 applied after every 10-ps. The systems were then allowed to equilibrate for 1 ns using the same conditions discussed above. Finally, the systems were subjected to a production run of 100-ns using Berendsen temperature coupling with NVT ensemble. As equilibration in density is 5

ACCEPTED MANUSCRIPT

T

attained, we used Berendsen to compute potential energies from a trajectory file. We also used NVT to tighten the shake tolerance, the particle-mesh-Ewald method (PME) tolerances and reduce the time step to get ideal energy conservation. This was important to compute potential energies from a trajectory file and tighten the shake tolerance and reduce the time step to get ideal energy conservation. The production run was performed with time step value of 2-fs and non-bounded interaction cut off 8.0 Å. Simulated trajectories were analyzed using CCPTRAJ [46] of AMBER. The snapshots taken from trajectories were visualized using UCSF Chimera [47] and Visual Molecular Dynamics (VMD) [48].

IP

2.5. RDF and AFD

PT

ED

M

AN

US

CR

Hydrogen bond interactions in addition to intermolecular close contacts between enzyme active site residues and that of inhibitor crucial for complex stabilization were scrutinized through a Perl written script used in VMD. This was vital in highlighting enzyme active residues that played significant role in recognizing and binding inhibitor over the period of simulation. The screened interactions were then used in RDF [22] and AFD [24]. RDF was applied for finding the probability and to quantify ligand atoms present in the vicinity of protein pocket. RDF assists in getting valuable insights about the conformation variations induced as a result of intermolecular interactions between the inhibitor and protein active residues. This was achieved by employing the PTRAJ module of AMBER14. To further express the coordination geometry of the ligand atoms with respect to protein in a three dimensional (3D) histogram, an indigenously engineered analytical tool termed as AFD was used. AFD is a novel analytical tool designed by the Computational Biology Group at National Center for Bioinformatics, Quaid-i-Azam University, Islamabad, Pakistan. AFD is very sensitive to local structural rearrangements and can recognize local displacements, providing in-depth informations about system stability. The output of the analysis is in the form of a 3D histogram, in which ligand atom coordinates are plotted on XY plane while protein atom was set as a reference [24].

CE

2.6. Binding free energy calculations

AC

In rational drug designing, structure-based virtual screening of drug libraries, combined with MD simulation and binding free energy computation is a widely recommended approach for identification of potent inhibitors against the intended biological target [49]. MMGBSA [25] is a popular technique for predicting estimated binding affinity of inhibitors at the expense of modest computational efforts. The net values of binding free energies are extracted by considering the difference between complex binding energy and that combined for receptor and ligand. This can be represented as below, ∆Gbinding energy = ∆Gcomplex– [∆Greceptor + ∆Gligand] Although quite successful in improving the findings of molecular docking, yet MMGBSA contains several questionable and crude approximations. Among these, the most disregard parameter is the use of implicit water model, which ignores the free energy of water molecules in the protein active 6

ACCEPTED MANUSCRIPT

T

site. This could result remarkably on the approximation of binding free energy. To circumvent the limitations of MMGBSA, Christopher et al [26-27] introduced a more sophisticated method of “WaterSwap” based on the explicit solvent model system. This method is based on the principle of swapping the ligand with an equal shape and volume of explicit water molecules present in the protein active cavity. A couple of simulation boxes are used in the method. The first box contains protein, ligand, and explicit water molecules while the second box is composed of only explicit water molecules. Both the boxes are placed in a heat bath to ensure the transfer of heat and entropy used for decoupling between the boxes.

CR

IP

WaterSwap was performed for 1000 iterations with sampling size for Monte Carlo simulation set to 1.6 × 109. Three algorithms for energy calculation like Bennets, Free energy perturbation (FEP), and Thermodynamic Integration (TI) were used. Agreement value of < 1 Kcal/mol among these, is regarded reasonable and indicates the highly stable nature of the complex.

US

3. Results and discussion

CE

PT

ED

M

AN

Current study presents the applications of computer-aided drug designing for exploring Ddl as a potential drug target against many intractable complications of A. baumannii. With the recently reported Ddl crystal structure, there is a dire need for understanding the conformational changes, which this enzyme undergoes in the catalytic mechanism of the drug-like natural inhibitors. This could provide meaningful insights, guiding the discovery of novel antibiotics against A. baumannii. The protein after retrieval from PDB was subjected to a preparation phase, where only one monomer of the structure was retained while the second monomer was removed (Fig.1). This was logical as both the monomers contain the same active cavity and proceeding further with one subunit could reduce the computational efforts. Afterward, the protein was minimized for total steps of 1500. Superimposing the minimized protein over the original structure revealed Root Mean Square Deviation (RMSD) value of 0.4 Å, which indicates the improved quality of the protein after minimization. 3.1. Active site prediction

AC

The mapping of protein active site for ligand binding is significant for a wide range of applications including de novo drug designing, functional site comparison, and molecular docking studies [50]. The protein active site was predicted through multiple ways. First, an online server, Metapocket [33] was used. Two potential sites were predicted (S-Table.1). The first cavity was large comprises 31 amino acids while the second cavity was composed of 15 amino acids. It was observed that the second cavity lies within the first predicted cavity as a sub-cavity. The cavity was further confirmed from DoGSiteScorer (S-Fig.2) [34] and literature search [28]. In order to select an amino acid atom coordinates for molecular docking through GOLD and AutoDock Vina, MSA approach was utilized (Fig.2) [35]. The basic idea behind this was to look for the residue in protein active site that is highly conserved among Ddl enzyme orthologues in bacterial pathogens. In this regard, OE2 atom of residue Glu247 was set as a site for docking of compounds. The residue is present in the central domain and is completely conserved among the various orthologues of Ddl 7

ACCEPTED MANUSCRIPT

ED

M

AN

US

CR

IP

T

enzymes. As the residue is highly conserved and vital for enzyme reaction mechanism, any point mutation of this residue will lead to loss of enzyme catalytic function [51, 52], indicating less chances of point mutation in this residue of the enzyme and hence lesser opportunity for generating resistance [53].

3.2. Molecular docking

PT

Fig.1. The three dimensional structure of Ddl monomer. The C-terminal, Central and N-terminal domains are shown in orange, purple and cyan, respectively.

AC

CE

Molecular docking phase was initiated with a virtual screening of Asinex antibacterial library. The library contains 4800 compounds at the time of retrieval. To proceed only with drug-like compounds, the library was screened for inhibitors that fulfill all the parameters of Lipinski rule of five. According to this rule, compounds will be considered drug-like if their molecular weight is less than 500 dalton (Da), log P score < 5, and hydrogen bond donors and acceptors < 5 and 10, respectively. This screening was accomplished using LigandScout 4.1, which revealed 1440 compounds that completely satisfy Lipinski rule of five. The compounds were minimized and subjected to molecular docking. In molecular docking, structure-based virtual screening of druglike compounds was performed to identify compounds having better affinity for the receptor enzyme active site. Inhibitors were first docked into enzyme active site using GOLD, followed by docking of the best sorted 10 compounds using AutoDock Vina. The top ten best inhibitors based on GOLD fitness score and AutoDock Vina binding energy along with druglike rules which they voilates are shown in Table.1. The correlation coefficient between GOLD fitness score and AutoDock Vina binding energy can be found in Fig.3. After comparative analysis, compound8

ACCEPTED MANUSCRIPT

PT

ED

M

AN

US

CR

IP

T

331(6-(4-((3- methoxyphenylsulfonamido) methyl) phenyl)-2-methylnicotinamide)) was found the best docked

AC

CE

Fig.2.Multiple sequence alignment of Ddl enzymes from different bacterial enzymes. (2187,Staphylococcus aureus), (3E5N, Xanthomonas oryzae), (3LWB,Mycobacterium tuberculosis), (3R23,Bacillus anthracis),(4FUO,Enterococcus faecalis), (5D8D,Acinetobacter baumannii), (3V4Z,Yersinia pestis), (4DGJ,Burkholderia xenovorans), (4EGO,Burkholderia ambifaria), (4EGQ,Burkholderia pseudomallei). inhibitor in both GOLD and AutoDock Vina. The GOLD fitness score and AutoDock Vina binding energy attained by the compound was 84.2 and -7.2 kcal/mol, respectively. In both tools, the compound was investigated to docked in the same position and interacts with almost same residues of the active site (Fig.4). Importantly, the compound follows all prominent rules of druglikenss i.e. Lipinski [30], Ghose [38], Muegge [38], Veber [38], and Egan rules [38] and increases the likelihood of surviving the well-documented high rates of attrition in drug discovery. In GOLD, the compound was found docked deep inside the cavity of the enzyme (Fig.5), with minor effect on the overall structure. It was observed the compound shields majority of the active cleft length present between the Central and C-domains of the protein. The methoxybenzene ring (anisole) of the compound was found extended towards N-terminal of the enzyme. These findings support the 9

ACCEPTED MANUSCRIPT

ED

M

AN

US

CR

IP

T

compound credibility in hindering the accessibility of the active site to the substrate and subsequent enzyme-substrate complex formation and final product formation. The preferred orientation of the compound was in such a way that formamide moiety of the 2methylnicotinamide interacts with Trp177 residue of the loop region from the Central domain. This results in tilting the 2,6-dimethylpyridine ring to make possible a strong hydrogen bond of 1.6 Å between nitrogen atom of 2,6-dimethylpyridine and hydrogen of Lys140 from the Central domain. The major contribution from the compound towards interaction with the enzyme active site came from sulfonamide group. This moiety interacts with residues from both the Central domain and the C-terminal domain. All the atoms of sulphonamide group tend to strongly interact with both these domains notably the Arg232 from C-terminal domain and Ser146 from the Central domain. The anisole ring was positioned in such a way to make possible interactions between its oxygen atom and Glu20 and Glu71 residue of N-terminal domain of the enzyme. In comparison to GOLD, AutoDock Vina generated pose binds totally in opposite direction. The formamide moiety of 2-methylnicotinamide ring of the inhibitor tends to move towards the N-terminal of the enzyme and favors to interact with this domain residues through its amino ethanol site chain. The hydrogen atom of pyridine ring of the inhibitor interacts with Ser146 oxygen through a strong hydrogen bond of strength 2.5 Å. The middle benzene ring of the inhibitor in contrast to GOLD move a bit up and rotate in clockwise direction. This movement allows the inhibitor to positioned the sulfonamide group upward. In AutoDock Vina, this pose resulted in no interaction between sulphonamide atoms and receptor residues. The behaviour could be assumed to continue for proper positioning of the inhibitor in the active cavity to achieve maximum interactions between the inhibitor atoms and receptor active site residues.

81

82

-3 -4

CE

-2

PT

-1

83

84

GOLD fitness score 85

86

87

88

89

y = -0.3157x + 18.654 R² = 0.7575

AC

AutoDock Vina binding energy (kcal/mol)

0

-5 -6 -7 -8 -9

-10

Fig.3. Correlation coefficient between GOLD fitness score and AutoDock Vina binding energy for the top ten best inhibitors. 10

ACCEPTED MANUSCRIPT

T P

I R

C S

U N

A

D E

M

T P

E C

C A

11

ACCEPTED MANUSCRIPT Fig.4. 2D representation of binding mode and interactions of the best characterized drug-like inhibitor in the active pocket of Ddl enzyme.

T P

I R

C S

U N

A

D E

M

T P

E C

C A

12

AC

CE

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

13

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

Fig.5. (A) Binding mode of compound-331 in active pocket of KdsA enzyme. (B) Closer view.

14

ACCEPTED MANUSCRIPT Table.1. Top ten best natural inhibitors shortlisted in the current study. GOLD fitness score

Compound

AutoDock Vina binding energy (kcal.mol-1)

Drugllikenss rule(s) violation

-9.2

Ghose rule 3 violations: MW>480, MR>130, #atoms>70

-8.5

Ghose rule 3 violations: MW>480, MR>130, #atoms>70

T P

I R

C S

U N

88.3

A

Compound-965

D E

M

T P

C A

E C

85.6

15

ACCEPTED MANUSCRIPT

Compound-1150

T P

I R

C S

85.1

-8.8

Ghose rule 1 violation: WLOGP>5.6 Muegge filter 1 violation: XLOGP3>5

-8.1

Ghose rule 2 violations: MR>130, #atoms>70

U N

A

Compound-1137

D E

M

T P

E C

C A

85.0

Compound-1131 16

ACCEPTED MANUSCRIPT

T P

-7.2

No violations

82.9

-7.6

No violations

84.2

I R

C S

U N

Compound-331

A

D E

T P

E C

M

C A

Compound-1130

17

ACCEPTED MANUSCRIPT

T P

I R

82.5

-6.9

No violations

C S

U N

A

Compound-355

D E

M

T P

C A

E C

Ghose rule 82.2

-7,4

1 violation: MW>480

18

ACCEPTED MANUSCRIPT Compound-1139

T P

I R

C S

82.1

-7.3

No violations

-81.4

-7.4

No violations

U N

A

D E

Compound-1133

M

T P

E C

C A

Compound-370 19

ACCEPTED MANUSCRIPT MW, molecular weight, MR, molecular refractivity, #atoms, number of atoms

T P

I R

C S

U N

A

D E

M

T P

E C

C A

20

ACCEPTED MANUSCRIPT 3.3. SwissADME and preADMET analysis

AC

CE

PT

ED

M

AN

US

CR

IP

T

The high absorption and well distribution of a drug in required time are important for its effective metabolism and action. In addition, toxicity that often overshadows ADME behaviour is another important consideration [54]. This is because majority of the drugs, which fails in clinical stages due to high toxicity proves very expensive and detrimental in the drug designing process [55]. Computational prediction of druglikness together with ADMET properties predication at the drug design stage could possibly provide an array of opportunities for accelerating lead compounds identification with predicted biological activities [55]. Druglikeness is the qualitative inspection of compound physicochemical or structural properties, to investigate the compound likelihood as an oral drug-like candidate [56]. Five filters were used in order to ensure the improved quality of the inhibitors for future lead optimization. The compound was unveiled to follow all the prominent druglike rules including Lipinski rule of five [30], Ghose rule [38], Veber filter [38], Egan filter [38], and Muegge filter [38].The molecular weight of compound is 419.5 g/mol, number of heavy atoms (29), number of aromatic heavy atoms (12), number of rotatable bonds (7), number of hydrogen bond acceptor and donor 7 and 3, respectively. The presence of rotatable bonds pointed the inhibitor to be a good adaptive inhibitor candidate. Adaptive inhibitors achieve serval alternative functionalities at variable regions of the targeted proteins and as such are regarded as to overcome point mutations that are responsible for decreasing inhibitor binding affinity [57]. The Topological Polar Surface Area (TPSA) of the compound is 122.92 Ų. TPSA calculates the presence of polar amino acids at surface of compound [58]. Lower TPSA values for a given compound is preferred as higher values diminishes compound membrane permeability and act as a substrate for p-glycoprotein [54]. This protein is a drug efflux protein and aids in lowery intracellular drug concentration [38, 54]. The ADME behaviour of the compound was evaluated from its pharmacokinetics. Better absorption of the compound illustrated higher absorption of the compound from the intestinal tract upon oral administration [54]. The carcinogenicity test in rats revealed the compound non-carcinogenic, while the in vitro Ames test in TA100 and TA1535 strains of Salmonella typhimurium declared the compound non-mutagen. The in vitro Human ether-a-go-go-related gene channel inhibition [55] was found ambiguous for the compound. The water solubility of the compound was investigated through three models of SwissADME server. According to Estimated Solubility (ESOL) model [59], the log S value for the compound is -3.38 and is placed in soluble category. Log S [60]; indicates compound solubility; less the log S score, higher will be the absorption and viseversa. The second model, which adapted from Ali et al, predicted the log S value of -3.94 with category of soluble. In third model, which proposed by SILICOS-IT (http://silicos-it.be.s3-website-eu-west-1.amazonaws.com/software/filterit/1.0.2/filter-it.html), the compound score log S value of -5.66 that placed the compound in moderate soluble category. The inhibitor molecule was found less skin permeate as indicated by the negative log Kp value (-7.6 cm/s). One of the concerns associated with the inhibitor is being a substrate for the P-glycoprotein. This is particularly important as the action of P-glycoprotein results in lowering the intracellular concentration of drug. Further, the compound was unraveled as non-inhibitor of cytochrome P450 superfamily isozymes. Solubility of compounds greatly 21

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

facilitates drug development process particularly the handling and formulationn phase [38]. The log Po/w value, which is a partition coefficient between n-octanol and water is an important descriptor for lipophilicity [61]. In SwissADME, the log Po/w is the arithmetic mean of values generated by five methods (iLOGP, XLOGP3, WLOGP, MLOGP, and SILICOS-IT). The value of log Po/w for the compound is 2.21. The bioavailability score of the compound is 0.55. The druglikness and pharmacokinetics of this compound were compared with D-cycloserine (a broad spectrum antibiotic used for the treatment of tuberculosis and inhibit Ddl and Alanine recemase enzyme). It was revealed that in contrast to the compound of our study that follows all drug rules, D-cycloserine violates two drug rules: Ghose and Muegge. D-cycloserine has low GastroIntestinal absorption when compared to our scrutinized compound and placed in low category. The bioavailability score of D-cycloserine was found the same predicted for our compound of interest. Lipophilicity of D-cycloserine was described as improved, and have enhanced solubility in water similar to the compound reported in the study. One of the major limitation in the use of Dcycloserine is its neurological side effects [14], as D-cycloserine penetrates into Central Nervous System (CNS) and results in several adverse effects. The screened compound in this study, on the other hand, could be assumed to have low CNS penetration due to its molecular weight > 400 Da and TPSA value of higher than > 90 Ų, thus can be less likely be related with CNS penetration associated effects. As the protein and compound contains ionisable groups, the knowledge of their pKa was essential [62-64]. Therefore, pKa value, which can be define as the negative log of the dissociation constant, calculation of the drug and the entire protein molecule was determined using online version of Chemicalize ChemAxon [65] and H++ software [66]. The pKa calculations results for the compound can be found in S-Fig.3 while titration curve for the protein can be seen in S-Fig.4. List of all computed pKs can be found S-Table.2 while residues that contribute the most to each pK shift is tabulated S-Table.3. The strongest acidic and basic pKa value of the compound was revealed 9.99 and 8.15, respectively. The strongest acidic pKa value was for the NH of sulphonamide while the strongest pKa was for NH2 group of 2-methylnicotinamide ring. The strongest acidic pKa value of NH of sulphonamide is due to resonance of electron pair towards the Sulphur after dissociation of the proton while the strongest basic pKa value of NH2 is due to high ability of accepting protons. Six microspecies of the compound was observed and there distribution at different pH can be seen in S-Fig.3. Overall, the predicted compound of this study have improved druglikenss and pharmacokinetics and can be targeted further for designing potent derivatives having improved pharmacokinetics and ADMET properties. 3.4. Molecular dynamics simulation Although several static structures of the Ddl enzyme have been proposed, alone or with ligand, yet the dynamics of the enzyme in aqueous milieu is still needed to decipher under physiological conditions. The dynamics of protein alone and in complex with the best inhibitor were investigated for a timescale of 100-ns.This leads to in-depth understanding of the structural adjustment of the protein adopted in the presence or absence of ligand and can potentially guide the discovery of novel antibiotics against the enzyme. Four statistical parameters; RMSD, Root Mean Square 22

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

Fluctuation (RMSF), Beta-factor (β-factor) and Radius of Gyration (Rg) were performed to unveil the protein stability over the course of simulation period (Fig.6). RMSD is the average measure of distances between backbone atoms of superimposed proteins. The relative fluctuations in the RMSD of Cα atoms of undocked protein was seen higher when compared to the docked, indicating the stability of enzyme-inhibitor complex. The mean RMSD value estimated for docked complex was 1.6 Å and the maximum RMSD value was 2.4 Å at 67th ns while the mean RMSD for undocked system was 1.9 Å. Intriguingly, two main deviations were noted in undocked system when compared to docked. The first deviation occurred from 0 to 27-ns with maximum RMSD of 2.3 Å (25th ns). The second deviation was observed from 62th ns till 100-ns with maximum RMSD of 2.4 Å. On the contrary, it was observed that studied inhibitor sustained its interactions and binding with Ddl enzyme with no aforesaid deviations. The convergence values for both the systems and to reveal the structural mobility for each residue of the system, RMSF values were monitored. It is quite evident that in bounded form the enzyme residues remained highly stable in contrast to unbounded form. The mean RMSF values for undocked and docked system were reported 1.08 Å and 0.90 Å, respectively. As depicted clearly in Fig.6.B, majority of the residues of docked protein have RMSF values less than 2 Å. This further demonstrates the highly stable nature of Ddl enzyme. In undocked protein, the residues from Ser2, Glu155 to Gln169 and Ala172 to Gly182 of the Central domain, Glu21 to Glu44, and Gly74 of the N-terminal domain, were found to have comparatively high flexibility than others but still have in stable range. In docked protein, residue Lys18, 47-54, 189-190, 197-200, 239-240, were reported to have comparatively high RMSF values than in undocked form. Binding of ligand modified the enzyme structure in such a way that it reduces the distances between the interacting active residues and move the other residues mentioned above bit far. The most important active pocket residues; Glu20, Glu71, Ser146, Lys140, Trp177, Arg232 in undocked system were revealed with an average value of 1.2 Å while in docked system their value was found 0.7 Å. The systems were further subjected to β-factor calculation. β-factor estimates flexibility and thermal stability of protein and its side chains with respect to time. A complete coherence was observed between β-factor values and RMSF values for the systems. The mean β-factor value for undocked and docked systems were determined as 38.2 Å and 26.5 Å, respectively. In last, equilibrium conformation of the systems was investigated for better understanding of the systems. This was achieved by calculating Rg, which estimates the compactness of protein structure. Variations in the Rg measures determines compactness of protein during simulation run. Higher Rg value implies less tight packing of protein amino acids whereas lower Rg value entails tight packing of protein amino acids. The mean Rg value computed for undocked system was 22.4 Å, with maximum value of 23.06 Å at 49thns. For docked system, the maximum value of 22 Å was observed at 48 ns with mean value of 22.3 Å. The Rg values reported herein for both the systems were analogous to the RMSD values and further validate the stability of the protein in docked form.

23

ACCEPTED MANUSCRIPT 3.5. RDF and AFD analysis

AC

CE

PT

ED

M

AN

US

CR

IP

T

A Tcl/Tk script was designed to screen protein residues involved in hydrogen bonding and intermolecular contacts at the start (first ten nanoseconds) and towards the end of simulation (last ten nanosecond) was utilized in the VMD software. The purpose of this was to shed light on vital residues of the enzyme active site, which played a significant role in recognizing, binding, and stabilizing the inhibitor over the entire course of simulation. The shortlisted residues were then subjected to RDF [22] and AFD [23] to complement their role in inhibitor binding. Analysis comprised of a correlation between enzyme active residues and ligand atoms on a XY plane in the form of 3D histogram. Only four interactions involving Lys176 and Trp177 residues of the pocket and inhibitor hydrogen and nitrogen atoms were selected for RDF and AFD, due their profound role in complex stability and structural integrity. Both the protein residues were revealed to efficiently trap the inhibitor by virtue of hydrogen bonds. RDF graphs were generated for the four aforesaid interactions and are illustrated in Fig.7.The first plot in the figure represents the interaction distance between Lys176-O and inhibitor-H13. At the start 10-ns of simulation, the largest peak appears at 1.71 Å having a g (r) value of 1.89 while towards the end 10-ns, the largest peak appears at 2.00 Å having a g (r) value of 1.30. The graph demonstrates that both at the start 10-ns and end 10-ns of simulation period, the oxygen atom of Lys176 and H13 of inhibitor spend majority of the time in close proximity of each other to exhibit the activity. For Lys176-O and inhibitor H14 atom, the largest peak at the start 10-ns of simulation was observed at 3.39 Å with a g (r) value of 0.4. Towards the end 10-ns of simulation, the peak shortened and broadened abit with the highest peak noted at 3.60 Å having a g (r) value of 0.36. The graph explains that towards the end 10-ns of simulation the distribution of inhibitor H14 atom around Lys176 oxygen is more dispersed and moved away about 0.2 Å from each other. The RDF for Lys176-O and inhibitor-N revealed that at the starting 10-ns time of simulation the largest peak is present at 2.84 Å having a g (r) value of 0.92, while towards the end 10-ns the largest peak was observed at 2.85 Å having a g (r) value of 0.56. Although towards the end 10-ns of simulation, the distribution of inhibitor atom is dispersed more than at start 10-ns of simulation, however, it was noted that atoms tend to come close to each other. The graphs generated for Trp177 nitrogen atom and inhibitor H13 atoms were found bit distorted at the start 10-ns and towards the end 10-ns of simulation. The maximum possibility of finding inhibitor H13 atom around Trp177 nitrogen was found at 3.34 Å with g (r) value of 0.29 while towards the end, the H13 atom of inhibitor was found at maximum of 3.64 Å with g (r) value of 0.28. Towards the end, the increase distance of 0.1 Å was observed. All the graphs infer that these four interactions held the ligand in the close proximity of the protein active cavity, and are contributing factors in formation of stable complex. The minor variations depict little adjustment of inhibitor in the protein active pocket for achieving the best orientation.

24

ACCEPTED MANUSCRIPT

T P

I R

C S

U N

A

D E

M

T P

E C

C A

Fig.6. RMSD (A), RMSF (B), Beta factor (C), and Radius of Gyration (D) for enzyme and enzyme-inhibitor complex. 25

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

Further, the RDF demonstrate the system composition in a correct way as indicated by the lower flexibility of the interaction and showed good agreement between the findings of docking and MD simulations. In simulation trajectories particles are positioned with respect to time, which are subjected to several different statistical parameters for explaining the molecular nature of the system under study. Among these, RDF gives insights about the local displacement of an atom with respect to a reference. RDF is less responsive to local structural reorganizations that are mainly because of the moving inhibitor in protein active cavity to achieve an energetically stable confirmation as the simulation advances. To get more understanding of local movements and structural changes, AFD was utilized. Generally, the pattern of AFD graph signifies that the inhibitor is retained in enzyme active site for the whole time of simulation. The contour plot for AFD analysis is illustrated in Fig.8. The AFD plots were observed with maximum values having little variation and small distribution area from top to down graph perspective that indicate highly stable nature of inhibitor in enzyme active pocket. In the AFD plot generated for oxygen atom of Lys176 and ligand H13 atom, it was observed that the inhibitor moved closer to protein active site both on X-axis and Y-axis towards the end 10-ns of simulation. Although, the maximum value is decreased towards the end 10-ns of simulation yet the density distribution of inhibitor atom is more refine. Area of distribution was revealed to decreased by 0.5 Å. A slight shift in inhibitor displacement was also observed that may be due to increased affinity for the protein. In case of oxygen of Lys176 and inhibitor H14 atom, the distribution area towards the end of simulation is slightly increased. The inhibitor moved both on X-axis and Y-axis towards the center of Lys176 oxygen atoms, which signifies the stability of this interaction towards the end 10-ns of simulation. The tilting behavior of the inhibitor was seen in all the distribution. For oxygen atom Lys176 and inhibitor N3, the inhibitor atom was revealed to come close in vicinity of Lys176 oxygen atom mediated by inhibitor rotation. Towards the end 10-ns of simulation the distribution area was more refined when compared to the start 10-ns of simulation, however, its magnitude decreases. In case of Trp177-N and inhibitor-H13, the distribution area had increased a little towards simulation end that indicate the inhibitor atom movements towards on X-axis plane. The magnitude of graph was also found decreased at the end of simulation. By comparison of RDF graphs there is no significant shift in the distance between the interacting partner but the AFD graphs represent a tilting behaviour of inhibitor towards the protein active site. This behavior of inhibitor can be seen with visual representation of snapshots of simulation at 1-ns and at 100-ns (Fig.9). From the figure side groups of the inhibitor are shifting more towards the protein. AFD highlighted these movement which RDF cannot represent using only distance for analysis. Overall these graphs interpret that through the entire simulation period, the inhibitor is present in the close vicinity of active pocket and remained bonded through multiple strong hydrogen bond interactions with protein active site residues.

26

AC

CE

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

27

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

Fig.7. Radial distribution function graph for: A. Lys-176-O-Lig-H13, B. Lys176-O-Lig-H14, C. Lys176-O-LigN3, D. Trp177-N-Lig-H13.

28

ACCEPTED MANUSCRIPT

T P

I R

C S

U N

A

D E

M

T P

E C

C A

29

ACCEPTED MANUSCRIPT Fig.8. Axial frequency distribution:A. Lys-176-O-Lig-H13, B. Lys176-O-Lig-H14, C. Lys176-O-Lig-N3, D. Trp177-N-Lig-H13.

T P

I R

C S

U N

A

D E

M

T P

E C

C A

30

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

Fig.9. Inhibitor movement from 0-ns to 100-ns.

M

3.6. MMGBSA based energy calculations

AC

CE

PT

ED

The accurate prediction of a drug affinity for its intended biological protein is considered as the “holy grails” in computer-aided drug designing [67]. In addition to the development of pharmaceutically active molecules with improved bioavailability and toxicity, the high affinity of a drug for its intended target holds promising outcomes in designing highly potent, selective, and eventually efficacious compounds. Though the presence of structural information aided the rationalization of target-ligand but still conformational changes of protein and ligand, the solvent effect, and the entropy-enthalpy compensation make this process a complex task [68]. Advances in computational power and theory especially in the last decade or so, the physics-driven computer simulation holds a great potential in predicting the accurate binding affinities of compound for their targeted biological macromolecules [69]. These methods naturally encounter the complicating effects generated because of discrete nature of entropy and solvent [67]. The estimated binding free energy calculations and intermolecular interactions between the compound331 and the enzyme were elucidated using MMGBSA [25]. Other elegant method for prediction binding free energies is Linear Interaction Energy (LIE) [70-72]. The comparison of MMGBSA with LIE and other binding free energy calculation methods have been done in several studies and was revealed to be highly depended on the tested system [73-79]. Therefore, it is difficult to predict which of the method is more accurate and elegant [73, 78].Generally, the accuracy of MMGBSA is comparable with LIE [78]. From expert point of view, MMGBSA is more preferred in many of the present studies thanks to its fast nature and often give better accuracy [78]. Compared to MMGBSA, LIE is highly system dependent and cannot readily be employed in studying molecules 31

ACCEPTED MANUSCRIPT without introducing explicit salt [78]. Altogether, MMGBSA was preferred in this present study. The calculations were performed on entire trajectories of MD simulation of enzyme-inhibitor complex. The estimated binding free energies obtained for the complex can be visualized in Table.2. Table.2. Estimated MMGBSA based binding free energy values for enzyme-inhibitor complex. Contribution

Estimated binding free energy values (kcal.mol-1) -50.30 -38.89 -89.2 41.83 50.8 -8.96 -47.36

US

CR

IP

T

∆Evdw ∆Eele ∆Egas ∆Gsolv,GB ∆GGB ∆Gsol-np Htot,GB

AC

CE

PT

ED

M

AN

The term entropy was excluded due to convergence problems in some cases where it failed to be calculated [23, 78]. The average of estimated binding free energies calculated for the complex was determined as -47.36 kcal/mol. Usually the calculated values of binding free energy are greater than the actual values and as such are referred as estimated or calculated values rather than absolute binding free energy in the manuscript [78, 80]. Major contributions towards complex stability came from van der Waals energy that dominates the overall calculated binding free energy. Van der Waals energy determined for the system was -50.30 kcal/mol. In contrast, non-favorable contributions arise from polar component of solvation free energy with value of 50.8 kcal/mol. The non-polar solvation free energy deduced from the system was -8.96 kcal/mol. In polar solvent, all molecular interactions are driven by polar and non-polar components, where non-favorable interactions and its associated energies are more promising. Electrostatic energies were also found dominant and aid significantly to the overall calculated binding energy of the system. The electrostatic energy calculated for the system was -38.89 kcal/mol. ∆Egas is gas-phase interaction energy and important because of its contribution from internal energy as well as from electrostatic and van der Waals interaction energies and can be represented as Egas = Eint +EELE +EVDW. The ∆Egas energy for the complex was -89.2 kcal.mol-1. Further understanding of the complex interactions at the atomic level was achieved by decomposing total energy of the system into each residue of the enzyme and inhibitor (Fig.10). Residues having estimated binding affinity of < -1 kcal/mol were referred as hotspot amino acids due to their major contributions in stabilizing the complex [81]. The following residues were found to have calculated binding energy of < -1 kcal/mol; Lys100 (-2 kcal/mol), Ile138 (-1 kcal/mol), Gly145 (-1 kcal/mol), Glu175 (-1.6 kcal/mol), Trp177 (-2 kcal/mol), Ile178 (-1 kcal/mol),Leu246 (-2 kcal/mol), and Asn249 (-1.3 kcal/mol). All these residues are important constituents of active pocket and have considerably low estimated binding energy. These residues are revealed to be present within or around the 32

ACCEPTED MANUSCRIPT enzyme active pocket and contribute substantially to hydrogen and hydrophobic interactions. It was also noticed that loss of these residues could result in loss of catalytic mechanism for substrate. Moreover, the calculated binding free energy of the active site residues that contribute significantly was decomposed into constituent like electrostatic, van der Waals, and polar solvation, and nonpolar solvation contribution (S-Table.4).The total calculated binding energy of complex, protein, inhibitor in 500 frames of simulation is shown in S-Fig.5.

T

3.7. WaterSwap based energy calculations

ED

M

AN

US

CR

IP

WaterSwap is a method of finding absolute binding free energy by overcoming the problems associated with cavitation and large value differences of double decoupling [26]. The use of explicit water model in WaterSwap means to include molecular details of protein-water-ligand, protein-water, and ligand-water interactions that are omitted in continum solvent methods like MMGBSA [26-27]. This is of particular importance in cases where water molecules act as bridging interactions between the protein and the ligand. WaterSwap has been successfully applied for mutants prediction of influenza neuraminidase enzyme that showed reduced drug affinity [82]. The binding free energy calculated through WaterSwap is presented in Table.3. It is evident from the findings of three algorithms that the agreement values is < 1 kcal.mol-1, indicating the high affinity of the inhibitor towards enzyme active pocket. For further understanding of the binding energy at protein residues level, the total energy was decomposed into residue level to shed light on the hotspot amino acids involved in inhibitor binding.It was revealed that active site residues:Lys176 and Trp177 have considerably low binding energy and played vital role in inhibitor recognition and binding throughout the period of simulation. The mean binding energy of Lys176 and Trp177 is -1.94 kcal.mol-1and -2.73 kcal.mol-1, respectively.

PT

Table.3. WaterSwap based absolute binding free energy calculation for enzyme-inhibitor complex.

CE

Algorithm

AC

Bennetts FEP TI

Binding free energy Kcal.mol-1 -19.01 -18.02 -18.99

33

ACCEPTED MANUSCRIPT

T P

I R

C S

U N

A

D E

M

T P

E C

C A

34

ACCEPTED MANUSCRIPT A266

Y269

D272

Q281Lig284 C275 I278

S2

T5

G8

A11

L14

K17

E20 V23

A263

-3.0

V260 H257

D26 Q29 L32

M254

L35

V251

S38

-8.0

V248

T P Q41

L245

A44

N242 E239

I R

-13.0

M236

-18.0

A227

-23.0

M

-28.0

D E

K212 T209

V56 D59 F62 L65

A

C218 Q215

T53

U N

A224 A221

R50

C S

I233 W230

P47

R68 E71 Q74 G77 E80

T P

L206 P203

N83 Y86

E C

Y200 R197 P194 Q191 L188 I185

T89 Q92

C A

A95 M98 V101 K104

E182

W107

T179

S110

K176

P113

M173 A170 Q167

P116 I119 K122

K164 A161

F158

A155

V152 M149

S146

G134 H143 K140 V137

A131

S128

D125

35

ACCEPTED MANUSCRIPT Fig.10. MMGBSA based binding free energy decomposition into each residue of the enzyme. Amino acids are reppresented by a single lettere code.

T P

I R

C S

U N

A

D E

M

T P

E C

C A

36

ACCEPTED MANUSCRIPT Conclusions

ED

M

AN

US

CR

IP

T

Among several targetable bacterial species, Ddl enzyme is one of the systematically searched and attracted candidate holding great promise for controlling the dissemination of A. baumannii resistant clones. The potential of natural antibacterial against DDI was revealed to be promising therapy and could serve as an alternative to D-cycloserine because of improved druglikeness and pharmacokinetics. Molecular docking unravels the deep binding of the inhibitors in the enzyme active cavity, validated by MD simulations as studied by the formation of highly stable complex between the enzyme and the best screened natural inhibitor over the course 100-ns of simulation. The active site core of the enzyme was found highly stable and showed great affinity towards the inhibitor by making strong hydrogen and hydrophobic interactions. RDF analysis demonstrated Lys176 and Trp177 as the main residues involved in stable interactions with the inhibitor while novel analytical AFD unveiled the retention of hydrogen between Lys176 and Trp177 and inhibitor throughout the simulation time. Estimated binding free energy based on MMGBSA depicted further, the exceptionally high stable value for the complex of -47.36 kcal.mol-1 dominated mainly by van der Waals interactions. The use of more sophisticated and recently introduced approach of WaterSwap discovered a more reliable agreement on the affinity of inhibitor towards enzyme active residues. The agreement value among Bennetts, FEP, TI algorithm is < 1 kcal.mol-1, signifying the highly stable nature of complex. Understanding the functionality of antibacterial enzymes in aspect of inhibitor binding and catalysis in addition to the knowledge of binding thermodynamics is vital in structure based medicinal chemistry research. Therefore, the outcomes of this study will provide significant directions for designing novel potent inhibitors with improved druglikeness, pharmacokinetics, and antibacterial activity that may serve as a treatment for A. baumannii related infections.

PT

Supplementary

S-Fig.1. The reaction mechanism of Ddl enzyme.

CE

S-Fig.2. DoGSiteScorer predicted binding cavity for Ddl enzyme.

AC

S-Fig.3. PKa calculations for the compound-331. S-Fig.4. Titration curve for the entire DDl molecule. S-Fig.5. MMGBSA free energy versus number of frames. S-Table.1. The potential ligand binding sites in Ddl enzyme predicted by Metapocket. S-Table.2. List of all computed pKs for the DDl enzyme. S-Table.3. Significant contributions to pK of each site from neighbouring groups. S-Table.4. MMGBSA free energy decomposition into most active residue of the enzyme.

37

ACCEPTED MANUSCRIPT Conflict of Interest The authors declare that they have no conflict of interest. Acknowledgment Authors are highly grateful to the Pakistan-United States Science and Technology Cooperation Program for granting the financial assistance.

T

References

A.M. Gonzalez-Villoria, V. Valverde-Garduno, Antibiotic-resistant Acinetobacter baumannii increasing success remains a challenge as a nosocomial pathogen, J. Pathog. 2016 (2016).

[2]

M.J. McConnell, L. Actis, J. Pachón, Acinetobacter baumannii: human infections, factors contributing to pathogenesis and animal models, FEMS Microbiol. Rev. 37 (2013) 130– 155.

[3]

C. for Disease Control, P. (US), Antibiotic resistance threats in the United States, 2013, Centres for Disease Control and Prevention, US Department of Health and Human Services, 2013.

[4]

C.R. Taitt, T. Leski, M.G. Stockelman, D.W. Craft, D. V Zurawski, B.C. Kirkup, G.J. Vora, Antimicrobial resistance determinants in Acinetobacter baumannii isolates taken from military treatment facilities, Antimicrob. Agents Chemother. (2013) AAC--01897.

[5]

W.H. Organization, others, Global priority list of antibiotic-resistant bacteria to guide research, discovery, and development of new antibiotics, Geneva World Heal. Organ. (2017).

[6]

D. Barton, Comprehensive natural products chemistry, Newnes, 1999.

[7]

T.D.H. Bugg, D. Braddick, C.G. Dowson, D.I. Roper, Bacterial cell wall assembly: still an attractive antibacterial target, Trends Biotechnol. 29 (2011) 167–173.

[8]

H. Barreteau, A. Kovač, A. Boniface, M. Sova, S. Gobec, D. Blanot, Cytoplasmic steps of peptidoglycan biosynthesis, FEMS Microbiol. Rev. 32 (2008) 168–207.

[9]

V. Škedelj, E. Arsovska, T. Tomašić, A. Kroflič, V. Hodnik, M. Hrast, M. Bešter-Rogač, G. Anderluh, S. Gobec, J. Bostock, others, 6-Arylpyrido [2, 3-d] pyrimidines as novel ATP-competitive inhibitors of bacterial D-alanine: D-alanine ligase, PLoS One. 7 (2012) e39922.

[10]

L.E. Zawadzke, T.D.H. Bugg, C.T. Walsh, Existence of two D-alanine: D-alanine ligases in Escherichia coli: cloning and sequencing of the ddlA gene and purification and characterization of the DdlA and DdlB enzymes, Biochemistry. 30 (1991) 1673–1682.

[11]

Y. Gholizadeh, M. Prévost, F. Van Bambeke, B. Casadewall, P.M. Tulkens, P. Courvalin, Sequencing of the ddl gene and modeling of the mutated D-alanine: D-alanine ligase in glycopeptide-dependent strains of Enterococcus faecium, Protein Sci. 10 (2001) 836–844.

AC

CE

PT

ED

M

AN

US

CR

IP

[1]

38

ACCEPTED MANUSCRIPT I. Tytgat, E. Colacino, P.M. Tulkens, J.H. Poupaert, M. Prévost, F. Van Bambeke, DDligases as a potential target for antibiotics: past, present and future, Curr. Med. Chem. 16 (2009) 2566–2580.

[13]

F.C. Neuhaus, J.L. Lynch, Studies on the inhibition of D-alanyl-D-alanine synthetase by the antibiotic D-cycloserine, Biochem. Biophys. Res. Commun. 8 (1962) 377–382.

[14]

P.-L. Lu, C.-F. Peng, J.-J. Hwang, Y.-H. Chen, Activity of twelve second-line antimicrobial agents against Mycobacterium tuberculosis in Taiwan, J. Chemother. 20 (2008) 202–207.

[15]

A. Kovač, V. Majce, R. Lenaršič, S. Bombek, J.M. Bostock, I. Chopra, S. Polanc, S. Gobec, Diazenedicarboxamides as inhibitors of D-alanine-D-alanine ligase (Ddl), Bioorg. Med. Chem. Lett. 17 (2007) 2047–2054.

[16]

D. Wu, Y. Kong, C. Han, J. Chen, L. Hu, H. Jiang, X. Shen, D-Alanine: D-alanine ligase as a new target for the flavonoids quercetin and apigenin, Int. J. Antimicrob. Agents. 32 (2008) 421–426.

[17]

G. Triola, S. Wetzel, B. Ellinger, M.A. Koch, K. Hübel, D. Rauh, H. Waldmann, ATP competitive inhibitors of D-alanine--D-alanine ligase based on protein kinase inhibitor scaffolds, Bioorg. Med. Chem. 17 (2009) 1079–1087.

[18]

A. Kovač, J. Konc, B. Vehar, J.M. Bostock, I. Chopra, D. Janežič, S. Gobec, Discovery of new inhibitors of D-alanine: D-alanine ligase by structure-based virtual screening, J. Med. Chem. 51 (2008) 7442–7448.

[19]

B. Vehar, M. Hrast, A. Kovač, J. Konc, K. Mariner, I. Chopra, A. O’Neill, D. Janežič, S. Gobec, Ellipticines and 9-acridinylamines as inhibitors of D-alanine: D-alanine ligase, Bioorg. Med. Chem. 19 (2011) 5137–5146.

[20]

G.E. Besong, J.M. Bostock, W. Stubbings, I. Chopra, D.I. Roper, A.J. Lloyd, C.W.G. Fishwick, A.P. Johnson, A De Novo Designed Inhibitor of D-Ala--D-Ala Ligase from E. coli, Angew. Chemie Int. Ed. 44 (2005) 6403–6406.

[21]

M. Sova, G. Čadež, S. Turk, V. Majce, S. Polanc, S. Batson, A.J. Lloyd, D.I. Roper, C.W.G. Fishwick, S. Gobec, Design and synthesis of new hydroxyethylamines as inhibitors of D-alanyl-D-lactate ligase (VanA) and D-alanyl-D-alanine ligase (DdlB), Bioorg. Med. Chem. Lett. 19 (2009) 1376–1379.

[22]

J. Donohue, Radial distribution functions of some structures of the polypeptide chain, Proc. Natl. Acad. Sci. 40 (1954) 377–381.

[23]

A. Abro, S.S. Azam, Binding free energy based analysis of arsenic (+ 3 oxidation state) methyltransferase with S-adenosylmethionine, J. Mol. Liq. 220 (2016) 375–382.

[24]

S. Ahmad, S. Raza, R. Uddin, S.S. Azam, Binding mode analysis, dynamic simulation and binding free energy calculations of the MurF ligase from Acinetobacter baumannii, J. Mol. Graph. Model. (2017).

[25]

B.R. Miller III, T.D. McGee Jr, J.M. Swails, N. Homeyer, H. Gohlke, A.E. Roitberg, MMPBSA. py: an efficient program for end-state free energy calculations, J. Chem. Theory

AC

CE

PT

ED

M

AN

US

CR

IP

T

[12]

39

ACCEPTED MANUSCRIPT Comput. 8 (2012) 3314–3321. C.J. Woods, M. Malaisree, S. Hannongbua, A.J. Mulholland, A water-swap reaction coordinate for the calculation of absolute protein--ligand binding free energies, J. Chem. Phys. 134 (2011) 02B611.

[27]

C.J. Woods, M. Malaisree, J. Michel, B. Long, S. McIntosh-Smith, A.J. Mulholland, Rapid decomposition and visualisation of protein--ligand binding free energies by residue and by water, Faraday Discuss. 169 (2014) 477–499.

[28]

K.-H. Huynh, M. Hong, C. Lee, H.-T. Tran, S.H. Lee, Y.-J. Ahn, S.-S. Cha, L.-W. Kang, The crystal structure of the D-alanine-D-alanine ligase from Acinetobacter baumannii suggests a flexible conformational change in the central domain before nucleotide binding, J. Microbiol. 53 (2015) 776–782.

[29]

V. Vukić, D. Hrnjez, S. Milanović, M. Iličić, K. Kanurić, E. Petri, Comparative Molecular Modeling and Docking Analysis of β-galactosidase Enzymes from Commercially Important Starter Cultures Used in the Dairy Industry, Food Biotechnol. 29 (2015) 248– 262.

[30]

C.A. Lipinski, Lead-and drug-like compounds: the rule-of-five revolution, Drug Discov. Today Technol. 1 (2004) 337–341.

[31]

G. Wolber, T. Langer, LigandScout: 3-D pharmacophores derived from protein-bound ligands and their use as virtual screening filters, J. Chem. Inf. Model. 45 (2005) 160–169.

[32]

T.A. Halgren, Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94, J. Comput. Chem. 17 (1996) 490–519.

[33]

B. Huang, MetaPocket: a meta approach to improve protein ligand binding site prediction, Omi. A J. Integr. Biol. 13 (2009) 325–330.

[34]

A. Volkamer, D. Kuhn, F. Rippmann, M. Rarey, DoGSiteScorer: a web server for automatic binding site prediction, analysis and druggability assessment, Bioinformatics. 28 (2012) 2074–2075.

[35]

S. Raza, G. Sanober, T. Rungrotmongkol, S.S. Azam, The vitality of swivel domain motion in performance of Enzyme I of phosphotransferase system; A comprehensive molecular dynamic study, J. Mol. Liq. (2017).

[36]

G. Jones, P. Willett, R.C. Glen, Molecular recognition of receptor sites using a genetic algorithm with a description of desolvation, J. Mol. Biol. 245 (1995) 43–53.

[37]

O. Trott, A.J. Olson, AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading, J. Comput. Chem. 31 (2010) 455–461.

[38]

A. Daina, O. Michielin, V. Zoete, SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules, Sci. Rep. 7 (2017) 42717.

[39]

S.K. Lee, G.S. Chang, I.H. Lee, J.E. Chung, K.Y. Sung, K.T. No, The PreADME: pc-based

AC

CE

PT

ED

M

AN

US

CR

IP

T

[26]

40

ACCEPTED MANUSCRIPT program for batch prediction of adme properties, EuroQSAR. 2004 (2004) 5–9. D.A. Case, V. Babin, J. Berryman, R.M. Betz, Q. Cai, D.S. Cerutti, T.E. Cheatham Iii, T.A. Darden, R.E. Duke, H. Gohlke, others, Amber 14, (2014).

[41]

P.K. Weiner, P.A. Kollman, AMBER: Assisted model building with energy refinement. A general program for modeling molecules and their interactions, J. Comput. Chem. 2 (1981) 287–303.

[42]

S. Andleeb, M.K. Rauf, S.S. Azam, A. Badshah, H. Sadaf, A. Raheel, M.N. Tahir, S. Raza, others, A one-pot multicomponent facile synthesis of dihydropyrimidin-2 (1 H)-thione derivatives using triphenylgermane as a catalyst and its binding pattern validation, RSC Adv. 6 (2016) 79651–79661.

[43]

M.G. Paterlini, D.M. Ferguson, Constant temperature simulations using the Langevin equation with velocity Verlet integration, Chem. Phys. 236 (1998) 243–252.

[44]

V. Kräutler, W.F. Van Gunsteren, P.H. Hünenberger, A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations, J. Comput. Chem. 22 (2001) 501–508.

[45]

H.C. Andersen, Molecular dynamics simulations at constant pressure and/or temperature, J. Chem. Phys. 72 (1980) 2384–2393.

[46]

D.R. Roe, T.E. Cheatham III, PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data, J. Chem. Theory Comput. 9 (2013) 3084–3095.

[47]

E.F. Pettersen, T.D. Goddard, C.C. Huang, G.S. Couch, D.M. Greenblatt, E.C. Meng, T.E. Ferrin, UCSF Chimera—a visualization system for exploratory research and analysis, J. Comput. Chem. 25 (2004) 1605–1612.

[48]

W. Humphrey, A. Dalke, K. Schulten, VMD-Visual Molecular Dynamics J Mol Graph 14: 33--38, (1996).

[49]

B.N. Dominy, Molecular recognition and binding free energy calculations in drug development, Curr. Pharm. Biotechnol. 9 (2008) 87–95.

[50]

J.A. Capra, R.A. Laskowski, J.M. Thornton, M. Singh, T.A. Funkhouser, Predicting protein ligand binding sites by combining evolutionary sequence conservation and 3D structure, PLoS Comput. Biol. 5 (2009) e1000585.

[51]

E.Y. Krynetski, J.D. Schuetz, A.J. Galpin, C.-H. Pui, M. V Relling, W.E. Evans, A single point mutation leading to loss of catalytic activity in human thiopurine Smethyltransferase, Proc. Natl. Acad. Sci. 92 (1995) 949–953.

AC

CE

PT

ED

M

AN

US

CR

IP

T

[40]

[52] V. Kairys, M.K. Gilson, V. Lather, C.A. Schiffer, M.X. Fernandes, Toward the design of mutation-resistant enzyme inhibitors: further evaluation of the substrate envelope hypothesis, Chem. Biol. Drug Des. 74 (2009) 234–245. [53]

N. Singh, M.P. Frushicheva, A. Warshel, Validating the vitality strategy for fighting drug resistance, Proteins Struct. Funct. Bioinforma. 80 (2012) 1110–1122.

[54]

C.M. Nisha, A. Kumar, P. Nair, N. Gupta, C. Silakari, T. Tripathi, A. Kumar, Molecular 41

ACCEPTED MANUSCRIPT Docking and In Silico ADMET Study Reveals Acylguanidine 7a as a Potential Inhibitor of β-Secretase, Adv. Bioinformatics. 2016 (2016). A. Wadood, M. Ghufran, S.F. Hassan, H. Khan, S.S. Azam, U. Rashid, In silico identification of promiscuous scaffolds as potential inhibitors of 1-deoxy-d-xylulose 5phosphate reductoisomerase for treatment of Falciparum malaria, Pharm. Biol. 55 (2017) 19–32.

[56]

G. Vistoli, A. Pedretti, B. Testa, Assessing drug-likeness--what are we missing?, Drug Discov. Today. 13 (2008) 285–294.

T

[55]

IP

[57] H. Ohtaka, E. Freire, Adaptive inhibitors of the HIV-1 protease, Prog. Biophys. Mol. Biol. 88 (2005) 193–208. P. Ertl, B. Rohde, P. Selzer, Fast calculation of molecular polar surface area as a sum of fragment-based contributions and its application to the prediction of drug transport properties, J. Med. Chem. 43 (2000) 3714–3717.

[59]

J.S. Delaney, ESOL: estimating aqueous solubility directly from molecular structure, J. Chem. Inf. Comput. Sci. 44 (2004) 1000–1005.

[60]

J. Ali, P. Camilleri, M.B. Brown, A.J. Hutt, S.B. Kirton, Revisiting the general solubility equation: in silico prediction of aqueous solubility incorporating the effect of topographical polar surface area, J. Chem. Inf. Model. 52 (2012) 420–428.

[61]

J.A. Arnott, S.L. Planey, The influence of lipophilicity in drug discovery and design, Expert Opin. Drug Discov. 7 (2012) 863–875.

[62]

Y.Y. Sham, Z.T. Chu, A. Warshel, Consistent calculations of p K a’s of ionizable residues in proteins: semi-microscopic and microscopic approaches, J. Phys. Chem. B. 101 (1997) 4458–4472.

[63]

R. Borštnar, M. Repič, S.C.L. Kamerlin, R. Vianello, J. Mavri, Computational Study of the p K a Values of Potential Catalytic Residues in the Active Site of Monoamine Oxidase B, J. Chem. Theory Comput. 8 (2012) 3864–3870.

[64]

M. Repič, M. Purg, R. Vianello, J. Mavri, Examining Electrostatic Preorganization in Monoamine Oxidases A and B by Structural Comparison and p K a Calculations, J. Phys. Chem. B. 118 (2014) 4326–4332.

[65]

M. Swain, Chemicalize. org, (2012).

[66]

R. Anandakrishnan, B. Aguilar, A. V Onufriev, H++ 3.0: automating p K prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations, Nucleic Acids Res. 40 (2012) W537--W541.

[67]

M. Aldeghi, A. Heifetz, M.J. Bodkin, S. Knapp, P.C. Biggin, Accurate calculation of the absolute free energy of binding for drug molecules, Chem. Sci. 7 (2016) 207–218.

[68]

C. Bissantz, B. Kuhn, M. Stahl, A medicinal chemist’s guide to molecular interactions, J. Med. Chem. 53 (2010) 5061–5084.

[69]

J. Wereszczynski, J.A. McCammon, Statistical mechanics and molecular dynamics in

AC

CE

PT

ED

M

AN

US

CR

[58]

42

ACCEPTED MANUSCRIPT evaluating thermodynamic properties of biomolecular recognition, Q. Rev. Biophys. 45 (2012) 1–25. [70]

J. Aqvist, J. Marelius, The linear interaction energy method for predicting ligand binding free energies, Comb. Chem. High Throughput Screen. 4 (2001) 613–626.

[71] H. Gutiérrez-de-Terán, J. Åqvist, Linear interaction energy: method and applications in drug design, in: Comput. Drug Discov. Des., Springer, 2012: pp. 305–323. W. Wang, J. Wang, P.A. Kollman, What determines the van der waals coefficient $β$ in the LIE (linear interaction energy) method to estimate binding free energies using molecular dynamics simulations?, Proteins Struct. Funct. Bioinforma. 34 (1999) 395–402.

[73]

S. Genheden, MM/GBSA and LIE estimates of host--guest affinities: dependence on charges and solvation model, J. Comput. Aided. Mol. Des. 25 (2011) 1085–1093.

[74]

S. Wong, R.E. Amaro, J.A. McCammon, MM-PBSA captures key role of intercalating water molecules at a protein- protein interface, J. Chem. Theory Comput. 5 (2009) 422– 429.

[75]

X. Barril, J.L. Gelpi, J.M. Lopez, M. Orozco, F.J. Luque, How accurate can molecular dynamics/linear response and Poisson--Boltzmann/solvent accessible surface calculations be for predicting relative binding affinities? Acetylcholinesterase huprine inhibitors as a test case, Theor. Chem. Acc. 106 (2001) 2–9.

[76]

P. Mikulskis, S. Genheden, P. Rydberg, L. Sandberg, L. Olsen, U. Ryde, Binding affinities in the SAMPL3 trypsin and host--guest blind tests estimated with the MM/PBSA and LIE methods, J. Comput. Aided. Mol. Des. 26 (2012) 527–541.

[77]

H.S. Muddana, C.D. Varnado, C.W. Bielawski, A.R. Urbach, L. Isaacs, M.T. Geballe, M.K. Gilson, Blind prediction of host--guest binding affinities: a new SAMPL3 challenge, J. Comput. Aided. Mol. Des. 26 (2012) 475–487.

[78]

S. Genheden, U. Ryde, The MM/PBSA and MM/GBSA methods to estimate ligandbinding affinities, Expert Opin. Drug Discov. 10 (2015) 449–461.

[79]

B. Kuhn, P.A. Kollman, Binding of a diverse set of ligands to avidin and streptavidin: an accurate quantitative prediction of their relative affinities by a combination of molecular mechanics and continuum solvent models, J. Med. Chem. 43 (2000) 3786–3791.

[80]

T. Hou, J. Wang, Y. Li, W. Wang, Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations, J. Chem. Inf. Model. 51 (2010) 69–82.

[81]

F.U. Haq, A. Abro, S. Raza, K.R. Liedl, S.S. Azam, Molecular dynamics simulation studies of novel β-lactamase inhibitor, J. Mol. Graph. Model. 74 (2017) 143–152.

[82]

C.J. Woods, M. Malaisree, B. Long, S. McIntosh-Smith, A.J. Mulholland, Computational assay of H7N9 influenza neuraminidase reveals R292K mutation reduces drug binding affinity, Sci. Rep. 3 (2013) 3561.

AC

CE

PT

ED

M

AN

US

CR

IP

T

[72]

43

AC

CE

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

44

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Graphical abstract

45

ACCEPTED MANUSCRIPT Highlights

CE

PT

ED

M

AN

US

CR

IP

T

The anti-Ddl activity of natural drug-like inhibitors is investigated. The screened inhibitor have better druglikenss compared to D-cycloserine. The docked complex is more stable than enzyme alone using MD simulation. Lys176 and Trp177 are hotspot amino acid from enzyme active site. MMGBSA and WaterSwap revealed high affinity of the inhibitor for Ddl active site.

AC

    

46