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:


/////////////////////////////////////////////////////////////////////////////

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].

-----------------------------------------------------------------------------

(c) V.Babin and C.Sagui 11/27/2009

/////////////////////////////////////////////////////////////////////////////