Amber Molecular Simulations Workshop:

1. Amber Molecular Simulations Workshop: Yeng-Tseng Wang Email: c00jsw00@gmail.com

2. Outline Basic: 1. Introductions: How to set-up calculations --Case study: DNA (01_dna) --Case study: small agent with protein receptor(02_HIV) --Case study: Ibalizumab (03_Ibalizumab) Advance: 2. Free energies calculations --Case study: small agent with protein receptor (02_HIV) --Case study: Ibalizumab (03_Ibalizumab) 3. Enhanced sampling --Case study: small agent with protein receptor 4. Membrane-Bound Proteins 5. Micelle 6. Force field : 4-Hydroxyl-Proline (PR4)

3. Advance 2. Free energy calculations --2.1 LIE methodsietraj --2.2 MMPBSA needed force constraints

4. Implicit solvent simulations: 2.1(SIE): HIV-1 protease

5. 2.1(SIE): HIV-1 protease Solvated interaction energies (SIE) method http://www2.bri.nrc.ca/ccb/pub/sietraj_main.php

6. Implicit solvent simulations: Sietraj 1. cd bin and source export.sh 2. cd sietraj 3. cpptraj –i cpptraj.gb.in 4. sh sie.sh Prediction Binding Affinity (△GP): -43.655 KJ/mol Experimental Binding Affinity (△GExp): -44.3 ~ -55.2 KJ/mol (-10.59~ -13.20 kcal/mol) 2.1(SIE): HIV-1 protease

7. 2.1(SIE): HIV-1 protease Explicit solvent simulations Analyzing LIE method: sietraj Experimental Binding Affinity (△GExp): -10.59~ -13.20 kcal/mol

8. Creating proteins structure: 1. GB and Tip3 water box simulations (the same as before) Reference to the files (gb.md.in gb.min.in polyAT_wat_md1.in polyAT_wat_md2.in polyAT_wat_min1.in polyAT_wat_min2.in 3O2D.sh 3O2D.TIP3.sh) 2. $BRIMM_ROOT/bin/sietraj -pt 3O2D.top -trj 3O2D.md.x -sf 1 -ef 100 -inc 10 -tr 1-6633 -lr 6634-9434 -o gb.sie.out -sie 3. $BRIMM_ROOT/bin/sietraj -ave gb.sie.out 2.1(MMPBSA): Ibalizumab

9. 2.2(MMPBSA): HIV-1 protease • The acronym MM-PBSA stands for Molecular Mechanics- Poisson Bolzmann Surface Area • The MM-PBSA approach represents the postprocessing method to evaluate free energies of binding or to calculate absolute free energies of molecules in solution. P.A. Kollman at al.,Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models, Acc. Chem. Res. 2000, 33, 889-897

10. ΔGBind A AB B Solvent Solvent 2.2(MMPBSA): HIV-1 protease ΔGBind = ΔGAB – ΔGA - ΔGB

11. 2.2(MMPBSA): HIV-1 protease ΔGBind = ΔGAB – ΔGA - ΔGB Approximate ΔGBindas ΔGBind ≈ GAB – GA – GB Where GX is the calculated average free energy ΔGX = EMM + GSolv – TSMM

12. 2.2(MMPBSA): HIV-1 protease ΔGX =EMM + GSolv – TSMM WhereEMM is the average molecular mechanical energy: EMM = Ebond + Eangle + Etors + Evdw + Eelec GSolvis the calculated solvation free energy – TSMM is the solute entropy, which can be estimated by using normal-mode analysis

13. 2.2(MMPBSA): HIV-1 protease ΔGX =EMM + GSolv – TSMM Could be easily calculated • In practice entropy contributions is usually neglected • Could be computationally expensive • Tend to have a large margin of error that introduces significant uncertainty in the result. ?

14. 2.2(MMPBSA): HIV-1 protease Compromise: MD trajectory One carries out a MD simulation in a periodic box with solvent Snapshots of representative structures Evaluate EMM , GSolv and –TSMM for every saved snapshot

15. ΔGBind A AB B Solvent Solvent 2.2(MMPBSA): HIV-1 protease • In general case, one carries out three independent MD simulations: for ligand, receptor, and complex • Single trajectory approach: one makes the approximation that no significant conformational changes occur upon binding so that the snapshots for all three species can be obtained from a single trajectory for a complex

