2.3. Molecular dynamics (MD) simulations of SARS-CoV-2 PLpro and inhibitors

The docked complex of SARS-CoV-2 PLpro with four ligands (TTT, PB5, CQ, HCQ) were used as starting conformation in Molecular Dynamics Simulations (MDS). For assessing the binding stability and molecular interaction, apo and holo structures (docked complexes) which include CQ, HCQ and reference inhibitor (TTT and PB5) were subjected for MD simulation studies. A total of five MDS systems were prepared and subjected for MDS with 50 ns production run (SARS-CoV-2 PLpro; SARS-CoV-2 PLpro_PB5; SARS-CoV-2 PLpro_TTT; SARS-CoV-2 PLpro_CQ and SARS-CoV-2 PLpro_HCQ). MD simulation was performed with GROMACS ver.2016.4 (Van Der Spoel, 2005) with Amber99SB force‐field (Lindorff-Larsen et al., 2010) and the protocol was essentially the same as mentioned in our previous reports (Patel et al., 2017, 2018). Briefly, ligand parameter and topology was prepared by ACPYPE (Sousa da Silva & Vranken, 2012), wherever required. All the MDS systems were solvated with Three-site water model (TIP3P) in dodecahedron box, maintaining a distance of 1 nm from edges of protein in all directions. The MDS systems were neutralised with an equal number of counter ions (Na⁺/Cl⁻) followed by energy minimisation using the steepest descent algorithm to remove any steric clashes and bad contacts with maximum force <1000 kJ mol − 1 nm − 1 (50,000 steps max). After the energy minimisation, equilibration with position restraint was carried out under NVT (constant number [N], constant volume [V] and constant temperature [T]) and NPT (constant number [N], constant pressure [P] and constant temperature [T]) ensemble for 1 ns each. For NVT equilibration, modified Berendsen thermostat algorithm (Berendsen et al., 1984) was used for maintaining the system at constant volume (100 ps) and at a constant temperature (300 K). Subsequently, NPT equilibration was performed at a constant pressure (1 bar) for 100 ps maintained by Parrinello‐Rahman barostat (Parrinello & Rahman, 1980). The Particle Mesh Ewald approximation was applied with 1 nm cut-off for calculating long‐range electrostatic interactions, computing coulomb & the Van der Waals interactions (Darden et al., 1993). The bond length was constrained using the LINCS algorithm (Hess et al., 1997). Finally, 50 ns production run was executed with default parameters, and the coordinates were saved after every 2 fs time frame. The MDS trajectories were visualised using VMD (Humphrey et al., 1996) and Chimera (Pettersen et al., 2004). For calculation of root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bonds (H‐bonds) etc. ‘gmx’ commands were used in GROMACS and the plotting tool GRACE was used for generating the plots (http://plasma-gate.weizmann.ac.il/Grace) (Patel et al., 2018). For molecular visualization and other analysis Discovery studio was utilized (BIOVIA 2020).