Dear community
I have a question about umbrella sampling. The question is not all about groamcs commands but rather about the method. I have amyloid in a membrane (I relaxed it before pulling it out), and my goal is to evaluate the free energy profile during amyloid extraction from the membrane. And this is a Martini coarse-grained system.
To obtain umbrella windows for sampling, I used SMD simulation.
Here are some of the parameters of the pulling simulation.
pull = yes
pull_ncoords = 1
pull_ngroups = 2
pull_group1_name = protein
pull_group2_name = membrane
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = N N Y
pull_coord1_groups = 1 2
pull_coord1_start = yes
pull_coord1_rate = 0.00005
pull_coord1_k = 200
The pulling takes ten ns. Attached is the pull_pullf.png, which shows the pull force profile.
Com_components shows the center of mass distance evaluation and the components of that vector (all calculations are done using the vmd tcl command window).
So, you can see that the protein drifts in the x and y directions instead of steadily moving out. My first question is:
- is that acceptable during free energy evaluation?
In other words, can we have two windows close to each other in the z-direction (which is the reaction coordinate) but not so close in other coordinates?
My second question is how to evaluate if the pulling parameters are correct.
So, for example, if I change pull_coord1_k to 1000 instead of 200, I will get the e following (see attached pul_pullf_1000.png)
After I created windows and started sampling, I obtained the following histogram
.
The FE profile is the
.
I think something is wrong with this profile because there is no energy barrier to the membrane. So, before increasing the sampling, I wondered if I performed the SMD correctly.
Best,
Vardan.
Dear @base_editor ,
I think it should be right to have the molecule move on the xy plane. If you think about it, you are pulling along z, that is, the trans-membrane axis. As such, you are implicitly assuming that the free energy profile is independent from the ‘position’ of the z axis, i.e., along with (x,y) specifically the axis goes through. Basically you are saying that the membrane is homogeneous/symmetrical along the xy plane for a given z, otherwise you can’t really obtain a profile that is true for the whole membrane. In principle, is there any difference if the molecule goes through a given (x,y), or through a position that is two nm nearby? Also, keep in mind that when you restrain a molecule along z for sampling a given window, then the molecule will necessarily diffuse along the xy plane, especially when the window is in the solution phase of your system. I would argue that actually this is something you do want, as this will better describe the window.
The force you use to pull the molecule and generate the starting configurations needs to be high enough for the molecule to go through and low enough that you are not actually making a hole/destroying the bilayer. I think there is no general value, as a rule of thumb I would check the pulling trajectory an see if how the bilayer is reacting to my pulling. For example, I was pulling water molecules through a bilayer and I remember a trajectory where a molecule stuck to the hydrophilic group of a cholesterol and basically pulled the lipid inside the membrane. This is clearly nonsensical, and I discarded the trajectory and started a new one for window generation. So, a good starting point for me is to take a look at the trajectories and see if anything weird is happening and/or you are damaging the bilayer structure.
Personally, the windows seem decently spaced. There is maybe an empty space around 4.25 nm which may be worth filling up. Regarding the free energy profile, I have no idea about what are good values to expect for your system. However, for sure it is impossible to see if the free energy has converged just by looking at a single free energy profile plot. For how long did you simulate each window? Did you combine the windows with WHAM? Also, what do you mean no energy barrier? You have a 200 kj/mol drob, which is enormous. Where is the membrane with respect to the reaction coordinate? Is it centered in zero?
Dear obZehn
Thanks for your reply. I understand the logic behind diffusion in the xy plane, and I agree.
About the force. I used force restraints on the phospholipid groups of the lipids in z directions. The reason is that I did not want the membrane to deform. And I know I later have to prove that those restraints don’t affect sampling. I plan to run the same simulations without restraints to see the effect.
The membrane is centered in the simulation box, and the reaction coordinate is the center of mass distance between the membrane and the protein. The membrane’s upper leaflet is around 2nm.
No barrier means that the free energy profile suggests that the protein will enter the membrane without difficulties, which is hard to believe.
Each window takes 50 ns, and I used WHAM afterward. I mainly took from the tutorials of the gromacs, based on Prof. Lemkul’s article “Assessing the Stability of Alzheimer’s Amyloid Protofibrils Using Molecular Dynamics” (and tutorials are also created by him or his group).
best,
base_editor.
