Free Energy Calculation For Biomolecules: Accuracy And Precision

Computer simulation has been widely used to study biochemical processes like solvation of small molecules, drug-target binding and protein folding. Free energy calculation is a key challenge in computational biochemistry. Successful free energy prediction depends on three basic components, that is, the suitable Hamiltonian, the reliable estimator and enough sampling. In this thesis, I will explore how to refine and improve these three components to fulfill precise and accurate free energy calculation.In chapter two, the semi-empirical linear interaction energy method is ap-plied to calculate the binding affinities for avidin and β-secret ase binding com-plexes. In molecular simulation, the polarized protein-specific charge (PPC) scheme instead of the conventional mean-field charge is employed to give a reli-able description of the Coulomb interaction between ligand and its environment. By comparison, a remarkable consistency between theoretical prediction and ex-perimental measurement has been observed. However, evident limitations exist in linear interaction method. For the avidin and β-secretase binding complexes, the fitted parameters are different from each other. Therefore, linear interaction method possesses no transferability of parameters. Furthermore, for systems without known experimental measurements, linear interaction energy method is also powerless. During the derivation some assumptions are taken, so this method is semi-empirical. Most importantly, no matter for conventional mean-field charges or PPC, the wave function distorsion energy which is very important in biosystems is neglected because of the usage of fixed charge model.In chapter three, an efficient approach that combines the linear scaling quan-tum mechanical method (EE-GMFCC) with conductor-like polarizable contin-uum model (CPCM), termed EE-GMFCC-CPCM, is developed for ab initio cal-culation of the electrostatic solvation energy of proteins. Compared with the previous MFCC-CPCM method proposed by Prof. Ye Mei in 2006, quantum mechanical calculation is applied to deal with short-range non-neighboring inter-actions replacing the classical treatment. Numerical studies are carried out for proteins up to 3837 atoms at the HF/6-31G* level. As compared to standard full system calculation, EE-GMFCC-CPCM shows clear improvement over the MFCC-CPCM method for both the total electrostatic solvation energy and its components (the polarized solute-solvent reaction field energy and wave function distorsion energy of the solute). For large proteins with 1000-4000 atoms, where the standard full system ab initio CPCM calculations are not affordable, the EE-GMFCC-CPCM gives larger relative wave function distorsion energies and weaker relative electrostatic solvation energies of proteins, as compared to the corresponding energies calculated by the Divide-and-Conquer Poisson-Boltzmann (D&C-PB) method. Notwithstanding, a high correlation between EE-GMFCC-CPCM and D&C-PB is observed. This study demonstrates that the linear-scaling EE-GMFCC-CPCM approach is an accurate and also efficient method for the cal-culation of electrostatic solvation energy of proteins. Some points to be noted: (1) In practical calculation of solvation free energy, the non-electrostatic contri-bution should also be considered, which is generally expected to be related with the solvent accessible surface of solute; (2) It is better to consider the ensemble average as the result rather than just using single conformation; (3) Detailed atomic information on the solvent is missing from the used implicit model, espe-cially for the first shell of solvent molecules which might be crucial to the solute property.In chapter four, the reweighting scheme is used to calculate the solvation free energy of organic molecules. Calculations of the free energy difference between two states at quantum mechanical (QM) level directly are generally prohibitively expensive, because some intermediate states are usually required in order to increase the overlap in phase space between two adjacent states. Fortunately, free energy is a state function, of which the difference is path-independent. Therefore, in a more practical way the intermediate states can be described by molecular mechanics (MM). The free energy difference between two adjacent states can be estimated using the Bennett Acceptance Ratio (BAR), which has been shown to give the minimum variance with fixed number of samples. To further reduce the computational expense by refraining from sampling at QM level for the end states, Boresch proposed a method termed Non-Boltzmann Bennett’s acceptance ratio (NBB), which combines BAR and energy reweighting. However, we showed that the most efficient way to calculate the free energy difference between two QM states is BAR+TP, which gives the minimum variance of the results. In this scheme, the free energy differences between MM states are estimated using BAR, and at both ends a thermodynamic perturbation from MM to QM states are applied. In addition, the QM expense in this scheme is only half of that in NBB. We also show that defining the biasing potential as the difference of the solute-solvent interaction energy, instead of the total energy, can converge to the calculated solvation free energies much faster but possibly to different values. It has also been discovered in this study that BLYP yields better results than MP2 and some later functionals such as B3LYP, M06-2X, and wB97X-D.In chapter five, the reweighting scheme introduced in chapter four is applied to calculate the absolute protein-ligand binding free energy. Firstly, the double decoupling method is applied to obtain the binding free energy on MM potential. Then, the QM correction of the protein-ligand electrostatic interaction energy is carried out by TP method. For L99A mutant of T4 lysozyme whose conformation changes little upon binding, its absolute binding free energies with benzene and phenol are predicted successfully after considering the QM correction. However, for the lectin-monosaccharide complex, evident discrepancy between theoretical prediction and experimental data is found regardless of whether or not the QM correction is considered. Lectin experiences large conformational changes during binding process, although QM correction can give a better description of the Hamiltonian, longer simulations or some enhanced sampling methods are required to collect enough conformations.