We have chosen S. aureus (Sa) and P. aeruginosa (Pa) as gram-positive and gram-negative bacterial representatives for computational analyses of DnaK protein binding. 3D structures of DnaK proteins, from the aforesaid species, were generated via homology modelling using MODELLER version 9.2481. DnaK has two conformations, namely, the open or ATP-bound and the closed or ADP-bound conformation65. Herein, we focused on the open conformation of DnaK, to identify potential competitive inhibitors of ATP to prevent the proper functioning of DnaK protein.

We have obtained the protein sequences of Sa and Pa DnaK from UniProtKB with accession IDs of Q2FXZ2 and A6VCL882, respectively. To search for suitable homology modelling templates, we utilized both NCBI BLASTp and the MODELLER in-built build_profile.py81,83. For Pa DnaK (PaD), the templates were full-length ATP-bound E. coli (Ec) DnaK protein structures (PDB ID: 5NRO, Chain: A, Query Coverage (QC): 94%, Percent Identity (PI): 79.50%, Resolution (R): 3.25 Å; PDB ID: 4JNE, Chain: A, QC: 94%, PI: 78.80%, R: 1.96 Å; and PDB ID: 4B9Q, Chain: A, QC: 94%, PI: 77.96%, R: 2.40 Å). For Sa DnaK (SaD), besides the afore-mentioned Ec DnaK (EcD) models, we selected one additional template, from Geobacillus kaustophilus DnaK protein (PDB ID: 2V7Y, Chain: A), due to the high percentage of sequence identity expected as per the gram-positive character of S. aureus and G. kaustophilus. As this template structure was in the closed conformation and we were only interested in the open conformation, only the Nucleotide Binding Domain (NBD, residues 1 to 350 in template model) which does not differ much in both conformations, were taken into consideration for homology modelling, and the remaining C-terminal residues modelling were guided by the Ec models to shape an open conformation. Therefore, the templates for SaD were (PDB ID: 2V7Y, Chain: A, Template Residues: 1–350, QC: 57%, PI: 83.19%, R: 2.37 Å; PDB ID: 5NRO, Chain: A, QC: 93%, PI: 56.19%, R: 3.25 Å; PDB ID: 4JNE, Chain: A, QC: 92%, PI: 55.54%, R: 1.96 Å; and PDB ID: 4B9Q, Chain: A, QC: 94%, PI: 55.43%, R: 2.40 Å). The template sequences were aligned with the target sequence for homology modelling via the built-in function of MODELLER (Fig. S6A). 5 homology models were generated for each protein of SaD and PaD, and the models with the lowest DOPE (discrete optimized protein energy) scores were selected for downstream virtual screening for both. We then validated the SaD and PaD homology models via Swiss-Model Structure Assessment and SAVES v5.0 servers84 (Fig. S6).

To validate the druggability of the ATP docking pocket, we have conducted ligand binding site prediction using P2Rank from PrankWeb server85. P2Rank predicts the chemical druggability on protein solvent-accessible surface via a non-templated machine learning approach. The ATP binding pocket was predicted to be druggable and ranked first in both cases of SaD and PaD (Table S9; Fig. S7). Thus, we further considered these pockets from the SaD and PaD complexes to be targeted for virtual screening.

We utilized the POAP pipeline86 for an in silico virtual screening of the chemical compounds obtained through different chromatographic separation. We have obtained the SMILES notations of these compounds, and generated their 3D models (in mol2 format) through the POAP Ligand Preparation pipeline. To this end, we utilized Chimera to generate physiological protonation states of the ligands, and PDBQT files were prepared87. We also carried out ligand optimizations via the POAP Ligand Preparation pipeline utilizing the MMFF94 force field which is being optimized for drug-like organic molecules and molecular docking88. Out of the 50 conformers, generated for each ligand through the Weighted Rotor Search approach, only the best conformers were retained. Finally, we have subjected the ligands to energy minimization for 5000 steps by the conjugate algorithm.

We have prepared the macromolecule receptors, of the SaD and PaD proteins, using AutoDockTools. We utilized AutoDock 4.2, aided by the POAP pipeline, for the virtual screening process89. For AutoDock parameters, we have set 100 generations of Lamarckian Genetic Algorithm for each protein–ligand complex. To fit in the previously predicted pocket, we adjusted the docking grids into squares of 24 Å with x, y, z coordinates of 17.647, 75.43, 27.766, and 18.069, 74.299, 28.532, for SaD and PaD, respectively. For the silicon-containing compound among the set of ligands, we separately carried out molecular docking with AD4.1_bound parameter file, obtained from AutoDock, wherein we added the parameters for silicon atoms (Rii = 4.3; eii = 0.402)90.

We have validated the docking methodology via redocking of experimentally confirmed and deposited structures with the reference ligand (Fig. S8). Therein, we have retrieved the ATP-bound E. coli DnaK crystallized structures from PDB (PDB ID: 4B9Q, CHAIN ID: A; PDB ID: 4JNE, CHAIN ID: A). The molecular docking search grids were squares of 24 Å, in x, y, z coordinates of 108.958, 73.922, 100.622, and 19.081, 76.887, 31.16, for 4B9Q and 4JNE respectively.

Using SwissADME91, we have carried out predictions on the pharmacological properties, encompassing pharmacokinetics, drug-likeness, and molecular information, for each chemical compound.

Ensuing virtual and pharmacological screenings, we rationally selected potential drug candidates to undergo molecular dynamics (MD) simulation via GROMACS version 2019.334. We have utilized the CHARMM36 force field of version July 2020, along with the TIP3P water model, for macromolecule processing92. We used Avogadro software for mol2 format conversion and complete protonation (protonation of non-polar atoms)93. We also used a Perl script, sort_mol2_bonds.pl, written by Justin Lemkul for bond order arrangements in ligand mol2 files. Then, we generated the topologies of the ligand models through the CGenFF server, and utilized a python script (cgenff_charmm2gmx.py) to convert topologies for CHARMM to GROMACS94. We carried out solvation in a dodecahedron box ranged 1.0 Å from the protein–ligand complex. The system was then ionized to achieve electrostatic neutralization. Subsequently, we subjected the system to energy minimization via the steepest descent algorithm until convergence at a maximum force of less than 1000 kJ mol−1 nm−1 (Fig. S9). Herein, we have monitored the potential energy shifts of the systems.

We carried out equilibration of the systems via NVT and NPT ensembles for 50,000 steps (100 ps), with temperature, pressure, and density shifts being monitored therein. Subsequently, we have carried out the production MD simulations for 5,000,000 steps (10 ns) to observe protein–ligand interactions. We computed the RMSD (Root Mean Square Deviation) values of ligands and receptors, number of hydrogen bonds between ligands and receptors, and ligand-receptor interaction energies (Coulombic interaction energies and Lennard–Jones energies) throughout the MD simulations. We have also computed thetotal interaction energies, and estimated the errorsvia error propagation by addition.

We utilized Matplotlib, a python library, to tabulate binding energies of all screened compounds31. We generated all 3D structural images using UCSF ChimeraX32, and 2D chemical structures using MarvinSketch33. Finally, we retrieved the figures for MD simulation analyses from the GROMACS in-built functions34.

Note: The content above has been extracted from a research article, so it may not display correctly.

Please log in to submit your questions online.
Your question will be posted on the Bio-101 website. We will send your questions to the authors of this protocol and Bio-protocol community members who are experienced with this method. you will be informed using the email address associated with your Bio-protocol account.

We use cookies on this site to enhance your user experience. By using our website, you are agreeing to allow the storage of cookies on your computer.