GROMACS version: 2022.2
GROMACS modification: Yes/No
Here post your question
Hi,
I would like to get a PMF of drug translocation over the bilayer. To do so I have performed Umbrella sampling simulations of a drug over the entire bilayer. First I did steered MD to collect the initial coordinates of my drug in the membrane, and afterward performed MD for all US windows. I have followed the procedure described in several papers and set a reaction coordinate as COM of drug and membrane. When I construct my PMF for each leaflet with wham method it looks fine. However, when I try to reconstruct that over the whole bilayer (both leaflets), it is messed up. I suspect that the reaction coordinate is set as a distance between COM of my groups, and when I try to join both leaflets that distance overlaps. Can anyone suggest how I could create one profile considering both membrane sides? Is it actually a problem with defining reaction coordinate or maybe the assembly itself in WHAM? Any hints are appreciated.
Here is my one-sided PMF.
My pull code:
;Pull code
pull = yes
pull_ncoords = 1
pull_ngroups = 2
pull_group2_name = drug
pull_group1_name = POPC
pull_coord1_type = umbrella
pull_coord1_geometry = distance
pull_coord1_dim = N N Y
pull-pbc-ref-prev-step-com=yes
pull_coord1_groups = 1 2
pull_coord1_start = yes
pull_coord1_rate = 0.0
pull_coord1_k = 1500
Is the membrane symmetric? If so, then it makes no physical difference on which side of the membrane your drug is, so you can just duplicate the raw PMF and add a minus sign to one copy’s X coordinate.
If you want to e.g. assess convergence (which is always a good idea), you can run WHAM twice, once per every half of the dataset, and again reverse the sign in one of them.
For future free energy calculations, you might want to use direction instead of distance, and supply a vector pointing in the Z direction (0 0 1).
Thanks a lot. I tried distance but the particle moved too fast. Now I adjusted the coord rate and it looks fine. The membrane is slightly asymmetric, so I will try to connect PMF of each leaflet.
One more thing, is it fine to use my previous pull files (pullx and pullf from distance geometry) and modify only the .tpr files (with direction) for the wham analysis of both leaflets? Maybe I should additionally -rerun the trajectories with new tpr? Or rather these calculations cannot be reliable anymore. Thanks.
If the membrane is asymmetric, then you can still do this in one shot - as you say, rerun the trajectories with a .tpr that will give you signed values, and perform reweighting according to the original US weights.
(In fact, reweighting is so general that you can perform it with any new collective variable/CV as long as it was sufficiently sampled in the original set of simulations.)
If you ran your US using the distance CV \xi with an external potential in i -th window V_i(\xi) = \frac12k_i (\xi - \xi_i^0)^2 and obtained an original PMF F(\xi), then every frame you got has an (unnormalized) statistical weight w = \exp(\beta(V(\xi) - F(\xi))). What you can then do is calculate this per-frame weight, bin the normalized weights along your new CV, and voila - you have the probability density along the new CV. (This procedure can be made binless with techniques like KDE, but that’s a technicality.)
Then you do standard Boltzmann inversion, i.e. F(\xi') = -k_BT \ln(p(\xi')), and you have your new corrected PMF. That’s the fully formal version, and it probably takes 10-15 lines of Python code to implement. Splitting the inputs based on the rerun will be an approximation, perhaps good enough if you don’t need a very precise description of the barrier region.
Milosz, thanks a lot for your time and help. You’re right, I’ll write a script to comprehensively reweight all windows.
You’re welcome! There’s a part of it that always feels like reinventing the wheel but then you get to learn the most from doing it on your own.
