Ligand Gaussian Accelerated Molecular Dynamics 2 (LiGaMD2): Improved Calculations of Ligand Binding Thermodynamics and Kinetics with Closed Protein Pocket

Abstract

Ligand binding thermodynamics and kinetics are critical parameters for drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics from molecular simulations due to limited simulation timescales. Protein dynamics, especially in the ligand binding pocket, often plays an important role in ligand binding. Based on our previously developed Ligand Gaussian accelerated molecular dynamics (LiGaMD), here we present LiGaMD2 in which a selective boost potential was applied to both the ligand and protein residues in the binding pocket to improve sampling of ligand binding and dissociation. To validate the performance of LiGaMD2, the T4 lysozyme (T4L) mutants with open and closed pockets bound by different ligands were chosen as model systems. LiGaMD2 could efficiently capture repetitive ligand dissociation and binding within microsecond simulations of all T4L systems. The obtained ligand binding kinetic rates and free energies agreed well with available experimental values and previous modeling results. Therefore, LiGaMD2 provides an improved approach to sample opening of closed protein pockets for ligand dissociation and binding, thereby allowing for efficient calculations of ligand binding thermodynamics and kinetics.

Figure 1.

Comparison of LiGaMD and LiGaMD2 simulations on the T4L mutant systems with buried and open binding pockets: Computational models of benzene binding to the T4L:L99A with a burred binding pocket (A) and T4L:F104A with an open binding pocket (D); Time courses of ligand root-mean-square deviation (RMSD) in T4L:L99A calculated from 49.2 ns LiGaMD (B) and LiGaMD2 (C) equilibration simulations, respectively. Time courses of ligand RMSD in F104A T4L calculated from 49.2 ns LiGaMD (E) and LiGaMD2 (F) equilibration simulations, respectively.

Figure 2.

LiGaMD2 simulations captured repetitive dissociation and binding of two different ligands (benzene and indole) to T4L mutants: (A-D) time courses of ligand heavy atom RMSDs relative to X-ray structures calculated from three independent 1 μs LiGaMD2 simulations of (A) benzene binding to T4L:L99A, (B) benzene binding to T4L:F104A, (C) benzene binding to the T4L:M102A and (D) indole binding to the T4L:L99A. (E-H) The corresponding PMF profiles of the ligand RMSDs averaged over three LiGaMD2 simulations of (E) benzene binding to T4L:L99A, (F) benzene binding to T4L:F104A, (G) benzene binding to the T4L:M102A and (H) indole binding to the T4L:L99A. Error bars are standard deviations of the free energy values calculated from three LiGaMD2 simulations.

Figure 3.

2D Free energy profiles and low-energy intermediate conformational states of ligand binding to the T4L mutants: (A-D) 2D PMF profiles regarding the ligand heavy atom RMSD and the pocket volume in LiGaMD2 simulations of (A) benzene binding to T4L:L99A, (B) benzene binding to T4L:F104A, (C) benzene binding to T4L:M102A, (D) indole binding to T4L:F99A. (E-H) Low-energy “Intermediate” (“I”) conformations (blue) as identified from the 2D PMF profiles of (E) benzene binding to T4L:L99A, (F) benzene binding to T4L:F104A, (G) benzene binding to T4L:M102A and (H) indole binding to T4L:L99A. X-ray structures of the ligand-bound complexes (“Bound”) are shown in green. The ligands are shown in balls and sticks, and the helix are shown in cartoon. The helix C, D, F and G are labeled as they show significant changes between the “Bound” and “Intermediate” conformational states.

Figure 4.

Pathways of ligand binding and dissociation in the T4L mutants. (A) Cartoon representation of the protein with helices labelled. Binding and dissociation pathways are denoted by the arrow lines. Number of binding (B) and dissociation (C) events through the different pathways captured by the LiGaMD2 simulations.