16. 2.2(MMPBSA): HIV-1 protease • Explicit solvent simulations • Analyzing • mmpbsa • tleap –s –f tleap.mmpbsa.in • sh mmpbsa.in • Experimental Binding Affinity (△GExp): -10.59~ -13.20 kcal/mol

17. Advance 3. Enhanced sampling --3.1 aMD --3.2 GaMD --3.3 remd

18. Enhanced sampling Potential Energies Potential Energies Enhanced-sampling Molecular Dynamics Normal Molecular Dynamics

19. Method • Replica exchange

20. Method • Accelerated Molecular Dynamics V*(r) = V(r), V(r) ≥ E V*(r) = V(r) + ΔV(r), V(r) < E where V(r) is the original potential, E is the reference energy, and V*(r) is the modified potential energy. △V(r) is boost potential energies.

21. Creating proteins structure: 1. GB and Tip3 water box simulations (the same as before) Refe 3.1(aMD): Ibalizumab a) EthreshP. Average total potential energy threshold. b) alphaP. Inverse strength boost factor for the total potential energy. c) EthreshD. Average dihedral energy threshold. d) alphaD. Inverse strength boost factor for the dihedral energy. a) EthreshP: E(tot)= -219899.7170 kcal mol-1 + (0.16kcal*mol-1 atom-1 * 90654 atoms) = -205349 kcal mol-1 b) alphaP: Alpha(tot)= (0.16kcal mol-1 atom-1 * 90654 atoms) = 14504 kcal mol-1 c) EthreshD: E(dih)=6321.83 kcal mol-1 + (4kcal mol-1 residue-1 * 27073 solute residues) = 114613 kcal mol-1 d) alphaD: Alpha(dih)=(1/5)*(4kcal mol-1 residues-1 * 58 solute residues) = 21658 kcal mol-

22. Method • Gaussian Accelerated Molecular Dynamics DOI: 10.1021/acs.jctc.5b00436(Yinglong Miao) Y. Miao and J. A. McCammon (2016), Graded activation and free energy landscapes of a muscarinic G protein-coupled receptor. Proc Natl Acad Sci U S A, 113 (43): 12162–12167. Miao, Y.; Feher, V. A.; McCammon, J. A., Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation. J. Chem. Theory Comput. 2015, 11, 3584-3595.

23. 3.1(Gamd): Ibalizumab GAMD for Mu/cell mem &cntrl imin = 0, irest = 0, ntx = 1, nstlim = 50000000, dt = 0.002, ntc = 2, ntf = 2, tol = 0.000001, iwrap = 1, ntb = 1, cut = 8.0, ntt = 3, temp0 = 310.0, tempi = 310.0, ntpr = 500000, ntwx = 500000, ntwr = 500, ntxo = 1, ig = -1, igamd = 3, iE = 1, irest_gamd = 0, ntcmd = 200000, nteb = 200000, ntave = 100000, sigma0P = 6.0, sigma0D = 6.0, /

24. Membrane-Bound Proteins&Micelle&other reference: http://ambermd.org/tutorials/advanced/tutorial16/ http://www.charmm-gui.org/ http://www.charmm-gui.org/?doc=input/membrane&time=1405600478 Yeng-Tseng Wang Email: c00jsw00@kmu.edu.tw

25. http://www.charmm-gui.org/?doc=input/membrane Rhodopsin (PDB ID: 1U19)

26. Currently Supported CHARMM-GUI Lipids with Lipid14*Can be built with VMD builder and processed with charmmlipid2amber.py

27. Currently Supported amber lipid molecules

29. Building amber files charmmlipid2amber.py is a robust residue and atom renaming and reordering script available with AmberTools 14 update 1. It is used here to convert the CHARMM-GUI PDB format structures into a PDB format that can be read with LEaP and Lipid14. If you do not have AmberTools 14 update 1 available, you can download the zip file from this web site. This can be placed in your Amber directory and unzipped there. (1) charmmlipid2amber.py : rename residues > python charmmlipid2amber.py -i step5_assembly.pdb -c charmmlipid2amber.csv -o amber.pdb > vi amber.pdb (CHL1 substitute CHL) (residue name)

