membrane simulation

_episodes/17-Preparing_and_Simulating_Bilayer.md

@@ -18,9 +18,11 @@ keypoints:
 
 This in-depth tutorial guides you through the process of building a bilayer, preparing simulation input files, minimizing energy, equilibrating the system, and running an equilibrium molecular dynamics simulation. You will need to follow a number of steps to complete the tutorial. If you are already familiar with some of these topics, you can skip them and focus on the ones you don't know about.
 
-### 1. A bilayer of lipids is formed by packing them together.
+### 1. Generating a bilayer by packing lipids together.
 
-It has been shown that the composition of the phospholipids in the erythrocyte membrane is asymmetric. Here is an article that reviews this topic in detail:
+The packmol-memgen program allows the creation of asymmetric bilayers with leaflets composed of different lipid species. 
+
+Bilayer asymmetry is a common feature of biological membranes. For example, the composition of the phospholipids in the erythrocyte membrane is asymmetric. Here is an article that reviews this topic in detail:
 [Interleaflet Coupling, Pinning, and Leaflet Asymmetry—Major Players in Plasma Membrane Nanodomain Formation](https://www.frontiersin.org/articles/10.3389/fcell.2016.00155/full). 
 
 Asymmetric lipid bilayer can be generated by the following submission script:
@@ -43,44 +45,42 @@ packmol-memgen \
 ~~~
 {: .language-bash}
 
-
 The job takes about 30 minutes. It generates several output files:  
 - AMBER topology         bilayer_only_lipids.top 
 - AMBER coordinates      bilayer_only_lipids.crd
 - PDB-formatted system   bilayer_only_lipids.pdb 
 
-By default parametrization uses LIPID17, ff14SB and TIP3P force fields.
+By default LIPID17, ff14SB and TIP3P force fields are used.
+
+![Figure: Initial bilayer]({{ page.root }}/fig/bilayer_1.png){: width="500" }
 
-Figure:
-Neutral: POPC green, POPE blue, PSM white
-Charged: POPS red
+Figure 1. Initial bilayer. Neutral lipids: POPC (green), POPE (blue), PSM (white). Charged lipids: POPS (red).
 
-### 2. Preparing the bilayer for simulation with a different force field. 
+### 2. Using AMBER force fields not available natively in packmol-memgen. 
 
-The file bilayer_only_lipids.pdb can be used to create topology files using any AMBER force field distributed with AmberTools. This allows to prepare a simulation system using a force field that is not provided with packmol-memgen. 
+Packmol-memgen offers several recent AMBER force fields that can be used to parameterize a bilayer. However, sometimes it is necessary to use force fields that are not readily available in packmol-memgen. The parameterization of a prepared bilayer with tleap directly is a way to achieve this. To use tleap you will need to provide an input file that is in pdb format compatible with AMBER. 
 
-However, if your system contains SPM this file will not work because packmol-memgen erroneoiusly adds TER record between the lipid headgroup (SPM) and second acyl tail (SA). This splits PSM molecule in 2 parts: PA+SPM and SA.  
+Even though packmol-memgen creates a pdb file (bilayer_only_lipids.pdb), you should use it with caution. It may not work properly with some molecules. PSM is an example. When saving pdb file containing this type of lipids packmol-memgen erroneously adds TER record between the lipid headgroup (SPM) and second acyl tail (SA). This splits PSM molecule in two: PA+SPM and SA.   
 
-For this reason we will, create a pdb file from AMBER topology using parmed and then use tleap to reparameterize it using a polarizable water model OPC3-Pol which is not yet available in packmol-memgen natively. 
+To illustrate the parameterization process, we will create a pdb file from AMBER topology and coordinate files, and then parameterize the bilayer using the polarizable water OPC3-pol not yet available natively in packmol-memgen.
+
+Begin by loading the ambertools/23 module and starting the Python interpreter:
 
 ~~~
 module load StdEnv/2020 gcc/9.3.0 openmpi/4.0.3 ambertools/23
 python
 ~~~
 {: .language-bash}
 
+In the python prompt enter the following commands:
 ~~~
 import parmed
 amber = parmed.load_file("bilayer_only_lipid.top", "bilayer_only_lipid.crd")
-amber.save("bilayer_only_lipid_fixed.pdb")
+amber.save("bilayer-lipid21-opc3pol.pdb")
 quit()
 ~~~
+{: .language-python}
 
-Alternatively you can fix the file bilayer_only_lipid.pdb manually using standard unix commands grep and awk. The following command removes these TER records:
-~~~
-grep -n  "H22 SPM " bilayer_only_lipid.pdb | cut -f1 -d: |  awk '{print $0 + 1}' | awk -F, 'NR==FNR { nums[$0]; next } !(FNR in nums)' - bilayer_only_lipid.pdb > bilayer_only_lipid_fixed.pdb
-~~~
-{: .language-bash}
 
 We only need to know one more thing before we can proceed: the size of the periodic box. It can be extracted from the file bilayer_only_lipid.top:
 
@@ -103,7 +103,7 @@ module load StdEnv/2020 gcc/9.3.0 openmpi/4.0.3 ambertools/23
 tleap -f <(cat << EOF
 source leaprc.water.opc3pol
 source leaprc.lipid21 
-sys=loadpdb bilayer_only_lipid_fixed.pdb
+sys=loadpdb bilayer-lipid21-opc3pol.pdb
 set sys box {103.5 103.7 99.7}
 saveamberparm sys bilayer.parm7 bilayer.rst7
 quit
@@ -113,6 +113,36 @@ EOF)
 
 Don't forget to replace {103.5 103.5 99.9} with your own numbers.
 
+
+> ## Use shell commands to delete the erroneous TER records
+>This exercise is optional for those who wish to improve their efficiency when working with text files. Using shell commands at an advanced level is the focus of the exercise. 
+>Using shell commands delete the erroneous TER records in the file bilayer_only_lipid.pdb.
+>
+>>## Solution
+>>The following command removes all erroneous TER records:
+>>
+>>~~~
+>>grep -n  "H22 SPM " bilayer_only_lipid.pdb | \
+>>  cut -f1 -d: | \
+>>  awk '{print $0 + 1}' | \
+>>  awk -F, 'NR==FNR { nums[$0]; next } !(FNR in nums)' \
+>>  - bilayer_only_lipid.pdb > bilayer_only_lipid_fixed.pdb
+>>~~~
+>>{: .language-bash}
+>>
+>>Let's break it down:
+>>-  `grep -n  "H22 SPM "`  finds all lines containing "H22 SPM " and 
+prints their line numbers followed by the line content.
+>>- `cut -f1 -d:` selects the first field containing the line number
+>>- `awk '{print $0 + 1}` increments this number, after this command we have a list of lie numbers we want to delete
+>>- `awk -F, 'NR==FNR { nums[$0]; next } !(FNR in nums)' - bilayer_only_lipid.pdb` takes the list of line numbers from the pipe and deletes lines with matching numbers from the file bilayer_only_lipid.pdb
+>>
+>>If this command is executed, it will remove the line following the last atom of each SPM residue, effectively deleting the erroneous TER records.
+>>
+> {: .solution}
+{: .challenge}
+
+
 ### 3. Energy minimization
 A system must be optimized in order to eliminate clashes and prepare it for molecular dynamics. 
 
@@ -121,6 +151,13 @@ To run energy minimization we need three files:
 - Initial coordinates      bilayer.rst7
 - Minimization input file  minimization.in  
 
+First, we create a directory for energy minimization job.
+
+~~~
+mkdir 1-minimization && cd 1-minimization
+~~~
+{: .language-bash}
+
 Here is an input file for minimization (minimization.in) that you can use.
 ~~~
 Minimization
@@ -135,42 +172,53 @@ Minimization
  ~~~
  {: .file-content}
 
- Using the script below, you can submit a minimization job.
- ~~~
- #!/bin/bash
-#SBATCH --mem-per-cpu=4000 --ntasks=10 --time=6:0:0 
+Using the script below, you can submit a minimization job.
+~~~
+#!/bin/bash
+#SBATCH --mem-per-cpu=4000 --ntasks=10 --time=1:0:0 
 
 module purge
-module load StdEnv/2020  gcc/9.3.0  cuda/11.4  openmpi/4.0.3 amber/20.12-20.15
+module load StdEnv/2020 gcc/9.3.0 openmpi/4.0.3 ambertools/23
 
-srun sander.MPI -O -i min.in \
+srun sander.MPI -O -i minimization.in \
 -o minimization.out \
 -p ../bilayer.parm7 \
 -c ../bilayer.rst7 \
 -r minimized.rst7 
  ~~~
 {: .language-bash}
 
-When the optimization job is completed, the file minimized.rst7 will contain the optimized coordinates that were produced.
+When the optimization job is completed, the file minimized.rst7 will contain the optimized coordinates that were produced. Minimization job completes in about 40 min.
 
 ### 4. Heating
+
+~~~
+cd ..
+mkdir 2-heating && cd 2-heating
+~~~
+{: .language-bash}
+
+Our optimized simulation system does not move yet, there is no velocities. Our next step is to heat it up to room temperature. The system can be heated in a variety of ways with different initial temperatures and heating rates. It's important to warm up the system gently, avoiding strong shocks that can create conformations that aren't physiological. A typical protocol is described in Amber tutorial [Relaxation of Explicit Water Systems](https://ambermd.org/tutorials/basic/tutorial13/index.php).
+
+We will start the system at the low temperature of 100K and heat it up to 300K over 20 ps of simulation time at constant pressure. At this step we use weak coupling thermostat and barostat with both time constants set to  5 ps.
+
 ~~~
 Heating 
  &cntrl
-  imin=0,         	! Molecular dynamics
-  ntx=1,          	! Read coordinates with no initial velocities
-  ntc=2,          	! SHAKE on for bonds with hydrogen
-  ntf=2,         	! No force evaluation for bonds with hydrogen
-  tol=0.0000001,  	! SHAKE tolerance
-  nstlim=20000,   	! Number of MD steps
-  ntt=1,          	! Berendsen thermostat
-  tempi=100,       	! Initial temperature
+  imin=0,           ! Molecular dynamics
+  ntx=1,            ! Read coordinates with no initial velocities
+  ntc=2,            ! SHAKE on for bonds with hydrogen
+  ntf=2,            ! No force evaluation for bonds with hydrogen
+  tol=0.0000001,    ! SHAKE tolerance
+  nstlim=20000,     ! Number of MD steps
+  ntt=1,            ! Berendsen thermostat
+  tempi=100,        ! Initial temperature
   temp0=300,        ! Target temperature
   tautp=5.0,        ! Time constant for heat bath coupling
   ntpr=100,         ! Write energy info every ntwr steps 
   ntwr=10000,       ! Write restart file every ntwr steps
-  ntwx=2000,      	! Write to trajectory file every ntwx steps
-  dt=0.001,       	! Timestep 
+  ntwx=2000,        ! Write to trajectory file every ntwx steps
+  dt=0.001,         ! Timestep 
   ntb=2,            ! Constant pressure periodic bc.
   barostat=1,       ! Berendsen barostat
   taup=5,           ! Pressure relaxation time
@@ -184,10 +232,10 @@ Heating
 Submission script
 ~~~
 #!/bin/bash
-#SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=6:0:0 
+#SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=1:0:0 
 
 module purge
-module load StdEnv/2020  gcc/9.3.0  cuda/11.4  openmpi/4.0.3 amber/20.12-20.15
+module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 amber/20.12-20.15
 
 pmemd.cuda -O -i heating.in \
 -o heating.out \
@@ -197,24 +245,29 @@ pmemd.cuda -O -i heating.in \
 ~~~
 {: .language-bash}
 
+A minute is all it takes for this job to be completed.
+
 ### 5. The first phase of the equilibration process.
+This second step of relaxation will maintain the system at 300 K over 100 ps of simulation time at a constant pressure, allowing the box density to relax. The initial density is too low because there are gaps between lipids an water and lipids are not perfectly packed. 
+As the simulation parameters are the same as on the previous step, you might wonder why we did not simply continue with the previous step?  It is necessary to restart because the box size changes too much from its original size, and a running simulation cannot handle such large changes.
+
 ~~~
 Equilibration 1 
  &cntrl
-  imin=0,         	! Molecular dynamics
+  imin=0,           ! Molecular dynamics
   irest=1,          ! Restart
-  ntx=5,          	! Read positions and velocities
-  ntc=2,          	! SHAKE on for bonds with hydrogen
-  ntf=2,         	! No force evaluation for bonds with hydrogen
-  tol=0.0000001,  	! SHAKE tolerance
-  nstlim=100000,   	! Number of MD steps
-  ntt=1,          	! Berendsen thermostat
+  ntx=5,            ! Read positions and velocities
+  ntc=2,            ! SHAKE on for bonds with hydrogen
+  ntf=2,            ! No force evaluation for bonds with hydrogen
+  tol=0.0000001,    ! SHAKE tolerance
+  nstlim=100000,    ! Number of MD steps
+  ntt=1,            ! Berendsen thermostat
   temp0=300,        ! Set temperature
   tautp=5.0,        ! Time constant for heat bath coupling
   ntpr=100,         ! Write energy info every ntwr steps 
   ntwr=10000,       ! Write restart file every ntwr steps
-  ntwx=2000,      	! Write to trajectory file every ntwx steps
-  dt=0.001,       	! Timestep (ps)
+  ntwx=2000,        ! Write to trajectory file every ntwx steps
+  dt=0.001,         ! Timestep (ps)
   ntb=2,            ! Constant pressure periodic bc.
   barostat=1,       ! Berendsen barostat
   taup=5,           ! Pressure relaxation time
@@ -228,37 +281,42 @@ Equilibration 1
 Submission script
 ~~~
 #!/bin/bash
-#SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=6:0:0 
+#SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=1:0:0 
 
 module purge
-module load StdEnv/2020  gcc/9.3.0  cuda/11.4  openmpi/4.0.3 amber/20.12-20.15
+module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 amber/20.12-20.15
 
 pmemd.cuda -O -i equilibration.in \
 -o eq-1.out \
 -p ../bilayer.parm7 \
--c heated.rst7 \
+-c ../2-heating/heated.rst7 \
 -r eq-1.rst7 \
 ~~~
 {: .language-bash}
 
+This job completes in 5 about minutes.
+
 ### 6. The second phase of the equilibration process.
+
+Now density is mostly relaxed and we can apply temperature and pressure coupling parameters appropriate for a production run.  The relaxation of a bilayer may take many nanoseconds.
+
 ~~~
 Equilibration 2 
  &cntrl
-  imin=0,         	! Molecular dynamics
+  imin=0,           ! Molecular dynamics
   irest=1           ! Restart
   ntx=5,            ! Read  positions and  velocities
-  ntc=2,          	! SHAKE on for bonds with hydrogen
-  ntf=2,         	! No force evaluation for bonds with hydrogen
-  tol=0.0000001,  	! SHAKE tolerance
-  nstlim=2500000,   ! Number of MD steps
+  ntc=2,            ! SHAKE on for bonds with hydrogen
+  ntf=2,            ! No force evaluation for bonds with hydrogen
+  tol=0.0000001,    ! SHAKE tolerance
+  nstlim=2000000,   ! Number of MD steps
   ntt=3,            ! Langevin thermostat 
   temp0=300,        ! Target temperature
   gamma_ln=2.0,     ! Collision frquency 
   ntpr=100,         ! Write energy info every ntwr steps 
   ntwr=10000,       ! Write restart file every ntwr steps
-  ntwx=10000,      	! Write to trajectory file every ntwx steps
-  dt=0.002,       	! Timestep (ps)
+  ntwx=10000,       ! Write to trajectory file every ntwx steps
+  dt=0.002,         ! Timestep (ps)
   ntb=2,            ! Constant pressure periodic bc.
   barostat=1,       ! Berendsen barostat
   taup=1.0,         ! Pressure relaxation time
@@ -275,12 +333,12 @@ Submission script
 #SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=6:0:0 
 
 module purge
-module load StdEnv/2020  gcc/9.3.0  cuda/11.4  openmpi/4.0.3 amber/20.12-20.15
+module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 amber/20.12-20.15
 
 pmemd.cuda -O -i equilibration-2.in \
 -o eq-2.out \
 -p ../bilayer.parm7 \
--c eq-1.rst7 \
+-c ../3-equilibration/eq-1.rst7 \
 -r eq-2.rst7 \
 ~~~
 {: .language-bash}