2.2. Molecular Dynamics Simulation

The molecular dynamics (MD) simulations were performed in Gromacs 2019, ref [32] and we considered OPLS-AA force field. OPLS-AA is an interesting force field that considers all atoms of a protein explicitly and is useful in the study of the organic and biomolecular systems [33]. The protein was situated inside a cubic box with a distance of 1.5 nm to the edges of the box and the explicit TIP4P [34] water model was considered. Na+ ions were added to neutralize the system. The energy minimization using steepest descent method of 200,000 steps with 0.001 nm initial step-size was employed.

Leap-frog algorithm for integrating Newton’s equations was employed to MD simulations. The bond lengths were constrained by the LINCS (LINear Constraint Solver) algorithm, ref [35] and periodic boundary conditions (PBC) in all directions, were used. Cut-off of 0.9 nm distance for the short-range interactions and Particle-Mesh Ewald (PME) [36] method for the long-range interactions were considered. V-rescale thermostat and Parrinello-Rahman barostat regulated the temperature (300 K) and pressure (1 bar).

We applied three MD steps: annealing, equilibrium, and production simulation. The first, simulated annealing procedure consisted of 20 annealing cycles in 500 ps with a high temperature of 400 K and a cooling temperature until 300 K by 10 ns. In the second step, the equilibration was performed in the NVT (number of molecules, volume, and temperature constant) canonical ensemble by a short time of 10 ns for integration, considering position restraint. They were restrained with harmonic forces and the water molecules were allowed to equilibrate around the protein. In the third step, we considered the MD in the isobaric-isothermal ensemble without position restraint during 200 ns was calculated. Furthermore, five systems with these temperatures: 300 K, 350 K, 400 K, 450 K, and 500 K, were considered at this point.

The analysis of Root-mean-square deviation (RMSD), Root-mean-square fluctuation (RMSF), radii of gyration (Rg), Solvent Accessibility Surface Area (SASA), and hydrogen bond (Hb), were obtained using the GROMACS tools and plotted in Gnuplot 5.0 software. We considered the average values of the last 20 ns of the MD simulation. Salt bridges of 5 frames of the trajectory (180 nm, 185 nm, 190 nm, 195 nm, and 200 nm) were analyzed by VMD (Visual Molecular Dynamics) [37] software. All graphics representation was recreated in Chimera UCSF 11.1.2 software [38].