Problematic SAMS simulation

[sgill2]

This is a followup to the meeting we had recently. During the meeting we had discussed one SAMS example in particular which @jchodera and @Lnaden noted had some equilibration/convergence/sampling issues, . Here's a picture from the notebook showing the lambda state sampling and the equilibration detection.

Visually inspecting the trajectories I can't find see anything particularly notable in the solvent simulation, but the complex simulation seems to have a change of binding mode occur when a ring in the complex flips. Where it occurs in the simulation corresponds roughly to where the change in the lambda switching occurs.

starting pose

new pose

What would be the best way to treat cases like this?

Here's the link to the Yank simulation: https://drive.google.com/file/d/1_nI8QofsnM59yNdIiDB4Aj5hhEx9Cf3T/view?usp=sharing

[jchodera]

Thanks for sharing the trajectory data!

I've taken a look at the problematic solvent phase, and noticed a few very slow timescale phenomena:

• The single Na+ counterion becomes associated with the carboxylate moiety for extended periods of time. If we can increase the salt concentration to better mimic assay conditions, this would reduce correlation times by providing additional charge screening and more ions to compete with the alchemically decoupled ion. You can do this with the ionic_strength option
• The tetrahydroquinoline (THQ) scaffold appears to be able to undergo syn/anti transitions at the aminobenzene substituted position that are likely responsible for some of the slow timescale behavior we see in the correlation time trace. We could potentially speed this up by also softening some ligand torsions by reducing the lambda_torsions in the protocols: section. This may also require passing alchemical_torsions=True to AlchemicalRegion, though @andrrizzi will have to document how to do this in the YANK docs.
• The THQ ring nitrogen also seems to have a slow nitrogen inversion timescale. Speeding this up might require softening angles as well, following a similar protocol to torsions above, but we'd have to be careful this doesn't introduce too much flexibility.

Some snapshots are shown below to illustrate these points.

[jchodera]

Visually inspecting the trajectories I can't find see anything particularly notable in the solvent simulation, but the complex simulation seems to have a change of binding mode occur when a ring in the complex flips. Where it occurs in the simulation corresponds roughly to where the change in the lambda switching occurs.

I think these kind of "change in binding mode" issues will probably be best addressed by experimenting with a combination of techniques:

• Using the RMSD restraint we're waiting on (align via complex, compute RMSD of ligand) to restrict the ligand to ~2 A heavy atom RMSD from the defined pose to evaluate the free energy of binding for the specified pose
• Using the existing RMSD restraint at intermediate lambda values to try to "pull" the ligand back into the bound pose at weak or noninteracting lambda values, in an attempt to reduce correlation times. This will be most effective when we have the ability to add a second restraint to release the noninteracting ligand into a Boresch or Harmonic restraint.
• Alchemically softening residues lining the binding site at intermediate lambda values once we have the ability to add a second restraint
• Experimenting with alternative forms of the periodic Boresch restraint where we specify which atoms in the scaffold and protein are used in the restraint

[sgill2]

I tried running with the same settings as above, except with a 50mM salt concentration using both the replica exchange and SAMS samplers. Compared to using only neutralizing salt in the replica exchange case, the correlation was much better between the predicted and experimental free energies (for both charged and neutral ligands). However, there still is something strange going in the SAMS case for the complex, which you can see from the image below. The left graph is a previous SAMS run which only included neutralizing salts, and the right is from the SAMS simulation with a 50mM ionic strength using the ionic_strength option. I'm only showing this pair for comparison, but this behavior was present for all the ligands I tested.

I'm not sure what is the source of the spurious energy jumps for the righthand side. Looking at the trajectory, there doesn't seem to be anything obvious to me that would be causing this, except that the binding mode of the ligand tends to differ from the starting pose at the point at those energy jumps. Maybe this is somehow related to the restraints (PeriodicTorsionBoresch)?

I would appreciate any insight you might have @jchodera @Lnaden Here's a link to the aforementioned Yank simulation: https://drive.google.com/file/d/1218fXhgX1NTaq5VLnadDidfkow6DCdeE/view?usp=sharing

[jchodera]

@sgill2 : The jumps in the top plot look rather concerning! This corresponds to such a huge change in likelihood that I am concerned there is some sort of major bug here.

@Lnaden : Can you help investigate in the morning?

[jchodera]

I've looked at the trajectory you provided, and I don't see anything out of the ordinary. There's some occasional transient ion pairing that occurs, but I'd be surprised if it leads to this large of a change in energy. I'm rendering a movie right now that I'll post soon.

@Lnaden : We may have to see if we can identify the frames where we see these jumps and try to figure out what is special about them. I don't have any clever ideas right now other than to try to get an energy breakdown by force component, but it may be easier to do that by running a new simulation with a version of YANK modified to store this.

@sgill2 : Can you post the input files you used here?

[jchodera]

@sgill2 : Can you post a ZIP archive of the input files you used here?