Incorporating Target Protein Structure Flexibility and Dynamics in Computational Drug Discovery Using Ensemble-Based Docking Analysis

Summary

Computational methods hold promises for expediting drug discovery, yet they frequently overlook the dynamic nature of protein structures. Here, we discuss ensemble-based docking analysis to indirectly incorporate protein flexibility, potentially improving the accuracy and reliability of drug discovery efforts.

Abstract

The drug discovery process is a rigorous, time-consuming, and expensive operation. The computational approach in drug discovery allows researchers to prioritize the most promising compounds for further testing, which would greatly reduce the required resources, leading to an increment of the overall efficiency in the drug discovery pipelines. Structure-based drug discovery is a common approach that requires the structural information of the target protein in a three-dimensional format. However, the current limitation of most computer-aided drug discovery strategies is their inability to introduce the flexibility and dynamics of the target protein structure during the ligand-protein docking simulation. While both induced fit docking and ensemble-based docking aim to address protein flexibility in the docking procedure, the latter can provide a more comprehensive view of dynamic protein behavior by incorporating multiple conformations throughout the simulation. In this report, we demonstrate and discuss the application of a technique called ensemble-based docking analysis that indirectly introduces the target protein structure's flexibility and dynamics in the molecular docking process. The protein and ligand selected for ensemble-based docking studies were lysozyme and Flovokawain B (FB), respectively. FB has been previously reported to have binding activity with lysozyme. A molecular dynamics (MD) simulation was performed on lysozyme in the presence of water, and the total energy, root-mean-square deviation (RMSD), and root-mean-square fluctuation (RMSF) were examined. Conformation clustering was generated based on several clustering cutoff values and was chosen for additional docking analysis with FB. Cluster no 2 gives the lowest binding energy at -29.37 kJ/mol. Molecular docking images were generated to anticipate the presence of binding forces. By incorporating the structural dynamics of the protein, the ensemble-based docking approach can better capture the range of possible binding scenarios, leading to more reliable predictions of binding outcomes.

Introduction

In computational drug discovery (CDD), techniques from computer science, chemistry, biology, and physics are integrated to explore vast chemical space, predict drug-target interactions, and optimize drug candidates with higher efficiency and lower costs compared to traditional experimental methods alone. It is a powerful approach that leverages computational methods and algorithms to accelerate the discovery and optimization of new therapeutic compounds^1,2,3. CDD has revolutionized drug discovery. However, there are limitations associated with protein three-dimensional (3D) structure dynamics that can affect the accuracy and reliability of computational predictions⁴. Protein 3D structures serve as templates in CDD for the design or optimize drug candidates based on the target protein-drug interactions. While X-ray crystallography models of protein structures provide valuable structural information about protein, it is essential to recognize the dynamic nature of protein structures and the limitations of static models^5,6,7. In addition, recent advancements in cryo-electron microscopy (cryo-EM) and computational predictions such as AlphaFold have also greatly expanded the availability of structural data in capturing the full spectrum of protein flexibility and dynamics^8,9,10,11,12.

Molecular dynamics (MD) simulations simulate the movement and interactions of atoms and molecules over time, providing insights into the dynamic behavior and flexibility of protein 3D structure^13,14. MD simulations are used to generate 3D structures of protein representing various conformational states, which serve as input for ensemble-based docking analysis. Through the sampling of diverse protein conformations, ensemble-based docking analysis accounts for the inherent flexibility and dynamics of biological targets, allowing for a more comprehensive exploration of ligand binding modes and interactions³.

