3.10. Molecular Dynamics

Molecular dynamics (MD) was performed with Gromacs 2022.3 [33,34] using the Charmm36-jul2021 force field [35,36]. Force field parameters for AZ8838, 1-PPA, and RO54464 were derived using cgenff and manually added to complexes’ topology [37,38]. Simulations of 1 µs were performed as previously described [39]. Briefly, models were solvated with the TIP3P water model in a rectangular box with a minimum distance of 1 nm between the protein complex and the border. Simulations were performed with explicit solvent content and 0.15 M of NaCl was added to simulate physiological ionic strength. System energy was minimized by 5000 steps of steepest descent energy minimization, with a tolerance of 1000 kJ mol⁻¹ nm⁻¹. Subsequently, a 200 ps NVT MD simulation was used to heat the system from 0 to 310 K and equilibrated to 1 atm during a 1 ns NPT simulation. Energy minimization and equilibration simulations were performed with 1000 kJ mol⁻¹ nM⁻² on all atoms; however, they were removed for the production run. The V-rescale thermostat was used to equilibrate the temperature, whereas the C-rescale barostat was used to control the pressure [40,41]. Newton’s equation of motion was integrated using a leapfrog algorithm with a 2 fs time step. The particle mesh Ewald (PME) method was used to compute the long-range electrostatic force [42,43]. Rotational and translational motions of the system were removed, and all bonds were constrained with the LINCS algorithm. Protein conformation variations were analyzed by estimating global values of Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) from the original structure with Gromacs built-in rms and rmsf tools, respectively. PCA on the motion of protein–ligands complexes was performed using Gromacs built-in covar, anaeig, and sham tools. Covariance matrices were obtained from Gromacs anaeig tool and converted into cross-correlation matrices. Finally, similarity between cross-correlation matrices was assessed with Spearman rho test.