General questions about RBFE

[AminaMnhr]

Hi everyone,

I am a beginner in using RBFE methods to compute free energies for series of congeneric ligands. I am a researcher and we would like to use RBFE for different kind of proteins. In most cases, we would have one crystal structure of one ligand, and the other ligands would be analogs that I would have to use to compute their binding free energies.

Question 1: To build a network of my ligands, should I do a radial network with one single ligand as the reference? Or should I do a cyclic network for example?

Question 2: I heard LOMAP is a popular tool for building the network. Do you know if then I have to make an environment and install openfe for using LOMAP? I just installed lomap2 with conda but it seems that to use the lomap tool as in the showcase (https://colab.research.google.com/github/OpenFreeEnergy/ExampleNotebooks/blob/main/showcase/openfe_showcase.ipynb#scrollTo=f5506f1a ), i need rather to install openfe package?

Also I might be wrong but is LOMAP the name of one kind of network for my ligands? Or does LOMAP can also do different networks?

Question 3: It seems that Konnector can also be used to generating the network between ligands. Would you suggest me to use the CyclicNetwork of Konnector or just the LOMAP tool as in the showcase?

Question 4: when I generate my transformation network with LOMAP or Konnector, how do I know afterwards which edges of ligands I have to compute to do my RBFE? Do I need to visualize the network and find myself which edges to do or is there a tool that gives me directly the edges? I can see that I can visualize the network created by LOMAP or Konnector, but I am still confused on what edges I have to compute for my series of ligands.

Question 5: Once I have computed the edges for my ligands and get the RBFE (I am not using openfe for this step but another method to compute RBFE, I would like to compare openfe and another method), so I have the ddG for each edge, how do I get the calculated dG value for each ligand? So to do the conversion ddG to dG ?

For some proteins, I might have more than one experimental value dG for the ligands in the serie. But for others proteins and in most cases I think, I might just have ONE experimental value dG for one ligand.

I heard of cinnabar but I am also confused on how I need to structure my data of ddG and to use the tools.

Question 6 : Refers to question 5. Does the ddG to dG conversion depends on my network of ligands? I do not know if I need somewhat to fit cinnabar with lomap/konnector or something...

Question 7 : In my RBFE calculations, there is the cycle closure error that I have to take in consideration to correct my ddG values. I believe this cycle closure depends on my network. Is there an algorithm to compute the cycle closure correction itself to my ddG values? And do you have to correct the ddG values before applying the tool (cinnabar?) to go from ddG to dG?

Thanks so much for your help!!

Amina

[IAlibay]

Question 1: To build a network of my ligands, should I do a radial network with one single ligand as the reference? Or should I do a cyclic network for example?

This is unfortunately a research question that no one really has an answer to right now. It depends on transformation difficulty and how much redundancy you want. If you're trying to be thorough, then having a network with a certain amount of redundancy is advised.

[IAlibay]

Question 2: I heard LOMAP is a popular tool for building the network. Do you know if then I have to make an environment and install openfe for using LOMAP?

It depends on what you want to achieve and how. You can use lomap on its own and it will also install gufe which means you'll be able to use the newer API. However, if you want to use it with the rest of the openfe toolkit (as demonstrated in the showcase), you'll need to install openfe too.

If you can maybe expand on your needs here we might be able to help provide a better answer.

[IAlibay]

Question 3: It seems that Konnector can also be used to generating the network between ligands. Would you suggest me to use the CyclicNetwork of Konnector or just the LOMAP tool as in the showcase?

This is hard to know, it will depend on each system. Right now, we are advocating for LOMAP networks as a default when you need redundancy (as part of our industry benchmarks), but this will change over time.

[IAlibay]

Question 4: when I generate my transformation network with LOMAP or Konnector, how do I know afterwards which edges of ligands I have to compute to do my RBFE? Do I need to visualize the network and find myself which edges to do or is there a tool that gives me directly the edges? I can see that I can visualize the network created by LOMAP or Konnector, but I am still confused on what edges I have to compute for my series of ligands.

When you create a LigandNetwork, you can then use it to create an AlchemicalNetwork for further simulations. You can iteration through the edges of the network to get each mapping. This is better explained in our documentation, for example here is a cookbook that shows how to convert from a LigandNetwork to an AlchemicalNetwork: https://docs.openfree.energy/en/stable/cookbook/create_alchemical_network.html

[AminaMnhr]

Hi IAlibay,

Thank you for your fast and useful answers!

To be more specific, I have already my own free energy protocol that compute edges of ligands. I might use openFE alchemical simulation as well for comparison later but not now. So I do not need to create an Alchemical Network, but just to compute an optimal network with my congeneric series of ligands using LOMAP, so that I can know which edges I can compute afterwards with my RBFE method. Most of the time, I have a bunch of ligands analogs to an experimental ligand, and I want to compute the ddG but I need to use LOMAP for knowing the optimal edges.

So far, I just need to use LOMAP to have a network of ligands and give me the edges I need to compute. Is there a way to just compute the LOMAP network and then have a Python dictionary or something that give me the edges for my ligands?

Once I have computed the edges with my RBFE method, I have the ddG values and I am wondering therefore how to rank my ligands (so I might need the dG values, hence my question 5/6/7 as well, if you have an idea that would be very kind!). For that, I have the calculated ddG, but I might just have one experimental dG.

Thanks.

[AminaMnhr]

Hi,

I am still wondering for LOMAP if there is a way to get the edges of a ligand network, for example as a Python dictionary or something that give me the edges for my ligands? I do not need Protocol or to create an Alchemical network as I am computing the relative ddG with another method protocol.

For converting ddG to dG using cinnabar, I managed to use it. I was wondering, your example cinnabar tutorial is using experimental data for the ligands of TYK2 (that defines the Experimental block for cinnabar.csv file). However for the future, I might just have one experimental data for my ligands and not all of them, as I compute the free energy method to find free energies of congeneric ligands we do not know yet the binding (in order to rank them afterwards).

Does cinnabar still work for computing absolute dG from relative ddG if I just have one experimental value for my congeneric serie of ligands?

[hannahbaumann]

Hi @AminaMnhr ,

Yes, you can generate the Lomap network and extract the names of the ligands in the edges. Would you also want to have the mapping or are you getting that elsewhere?
You can generate a Lomap network like in this example:
https://github.com/OpenFreeEnergy/ExampleNotebooks/blob/main/cookbook/generate_lomap_network.ipynb
You will have to specify an atom mapper and a scorer. In this example, the Lomap mapper and Lomap scorer are used.

Then you can go through the edges and get e.g. the ligand names like this:

edges = []
for e in ligand_network.edges:
    edge = (e.componentA.name, e.componentB.name)
    edges.append(edge)

You could also extract the mapping. The edges contain the SmallMoleculeComponents for the two endstates, the mapping and annotations, such as the score.

Regarding your cinnabar question: Yes, you can calculate MLE derived DG values without have any experimental values. Please check out this function to give you an idea on how to do it: https://github.com/OpenFreeEnergy/IndustryBenchmarks2024/blob/1cf4accc3263d1adfcb141fcfa3d28a470698f39/industry_benchmarks/utils/extras/extract_results.py#L76
You can then afterwards shift the calculated DG values (that are by default centered around zero) if you have the experimental DG for one of the ligands.

[AminaMnhr]

Hi Hannah,

Thank for your answer!!

1. For LOMAP: I am not sure I need a mapper in my case as my method does not use soft core potentials anymore. It is called ATM : https://www.compmolbiophysbc.org/atom-openmm/the-science-behind-atom
   But I will try your method to get the edges, I will let you know how it goes

2. For your cinnabar suggestion: do you mean I just need to use the function get_dg ? However how should be structured the input calc_data dictionary (key and values) to apply get_dg?

Can I still otherwise use this method instead of using get_dg : https://github.com/OpenFreeEnergy/ExampleNotebooks/blob/main/plotting_rbfes_with_cinnabar/PlottingFreeEnergiesUsingCinnabar.ipynb

but by just having one exprimental value in my cinnabar_input.csv

3. Lastly in this same tutorial : https://github.com/OpenFreeEnergy/ExampleNotebooks/blob/main/plotting_rbfes_with_cinnabar/PlottingFreeEnergiesUsingCinnabar.ipynb

How can we print the values of calculated absolute free energies after they have been shifted (not centered anymore in 0). Because the tutorial "Writing out the MLE derived absolute free energies" at the end seems to print the calculated_results still centered around 0

[hannahbaumann]

Re 2: Yes, you can also use that notebook! It also shows you an example of how to create the dictionary that is taken in in the get_dg function, in the notebook called calculated_data.
Re 3: Cinnabar currently shifts the calculated data by the difference of the mean of the predicted and the mean of the experimental DG which doesn't make sense when you only have one experimental value. Therefore, I think you would want to first get the calculated DG values centered around zero and then shift them by the difference between the exp and the calc DG value of the one ligand where you have the experimental value present.
Please let me know if you have any questions to that!

[AminaMnhr]

Thank you so much for your reply!

Re 3) Sorry, I meant in the case of computing a benchmark like TYK2 I have indeed all the experimental values already and not one. So I was wondering how can I print the .dat file with the calculated dG shifted? So not centered anymore on 0

