The VM2 method, recently developed by VeraChem LLC, calculates protein-ligand binding affinities. It has already been shown to predict free energies of ligand binding that correlate well with experiment for a number of different protein receptors.

It also provides a unique level of insight into the physical basis of the binding free energy in terms of enthalpic, entropic, and solvation contributions. The VM2 method therefore has outstanding potential as a basic research tool as well as for computer-aided discovery of new therapeutics.

If you are interested in obtaining VM2 contact us at vc@verachem.com.

Key Strengths

  • Yields true free energies of binding
  • Ability to rank ligand series in order of binding strength
  • Flexible active site and fully flexible ligand
  • Highly effective conformational searching
  • Accurate Poisson Boltzmann continuum solvation
  • Configurational entropy rigorously included
  • Detailed free energy decomposition to aid ligand design

Theory: 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 configurational integrals associated with low energy minima of the system. For example, the following 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.

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.


Potential Energy and Solvation Energy

Currently the potential energy and solvation terms required in the local configuration integrals are calculated classically. The potential energy U through standard force fields


and the solvation energy through continuum solvent models. For geometry optimizations and force constant calculations the Generalized Born method is used. The final solvent energies used (W) are corrected with the Poisson Boltzmann Surface Area (PBSA) method.


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 of the ligand inside the binding pocket. 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.

Mode distort-minimize 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. Again 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, and so on. 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.

Calculation of Local Configuration Integrals

The Harmonic approximation is used to calculate the local configuration integrals for each minimum produced.


To account for anharmonicty, which is often present for low energy modes, a numerical integral is calculated and replaces the analytic Harmonic one. An energy profile is determined along each low energy mode and integration is carried out by the trapezoid rule.


Summary of Energy Terms Accounted for by VM2


Importance of including Entropy in Binding Affinity Calculations

The following plot demonstrates a trend of entropy-enthalpy compensation in experimental binding data. 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.



Main Citations

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

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

Background Citations

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