Calculation of Protein Folding Thermodynamics using Molecular Dynamics Simulations

Despite impressive advances by AlphaFold2 in the field of computational biology, the protein folding problem remains an enigma to be solved. The continuous development of algorithms and methods to explore longer simulation timescales of biological systems, as well as the enhanced accuracy of potential functions (force fields and solvent models) have not yet led to significant progress in the calculation of the thermodynamics quantities associated to protein folding from first principles. Progress in this direction can help boost related fields such as protein engineering, drug design, or genetic interpretation, but the task seems not to have been addressed by the scientific community. Following an initial explorative study, we extend here the application of a Molecular Dynamics-based approach −with the most accurate force field/water model combination previously found (Charmm22-CMAP/Tip3p)− to computing the folding energetics of a set of two-state and three-state proteins that do or do not carry a bound cofactor. The proteins successfully computed are representative of the main protein structural classes, their sequences range from 84 to 169 residues, and their isoelectric points from 4.0 to 8.9. The devised approach enables accurate calculation of two essential magnitudes governing the stability of proteins −the changes in enthalpy and in heat capacity associated to protein unfolding−, which are used to obtain accurate values of the change in Gibbs free-energy, also known as the protein conformational stability. The method proves to be also suitable to obtain changes in stability due to changes in solution pH, or stability differences between a wild-type protein and a variant. The approach addresses the calculation by difference, a shortcut that avoids having to simulate the protein folding time, which is very often unfeasible computationally.