Accurate Prediction of Enzyme Thermostabilization with Rosetta Using AlphaFold Ensembles

Abstract

Thermostability enhancement is a fundamental aspect of protein engineering as a biocatalyst's half-life is key for its industrial and biotechnological application, particularly at high temperatures and under harsh conditions. Thermostability changes upon mutation originate from modifications of the free energy of unfolding (ΔG_u), making thermostabilization extremely challenging to predict with computational methods. In this contribution, we combine global conformational sampling with energy prediction using AlphaFold and Rosetta to develop a new computational protocol for the quantitative prediction of thermostability changes upon laboratory evolution of acyltransferase LovD and lipase LipA. We highlight how using an ensemble of protein conformations rather than a single three-dimensional model is mandatory for accurate thermostability predictions. By comparing our approaches with existing ones, we show that ensembles based on AlphaFold models provide more accurate and robust calculated thermostability trends than ensembles based solely on crystallographic structures as the latter introduce a strong distortion (scaffold bias) in computed thermostabilities. Eliminating this bias is critical for computer-guided enzyme design and evaluating the effect of multiple mutations on protein stability.

Figure 1

X-ray structures of LovD variants. Flexible loops encompassing residues 162–170 (in green) and 257–265 (in pink) and the 1–11 N-terminus (in yellow) are highlighted.

Figure 2

Representation of ΔT_Mvs ΔΔG_f,mut for LovD1-9 variants computed with respect to the wild-type (WT) as the average of the three top scoring decoys from 500 ddg runs starting from crystallographic structures. ρ = Pearson’s correlation coefficient. Dashed red line: linear regression.

Figure 3

Representation of ΔT_Mvs ΔΔG_f,mut for LovD1-9 variants computed with respect to the wild-type (WT) as the average of the 25 top scoring decoys from 3000 relax runs starting from crystallographic structures. ρ = Pearson’s correlation coefficient. Dashed red line: linear regression.

Figure 4

Representation of ΔT_Mvs ΔΔG_f,mut for LovD1-9 variants computed with respect to the wild-type (WT) from 100 AF runs averaging over the 25 top scoring decoys for each of the five AF models and all models combined. ρ = Pearson’s correlation coefficient. Dashed red line: linear regression.

Figure 5

Overlay of the 25 top scoring decoys obtained with the mAF-min method for LovD, LovD6, and LovD9 (AF models 1 and 4). The 162–170 loop is represented in green, and the 257–265 loop is represented in pink. RMSD = average root-mean-square deviation of loops computed on all heavy atoms of each ensemble.

Figure 6

Representation of ΔT_Mvs ΔΔG_f,mut for LipA variants computed with respect to the wild-type (WT) from 100 AF runs averaging over the 25 top scoring decoys for each of the five AF models and all models combined. ρ = Pearson’s correlation coefficient. Dashed red line: linear regression.