3.7. Molecular Dynamics Simulation of Protein–Ligand Complexes

MD simulations of the protein–ligand complexes were performed using GROMACS 2018 [107]. Their protein topologies were generated using the CHARMM36 all-atom force field. The ligand topologies were generated using the CHARMM force field via the CHARMM General Force Field (CGenFF) server (available at https://cgenff.umaryland.edu/). Complexes were generated from the ligands and protein topologies for each of the selected cases under study. Each complex was solvated with the transferable intermolecular potential with a 3 points (TIP3P) water model in a cubic box of size 1.0 nm and neutralized with Na and Cl ions. Energy minimization of each complex was conducted for 50,000 steps using the steepest descent algorithm. The ligands were restrained before the NVT and further using the NPT ensemble. Equilibration of each complex was performed for 100 ps apiece and the final MD simulation was conducted for 10 ns with time steps of 2 fs under the PME. Extended MD runs (in duplicates) were also conducted for 100 ns for the unbound proteins and selected protein–ligand complexes of NANPDB2403 and ZINC95486008. Duplicate MD runs were carried out by generating random seeds for the initial velocities of each run.

The root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (Rg) of the unbound proteins and selected protein–ligand complexes were determined. RMSD is a frequently used measure of the differences between the structures sampled during the simulation and the reference structure [114]. MD simulations require systems to be close to their equilibrium (native) conformation. The time trajectory of RMSD shows how a protein structure deviates from a reference structure as a function of time [114].

RMSF measures the movement of a subset of atoms concerning the average structure over the entire simulation. RMSF indicates the flexibility of different regions of a protein, which can be related to crystallographic B factors [114]. Residues contributing to the complex structural fluctuation can be assessed using this stability profile analysis. Higher RMSF values imply greater fluctuations. Greater amounts of structural fluctuations occur in regions known to be involved in ligand binding and catalysis, notably the catalytic loop regions [115]. Adaptive variation in flexibility lies principally in these regions of the sequence that influence the conformational stabilities of the protein–ligand complex [115].

The radius of gyration (Rg) assesses the changes in compactness of a protein–ligand complex over the simulation time. The loss of compactness affects the stability of the complex by introducing weak intermolecular bonds. When the Rg of a complex is relatively steady, the compactness of the protein–ligand complex is high, and the protein is folded well, whereas the Rg value changes over time if the protein unfolds [45,116].

G_mmpbsa [48] was also used to calculate the free binding energies of each complex over the 10 ns simulation period, utilizing frames in steps of 0.1 ns. The binding free energy contribution per residue was calculated using the MM/PBSA and the plots generated using R programming.