AMBER 11.0 Introductory Tutorial

John E. Kerrigan, Ph.D. The Cancer Institute of New Jersey Robert Wood Johnson Medical School 195 Little Albany Street New Brunswick, NJ 08903 USA (732) 235-4473 phone (732) 235-6267 fax [email protected]

1

Synopsis. Today you will learn how to use one of the most versatile molecular mechanics/dynamics modeling packages, AMBER 11. [1] The AMBER suite of programs was developed by the late Peter A Kollman and colleagues at the University of California San Francisco (UCSF). The Amber package was designed with the ability to address a wide variety of biomolecules including proteins and nucleic acids as well as small molecule drugs. See http://ambermd.org for more information. The following programs are central to the AMBER package: Leap – For preparation of the parameter and coordinate files. xLeap – X-window version with a GUI for structure editing. tLeap – Terminal version (no GUI) Antechamber – special program for preparation of parameter sets (“prep” files) for small molecules (drugs/ligands, lipids). Sander – The number cruncher (aka the “brains” of the Amber MD package) that performs the computations. Ptraj – Trajectory analysis tool. Research Problem. We will study the x-ray crystal structure of the oxytocin-neurophysin complex (1NPO.PDB). [2] We will focus on the oxytocin structure itself. Remember, the x-ray crystal structure has no hydrogen atoms or lone pairs. The AMBER Leap program will take care of adding hydrogens, etc. When working with a crystal structure for the first time, you must carefully review the data in the beginning of the PDB file (REMARKS provide important information!). For example, the structure may be missing side chains in some of its residues or it may even be missing entire residues! A very important information item in a PDB file of a protein is the SSBOND records, which designate all of the disulfide bonds found in the structure. You will need this information for processing your file in Leap. Note the key disulfide bond in the oxytocin structure to the left. We will use dynamics to study the impact of mutation on this structure. What do you think the impact of removing this disulfide bond will be to the structure? (Image created using VMD [3] and Raster3D [4])

Make a project directory and use it for this exercise. mkdir myproj

2

