Apparent convergence issues - Peptide Mutations

[amin-sagar]

Hello.

I tried doing a protein-peptide FEP (#682 )
I initially tried with the standard 11 replicas/lambda states.

This seems to suggest not enough overlap between the replicas, so I increased the number to 15 and then 21 using

protocol_settings.simulation_settings.n_replicas = 21
protocol_settings.lambda_settings.lambda_windows = 21

Now, I get

The timeseries plot seems to show exchanges but the transition matrix is mostly empty.
It doesn't seem like the transitions are 0 for the initial replicas.
Why could this be happening?
I would be really grateful for any suggestions.
Best,
Amin.

[IAlibay]

@amin-sagar I converted this to a discussion topic since it seems to be more related to a usability question.

Could you provide a little bit more details on your system?

For example:

• How long are you running the simulation?
• How large is the peptide?
• What does the peptide mapping look like (i.e. how many atoms are being perturbed)?
• Does the peptide undergo any large rearrangements during the simulation?

[amin-sagar]

Thanks @IAlibay
All the other values are default. 1ns Equilibration, 5ns Production

 'simulation_settings': {'equilibration_length': 1.0 <Unit('nanosecond')>,
  'production_length': 5.0 <Unit('nanosecond')>,
  'minimization_steps': 5000,
  'time_per_iteration': 1 <Unit('picosecond')>,
  'real_time_analysis_interval': 250 <Unit('picosecond')>,
  'early_termination_target_error': 0.0 <Unit('kilocalorie_per_mole')>,
  'real_time_analysis_minimum_time': 500 <Unit('picosecond')>,
  'sampler_method': 'repex',
  'sams_flatness_criteria': 'logZ-flatness',
  'sams_gamma0': 1.0,
  'n_replicas': 21},

The peptide is 17 residues long.
The perturbation is a Tryptophan to Alanine mutation. So the side chain of Tryptophan except c-beta disappears.
The following are the RMSD plots

I couldn't find a way to visualize the trajectories. Are there some instructions that I missed?

Aren't the matrix plots and timeseries two representations of the same data. i.e. if I see transitions in the timeseries, shouldn't I see them in the transition matrix as well? Or are they plotting something different?
Thanks again.
Amin.

[amin-sagar]

Update - Tried 31 replicas and lambda states.
This leads to nan error for 2 of the simulations and I get some results in one of the simulations.
The results look like this.

[RiesBen]

Hej @amin-sagar,
I see you have issues with the overlap in your Amino Acid changes. It is really interesting to see, that you seem to be able to generate a bit of overlap with adding the lambda windows! However did you consider going over intermediate Amino Acids? :)
So you could use intermediate Amino Acids like PHE/HIS/VAL in order to generate the overlap in your transformations:

Let me know what you think :)

P.s.: forgot the TYR->PHE->VAL->ALA Path in the figure

[RiesBen]

by the way, how do you build the Mappings? and are you able to share them, just for checking? :)

[hannahbaumann]

@amin-sagar Thank you for reporting about your experience, this is very interesting! Regarding your question about the difference between the replica exchange transition matrix and timeseries plots: You're right that those contain similar information. In the plots you showed where you used 21 lambda windows, I think that the transition probability of zero between e.g. lambda 0 and 1, is due to rounding. In the timeseries you can see that there are a few transitions (maybe ~10), but since every 1 ps there is a replica exchange attempt, 10 transitions over 5 ns (total of 5000 exchange attempts) is very little, so the transition probability is close to zero.

[amin-sagar]

@RiesBen @hannahbaumann
Thanks for the suggestions.
I am generating the mappings using the example script as follows

complexed_transformations = []
solvated_transformations = []

for (ligA_name, ligB_name) in connections:
    ligA = complexed[ligA_name]['ligand']
    ligB = complexed[ligB_name]['ligand']
    
    mapping = next(mapper.suggest_mappings(ligA, ligB))
    
    complexed_transformations.append(
        openfe.Transformation(stateA=complexed[ligA_name], 
                              stateB=complexed[ligB_name], 
                              mapping={'ligand': mapping},
                              protocol=protocol) 
    )
    solvated_transformations.append(
        openfe.Transformation(stateA=solvated[ligA_name], 
                              stateB=solvated[ligB_name], 
                              mapping={'ligand': mapping},
                              protocol=protocol) 
    )

Could other mappings work better?

@RiesBen I tried just the last leg of the your path i.e. Val to Ala. Here are the results
Default - 11 replicas

