Recent Developments in Free Energy Calculations for Drug Discovery

Abstract

The grand challenge in structure-based drug design is achieving accurate prediction of binding free energies. Molecular dynamics (MD) simulations enable modeling of conformational changes critical to the binding process, leading to calculation of thermodynamic quantities involved in estimation of binding affinities. With recent advancements in computing capability and predictive accuracy, MD based virtual screening has progressed from the domain of theoretical attempts to real application in drug development. Approaches including the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA), Linear Interaction Energy (LIE), and alchemical methods have been broadly applied to model molecular recognition for drug discovery and lead optimization. Here we review the varied methodology of these approaches, developments enhancing simulation efficiency and reliability, remaining challenges hindering predictive performance, and applications to problems in the fields of medicine and biochemistry.

Keywords: MM-PBSA; alchemical simulation; binding affinity; drug discovery; free energy simulation; lie; molecular dynamics.

FIGURE 1

Citation counts for each method over the past 20 years. The development and utilization of molecular simulation to guide drug discovery has grown dramatically in recent years. The MM-PBSA method, which balances simulation rigor, high speed, and minimal setup complexity to allow high throughput screening, has seen extensive application reaching over 2,000 citations in 2020. Steep computational costs and challenges in generalizing protocols to work on broad sets of protein-ligand systems have limited the usage of absolute alchemical and LIE based approaches.

FIGURE 2

MM-PBSA thermodynamic cycle. The binding free energy in aqueous environment is calculated as the difference between the sum of binding in vacuum and solvating the complex with solvating the receptor and ligand individually. The information necessary to complete this cycle can be obtained by decomposing a single trajectory into the ensemble desolvated receptor, ligand, and complex configurations, and computing the solvation free energies for each state with the Poisson-Boltzmann equation. Normal mode analysis can be performed to determine the contribution of entropy to the binding process.

FIGURE 3

LIE binding free energy calculation. The binding free energy is computed from force field energy estimates of the differences in van der Waals and electrostatic energies for the ligand bound to the protein and free in solvent environment. The system dependent LIE parameters α and β are empirically determined and used to scale the non-polar and coulombic interaction energies to have minimal error with respect to available experimental data. The final term γ acts as an optional offset parameter to further tune the model. LIE requires no post-processing and can be completed from a single trajectory.

FIGURE 4

Absolute alchemical simulation thermodynamic cycle. Two trajectories are completed to model the unbinding process. The simulations start from the complex of protein-ligand bound and end with receptor and unbound ligand (top track), and from ligand alone in solvent to ligand removed (bottom track). The ligand is transformed through a series of unphysical states to decouple electrostatic and van der Waals interactions with the surrounding environment to reach the final state where it no longer interacts with the initial system. The binding free energy prediction is the sum of the coulombic and non-polar energies involved in the transformation eliminating protein-ligand interactions. A restraint is typically included to prevent the ligand from exiting the active site while the binding interactions keeping the protein and ligand together are scaled off in order to aid convergence, this is corrected for with an additional transformation progressively turning on the restraints for the complex track and an analytical correction for the ligand track.