Pharmacophore-Based Virtual Screening and Molecular Dynamics Simulation for Identification of a Novel DNA Gyrase B Inhibitor with Benzoxazine Acetamide Scaffold.

DNA gyrase B is one of the enzyme targets for antimicrobial drug development, and its absence in mammals makes it a suitable target for the creation of safe antibacterial drugs. We identified six novel hits as DNA gyrase B inhibitors in the present study by employing 3D-pharmacophore structure-based virtual screening. The lead compounds complied with drug-likeness rules and lacked toxicity. Compound 4 (ZINC32858011) showed the highest inhibitory activity with an IC50 value of 6.3 ± 0.1 μM against the DNA gyrase enzyme. In contrast, the positive controls ciprofloxacin and novobiocin used in enzyme inhibition assay had IC50 values of 14.4 ± 0.2 and 12.4 ± 0.2 μM, respectively. The molecular docking of the six hits demonstrated that compounds 1, 2, 4, and 6 had suitable fitting modes inside the binding pocket. Molecular dynamics simulations were carried out for the six hits and the rmsd, rmsf, radius of gyration, and solvent accessible surface area parameters obtained from 100 ns molecular dynamics simulations for the six compounds complexed with a DNA gyrase B protein indicated that compound 4 (ZINC32858011) formed the most stable complex with DNA gyrase B. The binding free energy calculation with the MM-PBSA method suggested that the van der Waals interaction, followed by electrostatic force, played a significant role in the binding. Per-residue free binding energy decomposition showed that Ile78 contributed the most for the binding energy followed by Asn46, Asp49, Glu50, Asp73, Ile78, Pro79, Ala86, Ile90, Val120, Thr165, and Val167.

Introduction

Bacterial diseases have become one of the most serious global dangers in recent decades. The rise of pathogenic bacteria resistant to antibacterial medications in both community- and hospital-acquired diseases poses a severe danger to global health because currently existing drugs will no longer be effective[1] Gram-negative bacilli Escherichia coli and Salmonella enterica (subtype typhi), which can cause deadly diarrheal disorders, and Pseudomonas aeruginosa, which colonizes in immune compromised patients as those with cystic fibrosis, cancer, or AIDS are among the bacteria deemed a priority for new treatments.[2] Furthermore, despite the availability of effective antibiotics, the majority of clinical isolates of Gram-positive Staphylococcus aureus are resistant to a variety of medications. In addition, the mortality owing to S. aureus bacteria remains about 20–40%.[3,4] Drug resistance is becoming more common in this microorganism, posing a major concern, particularly in conditions such as diabetes and post-operative complications, which can result in amputations and death.[5] Depending on the target-based drug discovery technique, the bacterial DNA gyrase that is a kind of topoisomerase class II, is considered as a worthy target for the design and development of new potent antibacterial drugs.[6] This enzyme is heterodimer A2B2 formed of two subunits A and two subunits B,[7] where subunit A is involved in the breakage and reunion of DNA and subunit B displays ATPase activity.[8]

Among the main antibiotics that act as DNA gyrase inhibitors, three classes which have shown high efficiency. The first class is 6-fluoroquinolones such as norfloxacin (Figure A) which acts by binding to DNA gyrase A subunit interfering with DNA cleavage and reunion.[9] The 6-fluoroquinolone drugs are the only DNA gyrase inhibitors used in clinical tests; they bind to DNA gyrase A-DNA complex preventing the reunion of DNA strands,[10] consequently stopping the bacterial replication.[11] Quinolones are used mainly to treat urinary and respiratory infections; however, the development of serious side effects endeavor the researchers to explore novel DNA gyrase inhibitors.[12]

Figure 1

(A) Chemical structure of norfloxacin, (B) the chemical structures of drugs belonging to the aminocoumarin class.

The second inhibiting DNA gyrase includes coumarins such as novobiocin and clorobiocin.[13] They inhibit DNA gyrase B subunits. The structural features of aminocoumarins include 3-amino-4,7- dihydroxycoumarin moiety, l-noviosyl sugar, and an aromatic acyl component attached to the amino group, clorobiocin bears an additional group; a methyl pyrrole ring coupled to noviose sugar moiety (Figure B). Clorobiocin and novobiocin are natural compounds from the Streptomyces organisms inhibiting DNA gyrase B subunit by competitively binding to the ATP-binding site.[1] Clorobiocin and novobiocin act by impeding DNA gyrase B subunit by competitively binding to the ATP-binding site.[14] Clorobiocin (CBN) has greater potent antimicrobial activity but has never been used in clinics,[15] while novobiocin developed resistance to the bacterial strain and showed severe side effects, necessitating its withdrawn from the market.[16]

The third one is cyclothialidines (e.g., GR122222X).[17] which act by inhibiting DNA gyrase B subunit.[18] The analogs belonging to this class consist of five amino acids (two alanines, a cysteine, a serine, and a cis-3 hydroxyproline) attached to substituted resorcinol ring and dihydroxy-2,6-dimethylbenzoate. A lactone ring is formed between the resorcinol and the two amino acids cysteine and serine[19] (Figure ). Cyclothialidines showed higher potency and selectivity than coumarins, excellent broad in vitro antibacterial activity against gram-positive pathogens such as S. aureus, Streptococcus pyogenes, and Escherichia faecalis that could circumvent the issue of resistance to clinically used drugs but it lacked antibacterial activity in vivo due to their high lipophilicity.[19]

Figure 2

Chemical structures of drugs belonging to cyclothialidine class.

