Protein-peptide association kinetics beyond the seconds timescale from atomistic simulations

Abstract

Understanding and control of structures and rates involved in protein ligand binding are essential for drug design. Unfortunately, atomistic molecular dynamics (MD) simulations cannot directly sample the excessively long residence and rearrangement times of tightly binding complexes. Here we exploit the recently developed multi-ensemble Markov model framework to compute full protein-peptide kinetics of the oncoprotein fragment ^25-109Mdm2 and the nano-molar inhibitor peptide PMI. Using this system, we report, for the first time, direct estimates of kinetics beyond the seconds timescale using simulations of an all-atom MD model, with high accuracy and precision. These results only require explicit simulations on the sub-milliseconds timescale and are tested against existing mutagenesis data and our own experimental measurements of the dissociation and association rates. The full kinetic model reveals an overall downhill but rugged binding funnel with multiple pathways. The overall strong binding arises from a variety of conformations with different hydrophobic contact surfaces that interconvert on the milliseconds timescale.

Fig. 1

Metastable states and transition rates for the binding of PMI to Mdm2. The PMI peptide is colored according to (). States are represented by discs with areas proportional to the natural logarithm of the equilibrium probability. Arrows indicate transitions with rate constants of at least 1 ms⁻¹ in either direction. Numbers quantify transition rate constants in ms⁻¹ M⁻¹ for association events and in ms⁻¹ for all other transitions. The definition of the states is hierarchical: between top-level states 0 and 13, transitions happen on timescales of 10 µs or slower. States in the lower part of the figure are sub-states of top-level state 13. There, PMI transitions between different states in the main binding pocket of Mdm2 on timescales of microseconds or slower (only states with large probabilities are shown)

Fig. 2

Computational predictions and experimental validations of binding affinities and kinetics of PMI–Mdm2. a, b Maximum likelihood estimates of the binding free energy and residence time, respectively, as a function of the amount of data used. Diamonds mark maximum likelihood estimates, error bars indicate 95% confidence intervals. x% of total data means that x% of all biased data and x% of all unbiased data were used. c, d Estimates of binding free energy and residence time as a function of data composition. The fraction of biased simulation data is varied between 0 and 100% of all biased data while keeping the sum of the amount of biased and the amount of unbiased data constant at 502 μs. e Validation of the simulation model: binding free energy differences (ΔΔG) of PMI–Mdm2 upon alanine mutations of the PMI peptide, compared between the present simulations and experiments. Only biased simulation data was used and analyzed with MBAR. Error bars mark standard deviations (simulation error computed using bootstrap, see Supplementary Note 3.3). f Cyan: average and standard deviations of anisotropy time traces from three repeats of a binding competition experiment, Mdm2 (10 nM), pre-incubated with labeled PMI (10 nM), was mixed with unlabeled PMI (10 μM). A mono-exponential function (black) was fitted to the average time trace. Gray: control without unlabeled PMI

Fig. 3

Binding mechanism comprised by the 60% most probable pathways. Structures of metastable (on-pathway) intermediates are shown, labels are as in Fig. 1. Arrows indicate the direction and relative magnitude of the reactive flux from the dissociated state to the crystal-like bound state. PMI residues that form PMI–Mdm2 contacts with at least a probability of 0.5 in a given macro-state are shown as sticks

Fig. 4

Illustration of computing rare-event kinetics with TRAMMBAR using a model for protein ligand binding. a Potential energy surface and transition rates between five states (bound, pre-bound, two mis-bound states, dissociated). Arrow thickness is proportional to rate. b Probability of computing the binding free energy ΔG within 1k_BT accuracy of the exact value for a given amount of simulation data using MEMMs (TRAMMBAR estimator) or MSMs. The vertical bar indicates the mean-first-passage time (MFPT) for dissociation. c Probability of computing the dissociation rate within factor $\frac{1}{2}$ to 2 accuracy of the exact value