MORC16

MORC16 is a newly designed lipid molecule, this page will describe how to model and simulate this molecule especially in different pH environment. I started with charmm to model this molecule and later moved on to Amber for better performance for bilayer lipids and liposomes.

MOLC16_chemdraw

Molecular Modelling

Create molecule by gaussian viewer,and then optimize the molecule:(on medusa) %chk=morc16.chk %Nproc=26 # opt b3lyp/6-31g(d) pop=nbo geom=connectivity save the molecule as morc16-opt.pdb(also save  as morc16-opt.mol2 for swisspram)   
















About/Help

 Submit the morc16-opt.mol2 to SwissParam server: retrieve morc16-opt.zip after about one minute, unzip the file and get: morc16-opt.crd/itp/mol2/par/pdb/psf/rtf

convpdb.pl -segnames -out charmm22 morc16-opt.pdb |sed "s/LIG/MOR/g"|sed "s/PRO0/LIG/g">morc16.pdb

mkdir data toppar

cp morc16.pdb morc16-opt.psf ./data/

cp morc16-opt.rtf morc16-opt.par ./toppar/

edit morc16-opt.rtf and change the RESI from LIG to MOR

genPSF.pl -par xtop=./toppar/morc16-opt.rtf,xpar=./toppar/morc16-opt.par morc16.pdb>./data/morc16.psf

The charges are poorly defined in topology file, we update the charges in topology file by charge information from gaussian log file .

assign_charge.pl -nbo -log morc16.log -rtf ./toppar/morc16-opt.rtf

mdCHARMM.pl -par param=22,xtop=./toppar/morc16-opt.rtf,xpar=./toppar/morc16-opt.par,dynsteps=10000 -trajout morc16_equi.dcd -cmd tmp.cmd -final morc16_equi.pdb morc16.pdb

mdCHARMM.pl -par param=22,xtop=./toppar/morc16-opt.rtf,xpar=./toppar/morc16-opt.par,gb,dynsteps=10000 -trajout morc16_gb_equi.dcd -log morc16_gb_equi.log -cmd tmp.cmd -final morc16_gb_equi.pdb morc16.pdb

convpdb.pl -solvate -cubic -cutoff 8 morc16.pdb >morc16_box.pdb

genPSF.pl -par param=22,xtop=./toppar/morc16-opt.rtf,xpar=./toppar/morc16-opt.par morc16_box.pdb >morc16_box.psf

minCHARMM.pl -par param=22,xtop=./toppar/morc16-opt.rtf,xpar=./toppar/morc16-opt.par,dynsteps=10000 morc16_box.pdb >morc16_box_min.pdb

mpirun -np 4 mdCHARMM.pl -par param=22,xtop=./toppar/morc16-opt.rtf,xpar=./toppar/morc16-opt.par,dynsteps=10000 -trajout morc16_ex_equi.dcd -log morc16_ex_equi.log -cmd tmp.cmd -final morc16_ex_equi.pdb morc16_box_min.pdb

 

 

 

  VA GB EX
CHAMRM
AMBER

 

Above are trajectories of morc16 simulation in different environment. From the equilibration we can tell the two tails of the lipid are more twiested in amber than in charmm, and the two tails are more separated in implicit solvent than in explicit solvent model.

 

Constant pH simulation of MORC16

Protonated patch for MOR: add a hydrogen to morc16-opt.mol2,

%chk=morc16.chk %Nproc=26 # opt b3lyp/6-31g(d) pop=nbo geom=connectivity change charge from 0 1 to 1 1

morc16_h

cp ./toppar/morc16-opt.rtf ./toppar/tmp.rtf edit tmp.rtf,add ATOM HD HNR 0.0000 to the end of ATOM list, then run: assign_charge.pl -nbo -log morc16_h.log -rtf ./toppar/tmp.rtf(check sum of charges, sometime it is not integer!!!)

edit ~/softwares/c37a2/toppar/top_all22_prot.inp create PRES MOR1 1.000 copy GROUP and the rest from tmp.rtf to  top_all22_prot.inp and add "BOND N HD" to BOND section, and "IC HD N C2 C1 1.04 102.14 41.43 109.10 1.53" to IC section.

