Molecular Dynamics Simulations for Drug Binding

Further drug screening was carried out by force field based molecular dynamic (MD) simulations. The initial protein-drug complexes was from top score conformation Autodock Vina docking, the ligand was edited by pymol software [32 ] to make it in correct protonation state at pH 7. In this study, we selected 14 drug binding complexes for MD simulation, including Adenosine, Amenamevir, Amoxicillin, Azithromycin, Clofarabine, Fipronil, Gemcitabine, Nitisinone, Pralatrexate, Raltegravir, Romidepsin, Sofosbuvir, Teriflunomide and Vidarabine, respectively.
We also refined a pocket molecular dynamics simulation (pocket MD, S11B Fig) to facilitate the simulation process by only keeping the binding pocket region for simulation. Binding free energy calculation can be estimated by metadynamics simulations to explore whether protein-ligand will bind in solution. Metadynamics relies on addition of a bias potential to sample the free energy landscape along a specific collective variable of interest [33 (link),34 (link)]. Note that the binding free energy calculations from Metadynamics may only be suitable for detect the general trend of binding in virtual screening.
The pocket MD is same as the classical MD simulation, except that we only using the pocket region to reduce system size for simulation (S11B Fig), which is inspired by a previous dynamic undocking (DUck) method [35 (link)]. An in-house script was used to extract the pocket region of the protein (1nm within the binding ligand), the N terminal and C terminal ends were capped with the ACE and NHE terminals, respectively. We applied a position restrain to the ACE and NHE terminals to maintain the relative conformation of the pocket. MD simulation was carried out by Gromacs with AMBER-99SB force field [36 ,37 (link)]. The topology of ligand and the partial charges of ligand was generated by ACPYPE [38 (link)], which relies on Antechamber [39 (link)]. Firstly, we created a dodecahedron box and put the target-ligand complex at the center. A minimum distance from the protein to box edge was set to 1 nm. We filled the dodecahedron box with TIP3P water molecules [40 (link)], the counter ions were added to neutralize the total charge using the Gromacs program tool [41 (link)]. The long-range electrostatic interactions under the periodic boundary conditions was calculated with Particle Mesh Ewald approach [42 (link)]. A cutoff of 14 Å was used for van der Waals non-bonded interactions. Covalent bonds involving hydrogen atoms were constrained by applying the LINCS algorithm [43 (link)].
We performed the energy minimization steps with a step-size of 0.001ns, 100 ps simulation with isothermal-isovolumetric ensemble (NVT), and 10ns simulation with isothermal-isobaric ensemble (NPT) for water equilibrium. After that, a 100ns NPT production run (step size 2 fs) was carried out. The Parrinello-Rahman barostat and the modified Berendsen thermostat were used for simulation with a fixed temperature of 308 K and a pressure of 1 atm. RMSD and hydrogen bond number of the trajectory were calculated using Gromacs tools.
The simulation was continued using the metadynamics approach for exploring the free energy landscape. The interface coordination number of atoms of protein ligand complex was used as collective variable (CV). The protein-ligand interface coordination numbers correlate with the numbers of atom contact, and larger coordination number usually indicates that protein-ligand is in binding state.
The coordination number C is defined as follows by Plumed:
$C=\sum {}_{i\in A}\sum {}_{j\in B}{S}_{ij}$
and
$S_{ij}=\frac{1–{\left(\frac{r_{ij}–d_0}{r_0}\right)}^n}{1-{\left(\frac{r_{ij}–d_0}{r_0}\right)}^m}$
In the simulation, n was 6, m was 12, d₀ was 0 nm and r₀ was 0.5 nm. d₀ is a parameter of the switching function. r_ij is the distance between atom i and atom j. The degrees of contacts between two groups of atoms can be estimated by above function(1) [44 (link)]. Metadynamics simulation for each protein-ligand system was performed for 100 ns (except protein-Azithromycin, which was extended to 300ns in order to reach the 0 Coordination Number and achieve convergences). During the metadynamics simulation, Gaussian values were deposited every 1 ps with a height of 0.3 kJ/mol. The widths of the Gaussians were 5 for the coordination number. The free energy landscapes of the metadynamics simulations along the CV were generated by the Plumed program and plotted using Gnuplot [45 ].