Molecular dynamics (MD) simulations and MM/PBSA analysis

Molecular dynamics simulation of the WT and mutant protein structures was performed using GROMACS 5.1.4 version⁷⁷ and Linux 5.4 package. The GROMOS96 54a7⁷⁸ forcefield was selected as the force field for proteins and the ligand topologies were generated from the Automated Topology Builder version 3.0⁷⁹ (ATB) server. Due to the enhanced capacity of the backbone NH and CO groups to form hydrogen bonds with each other in the GROMOS96 54A7 parameter set, this force field reproduces the folding equilibria slightly better and can sample more 314-helical or hairpin conformations than the previous 53A6 or 45A3 force fields⁸⁰. Also, based on fitting to a large set of high-resolution crystal structures, the torsional angle terms were reparametrized in this parameter set⁸¹. The proteins and mutants-ligand complexes were solvated using simple point charge (SPC) water molecules in a rectangular box where every structure was placed in the center at least 1.0 nm from the box edges. Required number of Na+ and Cl− ions were added to make the simulation system electrically neutral. The salt concentrations were set to 0.15 mol/L in all the systems. The solvated systems were subjected to energy minimization for 5000 steps using the steepest descent method. Afterwards, three steps were conducted in the MD simulation: NVT (constant number of particles, volume, and temperature) series, NPT (constant number of particles, pressure, and temperature) series, and the production run. The NVT and the NPT series were conducted at a 300 K temperature and 1 atm pressure for the duration of 100 ps. V‐rescale and Parrinello‐Rahman were selected as the thermostat and barostat respectively of the performed simulation. Finally, 3 independent production runs of nine proteins (WT and mutants) and six protein–ligand complexes were performed at 300 K for a duration of 100 ns (nanoseconds) in a supercomputing system provided by the Bioinformatics Division of National Institute of Biotechnology (NIB), Bangladesh. Thereafter, a comparative analysis was performed between WT and mutants measuring root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA) and hydrogen bonds. Qtgrace program was used to represent all these analyses in the form of plots⁸². Further, the g_mmpbsa⁸³ package of GROMACS was used to calculate the MM/PBSA (Molecular Mechanics/Poisson Boltzmann Surface Area) binding free energies followed by final MD run to get a more detailed overview of the biomolecular interactions between the mutated proteins and ligand. The tool was tested on 37 structurally divergent HIV-1 protease inhibitor complexes by performing comparison of the calculated relative binding energy with the experimental binding free energies⁸³. Also, the results obtained using g_mmpbsa package were comparable to results obtained with the AMBER package in general within differences of 1–2 kcal/mol. Furthermore, the package can be used to approximate the energy contribution per residue to the binding energy and it has been used to identify the crucial residues for binding a range of inhibitors with HIV-1 protease⁸³. The free solvation energy (polar and nonpolar solvation energies) and potential energy (electrostatic and Van der Waals interactions) of each protein–ligand complex were analyzed to determine the total ΔG_bind. The binding energies were calculated using the following equation in this method:

Here, the ΔG_binding = the total binding energy of the protein–ligand complex, G_protein = the binding energy of free protein, and G_ligand = the binding energy of unbounded ligand.