pH simulation of just one monomer

copy add_h.inp to ./MORC16/ make sure the name as in red are "mor" and "morp", "if ?selresn eq mor patch morp(or whatever patch you defined mor1,mor2...) @segid @ires" then run:

charmm in=morc16 <add_h.inp

utput: morc16_h.pdb/psf/crd

copy water.crd to ./toppar/ copy solvate.inp to ./MORC16/ charmm in=morc16 <solvate.inp(change Xleng or Yleng or Zleng if you want to) outputs:morc16_h_waterbox.crd/str/prm/pdb morc16_h_solvated.psf/pdb/str/crd copy add_ions_hphmd.inp to ./MORC16/ charmm in=morc16 conc=0.100<add_ions_hphmd.inp outputs:morc16_h_ions.pdb/str/prm/crd copy minimize.inp to ./MORC16/ charmm in=morc16 <minimize.inp outputs:morc16_h.solv.ion.min.psf/pdb/crd/pbc the system should be looks like this:   Now we use this system to find parameter A and B for this molecule. copy phmd.in radius_gbsw.str to ./toppar/ edit phmd.in MOR 5.5 0.0 0.0 0.0 ATOM [protonated charge] [unprotonated charge] ...... HD [protonated charge] [unprotonated charge]  1.0 0.0 copy hphmd_derive.inp to ./MORC16/ charmm <hphmd_derive.inp extract theta and dE/dtheta value from lamb file and fit the data by xgrace. morc16_a_b replace A=-1119.99 B=0.937338 and barrier=2.0   For replica exchange: copy hphmd_phrex.inp to ./MORC16/ mpirun -np ? charmm < hphmd_phrex.inp copy hphmd_drive.inp to ./MORC16/ the fundamental change in this input script is add DERI key word to PHMD, and set mass 1.E30 BARR 0 vmd morc16_equi.pdb package require Orient namespace import Orient::orient set sel [atomselect top "all"] set I [draw principalaxes $sel] set A [orient $sel [lindex $I 2] {0 0 -1}](make sure the head face towards -z) $sel move $A set I [draw principalaxes $sel] set A [orient $sel [lindex $I 1] {0 1 0}] $sel move $A $sel writepdb morc16.pdb(replace old one) exit vmd genPSF.pl -par param=22 -xplor morc16.pdb>morc16_xplor.psf restart vmd, source lipid.tcl lipid -l morc16 -x 20 -y 20 outputs:lipid.pdb.lipid_charmm.psf,lipid_xplor.psf copy add_h.inp to ./MORC16/ make sure the name as in red are "mor" and "morp", "if ?selresn eq mor patch morp(or whatever patch you defined mor1,mor2...) @segid @ires" then run: charmm in=lipid <add_h.inp output: lipid_h.pdb/psf/crd copy water.crd to ./toppar/ copy solvate.inp to ./MORC16/ charmm in=lipid <solvate.inp(change Xleng or Yleng or Zleng if you want to) outputs:lipid_h_waterbox.crd/str/prm/pdb lipid_h_solvated.psf/pdb/str/crd copy add_ions_hphmd.inp to ./MORC16/ charmm in=lipid conc=0.100<add_ions_hphmd.inp outputs:lipid_h_ions.pdb/str/prm/crd copy minimize.inp to ./MORC16/ charmm in=lipid <minimize.inp outputs:lipid_h.solv.ion.min.psf/pdb/crd/pbc the system should be looks like this: morc16_lipid copy phmd.in radius_gbsw.str to ./toppar/ edit phmd.in MOR 5.5 -70 1.0 0.0 ATOM [protonated charge] [unprotonated charge] ...... HD [protonated charge] [unprotonated charge]  1.0 0.0 copy hphmd.inp to ./MORC16/ charmm <hphmd.inp outputs:lipid.ph?.dcd/rst/lamb/ene For replica exchange: copy hphmd_phrex.inp to ./MORC16/ mpirun -np ? charmm < hphmd_phrex.inp copy hphmd_drive.inp to ./MORC16/ the fundamental change in this input script is add DERI key word to PHMD, and set mass 1.E30 BARR 0 and add PHTEST num1 set 0.4 PHTEST num2 set 0.6 ..... PHTEST num6 set 1.4 in phmd.in, edit B=0 in output *.lamb file 10  0.4000 dE/dtheta  0.6000 dE/dtheta ..... extract those values for each  residue and calculate average theta values theta  dE/dtheta morc16_a_b  

 

