G ROMACS Tutorial
Umbrella Sampling
Justin Lemkul
Department of Biochemistry,Virginia Tech
This tutorial will guide the user through the process of setting up and running pulling simulations necessa
ry to calculate binding energy between two species. The tutorial assumes the user has already successfully completed the Lysozyme tutorial,some other tutorial,or is otherwise well-versed in basic G ROMACS simulation methods and topology organization.Special attention will be paid to the methods for properly building the system and settings for the pull code itself.
)is derived from the potential of mean force(PMF), The binding energy(∆G
bind
extracted from a series of umbrella sampling simulations.A series of initial configurations is generated,each corresponding to a location wherein the molecule of interest(generally referred to as a"ligand")is harmonically restrained at increasing center-of-mass(COM)distance from a reference molecule using an umbrella biasing potential.This restraint allows the ligand to sample the configurational space in a defined region along a reaction coordinate between it and its reference molecule or binding partner.The windows must allow for slight overlap of the ligand positions for proper reconstruction of the PMF curve. The steps for such a procedure(and the ones utilized in this tutorial)are as follows:
1Generate a series of configurations along a single degree of freedom (reaction coordinate)
2Extract frames from the trajectory in step1that correspond to the desired COM spacing
3Run umbrella sampling simulations on each configuration to restrain it within a window corresponding to the chosen COM distance
4Use the Weighted Histogram Analysis Method(WHAM)to extract the PMF and
calculate∆G
bind
The tutorial assumes that the reader is using G ROMACS version4.0.5or later,to take advantage of bug fixes pertinent to the pull code.The pull code was completely re-written after version3.3.3,such that none of the information contained herein (beyond the basic theory of the technique)is applicable to any G ROMACS version
prior to  4.0.
Start the Tutorial!
Tutorials Home
G ROMACS Tutorial
Step One:Prepare the Topology
Generating a molecular topology for an umbrella sampling simulation is just like any other simulation.Obtain the coordinate file of the structure of interest, and generate the topology from pdb2gmx.Some systems will require special ,protein-ligand complexes,membrane proteins,etc).For protein-ligand systems,please consult this tutorial,and for membrane proteins, I recommend my own tutorial on the topic.The principles of umbrella sampling are easily extendable to these systems,though we will consider only protein molecules in this tutorial.
The system we will consider here is the dissociation of a single peptide from the growing end of an Aβ
42
protofibril,and is based on simulations we recently
published.The structure file of the wild-type Aβ
42
protofibril used in those simulations,acetylated at the N-terminus of each chain,can be found here.The original PDB accession code is2BEG.
Run the structure through pdb2gmx:
pdb2gmx-f input.
Choose the G ROMOS9653A6parameter set,"None"for the N-termini,and"COO-"for the C-termini.Modify topol_B.itp to include the following lines(at the end of the file):
#ifdef POSRES_B
#include"posre_B.itp"
#endif
We will be using chain B as an immobile reference later on in the pulling simulations,hence the need to specially position-restrain this chain only,and none of the others.
PLEASE NOTE:Unfortunately,pdb2gmx does not properly generate the CA-C-O angle that must be written to the topology.Without these parameters,grompp will fail later.You have two options to fix this problem:
1Add a line for this angle in p entry:"CA C O ga_30"in the[angles] directive
2Manually add"ga_30"to the line in the[angles]directive of topol_*.itp where it is missing
Back: Introduction Next:Defining the
Unit Cell
G ROMACS Tutorial
Step Two:Define the Unit Cell
Defining the unit cell for a pulling simulation is not unlike defining the unit cell for any other simulation.There is,however,one major consideration.One must allow enough space in the pulling direc
tion to allow for a continuous pull without interacting with the periodic images of the system.That is,the minimum image convention must be continually satisfied,and as well,the pull distance must always be less than one-half the length of the box vector along which the pulling is being conducted.Why,you may ask?
G ROMACS calculates distances while simultaneously taking periodicity into account.This,if you have a 10-nm box,and you pull over a distance greater than 5.0nm,the periodic distance becomes the reference distance for the pulling,and this distance is actually less than 5.0nm!This fact will significantly affect results,since the distance you think you are pulling is not what is actually calculated.We will be pulling a total distance of    5.0nm in a 12.0-nm box,to avoid the complications described above.The center of mass of the protofibril will be placed at (3.280,2.181,2.4775)in a box of dimensions 6.560x 4.362x 12.Use editconf to place the protofibril at this location:
editconf - - -center 3.2802.1812.4775-box 6.5604.36212
You can visualize the location of the protofibril within the box using,for example,VMD.Load the structure in VMD and open the Tk console.Type the following (note that >indicates the Tk prompt,not something you actually type):
>package require pbctools
>pbc box
You should see something like the following in the VMD window:
Back:The
Topology Next:Solvation and Ions
G ROMACS Tutorial
Step Three:Adding Solvent and Ions
This step is conducted much like any other simulation.Refer to the Lysozyme tutorial for a more detailed description of what is going on here if you are unsure. First,we will add water with genbox:
---p
Next,we will add ions using genion,utilizing this.mdp file.We are going to be conducting these simulations in the presence of100mM NaCl,on top of neutralizing counterions:
grompp-f ions.-p-o ions.tpr
genion-s ions.tpr-o -p-pname NA+-nname CL--neutral -conc0.1
Back:Defining the Unit Cell Next:Energy Minimization and Equilibration
G ROMACS Tutorial
Step Four:Energy Minimization and Equilibration
The energy minimization and equilibration steps are going to be conducted just like any other protein-in-water system.Here,we will perform steepest descents minimization followed by NPT equilibration.The.mdp file for minimization can be found here,and the one for NPT equilibration can be found here.
Invoke grompp and mdrun,as usual:
grompp-f minim.mdp-c -p-o em.tpr
mdrun-v-deffnm em
grompp-f npt.-p-o npt.tpr
mdrun-deffnm npt
Because these procedures are time-consuming,they are likely best run on a cluster using MPI-enabled :
mpirun-np X mdrun_mpi-deffnm npt
In the above command,"X"represents the desired number of processors over which the parallel calculation is conducted.
Back:Adding Solvent and Ions Next:Generating Configurations
G ROMACS Tutorial
Step Five:Generating Configurations
To conduct umbrella sampling,one must generate a series of configurations along a reaction coordinate,ζ.Some of these configurations will serve as the starting configurations for the umbrella sampling windows,which are run in independent simulations.The figure below illustrates these principles.The top image illustrates the pulling simulation we will run now,conducted in order to generat
e a series of configurations along the reaction coordinate.These configurations are extracted after the simulation is complete(dashed arrows in between the top and middle images).The middle image corresponds to the independent simulations conducted within each sampling window,with the center of mass of the free peptide restrained in that window by an umbrella biasing potential.The bottom images shows the ideal result as a histogram of configurations,with neighboring windows overlapping such that a continuous energy function can later be derived from these simulations.
For this example,the reaction coordinate is the z-axis.To generate these configurations,we must pull peptide A away from the protofibril.We will pull over the course of500ps of MD,saving snapshots every1ps.This setup has been established based on trial-and-error to obtain a reasonable distribution of configurations.In other systems,it may be necessary to save configurations more often,or sufficient to save configurations less often.The idea is to capture enough configurations along the reaction coordinate to obtain regular spacing of the umbrella sampling windows,in terms of center-of-mass distance between peptides A and B,the latter of which is our reference group.
The.mdp file for this pulling can be found here.A brief explanation of the pulling options used is as follows:
;Pull code
pull=umbrella
reaction masspull_geometry=distance

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。