21 replicas/lambdas

31 replicas/lambdas

For V2A, I get ddg estimates for all three attempts for 11, 21 and 31 replicas. It's quite surprising that I get strongly negative ddg considering this Valine fits nicely into a hydrophobic pocket and the mutation to Alanine is experimentally known to severely reduce the affinity.

I am working on creating a system that I can share for reproducing the results.

Are there any other things I should try?

In general, when residue mutations are done by perses or FEP+, are the "intermediate" amino acids generated internally to help with convergence?

Thanks a lot for your suggestions.
Amin.

[RiesBen]

@amin-sagar

1. ddG sign: hm ... ok, if that is your expectation. Do you have any alternative results backing that?
2. The intermediates are not required, but they can be a way to make the transformation easier if the pertubation is to large. (there are examples of commercial software using this: https://www.cresset-group.com/about/news/flare-fep-speed/). For the restraint question I would link to @IAlibay. But I'm not sure why you think those restraints would help, unless you are using small peptides? Can you give us here more details on your system? otherwise it will be very difficult to understand your problem.
3. see here: [docs] we should point out some example(s) on how to read out the trajectories from our output files  #916

[amin-sagar]

Thanks @RiesBen

1. There is SPR data showing that this mutation leads to a non binding peptide.

2. I am using small peptides (13-20 residues). Sorry, I should have made it clearer. The problem is described in slightly more detailed manner in Protein-Protein or Protein-Peptide alchemistry? #682 . To summarize,

   • I am working with protein-peptide complexes where the peptide can be quite flexible.
   • The peptides can have non-standard amino acids, hence I want to use espaloma as it has been shown to perform well for both small molecules and proteins.
   • There is always a point mutation in a peptide and I have associated experimentally determined affinities (using SPR).
   • My aim to to use OpenFE to see if I get some agreement with the experimentally determined affinities.

3. Thanks. The scripts work except it seems I have to replace state with state_id.

4. I will prepare a dataset that I can share for reproducing the results.

Best,
Amin.

[RiesBen]

@amin-sagar
Cool :)

I think one very important step would be to check your atom mappings. You want to make sure, that as many atoms are mapped as possible. (goes in the similar direction as the restraints) so for the ALA-VAL mutation I expect that for the Peptide (or protein) all heavy atoms are mapped except 2 ( 2 atoms are added by the side chain modifications) or less. For TRP -> ALA this would lead to maximally 8 missing atoms in the mapping if I counted correctly.

In hybrid topology approaches like our RBFE - Protocol uses you have something like "natural" restraints, as most of the coordinates should be represented as one. :)

[amin-sagar]

@RiesBen
The atom mapping seems to be correct. The resolution of the saved file with lomap_mapping.draw_to_file is too low for large molecules like peptides, but I think I can see the following for the reverse modification i.e. Ala to Val. The blue circles are just barely visible on the valine but I think they are there.

I did the reverse modification to see if I get positive dg and I do.

But it could be just one of the cases where the prediction doesn't agree with experiment. If the protocol seems correct, then I will do a larger set to see how often the predictions agree with the experimental results.
For W2A, I see this. Apologies for the low resolution.

The alanine in this transformation has a red circle on one hydrogen.
In addition, there are a bunch of smaller red spots. Is some documentation about the meaning of different kinds of circles?

[RiesBen]

@amin-sagar
this sounds like a cool plan :)

for the coloring check this page out: https://docs.openfree.energy/en/stable/guide/setup/creating_atom_mappings_and_scores.html#visualising-mappings

bottom line - red = vanishing/appearing atom; blue = element changing atom

[amin-sagar]

Hello Everyone.
I am continuing with my attempts to use OpenFE for protein-peptide FEP.
There is an interesting recent paper which has a protein-peptide system which can be used as a test system.
https://www.biorxiv.org/content/10.1101/2024.04.22.590325v1.full
6sba is a nice example of a protein-peptide system.
Supplementary Table S3 has the experimentally determined binding affinities

The structure files are available at https://github.com/schrodinger/protein-fep-benchmark
I have started with 4 Ala mutations.
S8A
V9A
T20A
W21A

T20A is noted to be a hard case as it is not at the interface and probably changes affinity due to some major rearrangement of the peptide.
In any case, the results till now are as follows

Clearly, the simulations are not converged. I will try increasing the number of replicas and the length of simulations.
Again, any suggestion are most welcome and I will test them and report the results.

Best,
Amin.