Conformational equilibria of methyl-alpha-L-iduronic acid in explicit water

Pyranose ring puckering

Note:This is intentionally an advanced tutorial, showing one way of using of an approximate PMF that describes slow modes of a small molecule to setup a replica exchange simulation of that small molecule in explicit solvent. The outcome is twofold: (1) enchanced sampling of the solute degrees of freedom; and (2) accurate PMF of the slow modes involved into the setup. You will want to make reference to the following two items:

V. Babin and C. Sagui. Conformational free energies of methyl-alpha-L-iduronic and methyl-beta-D-glucuronic acids in water. J. Chem. Phys. 132, 104018 (2010). Here is a preprint, which describes what is being computed, and why and how;
The tar file (in gzipped form) that contains the inputs and outputs.

The aim of this document is to present a simulation protocol that can be used to efficiently sample puckering configurations of the methyl-$\alpha$-L-iduronic acid in the explicit solvent bath. The methodology is explained in this paper: V. Babin and C. Sagui. Conformational free energies of methyl-alpha-L-iduronic and methyl-beta-D-glucuronic acids in water. J. Chem. Phys. 132, 104018 (2010) The purpose of this document is to describe the technicalities. To make a long story short: the idea is to setup a replica-exchange scheme by means of biasing of the system along its slowest, difficult to sample, modes. In general case the latter are not so easy to identify and capture, yet for small compounds it is often possible to identify a suitable ``collective variable'' (or a ``reaction coordinate'') that provides a ``good enough'' description of the conformational landscape. It is then possible to use the corresponding potential of mean force (PMF), known also as ``Landau free energy'', to setup an efficient replica exchange simulation. For the molecule considered here the slow modes in question are related to the changes of the pyranose ring puckering. These conformational changes can be described in terms of the intraring dihedral angles O5--C1--C2--C3 and C3--C4--C5--O5 (see Fig.2 in the preprint). This is not the only possibility, obviously. As explained in the preprint, we proceed in stages (the step numbers used below coincide with the folder names in the archive available here.

0. Preliminaries

We assume that the reader is familiar with basic AMBER concepts and do not repeat them here [that is, we do not explain how to create the topology file, initial coordinates, etc].

This step is a short, regular molecular dynamics run at room temperature in implicit solvent. The goal here is to examine the magnitude and characteristic time scale of the thermal fluctuations of the O5--C1--C2--C3 and C3--C4--C5--O5 dihedrals. The MDIN file can be found in the ncsu-hremd-tutorial/0/input folder of the supplementary files and basically is nothing significantly beyond a typical GB/SA run. To get the dihedrals we could dump the trajectory and use ptraj afterwards, yet in our case the dihedrals are readily available through the ncsu_abmd subsystem ran in the ANALYSIS mode. More precisely, the text just after the &cntrl namelist instructs sander to compute the values of the dihedrals every 100 steps and dump them into the ``monitor.txt'' file. For the ncsu_abmd ran in the ANALYSIS mode the latter simply contains MD time (in picoseconds) in the first column followed the values of the collective variables in the subsequent columns.

Once the run is complete, one can visualize the time series of the dihedrals and conclude that the magnitude of the fluctuations is about 0.2 radian and the pertinent time is about 0.5 ps (these values are obviously very sensitive to the definitions; we used the old, proven by time ``stare-and-pick-a-value'' approach which consists of looking at the plots and trying to come up with some number).

1. Implicit solvent

As explained in the preprint, the goal of this step is to find an (approximate) potential of mean force (PMF, or Landau free energy) of the O5--C1--C2--C3 and C3--C4--C5--O5 angles. We could have skipped this part and proceeded with explicit water from the beginning, yet it is reasonable to expect a great degree of similarity between the PMFs computed in implicit versus explicit water models. The idea therefore is to get initial approximation using implicit water and use it as the starting point for more demanding explicit water computations.

We proceed in two stages: coarse "flooding", followed by finer "flooding". The "flooding" technically amounts to changing the "mode" parameter in the ncsu_abmd section of the MDIN file, and defining the scales: the time scale, that controls the speed of the biasing potential buildup, and the resolutions for the collective variables (defined separately for each of them). There are no strong, rigorous arguments behind the choice or the parameters. On the other hand, as we argue more in the preprint, the scales observed in a short equilibrium run appear to give a reasonable first guess.

The values of the dihedrals, along with the value of biasing potential, again get dumped to the "monitor.txt" file. Once all relevant minima are visited (which in this case are: two chairs and several boats), there is no point to continue the "coarse" simulation further and we proceed to the "finer" run. A larger value is used for the time scale in this sub-step and hence the dynamics is closer to equilibrium (disturbed less by the "growing" biasing potential). The outcome is (supposedly) more accuarte PMF (even more accurate PMF can be computed by performing an equilibrium simulation biased by the last non-equilibrium PMF; this is not necessary for our needs, however).

2. Explicit solvent

Our ultimate goal is to have a "good enough" PMF for the dihedrals in explicit water and use it for the replica exchange scheme. As we alluded earlier, it is not unreasonable to expect this "explicit PMF" to be similar to the one computed using an implicit solvent model. We therefore use the "fine" PMF from the previous step as our starting point for the explicit water computations.

We do "flooding" first in order to get closer to the true PMF. We use the values of the parameters from the fine run in implicit water. As before, we stop the simulation once everything relevant has been visited.

At that stage we hope to have a reasonable approximation of the true PMF. In order to assess the quaility of the approximation, we do a relatively brief (in terms of the statistics gathered) biased run at equilibrium (umbrella subdirectory in the data folder). The approximation appears to be good enough as the dihedrals sweep all relevant values (please see the preprint for the more detailed presentation).

3. (Hamiltonian) replica exchange in explicit solvent

Having the (approximate) PMF for the intraring dihedrals, it is easy to setup a (Hamiltonian) replica exchange simulation that exploits the PMF to facilitate the barriers crossing (much like the high-temperature replicas in regular parallel tempering method do).

There are two folders holding the 'equilibration' and 'production' files respectively. The latter is way too short as for the real "production" (it is 3 x 1ns while something like 3 x 100ns would be more appropriate).

We assume that the reader is familiar with doing replica-exchange simulations using sander: basically the 'groups' file is expected to list command line arguments to be parsed by different instances of sander; every instance (replica) has its own MDIN file (as well as MDOUT, MDCRD, etc). Compared to the previous step, the ncsu_abmd section is replaced with ncsu_bbmd; this triggers a slightly different code path that does actually perform the "exchange" part. The replica with rank 0 (first line in the 'groups' file) has some additional parameters in the beginning of the ncsu_bbmd section: exchange attempts frequency, name of the file to which the exchange statistics gets dumped occasionally, and the name of the file that stores the state of the pseudo-random generator used for the "exchange" part (mt19937) [if the latter does not exist, the generator can be initialized with a seed; please consult the manual for the full reference].