New theoretical techniques are being developed and characterized. These efforts are usually coupled with software development, and involve the systematic testing and evaluation of new ideas. This development is driven by current needs and interests. Comparative normal mode analysis at different model resolutions: A tool for multiscale modeling As part of the effort to develop generalized approaches to multi-scale simulations, new analysis techniques must be developed to interpret the results of multi-scale normal mode analysis calculations. Validating the results of these calculations is a difficult problem because both the number of modes and frequencies and the lengths of the individual mode vectors differ from those produced by all-atom or more traditional coarse-grained NMA calculations (e.g. via the Elastic Network Model). Various approaches to this problems such as measuring spatial extents of modes or mapping lower level modes on a higher level representation are being studied and benchmarked. pH Replica Exchange Method Internal ionizable groups are essential for biological proton transfer. Because of the coupling between protonation and structural rearrangements, prediction of pKa values of internal ionizable groups can be challenging. Molecular dynamics at constant pH (PHMD) allows for the coupling between ionization and protein dynamics and presents potentially an excellent tool for prediction of pKa values. We have developed a new method for calculation of pKa values and titration curves that is a pH replica exchange method based on constant pH molecular dynamics (PHREM). The performance and accuracy of this method and the original PHMD method were compared. The PHREM algorithm exhibits better sampling of protonation states than the PHMD method. In particular the pKa values converge faster and the accuracy of predicted Hill coefficients is significantly improved. Direct simulation of transition path of rare events In computational chemistry, the reaction mechanisms of rare events can be represented as a reaction path on the potential energy surface (PES) of the system under study. By collaborating with Dr. Milan Hodoscek, we have been working on the methods to directly simulate the rare events in transition path space rather than phase space. The direct reaction path simulation can generate an ensemble of transition paths, which can be subject to statistical mechanics analysis. In addition the simulations may also reveal multiple reaction pathways, which cannot be revealed easily by applying reaction path optimization. Evaluation of the accuracy of relative free energy simulations for drug design Free energy simulations are among the most accurate and general methodologies in the field of computational chemistry. They provide means to study diverse processes such as the binding affinities of ligands or novel drugs, enzymatic reactions, the solvation of organic molecules, as well as the effect of point mutations. However, most applications of free energy simulations do not provide reliable data on the accuracy of the method, since either a.) there is no reference data for an assessment of the quality of the method or b.) the free energy calculations were conducted after the experimental reference results became available leading to biased results. To evaluate the accuracy of several free energy methods, relative free energy calculations based on molecular dynamics simulations were combined with available experimental binding free energies to the predict unknown binding affinities of host-guest systems in the blind SAMPL3 competition. Our predictions showed very good agreement with experimental results (considering the size of the molecules), yielding root mean square errors of about 2.6 kcal/mol for seven host-guest systems ranking among the top 3 groups in this competition. We compared the performance of three different approaches: Bennett's Acceptance Ratio Method and Thermodynamic Integration based on both the trapezoidal and Simpson's rule. In addition, we also evaluated the influence of the protonation states of the amine groups of the guest molecules, as well as the the buffer concentration on the binding affinity. Crossing energy barriers with Self-guided Langevin Dynamics The performance of Self-Guided Langevin dynamics (SGLD) was evaluated in terms of speedup and deviations from the original ensemble for multiple benchmark systems, including the proline cis-trans isomerization (involving energy barriers of 14-20 kcal/mol) and five-atomic benchmark systems where analytical results are available. While normal MD simulations are able to cross 6-7 kcal/mol, our results demonstrate that SGLD with the right parameters can overcome energy barriers of up to 20 kcal/mol. Thus, the effective speed up simulations can be improved by several orders of magnitude, while causing only slight deviations from the original ensemble (up to 7%). For example, while the cis-trans isomerization of proline never occurs in a normal simulation of 0.5 microseconds, SGLD can lead to up to 10000 transitions in the same amount of time. Since cis-prolines are often an important factor for the formation of turns in the structure of proteins, this speedup might free the path to more efficient fold predictions. Calculating free energy differences based on Self-guided Langevin Dynamics Self-Guided Langevin Dynamics (SGLD) were combined with the Non-Boltzmann Bennett's Acceptance Ratio method (NBB) to make free energy calculations with the biased trajectories of SGLD simulations possible. We have tested this approach for five-atomic benchmark systems where normal molecular dynamics simulations fail to give correct results (being unable to cross the energy barriers of about 10 kcal/mol). While normal Langevin Dynamics in combination with Thermodynamic Integration and 41 lambda intermediate states results in deviations of 2.6 kcal/mol (or 20%) from the analytical results, SGLD with only two lambda states is able to yield the correct results with a deviation of just 0.07 kcal/mol (or 0.4%). We are also working on studying the statistical, both equilibrium and non-equilibrium, behavior of the SGLD algorithm. The ability to use SGLD to calculate statistics of a macromolecular system is dependent upon the nature of the algorithm's sampling. The two important aspects are the rate of sampling of unique states and the accuracy of the sampling statistics. We are using both numerical and analytical techniques to characterize these two properties of SGLD applied to simple test systems. The goal of this work is to give insight to proper application of SGLD in practice. Excited states methods The time-dependent density functional theory (TDDFT) is the method of choice for studying excited states in biological systems. We achieved an efficient parallel implementation of the analytical gradient for TDDFT. We are applying to the studies of a) organic dye molecules such as polycarbazole; b) firefly bioluminescence; c) flurorescence quenching of tryptophan in water and protein environments. Combining Conformational Space Annealing with Replica Exchange Method for Improved Conformational Sampling Reservoir Replica Exchange Method (R-REM) has been combined with the Conformational Space Annealing (CSA) method has been shown to be able to determine the global energy minimum of proteins efficiently and has been used in structure prediction successfully. CSA uses a genetic algorithm approach to generate a diverse set of conformations to determine the minimum energy structure. We have used conformations generated through CSA method to build a reservoir. Replica exchange was then performed where the top replica was seeded with the reservoir structures and fast convergence at every temperature is observed.

