Deriving initial guess for BARRIER

[marcel...@gmail.com]

Hi everyone,

I am running into an issue to guess an reasonable initial BARRIER values for my OPES_METAD protocol.  My system A is composed of a small ligand bound to a rather amorphous surface, while my system B is composed by a modified ligand bound to the same surface. My CV is a simple distance between the COM of the ligand and the COM of the surface.

I've tried a regular wtMETAD approach but I had to increase HEIGHT to 5kj/mol and BIASFACTOR to 20 to see a single unbinding event (most of my previous METAD approaches used smaller values). The ligand binds back to the surface and remains there for a very long time.

I've tried OPES_METAD_EXPLORE with BARRIER=25 to see if it improves my sampling and, while both ligands do unbind a few times, they often get stuck to the surface for very very long times.

So maybe my initial guess for BARRIER is wrong? Is there any trick or calculation I could use to make an informed guess for my BARRIER value?

Best,

--

Marcelo D. Polêto, Ph.D.

Research Associate

Department of Biotechnology

Universidade de São Paulo - Escola de Engenharia de Lorena (USP-EEL)

Estrada Municipal do Campinho, 100 - Campinho

Lorena, SP 12.602-810

He/him/his

[alx.th...@gmail.com]

Hi Marcelo,

I would say that your BARRIER might indeed be too low. Two reasons: first, your ligands might indeed have an affinity of more than 25 kJ/mol, and second, your force field and parameterization might overestimate binding affinity (even if you'd expect it to be ~25 kJ/mol).

Second idea is that your CV seems to be suboptimal to study just a ligand-surface aggregation. If we consider surface to be roughly the same, and binding events to different parts of the surface to represent the same macroscopic event, then we need a CV that is invariant to where on the surface you bind. Otherwise a binding event near the center (lower COM-COM distance) is one event, and a binding event more sideways (larger COM-COM distance) are entirely different events for your CV, you are sampling a new region of the space and you wait until it is filled with potential enough to simulate unbinding once again. In principle, with long enough sampling, you might get into a diffusive regime, but then you will face the challenge of how to interpret such a free energy profile. So, I would suggest here to think of a CV that better captures the essence of the process you are studying. Might want to play with projections on a specific axis, or minimal distances to groups of atoms, etc.

To your last question about a trick, if you have a good CV, I would suggest to run the most ordinary metadynamics (not even well-tempered). It will give you very crude and rough estimates, but fast. And because of how it works, you will almost always get this estimates to be higher that "true" values which is also good for OPES.

Hope it helps,

Alex

[marcel...@gmail.com]

Hi Alex,


Thanks for your suggestions. I will try the ordinary metadynamics and see what I get from it.

About the CV: from what I described, you would be totally right. I forgot to mention before, but my 'receptor' is a spherical nanoparticle made up by a few long polymers, so my simple distance is basically radial in relation to the receptor COM, so my understanding is that I'm not facing the issue you mentioned.

These challenging system uses a polarizable FF, and their non-polarizable counterpart worked just fine with wtMETAD (HEIGHT=1.0, BIASFACTOR=15). I got a roughly 15-20kJ/mol barrier from the non-polarizable metadynamics, so I assumed a BARRIER of 25 kJ/mol (even 50 kJ/mol) for the polarizable would do the trick, but the ligand unbinds once and remains stuck for long times once it binds back.

For reference and record purposes, this is what I am using:

#######################################
# Defining CVs

dmin: DISTANCES GROUPA=ligand GROUPB=center_com LOWEST
dmax: DISTANCES GROUPA=ligand GROUPB=center_com HIGHEST

#######################################
# Defining UPPER and LOWER walls

lwall: LOWER_WALLS ARG=dmin.lowest  AT=47 KAPPA=2000.0 EXP=2 EPS=1 OFFSET=0
uwall: UPPER_WALLS ARG=dmax.highest AT=85 KAPPA=2000.0 EXP=2 EPS=1 OFFSET=0

opes: OPES_METAD ...
   FILE=Kernels.data
   TEMP=298
   ARG=dmin.lowest
   PACE=1000
   BARRIER=40
   #SIGMA=0.5
   SIGMA_MIN=0.01
   STATE_RFILE=restart.data
   STATE_WFILE=state.data
   STATE_WSTRIDE=1000*10
   STORE_STATES
   NLIST
...

[alx.th...@gmail.com]

Hi Marcelo,

thank you for the clarification!

With this context, why are you not just using 

CV1: COM ATOMS=..atoms of the nanoparticle...

CV2: COM ATOMS=..atoms of the ligand...

CV: DISTANCE ATOMS=CV1,CV2

lwall: LOWER_WALLS ARG=CV ...

uwall: UPPER_WALLS ARG=CV ...

OPES_METAD ARG=CV ... ?

The reason for it is to make metad/opes always consistent in how forces are generated. In how you wrote it above, in each step, a CV is just one distance, and the potential generates forces applied only to atoms that define this one distance. Meaning, that the whole ligand should then adapt to only one atom being pushed.

If you use COM and COM, then the forces are distributed along all atoms that participate in the construction of such a CV. So each step the whole ligand is being pushed uniformly, and it is only the water that needs to adapt. Such a setup is less prone to result in some weird conformations.

Also, we might imagine that there exist multiple bound conformations that differ in how close ligand atoms can get to the COM of the nanoparticle. Then these are two minima on your 1D CV that need to be filled independently, and this also damages sampling efficiency. If you are primarily interested in "bound/unbound" level of resolution, then COM-COM distance would work better both as a driving CV and as a CV to do FES construction afterwards. You can then add a second dimension to reweight upon if you are interested in nuances.

One thing also to check is if your "stuck" situation correponds to times when you reached the minimal sigma. OPES tends to do it, and it is frustrating - it "thinks" it has sampled evething and it is time to decrease sigma, but in reality the system is far from diffusive behavior.

Best,

Alex