2.2. System Setup and MD Simulations

In the current investigation, ten independent systems were prepared and used for MD simulations. All the prepared systems were simulated using Amber 16 suite [31] with Amber ff03.r1 [32] force field, individually. Initially, the tleap program available in Ambertools 16 [31] was employed to build the model systems by adding the (i) the force field parameters; (ii) the hydrogen atoms; (iii) an appropriate number of counter-ions to neutralize the unbalanced charges; and (iv) an adequate number of transferable intramolecular potential three-point (TIP3P) water molecules [33]. Then, the system was placed in an octahedral box extended by 9 Å in every dimension from the atoms of the solute to reduce the high computational cost. Subsequently, a step-by-step equilibration was performed as follows. Initially, the system was subjected to energy minimization based on the following steps: First, all the Cα atoms were constrained with a harmonic force constant of 2 kcal/mol/Ų, and allowing the rest of the atoms, e.g., water molecules, counter-ions and amino acid side-chains, to move freely. 500 steps of steepest descent and 1000 steps of conjugate gradient methods were applied during the minimization. Then, the system was gradually heated from 0 to 300 K, followed by density equilibration with weak restraints on the complex for 50 ps, respectively. Next, the whole system was set free i.e., without any constraints, and treated with NPT ensemble for 500 ps, while the temperature was maintained at 300 K and the pressure at 1 atm using Langevin dynamics during the simulation. Finally, considering the fact that the system built using the homology models may require a certain time period to reach dynamic stability, a 100 ns production run was carried out separately on all ten mMcl1—PAP complexes. Furthermore, to obtain reliable BFEs comparable with experimental values for all ten mMcl1—PAP complexes, the MD simulations were extensively sampled by (i) obtaining an average snapshot from the final phase of the 100 ns trajectory and used as a starting structure, (ii) the production run was extended to 25 ns with ten repeats each, yielding a total of 3.85 µs; (iii) all subsequent analyses were performed using the trajectories obtained from the extended time period (i.e., 25 ns). During the simulations, 2 fs time step integration was implemented for the entire simulation. The SHAKE algorithm [34] was used to constrain all hydrogen atoms. The 8 Å cut-off was applied onto the short-range—electrostatic and van der Waals—interactions and monitored every step, and the particle-mesh Ewald [35,36] method was applied to monitor the long-range electrostatic interaction at every third step. Periodic boundary conditions were applied to all dimensions of the system.