Physics-based energy function models for protein stability prediction have been around for decades. The basic idea is straightforward: represent the protein's energy as a function of its coordinates, calculate how a substitution changes that energy, and use the energy difference as a proxy for the change in thermodynamic stability. The conceptual clarity of this approach is appealing, and for specific classes of mutations under specific conditions, it works well. The failure modes are equally specific, and knowing them matters when you are deciding which tool to use for a given prediction task.
This post is a practitioner's comparison of the main energy function approaches used in stability prediction, written from the perspective of building a hybrid model that incorporates physics-based features alongside learned sequence representations. I am not going to name specific commercial software or make claims about proprietary systems. The comparison is about the functional classes of energy model and where each one gains and loses accuracy in real stability prediction benchmarks.
Molecular mechanics force fields
Molecular mechanics (MM) force fields model bonded interactions (bond stretching, angle bending, torsion), van der Waals interactions (Lennard-Jones potential), and electrostatics (Coulomb with some solvation treatment) using a set of empirically derived parameters for each atom type. AMBER, CHARMM, and OPLS are the most widely used in protein modeling contexts. When combined with explicit solvent molecular dynamics and thermodynamic integration or free energy perturbation methods, they can produce high-accuracy predictions of relative stability changes for point mutations.
The accuracy comes at a steep computational cost. A rigorous delta-delta-G calculation using free energy perturbation for a single mutation in a protein of 300 residues requires nanosecond-to-microsecond MD simulations with explicit solvent. For a panel of 200 mutations, this is practical only on high-performance computing infrastructure with substantial allocation. For a screening task where the goal is to narrow 200 candidates to 20, the throughput bottleneck is severe.
MM force fields also depend on accurate starting coordinates. An error in the input structure, particularly in a flexible loop region or at a position near the mutation site, can cause the predicted energetics to reflect the structural error rather than the genuine mutation effect. For proteins where only a low-resolution crystal structure or a homology model is available, the coordinate uncertainty propagates into the energy prediction.
Statistical empirical potentials and knowledge-based energy functions
Knowledge-based energy functions derive their parameters from the statistical distribution of structural features observed in the PDB. Distance-dependent contact potentials, burial energy terms, solvation estimates based on side-chain accessible surface area, and hydrogen bond geometry potentials all fall into this category. Tools like Rosetta's REF2015 energy function are hybrid: they incorporate both physics-based terms and statistically derived components.
These approaches run orders of magnitude faster than explicit MM with free energy perturbation. Rosetta's single-point mutation delta-delta-G calculation on a static backbone typically runs in seconds to minutes per mutation on a modern workstation. The tradeoff is accuracy: statistical potentials are calibrated to the distribution of observed structures, which reflects evolutionary selection rather than pure physical chemistry. They perform well on mutations at well-packed core positions but have known failure modes for surface mutations, where solvation treatment is more critical, and for mutations that cause significant backbone rearrangement.
The backbone flexibility problem is a persistent challenge for all structure-based methods. Most rapid energy function evaluations use a fixed backbone: the input structure's backbone atoms stay in place, and only the mutated sidechain is repacked. For mutations that cause backbone adjustment (insertions, deletions, or substitutions at glycine or proline positions that strongly constrain backbone dihedral angles), fixed-backbone predictions are systematically less reliable. Flexible backbone protocols (e.g., backrub motions, kinematic closure) improve this but are computationally more expensive and introduce additional degrees of freedom.
Coarse-grained and residue-level potentials
For very fast screening, coarse-grained models that represent each residue as one or a few interaction sites (rather than all-atom representations) enable scanning of entire proteins in seconds. FoldX is perhaps the most widely used tool in this category for stability prediction. Its parameterization is specifically tuned to reproduce delta-delta-G values from the thermodynamic database literature, which makes it useful as a first-pass filter even when the absolute values are not reliable for all mutation types.
The accuracy of FoldX and similar tools on the ProThermDB benchmark is a reasonably good correlator with stability outcomes for buried hydrophobic core mutations but degrades for polar mutations, charged residue substitutions, and mutations that alter the local backbone strain. These are well-documented failure modes. The tool is fast enough and reasonably calibrated enough for its intended purpose, rapid candidate triage, and should not be used as a quantitative replacement for more rigorous methods at later stages of candidate evaluation.
Where hybrid models fit and why we built one
The physics-based approaches described above are interpretable and structurally grounded. Their limitations are largely orthogonal to the limitations of sequence-based statistical learning: physics-based tools fail when structure is uncertain or when mutations cause backbone rearrangement; sequence-based tools fail when alignment depth is low or when the relevant evolutionary signal is absent. A hybrid model that uses physics-derived energy terms as features alongside sequence-fitness learned representations is better positioned to cover the failure modes of either approach alone.
In our model architecture, physics-based terms contribute most strongly to the prediction at buried core positions in well-studied protein families where the structural context is reliable and alignment depth is adequate. Sequence-based features contribute most strongly at surface positions, for families with moderate alignment depth, and for mutation types where the physics-based terms are known to have systematic biases.
The empirical result on our internal benchmark across 47 protein families is that the hybrid approach outperforms either pure-physics or pure-sequence models on every protein family subset, though the margin varies. For the 12 protein families where alignment depth was highest and structure quality was best, the hybrid advantage over pure-physics models was modest (about 4 to 8% improvement in top-3 recall). For the 15 families with limited alignment depth or lower-quality structural models, the hybrid advantage was larger: pure-physics methods showed substantially reduced accuracy while the hybrid maintained more consistent performance.
Benchmarking practices and why they matter
Evaluating any stability prediction tool requires careful benchmark design. The most common error is using a benchmark that overlaps with the tool's training data. If a statistical potential was parameterized using the ProThermDB, testing it on ProThermDB variants is not an independent evaluation. This is a pervasive problem in the published literature. Tools that report strong performance on standard benchmarks should be interrogated for data leakage: was the benchmark dataset available when the parameters were set?
A second issue is selection bias in published benchmark datasets. The mutations in ProThermDB are not representative of the mutations you will actually be engineering. They are enriched for large-effect mutations (which get published more often), for well-studied protein families, and for mutations at positions that were specifically chosen for investigation because they were expected to have large effects. A tool calibrated on this distribution will tend to over-predict the magnitude of changes for real engineering mutations, which often involve smaller and noisier effects than the literature benchmark suggests.
The most meaningful benchmarking we do internally is prospective: we make predictions on a protein family before seeing the experimental results from a pilot partner, then compare predictions to measurements after the fact. That is the evaluation regime closest to actual use, and it reveals the failure modes that retrospective benchmarking systematically misses.