VM2: Free Energy Calculations by Mining Minima
The VM2 method, developed by VeraChem, calculates true free energies of binding using the mining minima method – a rigorous statistical mechanics based approach. VM2 can provide accurate protein-ligand binding affinities, and recent advances, such as parallel processor implementations for CPU and GPU clusters, as well as ports to Cloud platforms, now allow routine use of VM2 for rank-ordering of ligand series in pharmaceutical industry drug development settings.
Conventional drug lead optimization requires multiple rounds of intensive chemical synthesis and testing to make improvements in properties such as solubility while maintaining potency. Accurate predictions of protein-ligand binding affinities, therefore, provide the potential for tremendous savings with respect to cost and time in the drug development process.
If you think VM2 can help in your drug development efforts contact us at vc@verachem.com or submit a request to get VM2 below.
VM2: Fast, Accurate Protein-Ligand Binding Affinity Predictions
Benchmarking of VM2 compared to experiment shows fast, accurate prediction of protein-ligand binding affinities. The following is a comparison of VM2 Pearson correlation with experiment R with FEP+ and PMX/GROMACS for the Gapsys benchmark set:
VM2 correlation with experiment R compared with that from FEP+, Glide Dock, Glide Score, and Prime (MM/GBSA) for the Merck benchmark set:
VM2 wall time per ligand is 20-30 minutes on 12 CPUs/7-10 minutes on CPUs/GPUs. As VM2 calculations are entirely independent, given a modest sized CPU cluster, each entire benchmark can be run in several hours total. See this video presentation for further details on these VM2 benchmarks.
Additional VM2 Performance Examples
These plots are of VM2 calculated binding free energy predictions versus experimental binding affinities for the protein receptors P38A kinase, HIV-1 protease, PDE10A, and BRCA1, and corresponding inhibitor series ranging between 16 and 38 ligands. These data were generated both by VeraChem and academic users. An all-atom AMBER force field was used in each case with PBSA corrected solvation energies.
Results generated by industry partners for two protein receptors, one with a series of 170 ligands; the other with a seres of 50. An all-atom OPLS force field with PBSA corrected solvation energies was used.
The Importance of Entropy
This plot demonstrates entropy-enthalpy compensation found in experimental binding data – source BindingDB. That is, the stronger the enthalpy of binding the larger the entropy penalty paid. This is physically reasonable as a more strongly interacting ligand will be more immobilized as will the protein active site, resulting in a loss of entropy. However, the graph is not perfectly linear i.e. there is a width to the slope, and so enthalpy changes cannot be relied on exclusively to provide accurate ordering of binding affinities and entropy terms should also be included.
These two plots further demonstrate the need for entropy inclusion for accurate binding affinity prediction. On the left is the potential energy of binding from VM2 calculations is plotted against experimental binding
affinities i.e. entropy terms are excluded. On the right the entropy terms are included. System is PDE10a and 20 inhibitors. VM2 calculations used united atom CHARMM force field.
Theoretical Basis of VM2
Statistical thermodynamics of binding
VM2 is solidly grounded in a statistical thermodynamics approach to the prediction of binding affinities. VM2 aims to predict standard free energies of binding. Experiments can determine standard binding free energies through the equilibrium constant Kbind (see below), where the equilibrium constant itself is determined via the receptor, ligand, receptor-ligand concentrations along with the standard concentration.
VM2 calculates the standard free energy of binding via the standard chemical potentials of the receptor, ligand, and receptor-ligand complex.
VM2 uses the classical formulation of the partition function for calculation of each chemical potential.
Here the external degrees of freedom have been integrated out and C0 provides a correction to the standard state. Formally, the configuration integral must be determined over all space along the remaining internal degrees of freedom, the VM2 method approximates this using the concept of mining minima.
Mining minima
The mining minima approach replaces the configurational integral over all space with a sum over separate local configuration integrals associated with low energy minima of the system. For example, the figure below depicts a system where three low-lying energy minima have been found.
Determination of the local configurational integrals Zi allows a probability to be associated with each energy well, and this, in turn, allows a Boltzmann averaged energy
<U+W> to be determined and subtracted from the total free energy to give the system’s configurational entropy, useful when analyzing and interpreting predicted binding affinites.
Local configuration integrals
A Harmonic approximation together with a numeric correction for anharmonic low-energy modes, the Harmonic Approximation with Mode Scanning (HAMS) method, is used to calculate the local configuration integral
for each minimum produced. The harmonic approximation is applied via square roots of the Hessian matrix eigenvalues.
The Mode Scanning anharmonic corrections are calculated by numerical integration of the energy surfaces produced by distortion along low-energy internal-coordinate modes.
Molecular System Flexibility
Ligand molecules
All atoms in ligands and other small molecules (600 atoms or less) are always free to move unless the user chooses to apply constraints, such as harmonic or flat bottomed well potentials on particular atoms or groups of atoms. Such constraints can be useful for pre-generation of ligand R-group conformations, while maintaining a lead molecule’s scaffold conformation.
Protein molecules
All protein atoms can, in principle, be flexible in VM2 calculations. However, in practice, this is not required to achieve good binding affinity predictions, and it is computationally more efficient to partition the protein system into mobile, fixed, and excluded sets of atoms.
Implementation of VM2
Basic algorithm
The mining minima approach is iterative, with successive searches for low energy minima and calculation of a cumulative free energy. Once no new lower energy minima can be found the cumulative free energy stops changing and the calculation is converged.
Conformational searching
The VM2 method requires reliable determination of the low energy minima in a system. For this a robust conformational search method that can drive the system over energy barriers is needed so the system does not get stuck in higher energy local minima. VeraChem has implemented two types of search, which can be used on their own or in combination, one, a rigid body translation/rotation search, and two, a mode distort-minimize search.
Rigid body translation/rotations
The search can sample systematically or randomly. The translation/rotation is broken into small steps. After each step a few energy minimization steps are carried out with respect to the non-ligand atoms in the system allowing relaxation in response to the ligand distortion. After each completed translation/rotation distortion a full quasi-Newton geometry optimization is carried out to produce a minimum.
Click the following image to see a movie of random translation/rotations distortions of a ligand in a protein binding pocket (no minimization).
Mode distort-minimze search
Molecular torsional modes are calculated via diagonalization of matrix of energy 2nd-derivatives transformed into internal coordinates with all bond and angle rows and columns removed. As such, the modes are linear combinations of dihedral energy functions that include the full set of energy terms (valence, nonbonded, solvation). The system is then distorted along these modes. The distortion is broken up into small steps and atoms not involved in the mode allowed to relax at each step. After a complete distortion the whole system is energy minimized via a quasi-Newton geometry optimization.
The following illustrates the case of a small molecule mode-distortion and minimization.
Each time a new lowest energy minimum is found it is subsequently used as a basis for the next mode distort-minimization. The number of torsional modes available depends on the size of the system.
The search can cycle through the available modes in order (low to high energy) or at random. In addition, a linear combination of randomly chosen pairs of these modes can be used to define the distortion.
To further illustrate this approach consider the following six-atom system and its internal coordinates.
The total energy 2nd-derivative (Hessian) matrix with respect to internal coordinates (bond, angle, torsion coordinates) can be calculated
The dimensionality of the Hessian is reduced by only allowing rows and columns with dihedral angles; furthermore, only one dihedral angle per central bond is retained
This reduced size Hessian is then diagonalized. The resulting eigenvectors are each a linear combination of the two dihedral “basis functions” and each represents a torsional mode distortion driver. The initial conformation can then be distorted along each driver D with step size α
and as previously described a quasi-Newton geometry optimzation carried out after the distortion is complete.
The following movies depict various types of mode distortions.
Random-pair mode distortions.
Ligand focused mode distortions.
Protein focused mode distortions.
Cycles of distort-minimize that use the previous lowest energy minimum as a basis for the next distortion provide a mechanism for efficient conformational searching to find the lowest energy minima of the system.
Solvation Models
Generalized Born (GB)
For geometry optimizations and Hessian calculations (energy 2nd-derivatives) the Generalized Born solvation model is used.
Poisson-Boltzmann Surface-Area (PBSA)
The solvent energies (W) for minima found with the GB solvation model are corrected with the more accurate Poisson Boltzmann Surface Area (PBSA) method.
Explicit water molecules
For cases where particular water molecules in the binding site participate directly in the protein-ligand binding .
Underlying Energy Models
Classical force fields
Currently for protein-ligand VM2 free energy calculations only classical molecular molecular mechanics force field underlying energy models are available.
Quantum mechanics
Development work is underway to interface VM2 with quantum mechanical potentials for application to small molecular systems such as ligand molecules and host-guest systems e.g. ligand complexes with cyclodextrin.
VM2 applications
F. Molani, S. Webb, and A.E. Cho (2023) “Combining QM/MM Calculations with Classical Mining Minima to Predict Protein−Ligand Binding Free Energy” Journal of Chemical Theory and Computation. DOI:10.1021/acs.jcim.2c01637
G. C. Schröder, W. B. O’Dell, S. P. Webb, P. K. Agarwal and F. Meilleur (2022) “Capture of activated dioxygen intermediates at the copper-active site of a lytic polysaccharide monooxygenase.” Chemical Science, Edge Article, 13, 13303. DOI:10.1039/d2sc05031e
P. Xu, T. Sattasathuchana, E. Guidez, S. P. Webb, K. Montgomery, H. Yasini, I. F. M. Pedreira, and M. S. Gordon (2021) “Computation of host–guest binding free energies with a new quantum mechanics based mining minima algorithm.” Journal of Chemical Physics 154:104122. DOI:10.1063/5.0040759 PDF
W. Chen, X. Ren, C-E.A. Chang (2018). “Discovery of CDK8/CycC Ligands with a New Virtual Screening Tool.” ChemMedChem. DOI:10.1002/cmdc.201800559 Preprint
W. Chen, L. Sun, Z. Tang, Z.A. Ali, B.M. Wong, and C-E.A. Chang (2018). “An MM and QM Study of Biomimetic Catalysis of Diels-Alder Reactions Using Cyclodextrins.” Catalysts 8(12), 51. DOI:10.3390/catal8020051
W. You, Y. M. M. Huang, S. Kizhake, A. Natarajan and C-E. A. Chang (2016) “Characterization of Promiscuous Binding of Phosphor Ligands to Breast-Cancer-Gene 1 (BRCA1) C-Terminal (BRCT): Molecular Dynamics, Free Energy, Entropy and Inhibitor Design.” PLoS Comput Biol. 12: e1005057. DOI:10.1371/journal.pcbi.1005057
Huang, Y. M. M., W. Chen, M. J. Potter, and C.-E. Chang (2012). “Insights from Free-Energy Calculations: Protein Conformational Equilibrium, Driving Forces, and Ligand-Binding Modes.” Biophysical Journal 103(2): 342-351. DOI: 10.1016/j.bpj.2012.05.046
Chen, W., M. K. Gilson, S. P. Webb, and M. J. Potter (2010). “Modeling Protein-Ligand Binding by Mining Minima.” Journal of Chemical Theory and Computation 6(11): 3540-3557. DOI: 10.1021/ct100245n
VM2 background references
Chang, C.-E., M. J. Potter, and M. K. Gilson (2003). “Calculation of molecular configuration integrals.” J. Phys. Chem. B 107: 1048-1055. DOI: 10.1021/jp027149c
Chang, C.-E. and M. K. Gilson (2003). “Tork: Conformational analysis method for molecules and complexes.” Journal of Computational Chemistry 24(16): 1987-1998. DOI: 10.1002/jcc.10325
Chen, W., C.-E. Chang, and M. K. Gilson ( 2004). “Calculation of Cyclodextrin Binding Affinities: Energy, Entropy, and Implications for Drug Design ” Biophysical Journal 87: 3035-3049. DOI: 10.1529/biophysj.104.049494
Gilson, M. K., J. A. Given, B. L. Bush, and A. McCammon (1997). “The statistical-thermodynamic basis for computation of binding affinities: A critical review.” Biophys. J. 72: 1047-1069. DOI: 10.1016/S0006-3495(97)78756-3
Gilson, M. K., J. A. Given, and M. S. Head (1997). “A New Class of Models for Computing Receptor-Ligand Binding Affinities.” Chem. & Biol. 4: 87-92. DOI: 10.1016/S1074-5521(97)90251-9
Gilson, M. K., The Binding Data Base (http://www.bindingdb.org), Skaggs School of Pharmacy & Pharmaceutical Sciences (2013).
Head, M. S., J. A. Given, and M. K. Gilson (1997). “Mining Minima: Direct computation of conformational free energy.” J. Phys. Chem. 101: 1609-1618. DOI: 10.1021/jp963817g
Luo, R., L. David, and M. K. Gilson (2002). “Accelerated Poisson-Boltzmann calculations for static and dynamic systems.” Journal of Computational Chemistry 23(13): 1244-1253. DOI: 10.1002/jcc.10120
Potter, M. J. and M. K. Gilson (2002). “Coordinate Systems and the Calculation of Molecular Properties.” The journal of physical chemistry. A 106(3): 563-566. DOI: 10.1021/jp0135407
Qiu, D., P. S. Shenkin, F. P. Hollinger, and W. C. Still (1997). “The GB/SA continuum model for solvation. a fast analytical method for the calculation of approximate born radii.” J. Phys. Chem. 101(16): 3005-3014. DOI: 10.1021/jp961992r