cd myproj Download 1NPO.pdb from the protein data bank (http://www.rcsb.org) if you are extracting the oxytocin coordinates from the original pdb file or just download oxyt.pdb from the course webpage. Review the PDB file in a text editor (nedit, gedit, kedit, etc.). nedit 1NPO.pdb You will notice that the structure is a dimer of two complexes. You should also notice that in the REMARK section we find that REMARK 6 informs us that the first several residues and the last several residues (86-95) in neurophysin were not found or resolved in the structure. We will ignore these for now as they are not important. However, if you know that you will need missing residues, you may rebuild them using the SEQRES information in your PDB file and the jackal suite of programs (another tutorial, another day!). Using your text editor, remove chains A, C, and D along with any water molecules and ions from the PDB file. We only need one complex! Save your new file as oxyt.pdb. We will use Swiss PDB viewer (DeepView) to analyze the structure for missing side chains. spdbv oxyt_np.pdb Swiss PDB viewer picks up one residue with a missing side chain (59). It is always a good idea to review your PDB file using DeepView! DeepView automatically rebuilds the missing side chain for you! Save your new file. File > Layer… > File Save: oxyt_np.pdb (Click on Save and write over the old file; CAUTION: make certain that you are saving the file to your project directory or DeepView will save the file to your home directory by default!) File > Quit xLeap and tLeap – These programs perform the same function with the difference being that xLeap opens in an X-window interface and tLeap operates from a terminal prompt (no GUI). The principal function of these programs is to prepare the AMBER coordinate (inpcrd or crd for short) and topology (prmtop or just top) files. xleap –s –f leaprc.ff03 Make adjustments for use of the SPC/E water model. HOH=SPC WAT=SPC loadamberparams frcmod.spce The leaprc.ff03 file loads the parameters for the AMBER03 force field. [5] In the xleap Xwindow prompt or tleap prompt (both are “>”) type the following commands one at a time:

3

oxy = loadPdb oxyt.pdb The bond command creates bonds between atoms. The syntax is: unit_name.residue_number.atom_name. Create the disulfide bond using … bond oxy.1.SG oxy.6.SG check oxy

An illustration of the xLeap Window After using the check command you will notice several warnings about close contacts and possibly a warning about the overall charge of the unit if the unit is not neutral in overall formal charge. The overall charge is important, as we must neutralize this charge for particle-mesh ewald electrostatics later on. To look at the unit in a graphical display use the “edit” command. edit oxy

4

The following mouse buttons perform as indicated: LEFT MIDDLE RIGHT MIDDLE + RIGHT -

Selection Rotate Translate Zoom In and Out

The structure should be fine at this point. You should notice that Leap has added hydrogen atoms to the peptide structure. To exit the unit editor, use Unit > Close from the menu. Save a topology and coordinate file for in vacuo runs.

5

saveAmberParm oxy oxy_vac.prmtop oxy_vac.crd We will use the solvateOct command to solvate the structure in SPC/E water [6] using a truncated octahedron periodic box. solvateOct oxy SPCBOX 10.0 You have told Leap to solvate the unit in a truncated octahedral box using a spacing distance of 10.0 angstroms around the molecule. Ideally, you should set the spacing at no less than 8.5 Å (~ 3 water layers) to avoid periodicity artifacts. [7] For particle-mesh ewald electrostatics, [8, 9] our box side length must be > 2 X NB cutoff. We will use a 10.0 Å cutoff in our solvated system; therefore, our box side must be > 20 Å. Our box side length will be (2 X 9) + peptide dimension, which should easily be greater than 20 Å. The system must be neutral in terms of overall formal charge. Fortunately, our system is neutral as is. If this had not been the case then we would have used the addIons command to neutralize the charge (use Na+ to counter a negative charge or Cl- to counter a positive charge). The saveAmberParm command saves the parameter file (prmtop or top file) and the initial coordinate file (inpcrd or crd file). We did this before solvating the system so we could perform an in vacuo dynamics simulation for comparision to the solvated system. saveAmberParm oxy oxy.prmtop oxy.crd Quit Terminal Leap can be used to automate this process. The script below can be used to perform the same file preparation you did in xLeap. The oxy.scr script for tleap. source leaprc.ff03 HOH=SPC WAT=SPC loadamberparams frcmod.spce oxy = loadPdb oxyt.pdb bond oxy.1.SG oxy.6.SG check oxy saveAmberParm oxy oxy_vac.prmtop oxy_vac.crd solvateOct oxy SPCBOX 10.0 saveAmberParm oxy oxy.prmtop oxy.crd Quit The tleap program is very useful in situations where you might need to process a very large structure that requires a large solvent box and ions. Use the following command format for running tLeap. tleap –s –f oxy.scr

6

You‟re done with file preparation (hooray!). You will notice a file called “leap.log” in your working directory. Leap records every command and result (as well as commands issued by Leap under the hood) in a log file. Some of the information in this log file is required for publication (e.g. periodic box size; # of water molecules; # of ions if any, etc). Make a separate directory for your in vacuo dynamics runs and move your oxy_vac.top and oxy_vac.crd files there. mkdir in_vacuo mv oxy_vac.* in_vacuo Molecular Dynamics The in vacuo Model. This computation will be performed in two steps using an NVT ensemble. An energy minimization will be done followed by a production dynamics run. Step 1. Energy Minimization of the System. We perform a steepest descents energy minimization to relieve bad steric interactions that would otherwise cause problems with or dynamics runs followed up with conjugate gradient minimization to get closer to an energy minimum. Using the Sander program. The SANDER program is the number crunching juggernaut of the AMBER software package. SANDER will perform energy minimization, dynamics and NMR refinement calculations. You must specify an input file to tell SANDER what computations you want to perform and how you would like to perform those computations. Study the input file for energy minimization below. min_vac.in oxytocin: initial minimization prior to MD &cntrl imin = 1, maxcyc = 500, ncyc = 250, ntb = 0, igb = 0, cut = 12 /

General format of the sander command and its flags (The following is just an example, do not try to run it as it will crash!): sander –O –i in –o out –p prmtop –c inpcrd –r restrt [-ref refc –x mdcrd –v mdvel –e mden –inf mdinfo] Issue the following command to run your in vacuo energy minimization.

7

nohup sander –O –i min_vac.in –o min_vac.out –p oxy_vac.prmtop – c oxy_vac.crd –r oxy_vacmin.rst7 & We use “nohup” to prevent hang-ups. The ampersand at the end of the command runs the job in background. Monitor the progress using “tail –20 min_vac.out” command. The in vacuo minimization should proceed very rapidly (It only took 1 second on our 3.1 GHz linux workstation!). We run our in vacuo dynamics in the next step for 1000 ps (1 ns). Step 2. Molecular dynamics. md_vac.in oxytocin MD in-vacuo, 12 angstrom cut off &cntrl imin = 0, ntb = 0, igb = 0, ntpr = 100, ntwx = 500, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0, nstlim = 500000, dt = 0.002, cut = 12.0, ioutfm = 1 /

nohup sander –O –i md_vac.in –o md_vac.out –p oxy_vac.prmtop –c oxy_vacmin.rst7 –r oxy_vacmd.rst7 –x oxy_vacmd.nc –ref oxy_vacmin.rst7 –inf mdvac.info & This computation took about 14.0 minutes on a 3.1 GHz linux workstation. The stuff in the input files in a nutshell. There are many more settings/flags above than we can possibly explain here. So, for more in-depth info please consult the AMBER User‟s Manual. The settings important for minimization are highlighted and explained below. &cntrl and / – Most if not all of your instructions must appear in the “control” block (hence &cntrl). cut = nonbonded cutoff in angstroms. ntr = Flag used to perform position restraints (1 = on, 0 = off) imin = Flag to run energy minimization (if = 1 then perform minimization; if = 0 then perform molecular dynamics). macyc = maximum # of cycles ncyc = After ncyc cycles the minimization method will switch from steepest descents to conjugate gradient. ntmin = Flag for minimization method. (if = 0 then perform full conjugate gradient min with the first 10 cycles being steepest descent and every nonbonded pairlist update; if = 1 for ncyc cycles steepest descent is used then conjugate gradient is switched on [default]; if = 2 then only use steepest descent) dx0 = The initial step length dxm = The maximum step allowed drms = gradient convergence criterion 1.0E-4 kcal/mol•Å is the default

8

Molecular Dynamics in a Water Box This job will be accomplished in 4 steps: Step 1. Restrained Minimization – relieve bad vander Waals contacts in the surrounding solvent while keeping the solute (protein) restrained. Step 2. Unrestrained Minimization – Relieve bad contacts in the entire system. Step 3. Restrained Dynamics – Relax the solvent layers around the solute while gradually bringing the system temperature from 0 K to 300 K. Step 4. Production Run – Run the production dynamics at 300 K and 1 bar pressure. ****************** Step 1. - Restrained Minimization oxytocin: initial minimisation solvent + ions &cntrl imin = 1, maxcyc = 1000, ncyc = 500, ntb = 1, ntr = 1, cut = 10 / Hold the protein fixed 250.0 RES 1 9 END END

Hold the peptide fixed 250.0 (This is the force in kcal/mol used to restrain the atom positions.) RES 1 9 (Tells AMBER to apply this force to residue #‟s 1 to 9). nohup sander -O -i min1.in -o min1.out -p oxy.prmtop -c oxy.crd -r oxy_min1.rst7 -ref oxy.crd & Step 2. Unrestrained Minimization – oxytocin: initial minimisation whole system &cntrl imin = 1, maxcyc = 2000, ncyc = 750, ntb = 1, ntr = 0, cut = 10 /

9

nohup sander -O -i min2.in -o min2.out -p oxy.prmtop -c oxy_min1.rst7 -r oxy_min2.rst7 -ref oxy_min1.rst7 & Step 3. Position-Restrained Dynamics. This initial dynamics run is performed to relax the positions of the solvent molecules. In this dynamics run, we will keep the macromolecule atom positions restrained (not fixed, however). In a position-restrained run, we apply a force to the specified atoms to minimize their movement during the dynamics. The solvent we are using in our system, water, has a relaxation time of 10 ps; therefore we need to perform at least > 10 ps of position restrained dynamics to relax the water in our periodic box. Contents of md1.in oxytocin: 25ps MD with res on protein &cntrl imin = 0, irest = 0, ntx = 1, ntb = 1, cut = 10, ntr = 1, ntc = 2, ntf = 2, tempi = 0.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 12500, dt = 0.002, ntpr = 100, ntwx = 500, ntwr = 1000, ioutfm = 1 / Keep protein fixed with weak restraints 10.0 RES 1 9 END END

ntb = 1 Constant volume dynamics. imin = 0 Switch to indicate that we are running a dynamics. nstlim = # of steps limit. dt = 0.002 time step in ps (2 fs) temp0 = 300 reference temp (in degrees K) at which system is to be kept. tempi = 100 initial temperature (in degrees K) gamma_ln = 1 collision frequency in ps-1 when ntt = 3 (see Amber 8 manual). ntt = 3 temperature scaling switch (3 = use langevin dynamics) tautp = 0.1 Time constant for the heat bath (default = 1.0) smaller constant gives tighter coupling. vlimit = 20.0 used to avoid occasional instability in dynamics runs (velocity limit; 20.0 is the default). If any velocity component is > vlimit, then the component will be reduced to vlimit. comp = 44.6 unit of compressibility for the solvent (H2O) ntc = 2 Flag for the Shake algorithm (1 – No Shake is performed; 2 – bonds to hydrogen are constrained; 3 – all bonds are constrained). tol = #.##### relative geometric tolerance for coordinate resetting in shake.

10

ioutfm = 1 Use netcdf file format for trajectory file output. You will note that we used a smaller restraint force (10 kcal/mol). For dynamics, one only needs to use 5 to 10 kcal/mol restraint force when ntr = 1 (uses a harmonic potential to restrain coordinates to a reference frame; hence, the need to include reference coordinates with the –ref flag.). Larger restraint forces lead to instability in the shake algorithm with a 2 fs time step. Larger restraint force constants lead to higher frequency vibrations, which in turn lead to the instability. Excess motion away from the reference coordinates is not possible due to the steepness of the harmonic potential. Therefore, large restraint force constants are not necessary. nohup sander –O –i md1.in –o md1.out –p oxy.prmtop –c oxy_min2.rst7 –r oxy_md1.rst7 –x oxy_md1.nc –ref oxy_min2.rst7 – inf md1.info & Monitor the progress using “tail –40 md1.out”. Our run took approximately 29 minutes with no other processes running in background on a 3.1 GHz linux workstation. Step 4. The Production Run. This is where we do the actual molecular dynamics run. You will do a 1000 ps (1 ns) run. Contents of md2.in oxytocin: 1000ps MD &cntrl imin = 0, irest = 1, ntx = 5, ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0, cut = 10, ntr = 0, ntc = 2, ntf = 2, tempi = 300.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 500000, nscm = 500, dt = 0.002, ntpr = 100, ntwx = 500, ntwr = 1000, ioutfm = 1 /

ntb = 2 Constant pressure dynamics. ntp = 1 md with isotropic position scaling. taup = 2.0 pressure relaxation time in ps pres0 = 1 reference pressure in bar. nscm = remove center of mass motion (translational only) every 500 steps (1 ps) Issue the following command to run your dynamics simulation. nohup sander –O –i md2.in –o md2.out –p oxy.prmtop –c oxy_md1.rst7 –r oxy_md2.rst7 –x met_md.nc –ref oxy_md1.rst7 –inf md2.info &

11

Monitor the progress using “tail –40 md2.out”. Our run took approximately 2.8 hours to run on a SunFire 6900 server with the pmemd program (high performance parallel version of sander) using 16 cores (1.1 GHz UltraSparc IV). Analysis - Copy process_mdout.perl to your working project directory. Use the “which” command to determine the location of your perl interpreter (type “which perl”). On the SGI, the location (after #!) should be /usr/sbin/perl and on linux this location should be /usr/bin/perl. Check to make sure that the path in the perl script (1st line in the script) matches whatever path was given by the result of your which command. Make sure that the program is executable (use chmod 755 filename) where filename = process_mdout.perl. ./process_mdout.perl md1.out md2.out The resulting files are readable by the Grace program (http://plasma-gate.weizmann.ac.il/Grace/ ) or as space delimited files in Excel. The Grace program is available on UMDNJ SGI and linux machines in the computer labs. The Grace program is available for free download from the URL. We will take a look at summary.EPTOT (potential energy plot). xmgrace summary.EPTOT The potential energy fluctuates throughout the simulation. There is a jump in energy early on during the water equilibration (restrained MD in the first 25 ps) phase followed by a general trend toward lower energy, which is a good sign that the dynamics is leading toward lower energy conformations. Plot other files. Use Microsoft Excel to read in the file as a space delimited file or use Grace on the unix machines. summary.TEMP gives the temperature fluctuation with time. summary.PRES gives the pressure fluctuation with time. Trajectory Analysis. We will use the ptraj program. We will compute the protein backbone RMSD, a structure average, and the B-factor (temperature factor) for the backbone atoms. First we re-image the periodic box for proper visualization of the trajectory using VMD. Contents of re_image.in trajin oxy_md2.nc center :1-9 image center familiar trajout oxy_md2image.nc netcdf

12

Non imaged

Imaged

Contents of rms_avg.in trajin oxy_md2image.nc center :1-9 mass origin rms first mass out rms_oxy_md2.dat :1-9@N,C,CA average oxy_avg.rst rest reference oxy_avg.rst.-1 rms reference mass out rms_oxy_fitavg.dat :1-9@N,C,CA trajout oxy_md2fitavg.nc netcdf

3.5

Protein Backbone RMSD (Å)

3 2.5 2 in_vacuo water_box

1.5 1 0.5 0 0

200

400

600

800

1000

Time(ps)

Contents of bfactor.in trajin oxy_md2fitavg.nc atomicfluct out oxy_bkbone_Bfac.dat :1-9@N,C,CA byres bfactor

13

25

Temperature Factor

20

15

10

5

0 1

2

3

4

5

6

7

8

9

Residue No.

trajin – specifies trajectory file to process (input) center – centers the selection geometrically (:1-9 means residues 1 through 9; mass keyword specifies to use center of mass rather than geometric center) image – under periodic boundary conditions, forces molecules into the primary cell. rms – computed RMS fit to the first structure of the first trajectory read in. out – specifies name of output file @N, C, CA – Atom mask specifier (backbone atoms) (Note: The @ symbol is the atom specifier; alternately or in combination, you may use the colon : to specify residue ID # as well) For example if you only desired to examine the backbone atoms of residue #23, use :23@N,C,CA time – tells ptraj that each frame represents 1 ps Issue the following command to run ptraj. ptraj oxy.top re_image.in > re_image.log xmgrace rms_oxy_md2.dat The first steps involve centering the box followed by mass weighted RMS fit analysis, which is required for computation of the final average structure. The water molecule positions are highly variable in the average structure and will be essentially worthless. Remove the water molecules. grep –v ‘WAT’ oxy_avg.pdb > oxy_avg_nw.pdb Use rasmol to view the structure. rasmol oxy_avg_nw.pdb You will notice that the structure is crude. The hydrogen atom positions are distorted. You will need to refine the structure average using energy minimization. A NOTE ABOUT AVERAGE STRUCTURES: 14

Average structures are okay when you need a structure that is representative of a specific time range in the trajectory (e.g. a region of stability or equilibrium). Average structures are not physically meaningful taken “as is” and typically need refinement (a little energy minimization). If you are interested in conformation changes, review the relevant range of snapshots from the trajectory (do not use a structure average!). A few of the peptide residues will be labeled differently than IUPAC format (e.g. Histidine and Cysteine). These names must be corrected for use with other programs. When converting an amber restart file (rst file) to pdb with ambpdb for use with other modeling programs us the – bres flag to change the residue names to the standard Brookhaven format. For now we will convert the amber coordinate file to a PDB file that can be used with Leap as follows … ambpdb –aatm –p oxy.prmtop < oxy_avg.rst.-1 > oxy_avg.pdb grep –v ‘WAT’ oxy_avg.pdb > oxy_avg_nowat.pdb or use your favorite text editor to remove the „WAT‟ residues. This should give a PDB file with the proper set of hydrogen atom names for use with Leap. Create a new in vacuo crd and top files with Leap and use the in vacuo minimization input to setup the refinement of your average structure. Analysis of hydrogen bonds over the course of the trajectory. Use the following input file (hbond.in) to ptraj.

Illustration of hydrogen bonds to be analyzed. (Image created using VMD [3] and Raster3D [4]) trajin oxy_md2.nc

15

donor TYR O acceptor ASN N H acceptor CYX N H hbond distance 3.5 angle 120.0 includeself donor acceptor neighbor 2 series \ hbond

DONOR/ACCEPTOR – use to specify the donor acceptor atoms. DISTANCE – use to specify the cutoff distance in angstroms between the heavy atoms participating in the interaction. (3.0 is the default) INCLUDESELF – include intramolecular H-bond interactions if any. ANGLE – The H-bond angle cutoff (H-donor-acceptor) in degrees. (0 is the default) SERIES – Directs output of H-bond data summary to STDOUT. ptraj oxy.prmtop hbond.in > hbond.dat The statistical analysis in the hbond.dat file will be most important. Look for those H-bonds with a high % occupancy (> 60%). These are the more stable H-bonds. Throw out results with angles less than 120 degrees. The higher the % occupancy; the better. When analyzing H-bond data, it is best to establish reasonable guidelines for the distance and angle cutoffs. A recent paper by Chapman et al. provides a nice discussion of hydrogen bond criteria. [10] Torsion angle measurements of the peptide backbone. (torsions.in). trajin oxy_md2.nc dihedral phi_1 :1@C :2@N :2@CA :2@C out phi_1 dihedral psi_1 :1@N :1@CA :1@C :2@N out psi_1 dihedral omega_1 :1@CA :1@C :2@N :2@CA out omega_1 dihedral phi_2 :2@C :3@N :3@CA :3@C out phi_2 dihedral psi_2 :2@N :2@CA :2@C :3@N out psi_2 dihedral omega_2 :2@CA :2@C :3@N :3@CA out omega_2 dihedral phi_3 :3@C :4@N :4@CA :4@C out phi_3 dihedral psi_3 :3@N :3@CA :3@C :4@N out psi_3 dihedral omega_3 :3@CA :3@C :4@N :4@CA out omega_3 dihedral phi_4 :4@C :5@N :5@CA :5@C out phi_4 dihedral psi_4 :4@N :4@CA :4@C :5@N out psi_4 dihedral omega_4 :4@CA :4@C :5@N :5@CA out omega_4 dihedral phi_5 :5@C :6@N :6@CA :6@C out phi_5 dihedral psi_5 :5@N :5@CA :5@C :6@N out psi_5 dihedral omega_5 :5@CA :5@C :6@N :6@CA out omega_5

ptraj oxy.prmtop torsions.in > torsions.log The individual files can be plotted with xmgrace or Microsoft Excel to view the dihedral angle fluctuation with time in the simulation. How to view a trajectory movie with VMD. 16

A viewer utility that you may use on various platforms to view Amber coordinate files and trajectory files is Visual Molecular Dynamics (VMD; see http://www.ks.uiuc.edu/Research/vmd/ ). [3] In your unix shell, type “vmd” to start the VMD program (IMPORTANT: it is best to start the program from within your working project directory!). VMD Main: File > New Molecule > Molecule File Browser > Filename: Click on Browse… > Choose a molecule file: oxy.top > Click OK Molecule File Browser > Determine file type: parm7 > Click on Load Keep the Molecule File Browser window open! Molecule File Browser > Filename: Click on Browse… > Choose a molecule file: oxy_md2.mdcrd > Click OK Molecule File Browser > Determine file type: crdbox > Click on Load You will notice lots of water molecules. Use Graphics Representations to view just the protein. VMD Main: Graphics > Representations > Graphical Representations > Selected Atoms: protein Click Apply Use the left mouse button for x-y rotation and the middle mouse button for z rotation. Hit the s key on your keyboard and the mouse goes into scale (zoom in and out) mode. Hit the t button and the mouse will be in translation mode. Hit the r button to return to rotation mode. The AMBPDB Conversion Program How to use the ambpdb program to convert an amber restart coordinate file to a PDB file. For example: ambpdb –aatm –bres –p oxy.prmtop < oxy_md2.rst7 > oxy_md2.pdb The –aatm flag insures that the hydrogen atom names conform to PDB standard and the –bres flag insures that the residue names conform to PDB standard. Do not use the –bres flag if you will be using this pdb file as input to xLeap for another computation! We will use this pdb file for LSQ fit to the in vacuo structure. We will use a utility from the Gromacs software package [11] known as g_confrms to perform the LSQ fit. First, remove the water molecules from the PDB file you created. grep –v ‘WAT’ oxy_md2.pdb > oxy_md2_nw.pdb Next, convert your in vacuo structure to a pdb file. ambpdb –aatm –bres –p oxy_vac.prmtop < oxy_vacmd.rst7 > oxy_vacmd.pdb 17

g_confrms –f1 oxy_md2_nowat.pdb –f2 oxy_vacmd.pdb –o oxy_vac_wet_fit.pdb Select Group 1 (Protein) for both structures. The program will process and report a RMSD for the fit in nm (gromacs uses nm not Angstroms!). Our RMSD = 0.34226 nm = 3.4226 Å Use “rasmol –nmrpdb oxy_vac_wet_fit.pdb” to view the fit structure. Use Colours > Model and Display > Cartoon to view the two structure backbones coloured independently.

Legend: Blue = water box dynamics structure; Red = in vacuo dynamics result Using Procheck to check the stereochemical integrity of your structure. Procheck is a protein structure analysis program [12] available from (http://www.biochem.ucl.ac.uk/~roman/procheck/procheck.html ) that you should use to evaluate the results of your dynamics study. Take your structure, read it into SPDBV[13] or Sybyl (Tripos, Inc.) and remove the hydrogen atoms to prepare it for analysis. Make a separate directory for the analysis and move the PDB file into this directory. If you are using the DOS/Windows version of Procheck, copy the procheck.prm and procheck.dat files to this directory. Issue the command below to initiate the analysis. The “2.0” is the resolution in angstroms (if you were analyzing a crystal structure). Two angstroms is a nominal value used for the analysis of structures derived from NMR or modeling (or, for the case of models use the resolution of the starting crystal structure). Here is how to run procheck on two different systems. Linux/Unix Procheck: procheck filename.pdb 2.0 DOS/Windows Procheck: pro filename 2.0 Do not use the .pdb suffix for DOS procheck!

18

Procheck writes a number of files as output. A number of postscript files (plots) are generated [use ghostscript (http://www.cs.wisc.edu/~ghost/) or imagemagick (http://freealter.org/doc_distrib/ImageMagick-5.1.1/ ) to convert to PDF format or some other graphics format.]. We have ImageMagick on our Linux machines. Use imagemagick‟s “convert” program (e.g. “convert file.ps file.png” converts file.ps to a PNG file). The following is a summary of our results (Note: the sample data presented is from a simulation run with the TIP3P water model[14] and will be different from your data.). +---------->---------+ | | | oxy_md2_nw.pdb 2.0 9 residues | | | +| Ramachandran plot: 83.3% core 16.7% allow 0.0% gener 0.0% disall | | | *| All Ramachandrans: 1 labelled residues (out of 7) | | Chi1-chi2 plots: 0 labelled residues (out of 0) | +| Main-chain params: 4 better 0 inside 2 worse | | Side-chain params: 5 better 0 inside 0 worse | | | *| Residue properties: Max.deviation: 3.6 Bad contacts: 0 | *| Bond len/angle: 6.0 Morris et al class: 1 1 2

Based on an analysis of 118 structures of resolution of at least 2.0 Angstroms and R-factor no greater than 20.0 a good quality model would be expected to have over 90% in the most favoured regions [A,B,L]. *| G-factors Dihedrals: -1.11 | | | | M/c bond lengths: 85.7% within limits | | M/c bond angles: 61.9% within limits |

Covalent:

-2.14

Overall:

-1.71

14.3% highlighted 38.1% highlighted

G-factors provide a measure of how unusual, or out-of-the-ordinary, a property is. Values below -0.5* - unusual Values below -1.0** - highly unusual

19

Important note: The main-chain bond-lengths and bond angles are compared with the Engh & Huber (1991) ideal values derived from small-molecule data. Therefore, structures refined using different restraints may show apparently large deviations from normality. The main chain bond lengths and angles are highly unusual in our case because this is an unrefined structure (i.e. not energy minimized). Perform an energy minimization on your structure and compare! file_01.ps is the Ramachandran plot Use the convert utility (ImageMagick) to convert the postscript file (better yet, use Ghostview/Ghostscript) to an acceptable image file. To convert to jpeg use the following command: convert file_01.ps file_rama.jpeg Load the image into imageview. imgview file_rama.jpeg or use ImageMagick‟s display program display file_rama.jpeg

20

Note that at least four of our residues fall in the helix region (A) of the Ramachandran map. [15] The triangles are the glycine residues, which are special (they may appear in disallowed regions of the plot!). Some points may overlap and not be visible.

21

H H H H

N





H



H

O



O

H

H

N O

H

H H

H

Dihedral angles in a polypeptide chain. Adapted from Schulz and Schirmer.[15] Questions to consider: 1. Do you notice any clustering of hydrophobic side chains in this model?

2. In your dihedral angle plots (, , and ), which was the most variable? What region of the backbone was most flexible (analyze the  and  angles as well as the bfactor results)? What region was more rigid? What interactions contributed to this observed rigidity?

3. Comment on the differences observed between the in vacuo vs the water box results. Why does the in vacuo structure appear to be more compact?

4. Mutate one of the cysteine residues to an alanine. You will need to change the residue name of the remaining cysteine residue from CYX to CYS (Amber uses CYX to designate cysteine residues that are disulfide bonded and the normal CYS residue name to designate cysteine residues not involved in disulfide bonding.) in your PDB file (start with oxyt.pdb). One approach is to use the mutate function in DeepView (spdbv). Before loading the structure into

22

DeepView, change all of the CYX residue names to CYS (otherwise DeepView will be confused by CYX and change these to HETATM records!) using the Search/Replace function of any text editor. Mutate one of the CYS residues to ALA using the mutate button in DeepView, then use File > Save > Layer to save the new pdb file. After you save the new PDB file from DeepView, edit out the SPDBV lines in the file using any text editor. Perform a 250 ps water box simulation just as you did in the tutorial. Compare your mutated structure final snapshot with the structure snapshot (oxy_md2.rst7) of the wild type structure. What differences do you notice in comparing the backbones? Are there new stabilizing interactions (H-Bonds or Hydrophobic?)? What role does the solvent (water) play in the outcome?

23

Bibliography: 1.

2. 3. 4. 5.

6. 7.

8. 9. 10.

11. 12.

13. 14. 15.

Case, D.A., T.E. Cheatham III, T. Darden, H. Gohlke, H. Luo, K.M. Merz, A. Onufriev, C. Simmerling, B. Wang, and R. Woods, The Amber biomolecular simulation programs. J. Comput. Chem., 2005. 26: p. 1668-1688. Rose, J., C. Wu, C. Hsiao, E. Breslow, and B. Wang, Crystal structure of the neurophysin-oxytocin complex. Nature Struct. Biol., 1996. 3(2): p. 163-169. Humphrey, W., A. Dalke, and K. Schulten, VMD - Visual Molecular Dynamics. J. Molec. Graphics, 1996. 14.1: p. 33-38. Merritt, E.A. and D.J. Bacon, Raster3D: Photorealistic Molecular Graphics. Methods Enzymol., 1997. 277: p. 505-524. Duan, Y., C. Wu, S. Chowdhury, M. Lee, G. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. Wang, and P. Kollman, A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J. Comput. Chem., 2003. 24(16): p. 1999-2012. Berendsen, H.J., J. Grigera, and T. Straatsma, The missing term in effective pair potentials. J. Phys. Chem., 1987. 91: p. 6269-6271. Weber, W., P.H. Hünenberger, and J.A. McCammon, Molecular Dynamics Simulations of a Polyalanine Octapeptide under Ewald Boundary Conditions: Influence of Artificial Periodicity on Peptide Conformation. J. Phys. Chem. B, 2000. 104(15): p. 3668-3575. Darden, T., D. York, and L. Pedersen, Particle Mesh Ewald: An N-log(N) method for Ewald sums in large systems. J. Chem. Phys., 1993. 98: p. 10089-10092. Essmann, U., L. Perera, M.L. Berkowitz, T. Darden, H. Lee, and L. Pedersen, A smooth particle mesh ewald potential. J. Chem. Phys., 1995. 103: p. 8577-8592. Fabiola, F., R. Bertram, A. Korostelev, and M.S. Chapman, An Improved Hydrogen Bond Potential: Impact on Medium Resolution Protein Structures. Protein Sci., 2002. 11: p. 1415-1423. Lindahl, E., B. Hess, and D. van der Spoel, GROMACS 3.0: a package for molecular simulation and trajectory analysis. J. Mol. Model, 2001. 7: p. 306-317. Laskowski, R., M. MacArthur, D. Moss, and J. Thornton, PROCHECK: A program to check the stereochemical quality of protein structures. J. Appl. Cryst., 1993. 26: p. 283291. Guex, N. and M.C. Peitsch, SWISS-MODEL and the Swiss-Pdb Viewer: An Environment for Comparative Protein Modeling. Electrophoresis, 1997. 18: p. 2714-2723. Jorgensen, W., J. Chandrasekhar, J. Madura, and M. Klein, Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys., 1983. 79: p. 926-935. Schulz, G.E. and R.H. Schirmer, Principles of Protein Structure. Springer Advanced Texts in Chemistry, ed. C.R. Cantor. 1979, New York: Springer-Verlag. 314.

24