AWH convergence for a charged ligand permeating a bilayer

[Zahraibt7]

GROMACS version: 2023.2-plumed_2.10.0_dev
GROMACS modification: No
Hello everyone,

I am applying AWH to study the permeation of a charged ligand (with -2 charge) through a DOPC bilayer. My setup is largely based on the following discussion:

https://gromacs.bioexcel.eu/t/using-awh-for-bilayer-permeability/6727

Following that thread, I tested several combinations of AWH parameters (e.g., awh1-dim1-diffusion, awh1-error-init, pull_coord1_k, and related settings). I have attached a short summary of the results from several 600 ns simulations. My main difficulty is that the ligand often becomes trapped on one side of the membrane or within the bilayer, with very few complete back-and-forth transitions. Consequently, the PMF is still changing noticeably after 600 ns, and I am unsure whether this is expected for a charged ligand or indicates that my AWH setup still needs improvement.

Before investing in much longer or multi-walker simulations, I would appreciate some guidance on whether my current approach is reasonable and which parameters you would recommend adjusting first. In particular, comments from AWH developers or experienced users would be greatly appreciated.

Thank you very much for your time. The attached files

summarizes the parameter sets and corresponding results.

Bests

Zahra

Here is mdp settings:

define = -DPOSRES_DOPC_P -DPOSRES_FC_LIPID=1000

integrator = md

dt = 0.002

nsteps = 300000000

nstxout = 5000

nstvout = 50000

nstfout = 50000

nstcalcenergy = 100

nstenergy = 1000

nstlog = 1000

nstxout-compressed = 5000

cutoff-scheme = Verlet

nstlist = 20

rlist = 1.4

vdwtype = Cut-off

rvdw = 1.4

coulombtype = PME

rcoulomb = 1.4

tcoupl = Nose-Hoover

tc_grps = SOLU MEMB SOLV

tau_t = 0.5 0.5 0.5

ref_t = 323 323 323

pcoupl = C-rescale

pcoupltype = semiisotropic

tau_p = 5.0

compressibility = 4.5e-5 4.5e-5

ref_p = 1.0 1.0

constraints = h-bonds

constraint_algorithm = LINCS

continuation = yes

gen_vel = no

nstcomm = 100

comm_mode = linear

comm_grps = SOLU_MEMB SOLV

pull = yes

pull_ncoords = 1

pull_ngroups = 2

pull_group1_name = MEMB

pull_group1_pbcatom = 9610

pull_group2_name = SOLU

pull_group2_pbcatom = 29

pull_nstxout = 100

pull_nstfout = 100

pull_coord1_groups = 1 2

pull_pbc_ref_prev_step_com = yes

pull_coord1_type = external-potential

pull_coord1_potential_provider = awh

pull_coord1_geometry = direction

pull_coord1_vec = 0 0 1

pull_coord1_dim = N N Y

pull_coord1_rate = 0

pull_coord1_start = yes

pull_coord1_k = 1000

awh = yes

awh_nstout = 50000

awh_potential = convolved

awh_nbias = 1

awh_nstsample = 100

awh_nsamples_update = 10

awh1_growth = exp-linear

awh1_target = constant

awh1_user_data = no

awh1_ndim = 1

awh1_dim1_coord_provider = pull

awh1_dim1_coord_index = 1

awh1_dim1_start = -4.0

awh1_dim1_end = 4.0

awh1_equilibrate_histogram = yes

awh1_dim1_force_constant = 75000

awh1_dim1_diffusion = 2e-3

awh1_error_init = 20

[hess]

How high an energy barrier do you expect? You are getting number higher than 100 kJ/mol. This will be very difficult to converge.

The only issue I see is that you should use geometry=cylinder instead of direction. But I don’t think that that will solve your issues.

[hess]

Another thing is that you should not use multiple comm removal groups. This could cause issues.

[Zahraibt7]

Thank you very much for your helpful comments and suggestions. I will test the AWH simulations using geometry = cylinder, as you recommended. Based on my understanding of your second comment, I will also revise the COM-motion removal settings. My current .mdp file contains:

nstcomm = 100

comm_mode = linear

comm_grps = SOLU_MEMB SOLV

and for the next simulations I will use COM-motion removal for the whole system instead.