MORC16 with Amber

we start from morc16.pdb from the beginning of this page, 

antechamber -i morc16.pdb -fi pdb -o morc16.mol2 -fo mol2 -c bcc -s 2

Amber atomtypes can be found here.

check the parameters that are required.

parmchk2 -i morc16.mol2 -f mol2 -o morc16.frcmod

later we load morc16.frcmod to leap to add missing parameters.

MORC16 in vacuum 

source leaprc.ff14SB
source leaprc.gaff
source leaprc.tip3p
loadamberparams morc16.frcmod
MOR =loadmol2 morc16.mol2                                                       
saveoff MOR morc16.lib
saveamberparm MOR morc16.prmtop morc16.inpcrd
quit

tleap -f morc16_tleap.in

Minimization of morc16 in vacumm
&cntrl
  imin = 1,
  ncyc=10000,
  maxcyc = 20000,
  ntpr=100,
  ntb=0,
  igb=0,
  cut=12
/ 		heating and equilibrating in vacumm
&cntrl
   imin=0,irest=0, ntx=1,ntxo=1,
   ntpr=1000, ntwx=1000, nstlim=30000,
   dt=0.002, ntt=3, tempi=50.0,
   temp0=300, gamma_ln=1.0, ig=-1,                              ntb=0,igb=0,ntc=2, ntf=2, cut=12,
   ioutfm=1,ntrelax=100,
/
&wt type='TEMP0', istep=50.0,istep2=5000,value1=0.0,value2=300.0 /
&wt type='TEMP0', istep=5001,istep2=30000,value1=300.0,value2=300.0 /
&wt type='END' / 		production run in vacumm
&cntrl
  imin=0, irest=0, ntx=1,ntxo=1,
  ntpr=10000, ntwx=10000,ntpr=10000,nstlim=1000000,
  dt=0.002, ntt=3,
  temp0=300, gamma_ln=1.0, ig=-1,
  ntb=0,igb=0,ntc=2, ntf=2, cut=12,                                           ioutfm=1,ntrelax=100,
/

 

 

 

 

 

 

 

 

 

sander -O -i min.mdin -o morc16_va_min.out -p morc16_va.prmtop  -c morc16_va.inpcrd -r morc16_va_min.rst

sander -O -i heat.mdin -o morc16_va_heat.out -p morc16_va.prmtop -c morc16_va_min.rst -r morc16_va_heat.rst

sander -O -i md.mdin -o morc16_va_md.out -p morc16_va.prmtop -c morc16_va_heat.rst -r morc16_va_md.rst -x morc16_va_md.mdcrd

 

MORC16 in implicit solvent

Minimization of morc16 in vacumm
&cntrl
  imin = 1,
  ncyc=100000,
  maxcyc = 200000,
  ntpr=100,
  ntb=0,
  igb=1,                                           cut=12
/ 		heating and equilibrating in vacumm
&cntrl
   imin=0,irest=0, ntx=1,ntxo=1,
   ntpr=1000, ntwx=1000, nstlim=30000,
   dt=0.002, ntt=3, tempi=50.0,
   temp0=300, gamma_ln=1.0, ig=-1,
   ntb=0,igb=1,ntc=2, ntf=2, cut=12,                                                ioutfm=1,ntrelax=100,
