Constant-pH Molecular Dynamics Simulations for Large Biomolecular Systems

Abstract

An increasingly important endeavor is to develop computational strategies that enable molecular dynamics (MD) simulations of biomolecular systems with spontaneous changes in protonation states under conditions of constant pH. The present work describes our efforts to implement the powerful constant-pH MD simulation method, based on a hybrid nonequilibrium MD/Monte Carlo (neMD/MC) technique within the highly scalable program NAMD. The constant-pH hybrid neMD/MC method has several appealing features; it samples the correct semigrand canonical ensemble rigorously, the computational cost increases linearly with the number of titratable sites, and it is applicable to explicit solvent simulations. The present implementation of the constant-pH hybrid neMD/MC in NAMD is designed to handle a wide range of biomolecular systems with no constraints on the choice of force field. Furthermore, the sampling efficiency can be adaptively improved on-the-fly by adjusting algorithmic parameters during the simulation. Illustrative examples emphasizing medium- and large-scale applications on next-generation supercomputing architectures are provided.

Figure 1

The protonation state of a titratable system is completely defined by its occupancy vector λ, where each of its m elements is either a one or a zero depending on whether the given site, s, is or is not occupied, respectively. The protonation state of individual residues is determined by a small subset of the elements of λ such that multiple system states may contain the same residue state. The average of a given element of λ yields the protonated fraction for that site and corresponds to a microscopic pK_a. Multiple sites may be equivalent such that a macroscopic pK_a can be determined by grouping two or more sites together (e.g., the neutral states of carboxylate moieties). However, even non-equivalent sites can be grouped into macroscopic transitions, although in these cases the relationship between the two sets of pK_a values is not always straightforward (e.g., methyl imidazole).

Figure 2

The constant-pH MD algorithm consists of two part cycles in which standard equilibrium MD (blue and green solid lines) is performed followed by a driven nonequilibrium switch (orange dotted lines), which changes both the configuration and protonation state (arbitrarily labeled A and B). Detailed balance is restored after the nonequilibrium steps by applying a MC procedure in which the new configuration/state is accepted or rejected.

Figure 3

A switch move contains three main steps: 1) an exact MC sampling of auxilliary sidechain atoms, 2) neMD propagation of the coordinates and coupling constant λ as the original coordinates and state (blue spheres) decouple and the new coordinates and state (green spheres) become coupled, and 3) removal of the non-interacting atoms. If the neMD/MC move is rejected, then the simulation continues from the original coordinates.

Figure 4

Titration curves are easily computed for the reference dipeptides after initial parameterization and subsequent constant-pH simulations. Data points represent explicitly sampled pH values while lines represent the analytic curves predicted by UWHAM.

Figure 5

Translocation of a terminally blocked, titratable pentapeptide, AADAA, across a 1–palmitoyl–2–oleoyl–phosphatidylcholine bilayer was performed by restraining the system at various separations (panel B, top). The insets depict the protonation probability of the central aspartic–acid residue for different positions of pentapeptide along bilayer normal (panel A). Here, z is the Euclidian distance separating the center of mass of the pentapeptide from that of the membrane, projected onto the direction normal to the interface (i.e., z = 0 corresponds to the middle of the ~27 Å thick membrane). The dashed red line in panel B corresponds to the pK_a in bulk water.

Figure 6

Representative macroscopic titration curves (8 of 22 total) for SNase indicate a wide range of pK_a values, even amongst similar residues. Residues are colored by type and have different line patterns to denote the same residue in different environments (in ascending order as solid, dashed, and dotted lines).