3.5. Molecular Dynamic (MD) Simulations and Binding Free Energy

As a powerful tool to avail in silico, the stability of the above-repurposed systems, classical molecular dynamics (MD) simulations were carried out by using the PMEMD module [60] of the Amber20 package. Preliminary, all ligands were structurally optimized at the quantum mechanics (QM) level by applying the Hartree–Fock (HF) method with 6-31G** as implemented into the Gaussian09 program [61]. Afterwards, MM charges for each molecule were computed using the restrained electrostatic potential (RESP) method [62] by using the antechamber module of AmberTools20 [63]. The GAFF [64] and ff19SB [65] MM parameter sets were chosen for ligands and enzymes, respectively. Particularly, for 3MJ5 and 7BV2 systems, the Zinc AMBER force field (ZAFF) [66] was selected to simulate their metal centers. In addition, the SWISS-MODEL web server [67] was used to add missing amino acid residues into the 7BV2 system. The protonation state of the ionizable residues of all systems, at pH = 7, was calculated on PROPKA 3.1 version [68] and protons were added by the tleap module of AmberTools20 [63]. Each system was solvated by the TIP3P water model [69], extending 8 Å away from the solute atoms. Next, the appropriate counter-ions were added in order to neutralize the charges in the solvated system.

Each solvated system was energetically minimized to avoid bad atomic contacts using the PMEMD module. The minimization was carried out in four successive stages by applying 5000 steepest descent steps followed by 5000 steps of the conjugate gradient method, where restraints were removed during the process. After successful minimization, each system was slowly heated up to 300 K over 100 ps under constant volume, where the solute was restricted with the positional restraints of 50 Kcal/mol·Å2. Next, maintaining the same solute restraints, 200 ps of MD was performed at the NPT ensemble (p = 1 atm and T = 298 K). Then, the force constant of the restraints was slowly removed over the eight stages of equilibration; each stage was carried out for 100 ps under the NPT ensemble. Finally, a 200 ns unrestrained MD simulation (named “production stage”) was performed for each equilibrated system. The SHAKE algorithm [70] with the time step of 2 fs was applied for all the hydrogen bonds. The non-bonded cut-off set to 10 Å was used for the non-polar and polar interactions calculated using the Particle Mesh Ewald (PME) method [71]. The same MD protocol was used for all systems.

The structural analysis of all the simulated systems was evaluated by computing the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) of the backbone atoms (Cα, N, O, C). Particularly, the RMSD calculations for each ligand were computed by considering only their respective heavy atoms. The trajectories were aligned by the main-chain atoms of the average structures from production stages by using the CPPTRAJ module [72].

The CPPTRAJ module [72] was used to extract 10 ns (a total of 1000 representative snapshots) from the production stage of the MD simulations on each system to be selected for the binding free energy calculations using the MM/GBSA approach [46,47,48], which was implemented into the MMPBSA.py module [49] of AmberTools20. The main equations of the chosen approach can be summarized as follows:

where ΔGcomplex, ΔGreceptor and ΔGligand indicate the free energies of the complex, the receptor, and the ligand, respectively (Equation (1)). The ΔGbind is obtained from the gas-phase MM energy (ΔEMM), solvation energy (ΔGSOLV) and the entropic term (−TΔS) (Equation (2)). The ΔEMM includes the changes in the internal (bond, angles, and dihedral energies) (ΔEint), electrostatic (ΔEelec) and van der Waals (ΔEvdw) contributions (Equation (3)). Here, as a single-trajectory scheme was used for the binding free energy calculations, the ΔEint is equal to zero. The ΔGSOLV is the sum of the polar (ΔGGB) and non-polar (ΔGSA) energies for ΔGbind (Equation (4)). To reduce the computational cost, the entropic term (−TΔS) was not computed into the binding free energy calculations [47,48]. Furthermore, a per-residual decomposition analysis was included to avail the relative contribution of each amino acid residue [46]. This method has been frequently applied as a useful tool in SARS-CoV-2 drug design studies [35,40,41,73,74,75,76].