Speeding up electrostatic computations for molecular dynamics

Molecular dynamics (MD) simulations are routinely used to study the structure and function of biological molecules. However the accuracy and duration of these simulations are constrained by their computational costs, thus limiting the ability to accurately simulate systems of realistic sizes over biologically relevant time periods. The two most computationally demanding steps in these simulations are (1) determining the charge state of ionizable sites in biomolecules, which is a key input to the simulation, and (2) calculating long range electrostatic interactions during the simulation. Presented here are two novel methods, the direct interaction approximation (DIA) and the hierarchical charge partitioning (HCP) approximation, for speeding up each of these two computations.

The average charge state of ionizable sites in biomolecules can be calculated as the statistical average over all possible (2^N) microstates for a molecule, where N is the number of ionizable sites. In general this computation scales exponentially as O(N² 2^N). The DIA is an O(²) approximation for calculating the average charge state of ionizable sites. For each site, the DIA treats direct interactions (interactions involving the site of interest) exactly, while using an average value for indirect interactions (interactions not involving the site of interest). The DIA was tested on two problems. The computation of thermal average properties for the 2-D Ising model of ferromagnetism, and the average charge state of ionizable residues in biomolecules. Compared to the commonly used non-deterministic Monte Carlo method, for the same computational cost, the deterministic DIA was found to be at least as accurate, as measured by RMS error relative to the exact computation. Thus, the DIA may be a practical alternative to the Monte Carlo method for some problems.

In atomistic MD simulations, the computation of long range electrostatic interactions, scale as O(n²), where n is the number of atoms. For most biologically relevant timescales the simulations involve 10^12–16 simulation steps. Thus, the computational cost of long range interactions become the limiting factor in the size and duration of MD simulations. The HCP is an O(n log n) approximation for computing long range electrostatic interactions. The approximation is based on multiple levels of natural partitioning of biomolecular structures into a hierarchical set of components. For components that are far from the point of interest, the charge distribution for each component is approximated by a much smaller number of charges. For nearby components, the HCP uses the full set of atomic charges. For large structures the HCP can be several orders of magnitude faster than the exact pairwise O(n²) all-atom computation. For a representative set of structures, the accuracy of the HCP is comparable to the industry standard explicit solvent particle mesh Ewald (PME), and is in general more accurate than the spherical cutoff method. And, unlike the PME, the DIA can be easily extended to implicit solvent GB models. 50 ns implicit solvent simulations for a representative set of four biomolecules suggests that the HCP could be a practical alternative for implicit solvent simulations, and preferable to the cutoff based method. The HCP is available for general use in the open source MD software, NAB within AmberTools.