/
&wt type='TEMP0', istep=50.0,istep2=5000,value1=0.0,value2=300.0 /
&wt type='TEMP0', istep=5001,istep2=30000,value1=300.0,value2=300.0 /
&wt type='END' / 		production run in vacumm
&cntrl
  imin=0, irest=0, ntx=1,ntxo=1,
  ntpr=10000, ntwx=10000,ntpr=10000,nstlim=1000000,
  dt=0.002, ntt=3,
  temp0=300, gamma_ln=1.0, ig=-1,
  ntb=0,igb=1,ntc=2, ntf=2, cut=12,                                             ioutfm=1,ntrelax=100,
/
 

sander -O -i min.mdin -o morc16_gb_min.out -p morc16_gb.prmtop -c morc16_gb.inpcrd -r morc16_gb_min.rst

sander -O -i heat.mdin -o morc16_gb_heat.out -p morc16_gb.prmtop -c morc16_gb_min.rst -r morc16_gb_heat.rst

sander -O -i md.mdin -o morc16_gb_md.out -p morc16_gb.prmtop -c morc16_gb_heat.rst -r morc16_gb_md.rst -x morc16_gb_md.mdcrd

 

MORC16 in explicit solvent

source leaprc.ff14SB
source leaprc.gaff
source leaprc.tip3p
loadamberparams morc16.frcmod
loadoff MOR morc16.lib
MOR =loadmol2 morc16.mol2
solvatebox morc16 TIP3PBOX 8
addions morc16 NA 0
addions morc16 CL 0
set MOR box {48 42 26}
saveamberparm MOR morc16_box.prmtop morc16_box.inpcrd
quit
 

"O" "os" 0 1 131072 28 8 -0.417600edit morc16.lib

change charge from -0.417600 to -0.421600(+0.004) to make the total charge 0. 

vmd_box_dims.sh -i morc16_box.pdb -s water

47.390998840332032, 41.118000030517579, 25.479000091552734

add  "set MOR box {48 42 26} " into morc16_tleap.in, and run tleap again.

now we are ready to run some MD of MORC16.

mpirun -np 4 pmemd.MPI -O -i min.mdin -o morc16_ex_min.out -p morc16_ex.prmtop -c morc16_ex.inpcrd -r morc16_ex_min.rst

mpirun -np 4 pmemd.MPI -O -i heat.mdin -o morc16_ex_heat.out -p morc16_ex.prmtop -c morc16_ex_min.rst -r morc16_ex_heat.rst

mpirun -np 4 pmemd.MPI -O -i md.mdin -o morc16_ex_md.out -p morc16_ex.prmtop -c morc16_ex_heat.rst -r morc16_ex_md.rst -x morc16_ex_md.mdcrd

 

Self assembly of MORC16

 

 

MORC16 Micelle Modelling

packmol <micelle_packmol.inp

tleap -f morc16_micelle_tleap.in

mpirun -np 12 pmemd.MPI -O -i min.mdin -o micelle_min.out -p morc16_micelle.prmtop -c morc16_micelle_tleap.inpcrd -r micelle_min.rst

mpirun -np 12 pmemd.MPI -O -i heat.mdin -o micelle_heat.out -p morc16_micelle.prmtop -c micelle_min.rst -r micelle_heat.rst

pmemd.cuda -O -i md.mdin -o micelle_md.out -p morc16_micelle.prmtop -c micelle_heat.rst -r micelle_md.rst -x micelle_md.mdcrd

vmd:

source ~/toolset/lipid_render.tcl

start_lipid_rendering -lipid MOR

  0ns 10ns 20ns
micelle_50lipids
micelle_100lipids
micelle_150lipids
micelle_200lipids
micelle_250lipids
mix_250lipids

 

Steered Molecular Simulation of MORC16 micelle

We will perform SMD of 10 lipids one at a time, the pulling speed is 50A/1000000steps=25A/ns.

cpptraj -p morc16_micelle.prmtop -y micelle_md.rst -x micelle_md.pdb (vmd does not read netcdf version of amber rst file)

./create_SMDINPUTS.sh 10 1000000 morc16_micelle.prmtop micelle_md.pdb

./create_job.sh 10 micelle_md.rst 1

./job.1.sh

set cell [pbc set {110 110 110} -all]

