Improving the accuracy of molecular dynamics simulations: parameterisation of interaction potentials for small molecules

Abstract

The ability to accurately describe molecular systems with classical molecular dynamics (MD) is critically dependent on an understanding of the sources of error in the design, evaluation and analysis of the underlying model. One significant source of error is in the degree to which the equations and parameters used to represent atomic interaction can reproduce the experimental properties relevant to a particular application. While highly optimised and well-validated interaction parameters are available for common biomolecules (amino acids, sugars, lipids etc.), this is not the case for heterogeneous small molecules such as co-factors, substrates and drugs. The enormity of the chemical space covered by small molecules necessitates the development of automated approaches to parameterising and validating interaction parameters. Several tools are available that can be used to generate interaction parameters for small molecules compatible with a particular existing (bio)molecular force field. Despite their popularity, many are poorly validated, while some have been shown to be inappropriate for many applications. In other cases, validation studies have demonstrated reasonable performance, however, they also suggest significant improvements can be made. One such tool is the Automated Topology Builder (ATB, http://atb.uq.edu.au/) which provides interaction parameters for small molecules compatible with the GROMOS biomolecular force field.As part of this work, a fully automated protocol for the calculation of solvation energy by thermodynamic integration has been developed and applied to the large-scale validation of the ATB against experimental solvation data in water and hexane. It is shown that the largest errors are due to the Lennard-Jones parameters used for functional groups not present in common biomolecules. Other sources of error included the atomic charges assigned by the ATB for the united atom carbons and incompatibilities between GROMOS Lennard-Jones parameters and the ATB charge model.Significant sources of error were also identified in the evaluation of MD models, both for the calculation of solvation energy as well as for MD simulations of lipid membranes. In particular, the use of multiple-time-step algorithms as a time-saving technique. For the calculation of solvation energy by TI, the twin-range cutoff scheme—commonly used in simulations with the GROMOS force field—was found to introduce a systematic error of about 1 kJ/mol. While for lipid systems, various changes made to the integration algorithm of the GROMACS simulation code were found to significantly alter the properties of lipid membranes when simulated with identical force fields.Finally, a method for parameterising transferrable interaction potentials with respect to a wide range of target properties is presented. The method is based on an analysis of surfaces corresponding to the difference between calculated and target data as a function of alternative combinations of parameters. The consideration of surfaces in parameter space as opposed to local values or gradients leads to a better understanding of the relationships between the parameters being optimized and a given set of target data. This in turn enables for a range of target data from multiple molecules to be combined in a robust manner, and for the optimal region of parameter space to be trivially identified. The effectiveness of the approach is illustrated by using the method to refine the chlorine 6-12 Lennard-Jones parameters against experimental solvation energies in water and hexane as well as the density and heat of vaporization of the liquid at atmospheric pressure for a set of 10 aromatic-chloro compounds simultaneously. Single-step perturbation is used to efficiently calculate solvation energies for a wide range of parameter combinations. It was found that good agreement could be obtained with a single set of Lennard-Jones parameters for liquid density, heat of vapourisation and solvation energy in water, while the average solvation energy in hexane was a slight outlier, about 1.5 kJ/mol lower than the experimental values.