I have been trying to carryout simulation of cell penetration of a peptide by pulling it across a lipid bilayer. But as the peptide is moving into the bilayer the bilayer is completely getting distorted and getting break into single bilayer.
Can it be the pull rate which may affect the system or not restraining the lipid bilayer may cause a problem? What should I do to rectify it.? Please help me out.
The md_pull.mdp parameters i have attached with the file and to clarify my point i have attached an image too.
title = umbrella pulling MD
define = -DPOSRES
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000 ; 2 * 500000 = 1000 ps (1 ns)
dt = 0.002 ; 2 fs
; Output control
nstxout = 1000 ; save coordinates every 2 ps
nstvout = 1000 ; save velocities every 2 ps
nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps
nstenergy = 1000 ; save energies every 2 ps
nstlog = 1000 ; update log file every 2 ps
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cels
nstlist = 5 ; 10 fs
rlist = 1.2 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = Nose-Hoover ; More accurate thermostat
tc-grps = Protein DPPC SOL_CL ; three coupling groups - more accurate
tau_t = 0.5 0.5 0.5 ; time constant, in ps
ref_t = 310 310 310 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
; Pull code
pull = umbrella
pull_geometry = distance ; simple distance increase
pull_dim = N N Y
pull_start = yes ; define initial COM distance > 0
pull_ngroups = 1
pull_group0 = DPPC
pull_group1 = Protein
pull_rate1 = 0.1 ; 0.1 nm per ps = 100 nm per ns
pull_k1 = 1000 ; kJ mol^-1 nm^-2
pull_init1 = 0
pull_vec1 = 0 0 1
The image given below is the is the peptide in lipid bilayer at different interval of time. The image showed that at time t=0ps the peptide is at bilayer and water interface. But once it starts moving in the second lipid layer of bilayer started going away at t=60ps and t =90ps.
Hope so now image is clear..
What exactly does the image show? Different snapshots over time? If the structure at the top is any indication, it is very clear that the unit cell is far too small for doing pulling across the membrane. Restraining the lipids in any way can cause problems. Pulling a peptide across the membrane requires a very slow pull to allow the lipids to deform sensibly around the peptide. A pull rate of 0.1 is way too fast.
Doing umbrella sampling with a non-zero pulling rate is probably not what you want. To get a free energy profile of transfer of the ligand through your membrane, you should simulate several "windows", which differ in their equilibrium positions (different pull_init1 values). Evaluation by WHAM is a common option.
Even though I don't understand your image completely, it may also be a problem of visualization, rather then of simulation...
Hello,
The image given below is the is the peptide in lipid bilayer at different interval of time. The image showed that at time t=0ps the peptide is at bilayer and water interface. But once it starts moving in the second lipid layer of bilayer started going away at t=60ps and t =90ps.
Hope so now image is clear..
Michal,
Thanking you for your suggestion but can you tell me in slight detail reagrding this free energy and WHAM??
Umbrella sampling is covered very well in the literature. It is used extensively. As far as conducting such simulations in Gromacs, see my tutorial:
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html
The box is clearly too small to calculate a reliable free energy, and as Michal says, the actual umbrella sampling is done in separate simulations. What you're doing now is a steered MD protocol that generates configurations along a path that will later be sampled during US.
If you're going to try pulling a peptide across the membrane, this is a _MAJOR_ undertaking, which will require a huge amount of computer time to equilibrate. Check out the paper by Neale et al, Biophysical Journal, 2014, v106, ppL29-31
- Badges
- Science topic
- Similar topics
- Bioinformatics
- Bioinformatics and Computational Biology