Simulation of Reversible Protein-Protein Binding and Calculation of Binding Free Energies Using Perturbed Distance Restraints

Abstract

Virtually all biological processes depend on the interaction between proteins at some point. The correct prediction of biomolecular binding free-energies has many interesting applications in both basic and applied pharmaceutical research. While recent advances in the field of molecular dynamics (MD) simulations have proven the feasibility of the calculation of protein-protein binding free energies, the large conformational freedom of proteins and complex free energy landscapes of binding processes make such calculations a difficult task. Moreover, convergence and reversibility of resulting free-energy values remain poorly described. In this work, an easy-to-use, yet robust approach for the calculation of standard-state protein-protein binding free energies using perturbed distance restraints is described. In the binding process the conformations of the proteins were restrained, as suggested earlier. Two approaches to avoid end-state problems upon release of the conformational restraints were compared. The method was evaluated by practical application to a small model complex of ubiquitin and the very flexible ubiquitin-binding domain of human DNA polymerase ι (UBM2). All computed free energy differences were closely monitored for convergence, and the calculated binding free energies had a mean unsigned deviation of only 1.4 or 2.5 kJ·mol^-1 from experimental values. Statistical error estimates were in the order of thermal noise. We conclude that the presented method has promising potential for broad applicability to quantitatively describe protein-protein and various other kinds of complex formation.

Figure 1

Thermodynamic cycle of protein–protein binding used for the calculation of the binding free energy (ΔG_bind) in this work. Orange and blue cartoons represent the binding partners for which a binding free energy is calculated. Arrows indicate the direction of the considered reaction for which a free-energy difference is calculated. Restraints used to keep the molecules in a defined conformation (ENs) are indicated by circles filled with a cross. Distance restraints are indicated by black lines.

Figure 2

Scaling prefactors f(λ) for the potential energy functions of the different intermolecular distance restraints along λ used for the simulation of protein binding and unbinding. The potential energy functions used in eqs 2 and 3 are formulated as . The scaling factor for the intermolecular Cα–Cα distance restraints are indicated with a red line (HR, f(λ) = 4(1 – λ)³) and the linear scaling factor for the radial COM–COM distance restraint or linear COM–COM distance restraint in the z-component is indicated with a violet line (f(λ) = λ).

Figure 3

Examples for classifications of replica trajectories along the coupling parameter λ: round-trip (red line), full-trip (violet line), and no special classification (brown line).

Figure 4

Forward (black lines) and reverse (red lines) cumulative averages along the trajectories of the free-energy differences calculated with BAR for protein binding and unbinding (ΔG_unbind^res) in system RS of the wt UBM2 domain (A) and the P692A UBM2 domain (B), for the free-energy contribution of the ENs and specific intermolecular distance restraints in the bound state (ΔG_en,dr^b) of the wt complex (C) and the complex involving P692A UBM2 (D), for the free-energy contributions of the ENs in the unbound state (ΔG_en,1^u,ΔG_en,2^u) of wt UBM2 (E), P692A UBM2 up to 50 ns total simulation time (F), P692A UBM2 up to 100 ns total simulation time (H) and UBI (G). Simulations of the processes used to obtain (C–H) were performed using SCR potential energy functions. Note that in panel F the forward cumulative average falls off the scale of the plot.

Figure 5

Order parameters ⟨cos θ_αβ⟩, where θ_αβ is the angle between two vectors α and β, anchored to each protein. Three orthogonal vectors were defined in each protein to construct a system of orthogonal axes that were parallel in the bound proteins. (A) shows the development of all nine order parameters along λ for the simulation of protein binding and unbinding with the wt UBM2 domain and system setup RS. (B) shows the values of the same order parameters in the unbound state (λ = 1) for all performed simulations of protein binding and unbinding. (C) COM positions of wt UBM2 sampled around UBI shown in gray cartoon representation and (D) COM positions of UBI sampled around wt UBM2 shown in gray cartoon representation in the simulation of protein binding and unbinding with the wt UBM2 domain and system setup RS. In (C) and (D), the sampled COM positions are colored according to their λ point: λ = 0 (bound state) corresponds to blue, λ = 1 (unbound state) corresponds to red.

Figure 6

(A) Percentages of NOE upper bound violations of structural ensembles of trajectory windows with a simulation length of 1 ns for the unrestrained state of unbound wt UBM2 of the simulation employing SCR potential functions (red line), for the unrestrained state of unbound wt UBM2 of the simulation employing HR potential functions (violet line), and the restrained state of unbound wt UBM2 of the simulation employing SCR potential functions (brown line). (B) Secondary structure trajectory (blue: α-helical, red: β-strand, yellow: turn, green: bend, black: β-bridge, gray: 3-helix, white: coil) calculated with the DSSP algorithm for the unrestrained state of unbound wt UBM2 of the simulation employing SCR potential functions.