30. Building amber files (modifying amber.pdb) edit original PDB file: rename D-protonated HIS to HID rename E-protonated HIS to HIE rename dual-protonated HIS to HIP D-protonated means the H sits on ND; E-protonated (the normal case), it sits on NE; double-protonated means there is a H on both Nitrogen atoms in the ring. rename CYS involved in S-S bond into CYX switch OE and NE of ASN and GLN where suggested by WhatIf remove terminal OT1, OT2 atoms from peptide chains remove terminal P, O1P, O2P from DNA chains remove hetero atoms (except you have prep files for them) remove all hydrogen atoms.

31. Building amber files (2) Estimate Periodic Box Size > Open vmd on windows os or linux os >set all [atomselect top all] >measure minmax $all Or more charm-gui/step5_assembly.str • For PBC BOX  X ~= 81 / Y~= 81 / Z~=120 (3) Cut protein files and remove the hydrogens from amber.pdb (4) tleap >tleap -s –f tleap.in source leaprc.lipid14 source leaprc.lipid11 source leaprc.ff12SB loadamberparams frcmod.ionsjc_tip3p aa = loadpdb amber.pdb set aa box { 86 84 108 } setBox aa centers saveAmberParm aa amber.top amber.crd quit

32. Building amber files (5)Still error Leap added 2710 missing atoms according to residue templates: 28 Heavy 2682 H / lone pairs The file contained 82 atoms not in residue templates removing previous box.. Box dimensions: 90.311000 90.290000 118.840000 Checking Unit. WARNING: There is a bond of 7.746053 angstroms between: ------- .R<PA 349>.A<C12 44> and .R<PC 350>.A<C11 5> WARNING: There is a bond of 3.479261 angstroms between: ------- .R<PA 349>.A<C15 35> and .R<PA 349>.A<H5R 36> FATAL: Atom .R<ILE 48>.A<CD 20> does not have a type. FATAL: Atom .R<ILE 54>.A<CD 20> does not have a type. FATAL: Atom .R<HSD 65>.A<N 1> does not have a type. FATAL: Atom .R<HSD 65>.A<CA 2> does not have a type. (6) HSD change to HIE/ remove CD of ILE (7) Back to step 4

33. Simulations

34. Minimization (8) Minimization >pmemd -O -i 01_Min.in -o 01_Min.out -p amber.top -c amber.crd -r 01_Min.rst 01_Min.in Lipid minimization &cntrl imin=1, ! Minimize the initial structure maxcyc=10000, ! Maximum number of cycles for minimization ncyc=5000, ! Switch from steepest descent to conjugate gradient minimization after ncyc cycles ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntc=1, ! No SHAKE ntpr=50, ! Print to mdout every ntpr steps ntwr=2000, ! Write a restart file every ntwr steps cut=10.0, ! Nonbonded cutoff in Angstroms /

35. Heating (9) Heating-1a After the initial minimization, slowly heat the system to production temperature. Choosing a production temperature is an important choice in the lipid bilayer simulation. Lipid bilayers have experimentally-measured phase transition temperatures from highly ordered gel-like phases to liquid phases. One major problem in lipid bilayer simulations is accurately simulating the phase transition of lipids. >pmemd -O -i 02_Heat.in -o 02_heat.out -p amber.top -c 01_Min.rst -r 02_Heat.rst -ref 01_Min.rst -x 02_Heat.nc

36. Lipid heating 100K &cntrl  imin=0, ! Molecular dynamics  ntx=1, ! Positions read formatted with no initial velocities  irest=0, ! No restart  ntc=2, ! SHAKE on for bonds with hydrogen  ntf=2, ! No force evaluation for bonds with hydrogen  tol=0.0000001, ! SHAKE tolerance  nstlim=2500, ! Number of MD steps  ntt=3, ! Langevin dynamics  gamma_ln=1.0, ! Collision frequency for Langevin dynamics  ntr=1, ! Restrain atoms using a harmonic potential! (See the GROUP input below)  ig=-1, ! Random seed for Langevin dynamics  ntpr=100,  ntwr=10000,  ntwx=100, ! Write to trajectory file every ntwx steps  dt=0.002, ! Timestep (ps)  nmropt=1, ! NMR restraints will be read (See TEMP0 control below)  ntb=1,  ntp=0,  cut=10.0,  ioutfm=1, ! Write a binary (netcdf) trajectory  ntxo=2, ! Write binary restart files / &wt   type='TEMP0', ! Varies the target temperature TEMP0  istep1=0, ! Initial step  istep2=2500, ! Final step  value1=0.0, ! Initial temp0 (K)  value2=100.0 / ! final temp0 (K) &wt type='END' / ! End of varying conditionsHold lipid fixed10.0 ! Force constant (kcal/(mol Angstroms^2))RES 1 804 ! Choose residuesENDEND ! End GROUP input

