Moving a ligand through a beta barrel membrane protein using steered MD

GROMACS version:2018

Hi all, I would like to move a ligand through a beta barrel shaped membrane protein using steered MD. For this purpose, I have inserted the ligand at the top of the beta barrel protein as close as possible. Then using the protein and the ligand as two pull-groups I am planning to move the ligand through the cylindrical beta barrel protein from top opening to bottom end. However, if I use typical pull settings like:

pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = LIG
pull_group2_name = betab_CA
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = Y Y Y ; pull along z
pull_coord1_groups = 1 2 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2

I see the ligand moving further away from the protein and not going inside the barrel; however if I use a negative pulling-rate as: pull_coord1_rate = -0.01, then I see the ligand moving towards the protein COM.

However, as soon as the ligand comes very close to the protein COM, the simulation is suddenly stopped with an error message: “Pull reference distance for coordinate 1 (-0.000003) needs to be non-negative”, the ligand is unable to cross the COM of the protein and stuck there. I have checked that the protein-ligand COM-COM distance is approximately 1.2 nm, they are not even close to zero, I do not understand what is the negative distance that gromacs is reporting.

Is there any possibility that I can use steered MD and the protein-ligand COM-COM distance as the collective variable and move the ligand inside the beta-barrel protein (reverse pulling if possible). Any help will be much appreciated thank you.

You are using a distance (and in 3D, not along z as in the comment). What you want is pull-geometry=direction. A distance can not be less than zero. The simulation stops with an error because of that.

Thank you for your suggestions. I have tried with pull-coord = direction and pull-coord1-vec = 0 0 1, so to push the ligand through the betabarrel protein along the axis of the protein or z-direction. I have a concern regarding that process. If I use pull-coord1= direction, will the ligand be able interact with the wall residues of the beta-barrel protein or it will just move through the protein along the z-axis without interacting with the protein residues at all.

I have already run a steered MD simulation with the parameters:

pull-coord1-geometry = direction
pull-coord1-vec = 0 0 1 (for moving the protein along z-direction)
pull-coord1-dim = N N Y (in order to get the z-component)
pull-rate = 0.001 (unit nm per picosecond)

The ligand takes approximately 6-7 ns time to move from one end to the other. Then I extract the relevant configurations for a range of reaction-coordinate values. In this case I extract structures with z-distances ranging from -2.5 nm to +2.5 nm, with a gap of 0.1 nm (approximatel covering the distance from top opening to bottom exit of the beta barrel protein). Then I run umbrella sampling simulations for each of those configurations; but during the umbrella sampling, I use pull-coord1=distance, so that the ligand can explore the conformational space a bit more;

Problem with gmx WHAM: In order to get the protein-ligand COM-COM distances with negative and positive signs, I determine the z-component of the distance vector and put negative signs before distance values for configurations with negative z-components (i.e. these configurations ligand’s COM lies above the protein COM); then I generate a tpr file with negative reference values for pull-coord1-init and pull-coord1-geometry = direction; in this way I generate tprs where the initial distances are negative; If use the actual umbrella sampling tpr and the pullx (pull coordinate output) file with negative distances, then gmx wham prints NAN for free energies; however, generating a tpr with negative starting distances work; If I use plain gmx wham, then the corresponding PMF looks utterly weird non-converging with an energy barrier increasing upto 800 kJ mol-1; So, I have tried to use a cyclic pmf, as in our case the reaction-coordinate is probably periodic or cyclic even it is a distance. (PLEASE IF YOU CAN COMMENT ON THE USE OF CYCLIC PMFs FOR THIS PURPOSE). so, during gmx wham I use a flag -cycl, that gives me cyclic pmf which looks much nicer than acyclic or non-periodic pmf. However, I am still not sure if this is the right approach of calculating pmf for this particular system.

If you can please comment on the method employed for translocation of the ligand: 1) running steered MD with reaction-coordinate along the z-direction, 2) running umbrella sampling simulations with steered MD generated configurations but with reaction-coordinate = distance, for better sampling, 3) running gmx wham to get a cyclic PMF, thanks in advance.

You should never use geometry=distance.
With geometry=direction, the ligand is not at all restricted by the pull code (but of course hindered by interactions with the protein) along x and y.

Thank you Professor Hess for your suggestion. I will then use a much slower pulling rate along with pull-coord1-geometry = direction and pull-coord1-vec = 0 0 1; BTW, is there any particular reason you are referring to for not using pull-coord1-geometry=distance? Are you suggesting not to use geometry=distance for this case or in general

You should not use distance for your case. You want to pull a ligand through a protein. Here is no distance involved at all. You have the position along z as the reaction coordinate.
You can have a distance when pulling a ligand into a binding pocket or are investigating the PMF between two molecules.

Yes, now I clearly understand your point, thank you very much for your suggestions.