MD/FEP simulations

The software package Q [23] was used for the MD sampling under SBC with the following settings: atoms lying outside the simulation sphere were tightly constrained to their initial coordinates (200 kcal/mol/Ų force constant) and excluded from the calculation of the non-bonded interactions. In the boundary of the sphere, solvent atoms subject to polarization and radial restrains using the surface constrained all-atom solvent (SCAAS) [23, 54] model to mimic the properties of bulk water at the sphere surface. Long range electrostatics interactions beyond a 10 Å cut off were treated with the local reaction field method [55], except for the atoms undergoing the FEP transformation where no cutoff was applied. Solvent bonds and angles were constrained using the SHAKE algorithm [56].

The system was then subjected to 10 parallel MD replicate simulations, which only differed in their initially assigned random (Maxwell-Boltzman) velocities. Each simulation includes an equilibration scheme as follows, which is the default in QligFEP: (1) an initial phase of 31 ps where the simulation sphere is heated from 0.1 to 298 K and a positional restraint of 25 kcal/mol/Ų initially imposed on all solute heavy atoms slowly released; (2) a 100 ps unbiased and unrestrained equilibration. The subsequent production phase followed, where the sampling along a λ transformation pathway (defining the FEP transformation) was performed with the following parameters in the different datasets (unless indicated the contrary): the transformation was divided into 51 λ steps, evenly distributed (i.e. linear sampling), the corresponding MD sampling for each λ step being 10 ps using a 1 fs time step. Thus, each FEP transformation involved a total simulation time of 10 × (131 ps + (51 × 10 ps)) = 6.41 ns for each of the two legs in the thermodynamic cycle, resulting in a total sampling time of 12.82 ns per perturbation. For comparison, a typical perturbation in the commercial alternative FEP+ includes one simulation of 12 lambda windows, with 5 ns sampling per window giving 120 ns under a bigger system using PBC [30]. However, in cases where large changes are considered (i.e. more than 6 heavy atoms) the sampling above might be insufficient to achieve the desired convergence, which might be easily diagnosed by a high SEM or even simulation crashes. In such cases, one can increase the phase-space overlap by increasing the λ windows, or either using a more dense (non-linear) sampling around the initial/ending states, or just increase the sampling of each λ window to increase convergence. In this work, a different scheme was used for the larger perturbations involved in the datasets of the relative hydration free energies and the A_2AAR antagonist binding.

In the dual topology paradigm, each perturbation involves parallel growing (A) or annihilation (B) of the relative atoms of the two ligands A and B, and includes the perturbation to/from soft-core potentials as an intermediate step to ensure sufficient convergence [57]. Half-harmonic distance restraints of 2.0 kcal/mol/Ų were applied to maintain pairs of equivalent (non-dummy) atoms between the two ligands (A and B), within a window distance of 0.0–0.2 Å. Since this restraint is only applied on atom pairs, of which one of them will have the interactions fully turned off in each of the end states, the energetic term of the restraints cancels and no additional correction needs to be applied. The relative free energy (∆∆G_A-B) is calculated by closing the corresponding thermodynamic cycle, where the ∆G corresponding to each of the two legs was calculated with the BAR approach [45]. When the ligand mutation involves a change in the total charge of the sphere, a Born correction term was added to the calculated free energies to account for the polarization effect [58], estimated as:

where Q_I is the net charge of the solute, and r_Born is the radius of the cavity in the macroscopic medium, with a dielectric constant ε, in which the charge is embedded.

All calculations were performed using the OPLS-AA/M [33] force field with TIP3P waters, except for the CDK2 inhibitor set where results are reported with the three families of force fields implemented in QligFEP. Standard errors of the mean (SEM) are estimated from the individual replicate simulations by QligFEP, whereas errors reported for FEP+ are either BAR analytical error estimations or bootstrap estimated errors as reported in their original publications [49].