Biomodeling: membrane
Intro
We will do 3_membrane in the tutorial files
Hello, membrane!
In popc128_700ns, I uploaded a 700-ns simulation of a POPC bilayer system.
conf.pdb
topol.top
grompp.mdp
all.1ns.xtc: 700-ns trajectory with 1-ns time step.
You can generate your own membrane system using CHARMM-GUI web server.
topol.top
topol.top looks like this:
; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"
; bonded and non-bonded parameters for lipid
#include "../lipid21.ff/ffnonbonded.lipid21.itp"
#include "../lipid21.ff/ffbonded.lipid21.itp"
; moleculetype definition of POPC lipid
#include "../itps/popc.itp"
; Include water topology
#include "amber99sb-ildn.ff/tip3p.itp"
; Include topology for ions
#include "amber99sb-ildn.ff/ions.itp"
[ system ]
; Name
Protein
[ molecules ]
; Compound #mols
POPC 128
SOL 5736
grompp.mdp
The major difference between water box and membrane box is the pressure coupling part.
Membrane is not isotropic. Most natural choice for memrbane is semi-isotropic pressure coupling (Pxx = Pyy and Pzz are independent).
By setting Pxx = Pyy = 1 bar and Pzz = 1 bar, we can simulate zero surface tension.
; pressure coupling
Pcoupl = c-rescale
Pcoupltype = semiisotropic
nstpcouple = -1
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1.0
compressibility = 4.5e-5 4.5e-5
ref-p = 1 1
Spatial distribution of lipid and water
Using all.1ns.xtc trajectory file, we can characterize the spatial distributions of water and lipid molecules along the normal direction.
gmx density is the program for this purpose.
gmx density -f all.1ns.xtc -d Z -dens mass -sl 30 -ng 2
-d Z: density as a function of z. In other words, average over xy plane.
-dens mass: mass density.
-sl 30: the number of slices along the z axis.
-ng 2: the number of groups for analysis.
Select 2 groups to calculate density for:
Group 0 ( System) has 34360 elements
Group 1 ( Other) has 17152 elements
Group 2 ( POPC) has 17152 elements
Group 3 ( Water) has 17208 elements
Group 4 ( SOL) has 17208 elements
Group 5 ( non-Water) has 17152 elements
Select a group: 2
Selected 2: 'POPC'
Select a group: 4
Selected 4: 'SOL'
Plot density.xvg
gnuplot> pl 'density.xvg' u 1:2 w l title "POPC", 'density.xvg' u 1:3 w l title "water"
The membrane has a hydrophobic core of ~3 nm width that has no water molecules at all.
What is the chance that a water molecule crosses the membrane?
We can answer that by calculating ∆G using the umbrella sampling and WHAM.
Diffusion of water and lipid
Using all.1ns.xtc trajectory file, we can compare the diffusion coefficients of water and lipid molecules.
Obviously, water diffusion is much faster.
Homework: calculate the diffusion coefficients of water and lipid by using gmx msd.
Umbrella sampling of a water molecule as a function of the z-position
In the ∆G calculation for a water-water pair, the reaction coordinate was the distance between the water pair.
In this case, the reaction coordinate is the z-position of a water molecule with respect to the center of mass of the membrane.
Umbrella sampling
First, we need to generate the index file
We will just select the first water molecule. The selection will not change the results. All water molecules are equal!
There are 128 POPC molecules. Thus, the residue index of the first water molecules will be 129.
gmx make_ndx -f topol.tpr -o pull.ndx
ri 129
q
Caution: The initial distance between chosen water molecules must be smaller than half of the box. Otherwise, it violated PBC.
The groups will look like this:
0 System : 34360 atoms
1 Other : 17152 atoms
2 POPC : 17152 atoms
3 Water : 17208 atoms
4 SOL : 17208 atoms
5 non-Water : 17152 atoms
6 r_129 : 3 atoms
Perform simulation in each window
10 Å window
See PULL section in pull.10.mdp
; COM PULLING
pull = yes
; Cylinder radius for dynamic reaction force groups (nm)
pull-pbc-ref-prev-step-com = yes
pull-nstxout = 50
pull-nstfout = 50
pull-ngroups = 2
pull-ncoords = 1
pull-group1-name = POPC
pull-group1-pbcatom = 1
pull-group2-name = r_129
pull-coord1-type = umbrella
pull-coord1-geometry = direction
pull-coord1-vec = 0 0 1
pull-coord1-groups = 1 2
pull-coord1-dim = N N Y
pull-coord1-init = 1.0 ; the equilibrium length of spring is 10 Å
pull-coord1-k = 1000 ; spring constant in kJ/mol nm
Generate pull.10.tpr for gmx mdrun
gmx grompp -c 2.pdb -f pull.10.mdp -n pull.ndx -o pull.10.tpr -t state.cpt
-t state.cpt: use information from the check point file.
Perform MDRUN for 1 ns
gmx mdrun -s pull.10.tpr -nt 4 -gpu_id 0 -cpi -noappend -deffnm pull.10 -nsteps 500000
pullx files contain the length of spring
pullf files contain the force on the spring (F = -kx)
Repeat same things for other windows. For example, 0, 2, 4, 6, 8 , 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34 Å
It will take much longer than the water box. Adjust the simulation time considering your computer.