Understanding protein flexibility is essential, as it influences how drugs exert their biological effects, determines the location and orientation of binding sites, and affects binding kinetics, metabolism, and transport^15,16. Capturing this dynamic nature can significantly enhance the accuracy and reliability of docking predictions. In 1994, Kearsley et al. introduced a flexible docking technique, a framework that models the flexibility of both ligands and proteins. This approach allows the protein's conformation to adjust during docking, improving predictions of ligand-receptor interactions by accounting for structural flexibility¹⁷. Similarly, in 1999, Carlson et al. reported on ensemble docking, which applies flexible pharmacophore modeling to both static and dynamic models of HIV-1 integrase, further highlighting the importance of accounting for protein dynamics in docking studies¹⁸. Furthermore, Cavasotto et al. also reported the improvement of ligand docking accuracy by incorporating receptor flexibility into the docking process using normal mode analysis¹⁹. More recently, ensemble-based techniques^{20,21,22,23,24,25,26,27}have expanded drug discovery by identifying potential new ligand binding sites and providing more accurate estimates of free ligand-receptor binding energy. These advancements have been applied to targets such as quadruplex-duplex DNA²⁰, vascular endothelial growth factor 165 (VEGF-165)²¹, the SARS-CoV-2 target enzyme²², human liver cytochrome P450 enzymes²³, and anticancer proteins²⁴.

Flavokawain B (FB), classified as a flavonoid, has been documented to exhibit various pharmacological properties^28,29,30. Based on experimental and computational analysis, FB has been reported to form a stable complex with lysozyme (LYZ)³¹, a protein widely recognized for its antimicrobial activity and has also been identified as a ligand transporter^32,33,34. In this report, we further analyze the nature of the interaction of FB with LYZ using ensemble-based docking analysis to incorporate the impact of protein flexibility in FB-LYZ complex formation. The goal of this method is to provide researchers with a step-by-step, repeatable process for ensemble-based docking analysis. In addition, it is advisable for researchers to select the protein structures from the target organism for research.

Protocol

1. Flavokawain B ligand structure preparation

1. Open PubChem website. Search for Flavokawain B. Select Download and select Save for SDF 2D structure as Structure2D_COMPOUND_CID_5356121.MDL SDF.
2. Open Avogadro software. Move the SDF file into the empty space of the software. Click Yes.
   NOTE: Avogadro software will automatically build a three-dimensional (3D) geometry.
3. On the tab, click the icon, as shown in Figure 1. On the left menu tab, change the force field to MMFF94^35,36. Adjust the steps per update to 15. Set the algorithm to Steepest Descent and click Start.
   NOTE: Do not move any atoms while the rotation is in progress. Once finished, click Stop. MMFF94 is suitable for a broad range of small-molecule ligands (including both neutral and charged) and is widely used for energy minimization and molecular modeling³⁷. The Steepest Descent method is chosen for initial energy minimization because it is straightforward, computationally efficient, and effective for removing large structural strains and achieving a basic minimization state^38,39,40.
4. On the menu tab, go to File and click Save As. When a Save Molecule As folder appears, write the file name as ligand.pdb and click Save.

2. Lysozyme protein structure preparation

1. Open the RCSB protein data bank website. Search for lysozyme, code: 1LYZ. Click Download Files and select PDB Format. Save as 1lyz.pdb.
2. Open Chimera software. On the tab, click File and Open. In the folder, select 1lyz.pdb file.  
   NOTE: A 3D structure of lysozyme will be constructed using software.
3. On the tab, click Select > Residue > HOH. Then click Actions > Atom/Bonds > Delete. Click Select and Clear Selection.
   NOTE: This will remove any water from the protein.
4. On the tab, click Select > Chain > A. Then click Tools > Structure Editing > AddH. A parameter box will come out and leave it as a default. Click OK. Click Select and Clear Selection (Figure 2).
   NOTE: This will add hydrogen to the protein. The histidine residue is also protonated to create a neutral system. The software can perform pKa calculations by integrating PROPKA^41,42.
5. On the tab, then click Tools, Structure Editing, and Add Charge. A parameter box will come out; select Gasteiger. Click OK. NOTE: This will add partial charges to the protein.
6. To save the file, in the tab click File > Save PDB. Save as protein.pdb.
   NOTE: In this study, the structure of lysozyme was chosen to compare with previously reported result³¹. Readers are advised to carefully select protein structures from the target organism of their interest.