37. Heating (9) Heating-1b The second phase of heating slowly increases the temperature to the desired production temperature. The positions and velocities are read from the previous restart file. This time an anisotropic Berendsen weak-coupling barostat is used to also equilibrate the pressure in addition to the use of the Langevin thermostat to equilibrate the temperature. >pmemd -O -i 03_Heat2.in -o 03_Heat2.out -p amber.top -c 02_Heat.rst -r 03_Heat.rst -ref 02_Heat.rst -x 03_Heat.nc

38. Lipid heating 303K&cntrlimin=0,ntx=5, ! Positions and velocities read formattedirest=1, ! Restart calculationntc=2,ntf=2,tol=0.0000001,nstlim=50000, ! Number of MD stepsntt=3,gamma_ln=1.0, ntr=1,ig=-1,ntpr=100,ntwr=10000,ntwx=100,dt=0.002,nmropt=1,ntb=2, ! Constant pressure periodic boundary conditionsntp=2, ! Anisotropic pressure couplingtaup=2.0, ! Pressure relaxation time (ps)cut=10.0,ioutfm=1,ntxo=2,/&wttype='TEMP0',istep1=0,istep2=50000,value1=100.0,value2=303.0 /&wt type='END' /Hold lipid fixed10.0RES 1 384ENDEND

39. HOLD (10) hold In order to equilibrate the system's periodic boundary condition dimensions, it is necessary if you are using the GPU code pmemd.cuda, to run 5ns MD with a barostat. The system's dimensions and density must equilibrate before proceeding with production MD. Because the periodic boundary condition box dimensions are changing, it is necessary to increase the "skinnb" value and to restart the MD simulation after 500ns. This should avoid most "skinnb" errors. The reason for this is that, for performance reasons, the GPU code does not recalculate the non-bond list cells during a simulation. If these cells change size, due to the box changing size, by too much then it will cause the code to halt with an error related to skinnb. Once the system is equibrated the box size fluctuations should be small and so this should not be an issue during production.

40. HOLD >pmemd -O -i 04_Hold.in -o 04_Hold_1.out -p amber.top -c 03_Heat2.rst -r 04_Hold_1.rst -x 04_Hold_1.nc >pmemd -O -i 04_Hold.in -o 04_Hold_2.out -p amber.top -c 04_Hold_1.rst -r 04_Hold_2.rst -x 04_Hold_2.nc >pmemd -O -i 04_Hold.in -o 04_Hold_3.out -p amber.top -c 04_Hold_2.rst -r 04_Hold_3.rst -x 04_Hold_3.nc >pmemd. -O -i 04_Hold.in -o 04_Hold_4.out -p amber.top -c 04_Hold_3.rst -r 04_Hold_4.rst -x 04_Hold_4.nc >pmemd -O -i 04_Hold.in -o 04_Hold_5.out -p amber.top -c 04_Hold_4.rst -r 04_Hold_5.rst -x 04_Hold_5.nc >pmemd -O -i 04_Hold.in -o 04_Hold_6.out -p amber.top -c 04_Hold_5.rst -r 04_Hold_6.rst -x 04_Hold_6.nc >pmemd -O -i 04_Hold.in -o 04_Hold_7.out -p amber.top -c 04_Hold_6.rst -r 04_Hold_7.rst -x 04_Hold_7.nc >pmemd -O -i 04_Hold.in -o 04_Hold_8.out -p amber.top -c 04_Hold_7.rst -r 04_Hold_8.rst -x 04_Hold_8.nc >pmemd -O -i 04_Hold.in -o 04_Hold_9.out -p amber.top -c 04_Hold_8.rst -r 04_Hold_9.rst -x 04_Hold_9.nc >pmemd -O -i 04_Hold.in -o 04_Hold_10.out -p amber.top -c 04_Hold_9.rst -r 04_Hold_10.rst -x 04_Hold_10.nc