pbc box -color yellow -width 1

display resize 1000 1000

***text size 2.0 when labelling 

***quicksurf for micelle and VDW for residue 0

micelle_50
micelle_100
micelle_500
micelle_1000

for i in $(seq 1 10) ; do cp ./ASMD_$i/asmd.work.dat.1 ./asmd$i.work.dat.1; done

awk '{a[FNR]+=$4;b[FNR]++}END{for (i=1;i<=FNR;i++)print ((i-1)*0.005),a[i]/b[i]}' asmd*.work.dat.1 > average.work.dat.1

gnuplot> set style line 1 lc rgb '#0060ad' lt 1 lw 2 pi -1 ps 1.0
gnuplot> set style line 2 lc rgb '#dd181f' lt 9 lw 2 pi -1 ps 1.0
gnuplot> set style line 3 lc rgb '#29c524' lt 6 lw 2 pi -1 ps 1.0
gnuplot> set style line 4 lc rgb '#7D72F9' lt 7 lw 2 pi -1 ps 1.0

gnuplot>set ylabel "work(kj/mol)"
gnuplot>set xlabel "distance(A)"
gnuplot>set key right bottom

gnuplot>set title "PMF of Pulling One Lipid From Micelles With Different Size"

gnuplot> plot "./morc16_micelle_50/SMD/average.work.dat.1" ls 1 title '50 lipids micelle',\
>"./morc16_micelle_100/SMD/average.work.dat.1" ls 2 title '100 lipids micelle',\
>"./morc16_micelle_500/SMD/average.work.dat.1" ls 3 title '500 lipids micelle',\
>"./morc16_micelle_1000/SMD/average.work.dat.1" ls 4 title '1000 lipids micelle'

gnuplot> set terminal png size 1000,1000

gnuplot> set output 'morc16_micelle_smd.png'

 

The unit of force is kcal/mol A=69.4786 pN

The unit of work is kcal/molgnuplot>

 

 

Mixed atomistic solute and coarse grained water for morc16 micelle simulation

One challenge of the micelle simulation is the computing cost as the system goes bigger and bigger, the most time spent on simulation is nonbonded interaction, and if the the system is explicitly solvated, that means most of time of calucation is spent on nonboned interaction of bulk water molecules, which is not our interest, One way to reduce this cost is to reduce the degree of water molecules as the distance between water molecules and solute increases, there are some hybrid methods to deal with this situation, here we just replace TIP3P water model with WT4 coarse grained water model. 

A.Snapshot takenfrom a MD simulation showing the typical ordering of bulk water molecules

B.The positions of each of the oxygen atoms at the corners of the tetrahedra in A are now indicated with red beads.

C.Structural organization of WT4 in the bulk solution taken from a MD snapshot.

D.The ideal organization of a tetrahedral water cluster leads to the geometry of WT4.White and red beads (hydrogen-like, H WT4 , and
oxygen-like, O WT4 ) carry positive and negative partial charges of 0.41e, respectively.

Simulation Setup

As to the simulation protocol, the only thing need to change is solvate with "WT4BOX" instead of "TIP3PBOX", and source "/home/zwei/work/MD/MORC16/CGLA/leaprc.cgla"

tleaf -f morc16_micelle_tleap.in

mpirun -np 12 pmemd.MPI -O -i min.mdin -o micelle_min.out -p morc16_micelle.prmtop -c morc16_micelle_tleap.inpcrd -r micelle_min.rst

mpirun -np 12 pmemd.MPI -O -i heat.mdin -o micelle_heat.out -p morc16_micelle.prmtop -c micelle_min.rst -r micelle_heat.rst

pmemd.cuda -O -i md.mdin -o micelle_md.out -p morc16_micelle.prmtop -c micelle_heat.rst -r micelle_md.rst -x micelle_md.mdcrd

The simulation accelerated about ~6(11*3/4) times faster than TIP3P water model, 

./create_SMDINPUTS.sh 10 1000000 morc16_micelle.prmtop micelle_md.rst

./create_job.sh 10 micelle_md.rst 1

./job.1.sh