I was trying to simulate a 512-DPPC membrane after developing it from Tieleman's model of 128-DPPC bilayer. I removed the initial periodicity and the water molecules using trjconv as suggested by Justin Lemkul in his tutorial of membrane protein simulation (KALP15 in DPPC probably). I used genconf to place the 128-DPPC membrane sidewise in a 2(x) 2(y) 1(z) format.
The first picture i.e. of nice bilayer is after energy minimization using a steepest descent algorithm with a gradient of 100 kJ/mol/nm. It is only after NVT equilibration with a high position restraint applied on the DPPC residues, it is making a hole inside the membrane i.e. the bilayers tear apart.
I tried this several times changing the temperature (to change the membrane fluidity) and the gradient too (I primarily used 500 kJ/mol/nm and then reduced it to 100).
Any idea on how it can be solved? I see quite a few water molecules inside the bilayers which go out during equilibration. Also, the ions show a tendency to make a cluster at one corner of the membrane.
Please help...
I assume in the second figure, you have used unwrapped trajectories. Then,
to me, it seems like an in-plane pressure build-up problem. Have you tried applying a barostat, either isotropically or just in-plane? You have done an energy minimization first, which corresponds to minimum potential energy and equilibrium at 0 K; but when you apply a finite temperature, the membrane undergoes thermal expansion. With the fixed box size in an NVT simulation, it is reasonable to assume the membrane will "buckle" the way it seems it does.
Usually, separation of monolayer leaflets in a bilayer happens if the water density is less or a semi-isotropic pressure coupling is not used.
For your case, I would simply generate the initial structure of the system with charmm-gui and equilibrate it according to the mdp files provided with it. It seems to work fine for my system (similar to yours).
@Swapnil Wagle, I did not use the CHARMM-GUI just because I wanted to use the Berger lipid parameters which seems to be perfect for my work (as I found from the literatures).
@all,
Should I run the NPT first in order to apply the semiisotropic scheme??? I usually do the NVT first and then NPT. Is it wrong???
Sorry about the previous answer, it was a software issue.
Your problem may be a poorly packed membrane, so lipid density is low and there are not enough lipids to form a full lipid bilayer completely filling the box. Since the smallest surface of apolar molecules exposure to water should be a cylinder, the simulation is showing you this fact. Try increasing the lipid density at the time of membrane creation, by packing the lipids more closely (i.e with packmol, slow but reliable, you may find it here: http://m3g.iqm.unicamp.br/packmol/citation.shtml
). Or instead, you may use a NPT short simulation to adjust the box size, pressurize the membrane and squeeze water out of the transmembrane region, and then you can run the NVT followed by the NPT as you needed. Do not stick to protocols just for imitation, read at the principles behind each protocol.
Best wishes,
Rogelio
@Souparno: In that case, after energy minimisation I would start equilibrating the system directly with NPT with semi-isotropic P-coupling (without doing the NVT first). Also, while equilibration, I would stick to the Berendsen Barostat rather than Parrinello-Rahman.
The Berendsen Barostat needs frequent V-rescaling to produce a canonical ensemble, make sure to apply it.
Best wishes
Rogelio
- Badges
- Science topic
- Similar topics
- Bioinformatics
- Bioinformatics and Computational Biology