Protein thermostability calculations using alchemical free energy simulations

Abstract

Thermal stability of proteins is crucial for both biotechnological and therapeutic applications. Rational protein engineering therefore frequently aims at increasing thermal stability by introducing stabilizing mutations. The accurate prediction of the thermodynamic consequences caused by mutations, however, is highly challenging as thermal stability changes are caused by alterations in the free energy of folding. Growing computational power, however, increasingly allows us to use alchemical free energy simulations, such as free energy perturbation or thermodynamic integration, to calculate free energy differences with relatively high accuracy. In this article, we present an automated protocol for setting up alchemical free energy calculations for mutations of naturally occurring amino acids (except for proline) that allows an unprecedented, automated screening of large mutant libraries. To validate the developed protocol, we calculated thermodynamic stability differences for 109 mutations in the microbial Ribonuclease Barnase. The obtained quantitative agreement with experimental data illustrates the potential of the approach in protein engineering and design.

Figure 1

Hybrid residues. A database of hybrid residues has been compiled. A script replaces one or more residues in a structure file by hybrid residues that represent one amino acid in state A and another amino acid in state B.

Figure 2

Thermodynamic cycle. The unfolded state was modeled with capped GXG peptides (X = any amino acid). From the thermodynamic cycle, the folding free energy difference ΔΔG = ΔG₃-ΔG₂ between the wild-type protein and the mutant can be calculated via ΔG₁–ΔG₄.

Figure 3

Folding free energy differences for Barnase mutations. (Left) Scatter plot of experimental values versus calculated values. The two thin lines parallel to the diagonal line represent deviations of ±1 kcal/mol. (Right) Deviation from experimental values.

Figure 4

Assessment of calculation accuracy. (A) Alanine, non-Alanine, and Glycine mutations; (B) dependence on secondary structure; (C) dependence on packing properties; and (D) dependence on size.

Figure 5

Accuracy versus total simulation time. With only 20% of the computational effort, an accuracy of >65% of the predictions within ±1 kcal/mol is obtained. Convergence of the accuracy is apparently not reached at 26 ns per mutation.

Figure 6

Accuracy for mutations with a net charge change. (Left) The accuracy for mutations which result in a change of the net charge of the system is significantly worse than for neutral mutations. Only 52% of the calculated values are within ±1 kcal/mol of the experimental value, and six of the 25 mutations deviate by >10 kJ/mol (one data point is not shown). (Right) Dependence of the accuracy to the relative solvent-accessible surface area. Mutations at highly solvent-exposed positions are in favorable agreement with experimental data, whereas mutations at partly buried positions are badly predicted.