3. MD simulation of lysozyme in water

NOTE: The custom-built computers for molecular dynamics (MD) simulation used in this study are Intel CORE i7 11^th Gen for CPU, NVIDIA Geforce RTX 2060 for GPU, and DDR4 128 GB memory. The OS system is Ubuntu 22.04.4 LTS. The MD simulation is using GROMACS software. The GPU support is CUDA.

1. Download all required documents from this link: http://www.mdtutorials.com/gmx/lysozyme/01_pdb2gmx.html (GROMACS Tutorial, Lysozyme in Water)⁴³.
   NOTE: The documents needed are em.mdp, ions.mdp, md.mdp. npt.mdp, nvt.mdp, charmm36 ff and protein.pdb (from step 2.5). For system setup, em.mdp is for energy minimization, and ions.mdp is to add ions and neutralize the system. For equilibration, nvt.mdp is for temperature equilibration, while npt.mdp is for pressure equilibration. For production, md.mdp is for simulation. Charmm36 ff is for topology.
2. Right-click the empty space in the folder (the working directory) and click Open Terminal. Type gmx to open GROMACS software. NOTE: GROMACS does not have a graphical user interface; it is all from the command written at the terminal.
3. For the generation of protein topology, type gmx pdb2gmx -f protein.pdb -o protein.gro –ignh. A list of force fields is given, and type 1 for CHARMM all-atom force field⁴⁴, followed by type 1 for TIP3P for the water model.
   NOTE: The updated version, CHARMM36, is also available and can be downloaded from http://mackerell.umaryland.edu/charmm_ff.shtml⁴⁵.
4. Define a cubic box to cover the entire protein structure. Type gmx editconf -f complex.gro -o newbox.gro -c -d 1.0 -bt cubic. NOTE: The protein is placed at least 1 nm from the box edge. Ensure the water box size is large enough for protein, solvent, and ions. In addition, the water box should be adequate to prevent interactions between periodic images of the protein. The recommended minimum distance between the protein and box edge is at least 1.0–1.5 nm.
5. For solvent configuration⁴⁶, type gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro.
6. Add ions such as sodium and/or chloride by typing gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro.
   NOTE: Sodium and/or chloride ions are/is used for neutralizing the system.
7. Neutralize the system by typing gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral. Select 13 for the SOL group.
   NOTE: For lysozyme, 8 solute molecules have been replaced with chloride ions (Figure 3).
8. In the CHARMM 36 file, search for ions.itp folder. The abbreviation for chloride ion is CLA.
9. Open topol.top file, add CLA and number 8 as shown in Figure 4.
10. Open solv_ions.gro file and exchange all CL into CLA.
    NOTE: Make sure the solv_ions.gro and topol.top files have the same abbreviation of respective ions, or else a warning will come out. Renaming the ion name is in accordance with the CHARMM 36 ff file version. Since the abbreviation for chloride is already CL, there's no need to rename the ions in the older version. Nevertheless, always check for the ions abbreviation.
11. Relax the protein structure by steepest descent energy minimization with the maximum number of 50,000 steps. For this, type gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr followed by gmx mdrun -v -deffnm em.
    NOTE: The solvent and ions are equilibrated around the protein at two phases: (i) number of particles (N), system volume (V) and temperature (T) are constant (NVT), and (ii) number of particles (N), system pressure (P) and temperature (T) are constant (NPT). Temperature coupling (NVT phase) is set at 300 K based on a modified Berendsen thermostat, while pressure coupling (NPT phase) is set at 1 bar based on Parrinello-Rahman. The long-range electrostatic interactions and the nearest-neighbor search are automatically calculated using the particle-mesh Ewald (PME) method and Verlet algorithm, respectively.
12. For equilibrium of NVT, type gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr followed by gmx mdrun -v -deffnm nvt.
    NOTE: The time taken for temperature equilibration is around 5 min, depending on the CPU and GPU. For system equilibrium, the simulation time is at 100 ps.
