Robust Prediction of Relative Binding Energies for Protein-Protein Complex Mutations Using Free Energy Perturbation Calculations

Abstract

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein-protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how Protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.

Keywords: binding affinity prediction; free energy methods; in silico mutational screening; protein binding interface optimization; protein-protein interactions.

Figure 1:

Overview of systems and mutations in the benchmark dataset. (a) Ribbon representation of each model in the benchmark dataset, colored by chain, with Cα positions of mutated residues shown as dark grey spheres. Original PDB accession codes are indicated. (b) Number of mutations per system, colored by FEP perturbation type. (c) Distribution of experimental binding $\mathrm{\Delta }\mathrm{\Delta }G$ values across all systems, also colored by perturbation type. Median and mean experimental values are shown as dashed and dotted vertical lines, respectively. (d) Heat map showing coverage of amino acid mutation space. Due to the distribution of mutations in the underlying experimental data, mutations to ALA are visibly over-represented.

Figure 2:

FEP+ results versus experiment. Correlation plots of FEP+ predicted binding $\mathrm{\Delta }\mathrm{\Delta }G$ vs. experimental binding $\mathrm{\Delta }\mathrm{\Delta }G$ relative to WT for the benchmark dataset mutations. Results are shown for 10 ns (left) and 100 ns (right) FEP+ calculations, with both naïve (top) and FEP+ Groups (bottom) treatment of protonation states. Grey and light grey diagonal shaded areas represent regions of absolute error less than 1 and 2 kcal/mol, respectively. A least squares best fit line (bold) and relevant correlation statistics are indicated for each plot.

Figure 3:

Protein FEP+ Groups treatment example cases. Ribbon representations of either the wt input structure or a representative mutant trajectory frame, with the chain(s) for the binding component that includes the mutation shown in light grey, the binding partner chain(s) in tan, and key polar interactions as dashed yellow lines. FEP-calculated 100 ns binding $\mathrm{\Delta }\mathrm{\Delta }G$ values for mutation edges are shown along solid arrows, and FEP-based calculated pKas for titration edges in the bound (complex) and unbound (solvent) forms are shown along the horizontal dashed arrows. Experimental, naïve, and FEP groups-treated $\mathrm{\Delta }\mathrm{\Delta }G$ values are listed in kcal/mol. (a) For 1IAR A:W91D the calculated pKas for the mutant Asp91 carboxylate indicated a change in dominant protonation state between the unbound and bound forms at experimental pH. FEP+ Groups treatment resulted in reduced absolute error compared to the naive result. (b) The 6NRQ C:Q138D mutant Asp138 carboxylate had elevated FEP-calculated pKa values in both unbound and bound states, reflecting the largely hydrophobic residue environment, limited solvent access, and the conformational change required to form a salt bridge with D:K81. FEP+ Groups treatment again increased accuracy with respect to the experimental value.

Figure 4:

Classes of outlier cases with absolute errors greater than 2 kcal/mol. Structural images in ribbon representation, with key wt sidechains shown as thin sticks, colored by chain and atom type. Mutations are indicated above each panel and mutated residues are shown in ball-and-stick representation. (a) Salt bridge-breaking cases include trans (left two panels) and cis (second from right panel) interface salt bridges, and an internal salt bridge between the heavy and light chains of an antibody Fv region (right panel). (b) Two outlier cases introduced a buried charge. (c) A variety of sampling related cases. (d) Outliers resulting from mutation of buried aromatic residues to much smaller alanine. (e) One force field-related case involved a disrupted $\mathrm{\pi }$-hydrogen bond interaction.

Figure 5:

Automated classification and empirical correction of FEP+ outliers for the benchmark dataset. (a) Correlation plots 100 ns FEP+ Groups-treated results vs. experiment, before and after automated outlier flagging and empirical correction of charged outlier cases. Flagged outlier categories are indicated. The direction and magnitude of the applied correction for charged outlier cases is indicated by orange vertical tails. One case (1JRH I:N53A) which demonstrated substantial improvement in preliminary testing of a force field with explicit polarization is indicated by a green circle and a green vertical tail in the right plot. Diagonal shaded areas, least squares fit line, and statistics are shown as in Figure 2. (b) Three-class confusion matrices for the data from (a), with $\mathrm{\Delta }\mathrm{\Delta }\mathrm{G}$ values classified as either favorable ($\mathrm{\Delta }\mathrm{\Delta }\mathrm{G}\le -0.5\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$), unfavorable ($\mathrm{\Delta }\mathrm{\Delta }\mathrm{G}\ge 0.5$), or neutral $(-0.5<\mathrm{\Delta }\mathrm{\Delta }\mathrm{G}<0.5)$.

Figure 6:

Application of automated outlier flagging and correction to the case studies dataset. (a) Correlation plots and (b) three-class confusion matrices for the 100 ns FEP+ Groups-treated results vs. experiment, before and after automated outlier flagging and empirical correction of charged outlier cases as in Figure 5. The magnitude of the empirical correction applied to the flagged charged outlier cases in (a) was determined using the flagged cases from the benchmark dataset only.