Project Start
Project End
Budget Start
Budget End
Support Year
Fiscal Year
Total Cost
Indirect Cost
National Heart, Lung, and Blood Institute
Zip Code
Lee, Juyong; Lee, In-Ho; Joung, InSuk et al. (2017) Finding multiple reaction pathways via global optimization of action. Nat Commun 8:15443
Huang, Jing; Mei, Ye; König, Gerhard et al. (2017) An Estimation of Hybrid Quantum Mechanical Molecular Mechanical Polarization Energies for Small Molecules Using Polarizable Force-Field Approaches. J Chem Theory Comput 13:679-695
Wu, Xiongwu; Lee, Juyong; Brooks, Bernard R (2017) Origin of pKa Shifts of Internal Lysine Residues in SNase Studied Via Equal-Molar VMMS Simulations in Explicit Water. J Phys Chem B 121:3318-3330
Li, Ying; Li, Hui; Pickard 4th, Frank C et al. (2017) Machine Learning Force Field Parameters from Ab Initio Data. J Chem Theory Comput 13:4492-4503
Simmonett, Andrew C; Pickard 4th, Frank C; Ponder, Jay W et al. (2016) An empirical extrapolation scheme for efficient treatment of induced dipoles. J Chem Phys 145:164101
Pickard 4th, Frank C; König, Gerhard; Simmonett, Andrew C et al. (2016) An efficient protocol for obtaining accurate hydration free energies using quantum chemistry and reweighting from molecular dynamics simulations. Bioorg Med Chem 24:4988-4997
Jia, Xiangyu; Wang, Meiting; Shao, Yihan et al. (2016) Calculations of Solvation Free Energy through Energy Reweighting from Molecular Mechanics to Quantum Mechanics. J Chem Theory Comput 12:499-511
Pickard 4th, Frank C; König, Gerhard; Tofoleanu, Florentina et al. (2016) Blind prediction of distribution in the SAMPL5 challenge with QM based protomer and pK a corrections. J Comput Aided Mol Des 30:1087-1100
Dybeck, Eric C; König, Gerhard; Brooks, Bernard R et al. (2016) Comparison of Methods To Reweight from Classical Molecular Simulations to QM/MM Potentials. J Chem Theory Comput 12:1466-80
Wu, Xiongwu; Pickard 4th, Frank C; Brooks, Bernard R (2016) Isotropic periodic sum for multipole interactions and a vector relation for calculation of the Cartesian multipole tensor. J Chem Phys 145:164110

Showing the most recent 10 out of 53 publications