4. Is there a way with cinnabar_plotting.plot_DGs(fe.to_legacy_graph(), figsize=5, shift=shift.m), which is the function to plot the DGs to have the light grey zone at -2/2 kcal/mol and the dark grey at -1/1 ? Instead of having them at -1/1 and -0.5/0.5 as it is currently the case.

Thank you so much!

[IAlibay]

Re 3) Sorry, I meant in the case of computing a benchmark like TYK2 I have indeed all the experimental values already and not one. So I was wondering how can I print the .dat file with the calculated dG shifted? So not centered anymore on 0

Apologies I'm not I understand what you mean here. If you have the dG values centered around 0 and also your experimental values, then you can shift your calculated results by passing in the experimental mean. You can see an example of this in this issue: OpenFreeEnergy/cinnabar#42

Is there a way with cinnabar_plotting.plot_DGs(fe.to_legacy_graph(), figsize=5, shift=shift.m), which is the function to plot the DGs to have the light grey zone at -2/2 kcal/mol and the dark grey at -1/1 ? Instead of having them at -1/1 and -0.5/0.5 as it is currently the case.

Right now I don't think there's an easy way of doing this by passing arguments directly to the method. If you could raise an issue about this on the cinnabar repository that would be great!

One possible option (I haven't tried this and maybe it won't work), is to take the matplotlib Figure object that is returned by the method and redo the fill_between calls. It might be a bit messy to do though.

[AminaMnhr]

Hi,

I am using this function : https://github.com/OpenFreeEnergy/IndustryBenchmarks2024/blob/1cf4accc3263d1adfcb141fcfa3d28a470698f39/industry_benchmarks/utils/extras/extract_results.py#L76

To convert my ddG values to dG values. However my problem is I have some ligands in my set where I know the experimental values already, so how do I take account of those values when I convert ddG=>dG for each ligand with this function?

(I did not use your cinnabar tutorial methods because I have one or two experiemntal values and not ALL the experimental values for my set of ligands, so it won't work.)