This study targeted DNA gyrase B subunit to develop novel antimicrobial agents by applying 3D pharmacophore-based virtual screening techniques.[20,21] Commercial compounds were retrieved from the ZINC database[22] and screened through the constructed 3D pharmacophore followed by molecular docking using MOE software (https://www.chemcomp.com/). In the computational studies, novobiocin and clorobiocin were used as reference compounds.

Based on comparative analyses, six compounds were prioritized, namely, ZINC9592782 (compound 1), ZINC8216235 (compound 2), ZINC23908246 (compound 3), ZINC32858011 (compound 4), ZINC20791579 (compound 5), and ZINC5346121 (compound 6). The six compounds were consequently analyzed for drug likeness properties and were observed to be potential drug-like candidates. The six compounds were further subjected to enzymatic biological assay against DNA gyrase B, ZINC32858011 (compound 4) was identified with promising with IC50 6.3 ± 0.1 μM.

In order to have insights about the interactions of the selected hits with the gyrase B and their stability within the active site, molecular dynamics (MD) simulation was carried for the six hits and the reference compounds. Post MD analyses were performed with as root mean square deviation (rmsd), root mean square fluctuation (rmsf), solvent accessible surface area (SASA), radius of gyration (Rg), H-bonds analyses, and hydrophobic interactions tools. Binding free energies were calculated for the six hits and the reference compounds using the MM-PBSA method.

Results and Discussion

Virtual Screening Procedure

The virtual screening strategy consists of 3D pharmacophore-based virtual screening and molecular docking. First, the 3D pharmacophore was built upon interaction of DNA gyrase B subunit and the ligand; clorobiocin (PDB code: 1KZN). Second, the ZINC database containing 300,000 commercially available compounds was filtered with the constructed 3D pharmacophore. Accordingly, 1800 small molecules that fit the features of the pharmacophore were obtained. The resultant 1800 compounds were screened by using Lipinski rules[23] to afford 975 compounds. Next, the remaining 975 molecules were then used in docking with MOE. The native ligand (clorobiocin) and novobiocin were used as positive controls in the virtual screening process and docked with the obtained compounds. The top 100 compounds with the highest score were selected and inspected visually. Six compounds were purchased for further biological evaluation. Figure illustrates the overall process of virtual screening, demonstrating the efficacy of computational techniques used in this study at screening a considerable number of drug-like small molecules, to identify novel anti-DNA gyrase B inhibitors. The discovered molecules are likely to be promising hits with a novel scaffold for further developing potent anti-DNA gyrase B inhibitors.

Figure 3

Outline of the virtual screening workflow for the identification of novel DNA gyrase B inhibitors.

Pharmacophore Based Virtual Screening

The protein–ligand interaction fingerprint (PLIF) protocol embedded in MOE was used to generate the 3D pharmacophore model based on the receptor–ligand complex of the DNA gyrase B crystal structure (PDB code: 1KZN). A set of well-known DNA gyrase B inhibitors and inactive compounds were prepared as a decoy set to validate the constructed 3D pharmacophore. The receiver operating characteristic curve (ROC) was 0.9, which indicated the reliability of the constructed pharmacophore model. As shown in Figure S1, the generated pharmacophore model consisted of seven features: four hydrophobic groups, (Val43, Ala47, Val71, Ile78, Pro79, and Ile90 residues). One hydrogen bond acceptor (F3) corresponded to Arg136, one hydrogen bond donor (F4) corresponded to Asp73, and one hydrogen bond donor/acceptor (F5) corresponded to Asn46. Around 300,000 commercial compounds were filtered through the constructed 3D pharmacophore model, with 1800 compounds fitting at least five of the model’s features. It could be noticed that the 3D pharmacophore model reduced the large database (300,000 molecules) to 1800 compounds, each of which was defined by its binding interaction with DNA gyrase B indicating the significant of this technique in virtual screening and drug design.

Docking Based Virtual Screening

The crystal structure of DNA gyrase B in a complex with clorobiocin (PDB code: 1KZN) was used for docking studies. The docking protocol was evaluated by comparing the docked pose of clorobiocin and its original crystal structure. The results showed that clorobiocin exhibited similar interactions to its original crystal structure with a rmsd of 0.4 Å. Also, the docking results of novobiocin showed good interactions and low rmsd value, these results suggesting that the docking procedure was efficient. The 975 compounds that emerged from the 3D pharmacophore and Lipinski’s filtration were docked with MOE software against DNA gyrase B, clorobiocin and novobiocin were docked with the same docking protocol and showed docking scores of −4.4 and −6.3 kcal/mol, respectively. The docked molecules were ranked based on their binding score. Thus, a list of 100 compounds was generated using the binding modes and docking scores of docked compounds. The 100 compounds were visually screened for certain characteristics; for example, molecules with scaffold diversity, at least one hydrogen bond with the key residues, and no binding pose in collision with the protein were chosen. The final shortlist of the small molecules was investigated by structure search of SciFinder database[24] and pan assay interference compounds (PAINS). Six novel hits, namely, ZINC9592782 (compound 1), ZINC8216235 (compound 2), ZINC2390824 (compound 3), ZINC32858011 (compound 4), ZINC20791579 (compound 5), and ZINC5346121 (compound 6) were selected and purchased for biological screening (Figure ). In this study, the compound number will be used instead of ZINC ID.

Figure 4

Molecular structure of selected six hit compounds and the reference compounds (novobiocin and clorobiocin), produced by 3D-pharmacophore based virtual screening.

Drug Likeness

The results of physicochemical properties analysis of the six hits are listed in Table . The findings revealed that the six hits satisfied the rule of five and did not violate Lipinski’s rule, suggesting that the hits are likely to have a good permeation and absorption in the biological system. Log S parameter is related to the water solubility of the candidates and affects their absorption; all the six molecules showed 6– ≤ log S ≤ −4, indicating that they were moderately soluble. Compounds with a total polar surface area TPSA ≤ 140 Å2 are effective in permitting the cell membrane. The six lead compounds had TPSA ≤ 140 Å2, expect for compound 5, which had showed a slightly higher value (147.3 Å2), which had no significant effect on the absorption. The bioavailability of the six hits showed a score 0.55 and a good absorption from GI, demonstrating that the 6 compounds are likely to be promising drug candidates. Medicinal chemistry analyses revealed that all the hits were passed both PAINS filter[25] and a list of 105 fragments recognized by Brenk et al.,[26] except for lead 6, which showed a coumarin alert. The synthetic accessibility scores of the compounds indicated a range between 3.54 and 4.14, revealing that the six molecules were easily synthesized. The six leads were observed to be non-carcinogenic. The drug likeness results of the six molecules suggested that they were likely good drug candidates due to high GI absorption, good orally bioavailability, and lower toxicity, encouraging us to investigate their biological activity.

Table 1

Drug Likeness Properties of Selected Compoundsa

	CP1 ZINC9592782	CP2 ZINC8216235	CP3 ZINC23908246	CP4 ZINC32858011	CP5 ZINC20791579	CP6 ZINC5346121
molecular weight	449.5	405.5	436.5	473.5	476.6	424.4
lipophilicity (log P)	4.1	2.4	2.0	2.5	1.8	3.3
H bond donor	0	1	0	0	1	0
H bond acceptor	5	4	5	4	6	5
log S	–5.4	–4.8	–4.1	–4.6	–4.1	–4.8
TPSA Å2	124.6	120.1	131.2	140.8	147.3	124.6
bioavailability	0.55	0.55	0.55	0.55	0.55	0.55
GI absorption	high	high	high	high	high	high
PAINS	0 alert	0 alert	0 alert	0 alert	0 alert	0 alert
Brenk	0 alert	0 alert	0 alert	0 alert	0 alert	1 alert
synthetic accessibility score	4.14	4.05	3.92	4.14	4.11	3.54
carcinogenesis	non	non	non	non	non	non

Molecular weight (<500 Da), lipophilicity (log P <5), H bond donor (≤5), H bond acceptor (≤10), log S (water solubility on scale −6 ≤ log S ≤ −4 moderate soluble), Total polar surface area (TPSA ≤ 140 Å2), Pains; pan assay interference compounds alert, Brenk; 105 fragments identified by Brenk database, synthetic accessibility score on a scale of 1–10 (1 easy to 10 difficult).

In Vitro Screening of Hit Compounds for Inhibitory Activity

The inhibitory activity of the six compounds that came out from virtual screening was investigated against DNA gyrase isolated from E. coli. The inhibitory activities of the six hits and reference compounds, ciprofloxacin and novobiocin, are illustrated in Table . The most potent compound against E. coli DNA gyrase B was compound 4 with an IC50 value of 6.3 ± 0.1 μM compared to ciprofloxacin and novobiocin which showed IC50 values of 14.4 ± 0.2 and 12.4 ± 0.4 μM, respectively. These finding indicated that compound 4 was the most potent DNA gyrase inhibitor.

Table 2

Inhibitory Activity of Promising Compounds against DNA Gyrase B

compound ID	ZINC ID	E. coli DNA gyrase supercoiling IC50 μM
compound 1	ZINC9592782	17.6 ± 0.4
compound 2	ZINC8216235	13.5 ± 0.3
compound 3	ZINC23908246	59.7 ± 1.4
compound 4	ZINC32858011	6.3 ± 0.1
compound 5	ZINC20791579	28.1 ± 0.7
compound 6	ZINC5346121	16.2 ± 0.4
ciprofloxacin		14.5 ± 0.2
novobiocin		12.4 ± 0.2

Molecular Docking of the Hits

Investigation of DNA gyrase B the loop formed by the residues 98–118 of the protein closes the ATP binding pocket. This loop is disordered in all structures.[27] The loop between residues 73 and 85 are well ordered and integral to the ligands’ interaction. This observation can be attributed to the low electron density of the loop, resulting in more mobility.[28] The binding pocket is lined with several hydrophobic residues, including Val43, Ala47, Ile78, Pro79, Ile90, Met91, Val120, Val167, Arg136, and Asn46. These residues are conserved amongst DNA gyrase A, B, and type II topoisomerases.[28] Asn46 coordinates with the Mg+2 ion in the ATP binding pocket, which is crucial for ATP hydrolysis.[28] However, Arg 136, plays an essential role in determining the coumarin resistance, as a variety of mutations at Arg136 were reported.[28]

The ligands clorobiocin, and novobiocin binds to the ATP binding site at its entry (Figure A,B), clorobiocin formed H-bond between the lactone coumarin and the conserved residue Arg136, and between hydroxy group of noviose and Asn46 at distance 2.8 Å distance. The NH of clorobiocin’s pyrrole ring formed H-bond with C=O of Asp73, and another H-bond between Asp73 and clorobiocin’s C=O via water (Figure A, Table ). Besides the network of H-bonds, hydrophobic interactions with hydrophobic residues Val71 and Val167 were observed. The hydroxybenzoate isopropyl moiety of clorobiocin rested outside the binding site and wrapped around Pro79, weakening the hydrophobic interactions with Ile78 and Pro79. The methyl pyrrole ring was deeply embedded in the hydrophobic pocket (termed as pyrrole pocket), which consisted of hydrophobic residues Val43, Ala47, Val71, and Val167.

Figure 5

Refined docking model of (A) clorobiocin (violet stick), (B) novobiocin (brown stick), (C) compound 1, (D) compound 2, (E) compound 3, (F) compound 4, (G) compound 5 and (H) compound 6. Left panel represented the binding mode, and the right panel represented the interactions between the ligands and DNA gyrase B protein. DNA gyrase B protein is shown in a magenta cartoon representation and the active site residues are shown in a cyan stick, ligands represented as brown stick. Hydrogen bond represented as black dotted.

Table 3

Docking Results of the Promising Hits within DNA Gyrase B Active Site

compounds ID	interacting group	amino acids	length (Å)	Docking score kcal/mol
Clorobiocin	O of lactone ring	NH1-Arg136	H-bond (3.0)	–4.4
	OH of sugar	C=O Asn64	H-bond (2.8)
	N2 of pyrrole ring	C=O Asp73	H-bond (2.8)
		Val43, Ala47, Val71, Ile78, Pro79, Met95, Val120, Val167.	hydrophobic interactions.
Novobiocin	C=O benzamide	NH Ala47	H-bond (3.5)	–6.3
		C=O Val71	H-bond (3.2)
	OH of 4-hydroxy benzamide	NH Gly77	H-bond (2.8)
		C=O Thr165	H-bond (3.2)
		C=O Asp73
	C=O coumarin
	NH of benzamide
		Val43, Ala47, Val71, Ile78, Pro79, Met95,	H-bond via water mediated
	NH of benzamide	Val120, Val167.	hydrophobic interactions.
compound 1 ZINC9592782	C=O of oxoethyl	NH Asn46	H-bond (3.5)	–5.7
		Val43, Ala47, Val71, Ile78, Pro79, Met95, Val120, Val167.	hydrophobic interactions.
compound 2 ZINC8216235	N pyrazole	NH Asn46 Val43, Val71, Val167, Ile90 and Met91	H-bond (3.7)	–4.9
			hydrophobic interactions
compound 3 ZINC23908246	S of benzothiazole C=O	NH Asn46	H-bond (3.3)	–4.6
		NH Arg76	H-bond (2.6)
		Val43, Val71, Val167, Ile90 and Met91	hydrophobic interactions
compound 4 ZINC32858011	C=O	NH Asn46	H-bond via water mediated	–4.9
		Val43, Val71, Val167, Ile90 and Met91	hydrophobic interactions
compound 5 ZINC20791579	NH	C=O Asn64	H-bond (3.5)	–4.5
		NH Arg136	H-bond (3.8)
	NH	NH Arg76	H-bond (3.8)
		Val43, Val71, Val167, Ile90 and Met91	hydrophobic interactions
compound 6 ZINC5346121	O of chromen C=O	NH2 Asn46	H-bond (3.5)	–5.7
		NH Gly77	H-bond (2.5)
		Val43, Val71, Val167, Ile90 and Met91	hydrophobic interactions

Investigation of the binding site of docked novobiocin showed an extensive H-bonding network, involving a coumarin ring and hydroxy benzamide groups. For example, the OH of coumarin formed H-bond with NH Gly77 (2.8 Å), and NH of benzamide formed H-bonds with C=O Thr165 at distance 3.2 Å, and the side chain of Asp73 via water molecule. Further, H-bonds were observed between C=O of benzamide and NH Ala47 (3.5 Å) and between the OH of 4-hydroxybenzamide and the main chain of Val71 (3.2 Å) (Figure B, Table ). In addition to H bonds, there were hydrophobic interactions with hydrophobic pocket, the 3′-isopentenyl-4′-hydroxybenzoate moiety embedded in the pyrrole pocket revealed hydrophobic interactions with Val43, Val71, Val120, Val167, Ile90, and Met91. Pro79 is located beneath the plane of a coumarin ring, the alkyl group of Arg76 lies above, and the dimethyl group of the sugar wraps around Pro79, leading to docking score of −6.3 kcal/mol.

The docking pose of compound 1 showed H-bond between C=O and Asn46 at distance 3.5 Å (Figure C), indicating that the compound was entirely embedded in the hydrophobic pocket (Figure C, Table ). As illustrated in Figure , no portion of compound 1 was found outside the binding site. Oxoethyl isobenzofuran moiety revealed hydrophobic interactions with Val43, Val71, and Val167. The dimethoxy isoquinolin interacted hydrophobically with an alkyl part of Arg76, Ile78, Pro79, and Ile90 forming a salt bridge with Glu50, Asp73, Gly77, and Thr156. A thiophene group lying in the pocket’s entry showed a salt bridge with Asn46 (Figure C). The binding energy score of the complex was found to be −5.7 kcal/mol. The binding mode of compound 2 showed that N of pyrazole formed a H bond with NH Asn64 at a distance of 3.7 Å (Figure D), similar to compound 1, compound 2 surrounded by the binding pocket with no protrusion extending beyond the pocket. 1-Phenyl-1H-pyrazolo[3,4-d] pyrimidin-4(5H)-one moiety that is wholly enclosed in the hydrophobic pocket formed hydrophobic interactions with Val43, Val71, Val167, Ile90, and Met91. Octahydroisoquinolin moiety showed hydrophobic interactions with Ile78 and Pro79 (Figure D). Compound 2 formed a salt bridge with Asp73, Arg76, Gly77, and Arg136. The complex was defined with a binding score of −4.9 kcal/mol.

The interactions between compound 3 and DNA gyrase B protein showed two H bonds: one between S of benzothiazole and NH Asn64 and another between C=O and NH Arg76 at a distance of at 3.3 Å e of 2.6 Å, respectively (Figure E). Moreover, the residues Val43, Ala47, Val 71, Ile 78, Pro79, Val120, and Val167 are involved in hydrophobic interactions with piperazine benzothiazole moiety. However, dihydrobenzo[b][1,4] oxazin-4-yl ethenone protrudes beyond the binding pocket in opposite direction from clorobiocin, indicating hydrophobic interactions with Leu52 and Ala53 (Figure E). The docking score of the complex is −4.6 kcal/mol. The binding affinity of compound 4 toward the DNA gyrase B active site was characterized by a H bond with Asn46 at distance 4.2 Å, and the complex was stabilized with a −4.9 kcal/mol docking score (Figure F, Table ). The 1-phenyl-1H-pyrazole moiety is located in the pyrrole pocket, forming hydrophobic interactions with Val43, Ala47, Val71, Ile90, Met91, Val120, and Val176. The alkyl group of Arg67 was positioned above the plane benzo [1,4] oxazine, while Pro79 was located beneath it, comparable to the coumarin of novobiocin. A piperidine moiety extended slightly outside the binding pocket, interacting with Ile90, Met91, and Val167 hydrophobic residues.

The docking results indicated that the compound 5—DNA gyrase B complex formed three H bonds with Asn46, Arg76, and Arg136, at distance of 3.5, 3.8, and 3.8 Å, respectively (Figure G, Table ). Further, benzyl-5,6,7,8-tetrahydro-3-morpholino-2,7-naphthyridine moiety was located within the hydrophobic back pocket, forming hydrophobic interactions with Val43, Val71, Ile78, Pro79, Val120, and Val167 residues. Ile78 and Pro79 were involved in hydrophobic interactions with the morpholino moiety. 3-Morpholinopropylamino protruded outside the binding pocket, facing Ala47 and Asp49. The docking score of the complex was found to be −4.5 kcal/mol. According to the molecular docking results, the complex compound 6—DNA gyrase B was stabilized within the active site via two H bonds with Asn46 and Gly77 located at distances of 3.5 and 2.5 Å, respectively (Figure H, Table ). The 2-oxoethoxy-3,4,7-trimethyl-2H-chromen-2-one was embedded in the hydrophobic pocket, forming hydrophobic interactions with Val43, Val71, Ile90, Met91, Val120, and Val167 residues. The hydrophobic residues Ile78, Pro79, and alkyl part of Arg76 interacted hydrophobically with 4-fluorophenylpiperazine. The complex was stabilized with a docking score −5.7 kcal/mol.

MD Simulations

The system stability of all six inhibitors/DNA gyrase B complexes, clorobiocin and novobiocin/DNA gyrase B complexes was determined using rmsd of Cα protein and ligand structures after 100 ns simulation (Figure , Table S1). Clorobiocin and novobiocin complexes showed average steady rmsd (0.21 ± 0.02 and 0.19 ± 0.03 nm), respectively, whereas a non-liganded protein showed an average rmsd values of 0.18 ± 0.03 nm. Compound 1 and compound 2 complexes reached the equilibrium with an average of 0.23 ± 0.04 and 0.24 ± 0.03 nm, respectively.

Figure 6

(A) rmsd of the unliganded DNA gyrase B protein (Black), clorobiocin-complex (yellow), novobiocin-complex (green). The colours represented unliganded DNA gyrase B protein (black), clorobiocin-complex (yellow), novobiocin-complex (green) were kept in figures (B–G). (B) Compound 1-complex (red), (C) compound 2-complex (red), (D) compound 3-complex, (E) compound 4-complex (red), (F) compound 5-complex (G) compound 6-complex.

Compound 3/DNA gyrase B complex showed a rmsd deviation around 0.25 ± 0.04 nm, while compound 4, 5, and 6 DNA gyrase B complexes showed steady equilibrium with rmsd average 0.19 ± 0.02, 0.22 ± 0.03, and 0.23 ± 0.03 nm, respectively (Figure , Table S1).

Although the non-liganded protein had the lowest rmsd, the rmsd pattern was not smooth through the simulation, and the protein’s rmsd fluctuated. These fluctuations were stabilized with all novel inhibitors. A compound 4-protein complex was the most stable complex, with a rmsd pattern and value comparable to that of the reference compound novobiocin.

In order to understand the stability of the inhibitors within the active site of the DNA gyrase B protein, rmsd of ligands, clorobiocin and novobiocin was calculated (Figure S2). Novobiocin reached the equilibrium at 25 ns and continued the remaining dynamics steadily with a rmsd average 0.13 ± 0.07 nm (Figure S2A). The rmsd pattern of clorobiocin showed fluctuations from 30 to 70 nm with a rmsd average 0.20 ± 0.09 nm (Figure S2B). Compound 1 demonstrated a rmsd value at 0.1 nm from the dynamics beginning until 60 ns (Figure S2C). It then adopted second conformation at 0.2 nm until the end of the simulations. Compound 2 reached the equilibrium at 28 ns with a rmsd average 0.28 ± 0.06 nm, at 80 ns compound 2 adopted second conformation (Figure S2D). Compound 3 reached the equilibrium at 20 ns with rmsd 0.3 ns until 60 ns, then continued with rmsd average of 0.23 ns in the remaining dynamics (Figure S2E).

Compound 4 exhibited a stable rmsd pattern over the 100 ns simulation with a rmsd average 0.26 ± 0.04 nm. Compound 5 showed small fluctuations overall the dynamics with a rmsd value around 0.17 ± 0.05 ns. On the other hand, compound 6 showed a rmsd average 0.12 ± 0.09 nm. From 60 ns, compound 6 adopted a second conformation with a rmsd average 0.15 nm. It could be concluded that novobiocin and compound 4 were the most stable ligands within the binding pocket. Moreover, compound 4 formed the most stable complex.

rmsf analysis was calculated for the non-liganded protein, DNA-gyrase B/clorobiocin, novobiocin, and compounds 1–6 complexes (Figure ). Fluctuations were observed in loops and coils. The great amplitude of fluctuations was observed in residues 94–96, which are part of loop that closes the ATP binding site. Residues 83–85 represented a last section of the active site loop which is characteristic by its mobility. Further, residues 83–85 did not interact with the reference and/or lead compounds. Residues 174–177 belong to another loop that is located far away the active site. Generally, key residues Asn46, Arg76, Gly77, Ile78, Pro79, Ile 90, Met91, Thr165, and Val167 showed low fluctuations amplitude (Figure ). The mean rmsf values were found to be 0.12 ± 0.05, 0.09 ± 0.06, 0.10 ± 0.03, 0.11 ± 0.05, 0.11 ± 0.04, 0.12 ± 0.02, 0.09 ± 0.04, 0.10 ± 0.06, and 0.11 ± 0.05 nm for non-liganded protein, clorobiocin, novobiocin, compounds 1–6, respectively (Table S1). It was observed that the non-liganded protein showed the highest average rmsf value 0.12 nm. On the other hand, clorobiocin, and compound 4 complexes showed the lowest mean rmsf among all complexes with a value around 0.9 nm, while compound 3 complex demonstrated the highest fluctuation of 0.12 nm, followed by compound 1, 2, and 6 with average rmsf 0.11 nm. These results revealed that the most stable compound was compound 4 and the least stable one was compound 3.

Figure 7

rmsf of the unliganded DNA gyrase B protein (Black), clorobiocin-complex (yellow), novobiocin-complex (green). The colours represented unliganded DNA gyrase B protein, clorobiocin-complex (yellow), novobiocin-complex (green) were kept in figures (A–F). (A) Compound 1-complex (red), (B) compound 2-complex (red), (C) compound 3-complex, (D) compound 4-complex (red), (E) compound 5-complex, (F) compound 6-complex.

Rg and SASA Calculations

The compactness of the protein–ligand complexes was evaluated by calculating the radius of gyration (Rg). The Rg of apo protein was found to be 1.58 ± 0.09 nm (Figure S3A,B). For DNA-gyrase B complex with clorobiocin, novobiocin, compound 1–6, Rg was found to be 1.61 ± 0.11, 1.60 ± 0.15, 1.61 ± 0.09, 1.61 ± 0.16, 1.61 ± 0.19, 1.60 ± 0.09, 1.60 ± 0.16, and 1.60 ± 0.09 nm, respectively. The Rg value of 1.61 nm for all hit compounds showed that the binding of the inhibitors does not cause significant stress on the backbone of DNA gyrase B. The Rg results indicated that all the complexes were compact throughout the simulation, and all the hit-complexes showed Rg values such as the reference compounds. Table S1 displays the averaged Rg values.

The ability of water accessibility at the active site of DNA gyrase B with a non-liganded protein and inhibitor–protein complexes was calculated using SASA calculation (Figure S3C,D). The SASA averaged values are shown in Table S1. The obtained results showed that the SASA for non-liganded protein was 97.05 ± 0.8 Å2, and SASAs for reference compounds clorobiocin and novobiocin complexes were 96.98 ± 0.7, 95.89 ± 0.8 Å2. Hit complexes showed SASA values 99.03 ± 0.7, 98.16 ± 0.7, 98.95 ± 0.9, 94.59 ± 0.6, 95.9 ± 0.6, and 94.03 ± 0.7 Å2 for compounds 1–6 complexes, respectively. The water accessibility of protein and ligands 4, 5, and 6 complexes was reduced slightly to around 94–95 Å2 which was lower than the SASA of non-liganded protein and the reference compounds. These findings revealed that fewer inner residues were in contact with the solvent, and this may be the reason for low SASAs values. It was noticed that compound 3 complex showed the highest Rg and SASA values among the six complexes, which may be due to compound 3 adjusted conformational change during the simulations, with a slight deviation from the binding pocket.

Hydrogen Bond and Distance Analyses

For investigating the binding affinity of hit compounds toward DNA gyrase B, the trajectories of the complexes were analyzed, and H bonds between the ligands and protein were calculated over the 100 ns simulation and plotted (Figure S4). The total H bonds formed between the docked reference compounds novobiocin were three but only two H-bonds were stable over the simulation with 99% occupancy (Table S2). Clorobiocin exhibited a maximum of two H-bonds (50%), while only one H bond was stable and continuous (Figure S4, Table S2). Compound 1 and 5 showed that the highest H number between the all hit complexes, they formed 3 H-bonds with the residues. The H bonds of compound 1 were continued until 80 ns of the simulation with an occupancy 48.2%. Moreover, the H bonds of compound 5 were steady and continued until the end of simulation with an occupancy 69.2%. Compound 2 formed 2 H-bonds for first 40 ns for 18% of the simulation. Compound 3 showed the least H bond occupancy of 8.8%, while compound 4 showed the most stable H bond with 48.2% occupancy. Compound 6 formed 2 H-bond with the receptor with 33.5% occupancy (Figure S4, Table S2).

Further investigation was carried out by measuring the center-of-mass distance between the hits and the pocket residues over the 100 ns MD simulation (Figure S5), the calculations were performed on the residues located within 4 Å of the inhibitors. The average distance of clorobiocin and novobiocin from the pocket were 0.22 ± 0.08 and 0.21 ± 0.05 nm, respectively. Compound 3 demonstrated a distance around 0.24 ± 0.09 nm. Compound 5 had the smallest average distance (0.18 ± 0.06 nm), which may be attributed to its excellent hydrogen network with the receptor and its high H-bonds occupancy (69.2%). Similar to the reference compounds, the other compounds displayed an average distance of 0.20–0.21 nm and smooth figures over the 100 ns simulation. Therefore, it was suggested that all hits, exhibit significant binding modes within the active site. Average distances are presented in Table S2.

Protein–Ligand Interactions

The interactions between the protein and the ligands were monitored along 100 ns simulations. These interactions were categorized into hydrophobic, salt bridge, and hydrogen bond interactions. The previous section discussed the number of H bonds, the distance between ligands, and key residues within the pocket. This section will discuss the residues incorporated in H bonds with the hits and the occupancy of interactions. Figure S6 shows a schematic of the ligands interactions with the pocket residues. Clorobiocin, novobiocin, and the six hits demonstrated hydrophobic interactions with Ile78, Pro79, Ile90, and Val120 which were maintained in more than 90% of the simulations. Additional hydrophobic residues incorporated in the protein–ligands interaction. For example, Val43 displayed hydrophobic interactions with novobiocin, compounds 2, 4, 5, and 6 in more than 40% of the simulation. Ala86, Ala96, and Met91 are significantly involved in the hydrophobic interactions, with more than 60% occupancy with the inhibitors.

Residue Asn46 formed a salt bridge with the reference compounds and inhibitors in more than 70% of the simulation. Unlikely, compound 5 showed only 20% occupancy. Other crucial amino acids such as Glu50, Asp73, Arg76, Gly77, Arg136, and Thr16 were involved in the salt bridge, with an occupancy of more than 70% (Figures S6 and S7)

Asn46, Asp73, and Ser85 were involved in H bonds with Novobiocin, with an occupancy of 18, 99 and 60%, respectively (Figures S6 and S7). Clorobiocin formed H bonds with Glu50 and Arg136 in 40 and 18% of the MD simulation, respectively (Figures S6 and S7). Asn46 was involved in H bonds with compound 1, 2, 4, and 5 with an occupancy 25, 18, 10, and 44%, respectively (Figures S6 and S7). Thr165 formed H bonds with compound 4 for 42% of the simulation, besides establishing H-bonds with compound 1 and 6. However, its occupancy was less than 15% (Figures S6 and S7). Compound 5 stabilized within the active site through hydrogen bonding with Asn46 and Asp49, with an occupancy of 40%. Moreover, compound 6 informed H bonds with Arg67, Gly77, and Thr165 with an occupancy of 20% (Figures S6 and S7)

MD analyses have demonstrated the following findings, the hits were stable within the active site expect compound 3 which was deviated from the pocket. The H-bonds analyses showed that compound 5 formed the highest number of H-bonds while compound 4 formed the most stable and steady H bonds, protein–ligand interactions indicated that the most hits were stabilized within the active site through hydrophobic interactions and salt bridge. However, the key residues involved in H-bonds in low occupancy.

Analysis of Binding Free Energy

The binding free energies for protein-reference and protein-hit complexes were calculated for the last 60 ns of MD trajectories. Six hits and references complexes showed negative binding free energies (Table ), indicating that all the complexes were stable. Clorobiocin demonstrated the highest free binding energy (−288.1 ± 1.6 kJ/mol), followed by Novobiocin (−133.4 ± 0.7 kJ/mol). Compound 3 showed the lowest free binding energy (−15.8 ± 2.1 kJ/mol). Compounds 1, 2, 4, 5, and 6 showed free binding energies of −111.2 ± 1.3, −121.1 ± 1.3, −122.7 ± 0.6, −131.7 ± 1.1, and −130.15 ± 0.6 kJ/mol, respectively.

Table 4

Energy Contribution of the 6 Hits, and Reference Compounds to the Total Free Energy (kJ/mol)

inhibitors	CPs no	ΔEvdw	ΔEele	ΔGpol	ΔGnp	ΔGbind
Clorobiocin		–206.1 ± 0.6	–225.5 ± 0.7	162.3 ± 1.6	–18.7 ± 0.1	–288.1 ± 1.6
Novobiocin		–206.8 ± 0.6	–46.4 ± 0.4	143.0 ± 0.7	–23.2 ± 0.1	–133.4 ± 0.7
ZINC9592782	CP 1	–163.8 ± 0.5	–17.2 ± 1.5	84.3 ± 0.9	–14.5 ± 0.1	–111.2 ± 1.3
ZINC8216235	CP 2	–188.0 ± 0.5	–18.3 ± 0.4	104.0 ± 0.6	–18.7 ± 0.1	–121.1 ± 1.3
ZINC23908246	CP 3	–79.9 ± 1.6	–1.0 ± 0.3	73.8 ± 1.6	–9.7 ± 0.1	–15.8 ± 2.1
ZINC32858011	CP 4	–204.5 ± 0.5	–13.1 ± 0.5	115.5 ± 0.9	–20.6 ± 0.2	–122.7 ± 0.6
ZINC20791579	CP 5	–234.2 ± 1.1	–263.1 ± 1.4	388.7 ± 1.7	–23.2 ± 0.1	–131.7 ± 1.1
ZINC5346121	CP 6	–204.4 ± 0.5	–48.2 ± 0.3	140.9 ± 0.5	–18.3 ± 0.2	–130.1 ± 0.6

For all complexes, the contributions of van der Waals (ΔEvdw), electrostatic (ΔEele), and non-polar solvation (ΔGnp) were attractive for protein–ligand interactions. However, the (ΔGnp) contributed with low values to the total free binding energy, due to the shielding of inhibitors from the solvent. The (ΔEele) contribution was negative for all the systems with low to medium values. Unlikely, (ΔEele) for clorobiocin and compound 5 showed a high contribution with (−225.5 ± 0.7 kJ/mol) and (−263.1 ± 1.4 kJ/mol), respectively. The polar free solvation energy (ΔGpol) was positive for all the systems, indicating the unfavorable bindings of inhibitors with the enzyme. Moreover, van der Waals interaction represented the main bonds for all inhibitors, while electrostatic interactions represented the main bonds for Clorobiocin and compound 5. Despite of the promising inhibition results of the most hit compounds, they demonstrated weaker free binding energy than the reference compounds.

Protein Inhibitors Affinity Relationship

To identify the hot spot residues involved the binding with the inhibitors, binding free energies of all systems expect compound 3 complex were decomposed into residues applying mmpbsa.py program. The obtained results are shown in Figure , left panel, whereas the last pose orientation of each inhibitor in MD simulation is presented in the right panel. The residues that demonstrated ΔGbind ≥ −1.0 kJ/mol were discussed. The finding was as follow: clorobiocin and novobiocin revealed nine hot-spot residues; Asp49, Glu50, Asp73, Ile78, Pro79, Ile90, Val120, The165, and Val167 (Figure A,B, respectively). The hydrophobic residue Ile78 contributed to all the systems with ΔGbind per residue ≥ −8.0 kJ/mol. Ile78 and Ile90 interacted significantly with compound 1, with an energy value of −8.2 kJ/mol, while Asp49 and Thr165 contributed to the total energy value of around −2.2 kJ/mol (Figure C). While Arg76 and Thr165 contributed to binding free energy of compound 2 (Figure D). Asn46 and Thr156 were common residues that contributed to binding free energy of compounds 4, 5, and 6 with energy values from −2.2 to −4.2 kJ/mol (Figure E–G respectively). Previous findings demonstrated that Asn46, Asp49, Arg76, and Thr165 interacted through H-bonds or salt bridge explaining the contribution these residues in the free binding energy. Hydrophobic residues Ile90, Val120, and Val176 interacted with all hits and contributed to the free binding energy (Figure ). Ala86 contributed to the free binding energy of compounds 4 and 5 with energy values −5.00 and −2.5 kJ/mol, respectively. Per residue decomposition showed that Asn46, Asp49, Glu50, Asp73, Ile78, Pro79, Ala86, Ile90, Val120, Thr165, and Val167 are the hot spot residues.

Figure 8

ΔGbind per residue of inhibitors-complex (left). The contributing residues involved in ligands binding in the last MD snapshot (right). Ligands shown as brown stick, protein represented as magenta cartoon, residues showed as cyan stick. (A) clorobiocin, (B) novobiocin, (C) compound 1, (D) compound 2. (E) Compound 4, (F) compound 5, (G) compound 6.

Conclusions and Perspectives

This study aimed to discover novel DNA gyrase B inhibitors using virtual screening techniques and biological methods. The virtual screening process resulted in six leads, demonstrate excellent binding score and interactions. Moreover, they demonstrated drug likeness properties and non-carcinogenic. The biological testing of their inhibitory effects on DNA gyrase B revealed that compound 4 showed a promising IC50 value 6.3 ± 0.1 μM which is better than clorobiocin and novobiocin. The molecular docking results showed that the six hits formed H-bonds with key residues and hydrophobic interactions with the hydrophobic residues. However, their binding modes revealed that compounds 1, 2, 4, and 6 occupied entirely the binding pocket, particularly the hydrophobic back pocket. Unlikely, compound 3 and 5 protruded outside the pocket. These findings are in agreement with biological results. The MD analyses revealed that the stability of the hits within the binding pocket. Compound 4 complex showed stability better than clorobiocin, with smooth rmsd pattern along the MD simulation. Also, the rmsf average of compound 4 was lower than the non-liganded protein and novobiocin while like clorobiocin. Consequently, the six hit complexes were found to have Rg like the reference compounds. Then, SASA demonstrated that compound 4 and 6 complexes have rmsf values lower than the non-liganded protein and references compounds. These results ensure that compound 4 stabilized the DNA gyrase B protein. Hydrogen bond and hydrophobic interaction analyses between the hits and the protein revealed the involvement of key residues, namely, Asn46, Asp49, Arg76, Gly77, Arg136, and Thr165 in H-bonds interactions. Furthermore, Val43, Ile78, Pro79, Ala86, Ile90, Met91, Ala96, Val120, and Val167 were involved in hydrophobic interactions. The binding free energy MM-PBSA revealed that the van der Waals interactions were the most critical force for the binding of all inhibitors, while electrostatic interactions contributed significantly for the binding of clorobiocin and compound 5. These findings provide an excellent example of employing 3D pharmacophore virtual screening as practical approach to discover novel anti-microbial agents. Finally, designing derivatives for compound 4, demonstrating more H bonds is required and minimum inhibitory concentration for the active compound is still needed.

Experimental Methods

Protein and Database Preparation

The crystal structure of DNA gyrase B in complex with clorobiocin with resolution 2.30 Å (PDB: 1KZN) was retrieved from Protein Data Bank (www.rcsb.org), the protein was protonated at pH 7.0 using the PDB2PQR server[29] (http://nbcr-222.ucsd.edu/pdb2.pqr2.0.0/) and the incomplete residues were fixed using Swiss-PDB Viewer v4.1.0[30] and UCSF Chimera.[31] The hydrogen atoms were added, water molecules were removed, and the energy was minimized with MOE software.

A data set of 300,000 compounds were downloaded from ZINC15 database and saved as mol2, the ligands were protonated, and the partial charges were assigned by MOE, the forcefield MMFF94x is used for parameterized the small molecules to set the partial charges of the ligands. These partial charges are dependent on the bond–charge increments. The energy was minimized by MMFF94x force field until rmsd gradient of 0.05 kcal/mol Å was obtained and the ligands converted from mol2 format to mdb format to be filtered by the pharmacophore model.

3D-Pharmacophore Modelling

The prepared protein (PDB code: 1KZN) was used for generating a 3D pharmacophore model utilizing the PLIF protocol in MOE 2008. The pharmacophore models were constructed upon the interactions of the native ligand clorobiocin and the protein. The obtained ones were validated by using a decoy set (200 compounds) collected from the Database of Useful Decoys (DUDe: http://dude.docking.org) and a set of known active compounds (25 compounds). The pharmacophore model showed a ROC 0.9 was selected to screen the prepared database. The compounds that passed the pharmacophore filtration were screened with Lipinski’s rule.

Molecular Docking

Molecular docking was carried out using MOE 2008, the prepared protein (PDB: 1KZN) was defined as a receptor, and the active site was determined using the center of the native ligand clorobiocin (X = 21.08, Y = 30.5, Z = 34.9). The docking protocol was performed using Dock workflow that divided into four stages. The first is a conformational analysis, in which DOCK generates conformations for each 3D single ligand conformation, followed by a systematic search involving all combinations of angles on a grid. The second stage is the placement phase, during which dock generates poses from ligand conformations using the triangle matcher method. The third phase was the scoring, which was conducted using the London dG scoring method. Finally, in the fourth refinement phase, the force field refinement scheme was applied, the energy of the system was minimized using MMFF94x forcefield, the charges of all atoms were assigned using the MMFF94x forcefield. The accuracy of the docking protocol was evaluated by docking the native ligand clorobiocin. The results revealed that the re-docked pose of clorobiocin had a rmsd value of 0.4 Å against the original pose of the co-crystalized ligand indicating the reliability of the docking protocol. The ligands and reference compounds (novobiocin and clorobiocin) were screened against the prepared protein (PDB: 1KZN) employing the reliable docking protocol. For each ligand, 10 docked poses were produced; the pose with the highest binding score was retained, then all different poses were ranked according to their binding scores. Compounds with a higher binding score than clorobiocin (native ligand) were selected. The compounds that demonstrated favorable interactions with the protein were shortlisted based on visual inspection. The shortlisted molecules were investigated for PAINS using the online PAINS filters at ZINC (docking.org).

Drug Likeness

The pharmacokinetics and drug-likeness of the selected hits were evaluated using the SwissADME server.[32] Carcinogenicity of hit compounds were predicted by the AdmetSAR 2.0.[33] The Lipinski rule was applied to predict the drug-likeness properties of the six hits; the rule of five stated that the compounds with MW ≤ 500 Da, H-bond donors ≤ 5, H-bond acceptors ≤ 10, log P ≤ 5 are likely to be a successful drug candidate.[23]

Inhibition Assay for DNA Gyrase

All six compounds were subjected to an in vitro enzyme inhibition assay, along with novobiocin and ciprofloxacin as positive controls. The assay was performed in the confirmatory diagnostic unit, Vacsera, Egypt. Compounds were profiled against DNA gyrase using the E. coli DNA gyrase kit.[34] The assay was conducted on streptavidin-coated 96-well microplates, and a wash buffer (20 mM Tris-HCl, pH 7.6), BSA, 0.01% (w/v), 137 mM NaCl, and Tween 20 0.05% (v/v) was used to hydrate the microplates. The biotinylated oligonucleotide was immobilized at room temperature for 5 min and the excess was washed with the buffer. The assay was carried out by adding 0.5 μg of relaxed pBR 322 as the substrate to 1 U DNA gyrase E coli in 30 μL reaction volume, followed by adding the tested compounds and positive reference compounds (novobiocin and ciprofloxacin) as 100× stock at a final DMSO concentration of 1%. The reaction mixtures were incubated at 37 °C for 30 min, and the TF buffer was added to the 96-well microplates and incubated at room temperature for 30 min. After washing, the unbound plasmid with TF buffer and terminating the enzyme reactions, the T10 buffer [1 mM EDTA and 10 mM Tris-HCl (pH 8)] was added and stained for 20 min. The fluorescence was then read and IC50 of the compounds was calculated using GraphPad Prism software.

MD simulations were performed for the docked hit pose complexes that showed high binding affinity scores using GROMACS 5.1.2 software.[35] The pdb2gmx program and Amber99-SB-ildn were used to generate the topology file of gyrase B protein and the forcefield, respectively. The ACPYPE[36] program was used to build the topology parameters of the six inhibitors, clorobiocin and novobiocin applying general Amber Force Field GAFF force field which is designed to be compatible with most organic molecules.[37] Parameterizations were first performed on bond angles and bond lengths and bond angles, then torsional angle parameters were parametrized.[38,39] All systems were solvated with a water model in a cubic periodic box with a minimum margin of 1 nm and the box was filled with TIP3P water.[40] The water molecules were replaced by chloride and sodium ions, resulting in a neutral (generally uncharged) system with an ionic strength of 0.15 M (physiological condition). After minimizing the energy in each system using steepest descent and conjugate gradients, the systems were heated gradually to 300 K under NVT condition for 100 ps and then equilibrated with NPT conditions for 100 ps. The particle mesh Ewald method was applied to define the long-range electrostatic forces, with a cut-off of 1.2 nm for Van der Waals and columbic interactions. The bonds and angles were kept constrained by using LINC algorithm.[41] Finally, MD were run for 100 ns on each system using md run program. The temperature and pressure of all systems were kept constant using the velocity-rescale algorithm and the Parrinello-Rahman algorithm during the simulation. The velocities, energies, and trajectories of all systems were updated at the time interval of 10 ps. All system trajectories were visualized using the visualization MD software package[42] and Pymol software.[43] The rmsd of Cα-atom deviations of the protein and ligand structures, and the rmsf of each residue and number of hydrogen bonds were computed using g rms, g rmsf, and gmx hbond tools. The protein compactness in all complexes was measured using radius of gyrations (Rg), and the electrostatic contributions of molecular solvation were calculated with the SASA using gmx gyrate and gmx sasa tools of GROMCS.

MM-PBSA Binding Free Energy Calculation

The binding free energy of all protein–ligand complexes was calculated by the MM-PBSA method[44] with the g_mmpbsa tool.[45] The binding free energy of each complex was computed using the script MmPbSaStat.py.[46] The snapshots were collected every 100 ps, starting from 40 ns giving 600 frames. The binding free energy estimated the electrostatic energy, van der Waals energy, polar solvation, and non-polar solvation energies based on the SASA model. The parameters which were used for calculating free energies were (inp = 2), cavity_surften = 0.037 and cavity_offset = −0.569, fillratio = 4, and scale = 2.0. The ionic strength (istrng) was 0.1 and the external and internal dielectric constants were (exdi = 80, indi = 1.0) respectively. The number of iterations was 1000 to perform of the linear PB equation (linit). The MmPbSaDecomp.py script was used to calculate the final contribution energy of each residue to the free binding energy of each complex.