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.

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.image_2

VM2 calculates the standard free energy of binding via the standard chemical potentials of the receptor, ligand, and receptor-ligand complex.


image_3VM2 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.

image_6

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).

Random translation/rotation distortions in the binding pocket.

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.image_9Each 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.

image_10The total energy 2nd-derivative (Hessian) matrix with respect to internal coordinates (bond, angle, torsion coordinates) can be calculated

image_11

 

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

 

image_12

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 α

image_13

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.

Pair-mode distortions.
Ligand focused mode distortions.
Receptor focused mode distortions.

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 references

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

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