In this example, we show how to set up and run an NMR refinement for a 10-base-pair DNA duplex, starting from LEaP through the full simulated annealing. Unlike the examples above, this uses explicit solvent, which should be an increasingly popular way of carrying out refinements. The files here were created by Vickie Tsui, based on work she has done in collaboration with Jarrod Smith and Dave Case at TSRI. Note that one "disadvantage" of explicit water simulations is that a significant effort must be made to prepare and equilibrate the sample. These steps would likely only be used after a conventional (vacuum or generalized Born) refinement had been done to get reasonable starting structures.
We start with a script that can be fed to tleap or xleap:
+---------------------------------------------------------------------------------+ | Running LEaP | +---------------------------------------------------------------------------------+ |gcg=loadpdb gcg.pdb (read in starting structure) | |addions gcg Na+ 0 (neutralize system with sodium counterions) | |solvatebox gcg WATBOX 216 9.0 0.8 (put in explicit water) | |saveamberparm gcg prmtop gcg_start.x (save topology and AMBER coordinate files) | +---------------------------------------------------------------------------------+
Note that the input file, gcg.pdb must (of course) have the proper
sequence, but it does not need to have hydrogens, etc. It should be
possible to use a file directly from Brookhaven: LEaP will take care of
conversions to Amber nomenclature, adding hydrogens, etc. As always,
though, consider carefully any errors or warnings that LEaP gives.
Next we carry out a sander minimization, optimizing just the waters
and counterions, and keeping the solute and fixed:
+---------------------------------------------------------------------------------+ | Minimize water and ions | +---------------------------------------------------------------------------------+ | | | Minimization protocol for water molecules and conterions only | | | | &cntrl | | imin=1, (do minimization) | | ncyc=50, maxcyc=1000, (50 steps steepest descent, then conjugate gradient) | | nmropt=1, (distance and angle restraints) | | ntc=2, tol=0.000001, (SHAKE) | | cut=9.0, (cutoff for long-ranged interactions) | | ntpr=100, (periodic printout) | | ntb=1, (constant volume) | | ntr=1, (MD with restraint of atoms specified) | | &end | | &wt type='REST', istep1=0,istep2=1000,value1=1.0, (NMR restraints) | | value2=1.0, &end | | &wt type='END' &end | |LISTOUT=POUT | |DISANG=RST | |Group input for restrained atoms (Harmonic restraints on solute coordinates) | | 10.0 | |RES 1 20 | |END | |END | +---------------------------------------------------------------------------------+
The NMR-based restraints would be in the RST file. We do not use PME
here, since this is just the beginning of an equilibration. The next step
illustrates an initial MD run:
+-----------------------------------------------------------------+ | Initial equilibration | +-----------------------------------------------------------------+ | | | Heating up from 0 K to 300 K | | | | &cntrl | | imin=0, (do MD) | | nmropt=1, | | ntx=1, irest=0, (format of input) | | cut=9.0, | | ntt=1, temp0=300.0,tautp=1.0, (temperature specifications) | | ntb=1, ntc=2, tol=0.000001, | | ntpr=100, ntwr=99999, | | ntr=1, | | nstlim=30000, (30000 steps of MD) | | &end | | &wt type='REST', istep1=0,istep2=1000,value1=1.0, | | value2=1.0, &end | | &wt type='END' &end | |LISTOUT=POUT | |DISANG=RST | |Group input for restrained atoms | | 5.0 | |RES 1 20 | |END | |END | +-----------------------------------------------------------------+
Next we continue the equilibration with constant pressure, to get the density to be about right. Note that we are still keeping the starting structure of the DNA nearly fixed, so this still is part of the equilibration period fo the water
+--------------------------------------------------------+ | Constant pressure equilibration | +--------------------------------------------------------+ | | | Constant pressure simulation | | | | &cntrl | | imin=0, nmropt=1, ntx=7, irest=1, | | scee=1.2, cut=9.0, | | ntt=1, tempi=300., temp0=300.0, | | ntb=2, (constant pressure) | | ntp=2, (pressure regulation) | | ntc=2, tol=0.000001, | | ntpr=100, ntwr=99999, | | ntcm=1, nscm=1000, ndfmin=-3, | | ntr=1, | | nstlim=20000, | | &end | | &wt type='REST', istep1=0,istep2=1000,value1=1.0, | | value2=1.0, &end | | &wt type='END' &end | |LISTOUT=POUT | |DISANG=RST | |Group input for restrained atoms | | 5.0 | |RES 1 20 | |END | |END | +--------------------------------------------------------+
Here is a final equilibration step for the water and counterions:
+-----------------------------------------------------------------------------+ | More constant pressure equilibration | +-----------------------------------------------------------------------------+ | | | Constant pressure with larger tautp -- make sure the density has stabilized | | | | &cntrl | | imin=0, nmropt=1, irest=1, ntx=7, | | scee=1.2, cut=9.0, | | ntt=1, tempi=300., temp0=300., | | ntb=2, ntp=2, tautp=2.0, (specify tautp) | | ntc=2, tol=0.000001, | | ntpr=100, ntwr=999999, | | ntwx=1000, (start recording trajectory every 1ps) | | ntr=1, | | nstlim=40000, | | &end | | &wt type='REST', istep1=0,istep2=1000,value1=1.0, | | value2=1.0, &end | | &wt type='END' &end | |LISTOUT=POUT | |DISANG=RST | |Group input for restrained atoms | | 5.0 | |RES 1 20 | |END | |END | +-----------------------------------------------------------------------------+
Next, we return to constant volume MD (although this is mostly a matter of taste--at constant volume, energy conservation can be monitored). We also reduce the weight of the constraints keeping the DNA tethered to the intial conditions.
+--------------------------------------------------------+ | Final equilibration | +--------------------------------------------------------+ | | | Decrease harmonic restraint of solute | | | | &cntrl | | imin=0, nmropt=1, irest=1, ntx=7, | | scee=1.2, cut=9.0, | | ntt=1, tempi=300., temp0=300., | | ntb=1, tautp=2.0, | | ntc=2, tol=0.000001, | | ntpr=100, ntwr=999999, ntwx=1000, | | ntr=1, | | nstlim=50000, | | &end | | &wt type='REST', istep1=0,istep2=1000,value1=1.0, | | value2=1.0, &end | | &wt type='END' &end | |LISTOUT=POUT | |DISANG=RST | |Group input for restrained atoms | | 0.5 (lowered force constant) | |RES 1 20 | |END | |END | +--------------------------------------------------------+
Finally, we are ready for a "production" run, with no artificial constraints. This step could also include some (mild) heating and cooling for simulated annealing, as in the example above.
+---------------------------------------------------------+ | Production dynamics | +---------------------------------------------------------+ | | | Production run | | | | &cntrl | | imin=0, nmropt=1, irest=1, ntx=7, | | scee=1.2, cut=9.0, | | ntt=1, tempi=300., temp0=300., | | ntb=1, tautp=2.0, | | ntc=2, tol=0.000001, | | ntpr=100, ntwr=999999, ntwx=1000, | | nstlim=200000, | | &end | | &wt type='REST', istep1=0,istep2=1000,value1=1.0, | | value2=1.0, &end | | &wt type='END' &end | |LISTOUT=POUT | |DISANG=RST | +---------------------------------------------------------+
The techniques illustrated above are really fairly general, and illustrate one (of many) ways to prepare a system for solvated simulation. The ideas apply equally well to an NMR refinement as to a "regular" simulation. A key point is that it takes a long time to equilibrate systems, and that it is typical (but not required) to equilibrate the solvent first, and then the whole system: this decreases that chances that bad initial solvent coordinates will disrupt the solute structure at early stages of the process.