In the meantime, while waiting for feedback, I completed another 1 μs AWH simulation using a stronger force constant and have attached the updated results. The ligand still does not move frequently back and forth across the membrane and instead remains trapped on one side for a considerable part of the simulation.

Regarding the high free-energy barrier, I was wondering whether at least part of it could be expected from the physical nature of this system. The ligand carries a −2 charge and contains a cyclic phosphate group, so I would expect both a significant electrostatic penalty inside the low-dielectric membrane and a considerable desolvation penalty when moving from bulk water into the membrane core. The trajectories also show that the ligand remains strongly hydrated, suggesting that disruption of its hydration shell may contribute substantially to the permeation free energy. However, since this is my first AWH study, I am not sure whether these effects alone could explain such a high barrier or whether the results still indicate that something in my setup should be improved. I will post the results from the revised setup as soon as they are available. Thank you again for your valuable guidance.

[hess]

I was not suggesting that the high barrier is unreasonable. I think it is expected. My point is that it is very difficult to converge sampling of a process with such a high barrier. My guess is that this is not feasible within the simulation times that can be reached.

[Zahraibt7]

Yes, I agree that sampling convergence is the main challenge for such a high-barrier process, rather than the barrier itself necessarily being unreasonable. Following your suggestions, I started a new AWH simulation with a modified setup, and after approximately 150 ns I am already observing changes in the sampling behavior.

During these tests, I found that the reaction coordinate occasionally exceeded the current AWH interval of −4.5 to 4.5 nm, reaching approximately −4.58 nm and causing the simulation to stop because the coordinate moved too far outside the defined range. I am now adjusting the AWH and cylindrical pulling parameters and will first perform shorter test simulations to identify a stable setup before starting longer production runs. In particular, I am testing the cylinder radius and PBC-reference settings, together with moderately reduced values of the AWH force constant and diffusion estimate. I will also examine whether a slightly wider reaction-coordinate interval and different initial ligand positions improve stability and sampling.

In parallel, I am considering expanding the study to a small set of related ligands with different biological activities. The current ligand has micromolar activity, so including ligands with nanomolar activity may provide a broader comparison and help assess whether the calculated permeation free energies show any relationship with biological activity.

I will update the thread once the short stability tests and subsequent longer simulations with the revised setup are completed. Thank you again for your helpful suggestions.

this is the modified .mdp

define = -DPOSRES_DOPC_P -DPOSRES_FC_LIPID=1000
integrator = md
dt = 0.002
nsteps = 25000000

nstxout = 5000
nstvout = 50000
nstfout = 50000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
nstxout-compressed = 5000

cutoff-scheme = Verlet
nstlist = 20
rlist = 1.4

vdwtype = Cut-off
rvdw = 1.4
coulombtype = PME
rcoulomb = 1.4

tcoupl = Nose-Hoover
tc_grps = SOLU MEMB SOLV
tau_t = 0.5 0.5 0.5
ref_t = 323 323 323

pcoupl = C-rescale
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0

constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
gen_vel = no

pull = yes
pull_ncoords = 1
pull_ngroups = 2

pull_group1_name = MEMB
pull_group1_pbcatom = 9610
pull_group2_name = SOLU
pull_group2_pbcatom = 29

pull_nstxout = 100
pull_nstfout = 100

pull_coord1_groups = 1 2
pull_pbc_ref_prev_step_com = no

pull_coord1_type = external-potential
pull_coord1_potential_provider = awh
pull_coord1_geometry = cylinder
pull_cylinder_r = 1.5
pull_coord1_vec = 0 0 1
pull_coord1_dim = N N Y
pull_coord1_rate = 0
pull_coord1_start = yes
pull_coord1_k = 1000

awh = yes
awh_nstout = 50000
awh_potential = convolved
awh_nbias = 1
awh_nstsample = 100
awh_nsamples_update = 10

awh1_growth = exp-linear
awh1_target = constant
awh1_user_data = no

awh1_ndim = 1
awh1_dim1_coord_provider = pull
awh1_dim1_coord_index = 1
awh1_dim1_start = -4.5
awh1_dim1_end = 4.5
awh1_equilibrate_histogram = yes
awh1_dim1_force_constant = 50000
awh1_dim1_diffusion = 5e-3
awh1_error_init = 10