13. For equilibrium of NPT, type gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr followed by gmx mdrun -v -deffnm npt.
    NOTE: The time taken for pressure equilibration is around 5 min, depending on the CPU and GPU. For system equilibrium, the simulation time is at 100 ps.
14. Have the system go through equilibration runs for 1 ns followed by the production run using the mdrun function with a duration of 100 ns. Open md.mdp file. At n steps, change to 50000000; 2 ´ 50000000 = 100000 ps (100 ns).
15. For the production of molecular dynamics, type gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr followed by gmx mdrun -v -deffnm md.
    NOTE: The time for this is around ~2 days. The trajectory frames generated from the MD simulation is captured at 10 ps intervals and utilized for root mean square deviation (RMSD)-based clustering analysis.

4. RMSD-based clustering analysis

1. After the MD simulation, continue with the command analysis. Protein will diffuse through the unit cell and may appear "broken" or "jump" across to the other side of the box. Type gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -center -pbc mol -ur compact. Type 1 for centering the protein and 0 for the system's output. Type vmd em.gro to visualize the protein.
2. For total energy analysis, type gmx energy -f md.edr -o totalenergy.xvg, and continue to type 14 for total energy.
   NOTE: Monitoring total energy is fundamental for ensuring the reliability and correctness of MD simulations^47,48.
3. For RMSD analysis, type gmx rms -s md.tpr -f md_center.xtc -o rmsd.xvg -tu ns. Continue to type 3 for the C-alpha of the least squares fit and 3 for the C-alpha of the RMSD calculation (Figure 5).
4. For RMSF analysis, type gmx rmsf -s md.tpr -f md_center.xtc -o rmsf.xvg –res. Continue to type 1 for protein.
5. For using grace, type this command xmgrace totalenergy.xvg. Adjust the axes by double-clicking on the graph box line. A Grace: Axes box will appear. Click Accept.
6. On the tab menu, click File and Print_setup. A Grace: Device setup will appear. Change the device from PostScript to JPEG. Click Accept. Now click Print on the file menu. Save as totalenergy.jpg.
7. Repeat steps 4.5 and 4.6 for xmgrace rmsd.xvg and xmgrace rmsfxvg.
8. For cluster analysis, type gmx cluster –s md.tpr –f md_center.xtc –g cluster.log –sz cluster-size.xvg –clid clus-id.xvg –cl cluster.pdb –cutoff 1.0. Continue to type 1 (protein group) to calculate the least squares fit and RMSD and type 1 (protein group) for system output.
   NOTE: Adjust the cutoff value according to the clustering result. Command -cl represents outputs of average for each cluster.
9. Open cluster-size.xvg.
   1. Based on the information provided, increase the RMSD cutoff value if the number of clusters is low or decrease the RMSD cutoff value if the number of clusters is high.
10. Repeat step 4.8 with different cutoff values.
    NOTE: Clustering with different cutoff values till an optimal RMSD cutoff value has been determined based on the following criteria: (1) The total number of clusters should be limited to fewer than 30. (2) Minimizing the presence of clusters with only one member is preferred. (3) It is desirable for over 90% of the trajectory to be represented in less than 10 clusters.
11. Open Chimera software and search for cluster.pdb.
    NOTE: cluster.pdb contains an average for each cluster group.
