normal md simulation code

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 8 months ago #1097 by dujiangfeng
md simulation code was created by dujiangfeng
Dear Everyone,

I am doing a martini simulation of my protein (159 AA) and membrane (DOPS/DOPC). The protein is supposed to bind to membrane because of hydrophobic loop and electrostatics interaction.

However, When I am running it by the recommended parameters here (shifted cutoff for both coulomb and vdw), the protein seems never to move the membrane. But if I use PME for coulomb (PME is a not reasonable choice for normal waters) and shifted cutoff for vdw, the protein went to membrane quickly, which actually is what I want to see. Then here a question comes: How dangerous it is if I pick up the "wrong" parameters for simulation? It anyway can generate the binding of protein and membrane.

Moreover, I don't know why the cg model which was created by martinize.py is very unstable during the simulation. The conformation is almost broken by MD.

I tried six ways, but neither was stable.
1, atomistic model from PDB database --> martinize.py --> cg.model --> it is broken during MD.
2, atomistic model from PDB --> EM & MD --> martinize.py --> cg.model --> broken during MD.
3, atomistic model from PDB --> martinize.py --> cg.model --> EM --> broken during MD.
4, atomistic model from PDB --> EM & MD --> martinize.py --> cg.model --> EM --> broken during MD.
5, atomistic model from PDB --> EM & MD --> martinize.py + elastic network --> cg.model --> EM --> broken during MD.
6, atomistic model from PDB --> EM & MD --> martinize.py + elastic network --> cg.model --> broken during MD.

I appreciate any reply about my topic.

Thank you very much,

Jiangfeng.

Please Log in or Create an account to join the conversation.

More
11 years 8 months ago #1098 by xavier
Replied by xavier on topic md simulation code
The deformation of the protein using regular Martini is expected, especially if the protein is large and a complex mix of secondary structure elements. B-sheet structures are extremely unstable.

An alternative to regular Martini for large protein that deform substantially is the ElNeDyn approach, which will define an elastic network based on the Calpha atoms of the protein. It is as an option of the Martinize script.

For the PME/Cutoff difference you observe, this might be the case if you have the protein far from the bilayer. The use of PME is not recommended as we have observed strange behavior for highly charged systems. We are trying to determine the exact cause of the behavior but in the mean time we recommend the use of cufoff.

Have you tried the polarizable water model? if your protein/bilayer are charged it might help improve the behavior of the system.

I hope this helps.
XAvier.

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 8 months ago #1100 by dujiangfeng
Replied by dujiangfeng on topic md simulation code
Hi Xavier,

Thank you for your reply. In fact, the proteins is a beta-barrel(comprised of beta-sheet), which is probably the reason deforming the protein by coarse graining it. Actually, I have also tried the ElNeDyn approach as two of my six different ways to make the cg model. Unfortunately, the protein deformed still. I don't know if I can increase the rubber force directly from itp file to stable the protein conformation?


It is right that the distance between my protein and membrane is pretty big as ~3 nm. If I really want to keep the distance as such big, what else should I do since the shifted cutoff of coulomb and vdw didn't drive the protein to membrane? Meanwhile, I am going to try the polarized water system.

Thank you again,

Jiang.

Please Log in or Create an account to join the conversation.

More
11 years 8 months ago #1101 by djurre
Replied by djurre on topic md simulation code
Does your protein have more then one chain (with different secondary structures)? In that case it might be that the topology created by martinize is not correct and you want should download the fixed version (martinize.py 2.1) which you can find on the website.

It is really strange the protein deforms even when using Elnedyn. What do you mean by "the protein is broken"? Are the different strands of the barrel falling apart? Increasing the force constant won't help (if you use the 'normal' value already).

I agree with Xavier it would be best not to use PME. If you want to keep the distance big, you just have to wait long, eventually it will diffuse to close to the membrane.

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 8 months ago #1103 by dujiangfeng
Replied by dujiangfeng on topic md simulation code
Hi Djurre,