41. Lipid production 303K 500ps&cntrlimin=0,ntx=5,irest=1, ntc=2,ntf=2,tol=0.0000001,nstlim=250000,ntt=3,gamma_ln=1.0,temp0=303.0,ntpr=5000,ntwr=5000,ntwx=5000,dt=0.002,ig=-1, ntb=2,ntp=2,cut=10.0,ioutfm=1,ntxo=2,//&ewaldskinnb=5, ! Increase skinnb to avoid skinnb errors/

42. Production Run (11) Production run Actual production molecular dynamics can be run with the input file shown below. Temperature is controlled here using Langevin dynamics while pressure is controlled using the anisotropic Monte Carlo barostat included in Amber 14. >pmemd.cuda -O -i 05_Prod.in -o 05_Prod.out -p amber.top -c 04_Hold_10.rst -r 05_Prod.rst -x 05_Prod.nc -inf 05_Prod.mdinfo

43. Lipid production 303K 125ns &cntrl  imin=0, ! Molecular dynamics  ntx=5, ! Positions and velocities read formatted  irest=1, ! Restart calculation  ntc=2, ! SHAKE on for bonds with hydrogen  ntf=2, ! No force evaluation for bonds with hydrogen  tol=0.0000001, ! SHAKE tolerance  nstlim=62500000, ! Number of MD steps  ntt=3, ! Langevin dynamics  gamma_ln=1.0, ! Collision frequency for Langevin dyn.  temp0=303.0, ! Simulation temperature (K)  ntpr=5000, ! Print to mdout every ntpr steps  ntwr=500000, ! Write a restart file every ntwr steps  ntwx=5000, ! Write to trajectory file every ntwc steps  dt=0.002, ! Timestep (ps)  ig=-1, ! Random seed for Langevin dynamics  ntb=2, ! Constant pressure periodic boundary conditions  ntp=2, ! Anisotropic pressure coupling  cut=10.0, ! Nonbonded cutoff (Angstroms)  ioutfm=1, ! Write binary NetCDF trajectory  ntxo=2, ! Write binary restart file  barostat=2, ! Use Monte Carlo barostat (Amber 14) /

44. Analysis (11) Analysis Area per Lipid: Area per lipid is the average area that a single phospholipid occupies in an interface and is usually reported in Angstroms squared. It can be calculated with: Area per lipid = (box X dimension) * (box Y dimension) / (number of phospholipids per layer) >sh area_per_lipid.sh box_dimension.cpptraj # Read in trajectory file trajin 04_Hold_7.nc trajin 05_Prod.nc # Write the vector named "test" (selecting all atoms) # of box coordinates out to a file named "vector.dat" vector test box out vector.dat

45. Analysis (11) Analysis Electron-Density Profile : >cpptraj -i density.cpptraj -p amber.top density.cpptraj # NOTE # Your frames will be different depending on the part of # the trajectory you choose to analyze # Load the trajectory file frames 2500-7500 inclusive trajin 03_Heat.nc 1 50 1 # Center the lipid bilayer tail residues at zero coordinates center ":PC |:PA | :OL" origin # Image all water molecules to the zero coordinates using the # center of mass of the molecules image origin center ":WAT" # Calculate the electron density using 0.25 Angstrom slices. # Write out to the file electron_density.dat density out electron_density.dat electron delta 0.25 "*"

46. Micelle

47. http://www.charmm-gui.org/?doc=input/micelle

48. Build (1) Optimize the initial structure >cd charmm-gui/namd/ >namd2 step6.1_equilibration.inp Then Save pdb file with vmd software=> afternamd.pdb

49. Build (2) chamber >chamber -top top_all36_lipid.rtf -param par_all36_lipid.prm -xpsf step5_assembly.xplor.psf -crd afternamd.pdb -p amber.top -inpcrd amber.crd -nocmap -tip3_flex -box 82.0000000 82.0000000 82.0000000 (3) antechamber > antechamber -fi mol2 -fo prepi -i sds.mol2 -o sds.prepi -c bcc -j 4 -at gaff -nc -1 >parmchk -i sds.prepi -o sds.frcmod -f prepi >vi sds.prepi (MOL -> SDS) >tleap -s -f tleap.in (3) simulations > The same as Membrane-Bound Proteins

50. Analysis (4) radius of gyration (Approximate sphere) >ptraj amber.top ptraj.in (4) radius of gyration (Approximate sphere) Chamber for Membrane-Bound Proteins ??