12. Click Presents and Publication 1 (silhouette, rounded ribbon).
13. Then, go to File > Save Image > Save.
14. Click Select > Chain > (no ID) > cluster.pdb (#1). Click Select > Invert (all models). Click Actions > Atoms/Bonds > delete.
    NOTE: This will remove all group clusters except cluster 1.
15. Then, go to File > save PDB > save. Save as cluster1.pdb.
16. Repeat steps 4.11–4.15 for different clusters. Save as cluster2.pdb, cluster3.pdb and cluster4.pdb.
    NOTE: Since the total number of the top 4 clusters is more than 90% of the total trajectories, each representative of the top 4 clusters is extracted and subjected to molecular docking analysis using Chimera (See Figure 6).

5. Ensemble-based docking

1. Double-click Autodock Tools software.
   NOTE: For docking, this study used AutoDock and AutoDock tools software^13,49,50,51.
2. Place files cluster1.pdb and ligand.pdb in a new folder.
3. On the menu, click File > Preferences > Set. A Set User Preferences box will pop out. Copy the address of the "new folder" as a text. Paste the address at the Startup Directory in the Set User Preferences box. Click Set.
   NOTE: This is an important step if using Windows OS.
4. Click the blue folder image. A read molecule folder will pop out. Select cluster1.pdb.
   NOTE: Autodock tools will read the protein molecular structure.
5. Click Edit > Charges > Add Kollman Charges. Then click OK. Click Edit > Hydrogens > Merge Non-Polar.
   NOTE: Kollman charges are added to the protein.
6. Click Grid > Macromolecules > Choose. A Choose Macromolecules box will pop out. Select cluster1 and click Select Molecules. Click OK. A Modified AutoDock4 Macromolecule file will appear. Save as cluster1.pdbqt.
7. Empty the workspace by clicking Edit > Delete > Delete All Molecules. Click Continue.
8. Click Ligand > input > Open.  A Ligand file for Autodock4 folder will appear. Choose All Files, select ligand.pdb, and click Open. Click OK.
   NOTE: The setup ligand will include the incorporation of Gasteiger charges and the merging of non-polar hydrogen.
9. Click Ligand > Torsion Tree > Detect root.
10. Click Ligand > Output > Save as PDBQT. A Formatted Autotors Molecules folder will appear. Save as ligand.pdbqt.
11. Empty the workspace by clicking Edit > Delete > Delete All Molecules. Click Continue.
    NOTE: Follow the same procedure as described in step 5.6.
12. Click Grid > Macromolecules > Open. Select cluster1.pdbqt and click Open. Click Yes. Click Ok.
13. Click Grid > Set Map Types > Open Ligand. Select ligand.pdbqt and click Open.
    NOTE: In the workspace, there are proteins and ligands.
14. Click Grid and Grid Box. A Grid Options box will appear. At Number of points in x-dimension, adjust the parameter to 120, and set Number of points in y-dimension to 120, set Number of points in z-dimension to 120. Set spacing (angstrom) to 0.375. Leave Center Grid Box parameters as default. Click File and Close Saving Current.
    NOTE: The grid box is covering the whole protein, which means it's a blind docking.
15. Click Grid > Output > Save GPF. A Grid Parameter Output file will appear. In the filename, type grid.gpf and click Save.
16. Click Run and Run AutoGrid. A Run AutoGrid box will appear. At Parameter Filename tab, click Browse. An Autogrid Parameter file will appear. Select grid.gpf. Click Open. At Program Pathname > Browse. An  autogrid4 file will appear. Search for autogrid4.exe and click Open > Launch.
    NOTE: Ensure that the folder name does not contain any spaces, as this may result in an error during running. The folder autogrid4.exe can be installed from https://autodock.scripps.edu/download-autodock4/.
17. Click Docking > Macromolecules > Set Rigid Filenames. A PDBQT Macromolecules file will appear. Select cluster1.pdbqt and click Open.
18. Click Docking > Ligand > Choose. A Choose Ligands box appears. Select the ligand and click Select Ligand. An AutoDpf4 Ligand Parameter box will appear. Click Accept.
19. Click Docking > Search Parameter > Genetic Algorithm. A Genetic Algorithm Parameters box will appear. Change the number of GA Runs to 100. Leave the remaining parameters as default. Click Accept.
20. Click Docking > Output > Lamarckian GA(4.2). An Autodock4.2 GALS Docking Parameter Output file will appear. For the file name, type docking.dpf and click Save.
21. Click Run > Run AutoDock. A Run Autodock box will appear. At Parameter Filename click Browse. A autodock4 Parameter file appears. Select docking.dpf. Click Open. At Program Pathname click Browse. An autodock4 file will appear. Search for autodock4.exe and click Open. Click Launch.
    NOTE: Ensure that the folder name contains no spaces, as this may result in an error during running. Time taken is ~30 min. The folder autodock4.exe can be installed using the link: https://autodock.scripps.edu/download-autodock4/.
22. Delete all molecules as described in step 5.7.
23. Repeat steps from 5.1–5.22 for cluster2.pdb, cluster3.pdb and cluster4.pdb.
    NOTE: All 4 docking groups will be used for analysis.

6. Ensemble-based docking analysis

1. Use AutoDock Tools software to continue analysis. Click Analyze > Docking > Open. A docking log file will appear. Select docking.dlg and click Open. Then, click OK.
2. Click Analyze > Macromolecules > Open.
3. Click Analyze > Conformations > Play, ranked by energy. A ligand box will appear.
4. At the new folder, open docking.dlg folder using a notepad. Search for cluster analysis of conformation. From the information given, search the conformation run (stated as Run), which has the lowest binding energy.
   NOTE: Since there are 100 conformations run, only one has the strongest binding affinity between ligand and protein. The term for lowest binding energy is defined as having a stronger binding affinity^52,53,54.
5. From step 6.3 (ligand box), input the conformation run (from step 6.4), and press Enter.
   NOTE: The ligand will position itselft accordingly inside the protein.
6. At the ligand box, click the button open panel to change play option. A Set Play Options box will appear. Click Write Complex. A write complex of receptor folder will appear. Save as complex.pdb.
   NOTE: Since there are 4 cluster groups, complex.pdb file can be named accordingly, such as complex1.pdb for cluster group 1.
7. Open Chimera software and search for complex.pdb (as in step 6.6).
8. Click Presets > Interactive 1 (ribbons) > Publication 1 (silhouette, rounded ribbon).
9. Then go to File and click Save Image.
   NOTE: The image saved is the ribbon structure of the FB-LYZ complex.
10. Click Presets > Interactive 3 (hydrophobicity surface) > Publication 1 (silhouette, rounded ribbon).
    NOTE: The image saved is the ballon structure of the FB-LYZ complex.
11. Open Discovery Studio software.
12. Drag the complex.pdb into the empty space of the software.
13. Click on  the Tools tab and select Show 2D Diagram.
    NOTE: The software will automatically generate the 2D structure of the complex (Figure 7).
14. Click File and Save as. Write the file name and click Save.

Results

The chemical structure and the 3D structural representation of FB after optimization are shown in Figure 8A. Figure 8B shows the 3D structure of lysozyme with pdb code 1LYZ at the initial state prior to MD simulation. In order to study the dynamics and flexibility of lysozyme 3D structure, an MD simulation for 100 ns was conducted. The total energy of the protein structure was stable during the simulation, as shown in

Discussion

Ensemble-based application in computational drug discovery involves utilizing numerous conformational ensembles derived from crystal structures, nuclear magnetic resonance (NMR) studies, or molecular dynamics simulations²⁶. The exploitation of multiple computational conformations would allow the analysis to incorporate the dynamics and flexibility of a protein structure, leading to the enhancement of accuracy and reliability in identifying potential drug candidates.

To ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
AutoDock	The Scripps Research Institute, US	Version 4.2.6
AutoDock Tools	The Scripps Research Institute, US	Version 1.5.6
Avogadro	Geoffrey R Hutchison, Department of Chemistry, University of Pittsburgh, Pittsburgh, USA	Version 1.95
Discovery Studio	Dassault Systèmes, Massachusetts, US	Version 2021
GROMACs	University of Groningen Royal Institute of Technology Uppsala University, Sweden	Version 2023
UCSF Chimera	Resource for Biocomputing, Visualization, and Informatics University of California	Version 1.16