djurre wrote: Does your protein have more then one chain (with different secondary structures)? In that case it might be that the topology created by martinize is not correct and you want should download the fixed version (martinize.py 2.1) which you can find on the website.


No, my protein is only in one chain, which is comprised by beta-sheets and loops. I will try the new version now.

djurre wrote: It is really strange the protein deforms even when using Elnedyn. What do you mean by "the protein is broken"? Are the different strands of the barrel falling apart? Increasing the force constant won't help (if you use the 'normal' value already).



Sorry, "the protein is broken" is not a correct description of my situation. During the simulation, the conformation of the protein has been deforming from a structured form to a round-shape form. So I couldn't see where is the strand/loop region.

djurre wrote: I agree with Xavier it would be best not to use PME. If you want to keep the distance big, you just have to wait long, eventually it will diffuse to close to the membrane.


I had tried to perform several long simulation (250 nm) systems (there are identical) with cutoff. Well, the proteins performed differently in different systems. Some proteins went to opposite direction of membrane, while some proteins rotated themselves in the solvent. All of them seems no interest to membrane. The only common point of those systems is the deforming of the protein. :(


Thank you!

Jiang.

Please Log in or Create an account to join the conversation.

More
11 years 8 months ago #1105 by xavier
Replied by xavier on topic md simulation code

dujiangfeng wrote: Hi Djurre,

I had tried to perform several long simulation (250 nm) systems (there are identical) with cutoff. Well, the proteins performed differently in different systems. Some proteins went to opposite direction of membrane, while some proteins rotated themselves in the solvent. All of them seems no interest to membrane. The only common point of those systems is the deforming of the protein. :(

Jiang.


You probably mean 250 ns. This is not very long ... a few things you could try:
1- First the protein with ElNeDyn SHOULD not deform. If it does it means that either something funky happened during the construction of the topology or as Djurre mention you system is not as simple as a protein and the script did not understand it.
Please, precise what you mean by deformation and if the "protein" you simulate is a single chain protein or a multiple chain ...
2- Once you have a reasonable topology you could use a conformation of the system where the protein is in contact with the lipids (from a PME simulation) and switch to a cutoff scheme. This would tell you if the cutoff scheme can stabilize the protein/membrane interaction.

First make sure your topology is correct.

XAvier.

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 8 months ago #1106 by dujiangfeng
Replied by dujiangfeng on topic md simulation code
Dear XAvier,

I am sorry I made a mistake when I was replying Djurre's Post yesterday by quoting his comments. I had answered his questions after every paragraph of his post, so some of my sentences were in the quotation also. I just edited it.

XAvier wrote: You probably mean 250 ns. This is not very long ... a few things you could try:


Sorry again. It is 250 ns.

XAvier wrote: 1- First the protein with ElNeDyn SHOULD not deform. If it does it means that either something funky happened during the construction of the topology or as Djurre mention you system is not as simple as a protein and the script did not understand it.


The script works find to process my protein with ElNeDyn approach, and there is no error reported by the script. My protein is a beta-barrel, and it is comprised of beta-sheet and loops, which might be the problem for the script, is it?

XAvier wrote: Please, precise what you mean by deformation and if the "protein" you simulate is a single chain protein or a multiple chain ...


the deformation of my protein means that the original structured beta barrel was being changed into a round-shaped object. Well, my protein is in one chain.

XAvier wrote: 2- Once you have a reasonable topology you could use a conformation of the system where the protein is in contact with the lipids (from a PME simulation) and switch to a cutoff scheme. This would tell you if the cutoff scheme can stabilize the protein/membrane interaction.


Do you mean, at the beginning of the simulation, I may use PME to make protein in contact with lipids and then switch the simulation code to the shifted cutoff to continue the simulation for a long time, right?

XAvier wrote: First make sure your topology is correct.


The topology seems ok since there is no error at all.

Thank you very much for your reply,

Jiang.

Please Log in or Create an account to join the conversation.

More
11 years 8 months ago #1107 by xavier
Replied by xavier on topic md simulation code
The ElNeDyn topology does not allow the protein to deform more than according to atomic fluctuations observed in atomistic simulations. Thus the protein should not have an RMSD higher than 0.2 nm from the experimental structure.

The script might not detect the problem but it is sure you have one from the behavior you describe! a B-barrel going to a blob is not reasonable.

Yes, I meant starting with PME to generate contact and then turn on the cutoff scheme to test if in that case the protein does prefer to solvate in water.

Note that B-barrel proteins are often trans-membrane proteins, is this the case of yours?

You need to rebuild your topology checking that the atomistic structure you give is correct (no atoms are missing) and the RMSD should not go higher than 0.2 nm. If the problem persist get back to us with the atomistic pdb and CG topology. But please try to fix the problem yourself first. Notably edit the CG topology and check that all the parameters make sense.

XAvier.

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 7 months ago #1114 by dujiangfeng
Replied by dujiangfeng on topic md simulation code
Dear XAvier,

xavier wrote: The ElNeDyn topology does not allow the protein to deform more than according to atomic fluctuations observed in atomistic simulations. Thus the protein should not have an RMSD higher than 0.2 nm from the experimental structure.


Here if I am understood correctly, ElNeDyn restricts each bead of cg model within an RMSD less than 0.2nm from its initial structure, right?

xavier wrote: The script might not detect the problem but it is sure you have one from the behavior you describe! a B-barrel going to a blob is not reasonable.


Since I updated my DSSP database, and also I used another x-ray structure to make cg model with ElNeDyn, the structure has been performing fine. Command "g_rms -s pro1_psys1iqd.tpr -f pro1_psys1iqd.trr -o rms_psys1iqd.xvg" shows the protein RMSD from 0 to 0.25 nm (keep stable there). However, the previous x-ray structure I used before is almost same with the one I used now by Ca superposition.

xavier wrote: Yes, I meant starting with PME to generate contact and then turn on the cutoff scheme to test if in that case the protein does prefer to solvate in water.


When I used PME in normal water MD, the protein came to membrane quickly (300 ps moving 3 nm), then when I changed into shifted cutoff scheme, the protein kept there during the whole cutoff simulation (200 ns).
When I used PME in the polarized water MD, the protein kept in the water but didn't come to membrane (it was only 20 ns so far, maybe it is too short).

xavier wrote: Note that B-barrel proteins are often trans-membrane proteins, is this the case of yours?


The protein is supposed to bind to membrane but not totally go into membrane since it is just one part of the huge coagulant factor VIII.

xavier wrote: You need to rebuild your topology checking that the atomistic structure you give is correct (no atoms are missing) and the RMSD should not go higher than 0.2 nm. If the problem persist get back to us with the atomistic pdb and CG topology. But please try to fix the problem yourself first. Notably edit the CG topology and check that all the parameters make sense.

XAvier.


The new x-ray structure is good for cg model, I suppose.


Thank you,

Jiang.

Please Log in or Create an account to join the conversation.

More
11 years 7 months ago #1117 by xavier
Replied by xavier on topic md simulation code

dujiangfeng wrote: Dear XAvier,

xavier wrote: The ElNeDyn topology does not allow the protein to deform more than according to atomic fluctuations observed in atomistic simulations. Thus the protein should not have an RMSD higher than 0.2 nm from the experimental structure.


Here if I am understood correctly, ElNeDyn restricts each bead of cg model within an RMSD less than 0.2nm from its initial structure, right?

ElNeDyn defines an elastic network between the backbone beads (Calpha) and let the proteins fluctuate around the initial structure leaving enough flexibility to reproduce the atomistic behavior. Please have a look at the paper. (periole et JCTC-2009-5, 2531–2543).

dujiangfeng wrote:

xavier wrote: The script might not detect the problem but it is sure you have one from the behavior you describe! a B-barrel going to a blob is not reasonable.


Since I updated my DSSP database, and also I used another x-ray structure to make cg model with ElNeDyn, the structure has been performing fine. Command "g_rms -s pro1_psys1iqd.tpr -f pro1_psys1iqd.trr -o rms_psys1iqd.xvg" shows the protein RMSD from 0 to 0.25 nm (keep stable there). However, the previous x-ray structure I used before is almost same with the one I used now by Ca superposition.

That sounds much more reasonable. It is not all about the Calpha position but also how the chains are defined, how you clean the file, if the side chains are complete ... a lot of things can go wrong ...

dujiangfeng wrote:

xavier wrote: Yes, I meant starting with PME to generate contact and then turn on the cutoff scheme to test if in that case the protein does prefer to solvate in water.


When I used PME in normal water MD, the protein came to membrane quickly (300 ps moving 3 nm), then when I changed into shifted cutoff scheme, the protein kept there during the whole cutoff simulation (200 ns).
When I used PME in the polarized water MD, the protein kept in the water but didn't come to membrane (it was only 20 ns so far, maybe it is too short).

20 ns is short ... not that if your protein is binding the membrane it is probably charged and should interact with charged lipids ... do you have charged lipids? did you include the counter ions to neutralize the system? Martini polarizable water model should be better in modeling systems where the electrostatic is important, it give a better mimic of the relative dielectric properties of the aqueous and membrane environments ... You should think about these problems and again 20 ns is nothing for these type of phenomena.

dujiangfeng wrote:

xavier wrote: Note that B-barrel proteins are often trans-membrane proteins, is this the case of yours?


The protein is supposed to bind to membrane but not totally go into membrane since it is just one part of the huge coagulant factor VIII.

xavier wrote: You need to rebuild your topology checking that the atomistic structure you give is correct (no atoms are missing) and the RMSD should not go higher than 0.2 nm. If the problem persist get back to us with the atomistic pdb and CG topology. But please try to fix the problem yourself first. Notably edit the CG topology and check that all the parameters make sense.

XAvier.


The new x-ray structure is good for cg model, I suppose.


Thank you,

Jiang.

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 7 months ago #1123 by dujiangfeng
Replied by dujiangfeng on topic md simulation code

xavier wrote: 20 ns is short ... not that if your protein is binding the membrane it is probably charged and should interact with charged lipids ... do you have charged lipids? did you include the counter ions to neutralize the system? Martini polarizable water model should be better in modeling systems where the electrostatic is important, it give a better mimic of the relative dielectric properties of the aqueous and membrane environments ... You should think about these problems and again 20 ns is nothing for these type of phenomena.


I got 120 ns simulation (shifted, polarized water) so far. The protein went to membrane at 50 ns (mostly randomly) and then went away again.

I also got 40 ns simulation with PME, polarized water so far (it is still running). The protein didn't move anymore. But when I checked by g_mindist, the distance between protein and membrane decreased (a little but consecutively).

In my system, the membrane is negative charged, while the protein has positive charges. However, at this two trails, I didn't add counter ions to neutralize the system (I found the counter ions sometimes prevent the binding ability of the protein to membrane (maybe this is a wrong conclusion)).

What I am going to do are:
1. The simulation with PME should be continued.

2. Make another neutralized system.

Thank you for your help,

Jiang.

Please Log in or Create an account to join the conversation.

More
11 years 7 months ago #1125 by floris
Replied by floris on topic md simulation code
Hello Duijangfeng,

it is always important to add ions.
This can have a big influence on your results, as you probably can imagine. I also would recommend, as you said, to increase the length of your simulations, as many interactions take much more time.

Floris

Please Log in or Create an account to join the conversation.

More
11 years 7 months ago #1126 by xavier
Replied by xavier on topic md simulation code
As a matter of fact the PME equations are wrong for a non-neutral system. So the results should not be trusted. Then 20/120 ns is not relevant. The phenomena you are looking at are on the microsecond time scale so that what you should do!

But before you are investing a significant amount of computer time, make sure to generate a good topology. Notably find out if the "slow" motion of the protein is realistic ... it might be but it should be moving around ...

Please Log in or Create an account to join the conversation.

More
11 years 7 months ago #1128 by djurre
Replied by djurre on topic md simulation code

xavier wrote: As a matter of fact the PME equations are wrong for a non-neutral system. So the results should not be trusted.


I don't think that is completely correct. The way they are implemented in Gromacs there is a homogeneous background charge added when simulating net-charged system with PME. See for example this discussion on the gromacs mailing list: lists.gromacs.org/pipermail/gmx-users/20...eptember/017151.html

So I think in your case it might not be a huge problem. But it never hurts to add ions (maybe not just counter ions, but a physiological ions concentration?).

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 7 months ago #1129 by dujiangfeng
Replied by dujiangfeng on topic md simulation code
Hi Floris and XAvier,

The neutralized system is ready to run now. Here I encountered another question: should I use the minimum counter ions (15 NA+ is enough) to neutralize the system or neutralize it in the physiological concentration (the ions are too much (few hundreds of NA+ and CL-) in the latter situation)?

When I run the system by PME, the polar waters were linked by a list of horizontal and vertical lines (matrix) visualized by VMD. I don't know those lines indicate the coulomb forces or errors.

XAvier, you mentioned to have a good topology. What do you mean "topology" here? Is it that the position of protein, PW, lipids, ions in the initial system? or the *itp and *top files? Actually, I also think this is a very important step for the simulation, however I don't know how to check the quality of them.

Thank you,

Jiang.

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 7 months ago #1130 by dujiangfeng
Replied by dujiangfeng on topic md simulation code

djurre wrote: So I think in your case it might not be a huge problem. But it never hurts to add ions (maybe not just counter ions, but a physiological ions concentration?).


Hi Djurre,
Yes, I also want to add ions in physiological concentration (0.15mol/L), then I have to add 65 NA+ and 50 CL- around my protein.

At this time, I only added 15 NA+ in the system. The simulation made the polar waters linked by horizontal and vertical lines which are not observed in the system without ions. :(.

Jiang.

Please Log in or Create an account to join the conversation.

More
11 years 7 months ago #1131 by xavier
Replied by xavier on topic md simulation code

djurre wrote:

xavier wrote: As a matter of fact the PME equations are wrong for a non-neutral system. So the results should not be trusted.


I don't think that is completely correct. The way they are implemented in Gromacs there is a homogeneous background charge added when simulating net-charged system with PME. See for example this discussion on the gromacs mailing list: lists.gromacs.org/pipermail/gmx-users/20...eptember/017151.html

So I think in your case it might not be a huge problem. But it never hurts to add ions (maybe not just counter ions, but a physiological ions concentration?).


mmmmmm ... Djurre might be right but I can not be totally sure. I sure should look into this.

Please Log in or Create an account to join the conversation.

More
11 years 7 months ago #1132 by xavier
Replied by xavier on topic md simulation code

dujiangfeng wrote: Hi Floris and XAvier,

The neutralized system is ready to run now. Here I encountered another question: should I use the minimum counter ions (15 NA+ is enough) to neutralize the system or neutralize it in the physiological concentration (the ions are too much (few hundreds of NA+ and CL-) in the latter situation)?

This is up to you. These are two different issues. One has to do with the fact that a charged system is not realistic (unless specific cases) and might generate artifacts in your simulations (dipoles). The second is that at physiological conditions you have a certain amount of ions in solution and they will be very important in most processes as it has been shown experiemtnaly. This is especially important for systems where electrostatic interactions are important.

dujiangfeng wrote: When I run the system by PME, the polar waters were linked by a list of horizontal and vertical lines (matrix) visualized by VMD. I don't know those lines indicate the coulomb forces or errors.

This fine. It is visualization issue. When only a part of a water molecule goes across the periodic boundary conditions VMD keeps them linked.

dujiangfeng wrote: XAvier, you mentioned to have a good topology. What do you mean "topology" here? Is it that the position of protein, PW, lipids, ions in the initial system? or the *itp and *top files? Actually, I also think this is a very important step for the simulation, however I don't know how to check the quality of them.

Here by topology I mean the itp file describing the protein. You have to make sure the pdb file you use is complete ... martinize will then generate a "good" toplogy :))

dujiangfeng wrote: Thank you,

Jiang.

Please Log in or Create an account to join the conversation.

  • dujiangfeng
  • dujiangfeng's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 7 months ago #1133 by dujiangfeng
Replied by dujiangfeng on topic md simulation code
1 My neutralized protein-membrane system with polarized water (using PME) has been running 240 ns since now. The protein is dynamic but no motion (mindist between protein and membrane didn't change, but it should decrease), even no any ion particle can move more than 0.3 nm (NA+ should go to membrane also). I didn't restrict any thing during the simulation. the total energy increased and reached to plateau, the same as potential energy and enthalpy. The LJ energy is stable always while coulomb energy kept increasing weakly. I don't have experience for running a long simulation, so I don't in this case, Should I think about another simulation code or extend the simulation time. However, in my view, I don't see the point to continue it. What is your opinion?
I posted some of my simulation code here, see if they are all reasonable:
rlist=1.2; coulombtype=PME; vdwtype=shift; rcoulomb=1.2; rvdw=1.2; fourierspacing=0.12;
pme_order = 4; ewald_rtol = 1e-5; optimize_fft = yes; nstlist = 10.


By the way, in normal martini water model system, the protein went to membrane directly and quickly (300 ps) when using PME.


2 My another neutralized protein-membrane system with polarized water (using SHIFT cutoff) has been running ~100 ns since now. During the simulation, the protein not only moves but also rotates. The trajectory of the protein motion seems not explainable.

Since the initial distance between protein and membrane is around 3.5 nm, there is no any interaction between them via "shift cutoff 1.2" calculation, right?


Thanks a million,

Jiang.

Please Log in or Create an account to join the conversation.

More
11 years 7 months ago #1137 by xavier
Replied by xavier on topic md simulation code
Ok, what you describe is a bit funky ...

1- the options you show are reasonable. The protein only moving little during 240 ns is a bit strange. You may have used position restrains without noticing! Check out your mdp file and make sure to use the tpr file you think you are using. It is recommended to use different "root" names for the different stages of the simulation. Such as "mini" for minimization, "equil" for equilibration, "run" for the actual run ...

Did you by any chance remove the COM motion of the protein separately from the solvent?

2- you are correct, if the distance is 3.5 nm then the protein does not seen the bilayer surface.

dujiangfeng wrote: 1 My neutralized protein-membrane system with polarized water (using PME) has been running 240 ns since now. The protein is dynamic but no motion (mindist between protein and membrane didn't change, but it should decrease), even no any ion particle can move more than 0.3 nm (NA+ should go to membrane also). I didn't restrict any thing during the simulation. the total energy increased and reached to plateau, the same as potential energy and enthalpy. The LJ energy is stable always while coulomb energy kept increasing weakly. I don't have experience for running a long simulation, so I don't in this case, Should I think about another simulation code or extend the simulation time. However, in my view, I don't see the point to continue it. What is your opinion?
I posted some of my simulation code here, see if they are all reasonable:
rlist=1.2; coulombtype=PME; vdwtype=shift; rcoulomb=1.2; rvdw=1.2; fourierspacing=0.12;
pme_order = 4; ewald_rtol = 1e-5; optimize_fft = yes; nstlist = 10.


By the way, in normal martini water model system, the protein went to membrane directly and quickly (300 ps) when using PME.


2 My another neutralized protein-membrane system with polarized water (using SHIFT cutoff) has been running ~100 ns since now. During the simulation, the protein not only moves but also rotates. The trajectory of the protein motion seems not explainable.

Since the initial distance between protein and membrane is around 3.5 nm, there is no any interaction between them via "shift cutoff 1.2" calculation, right?


Thanks a million,

Jiang.

Please Log in or Create an account to join the conversation.

Time to create page: 0.129 seconds