TDC Complex Energy vs TDC Delta Energy

[Isaac Silverman]

Hi,

I have performed GB energy calculations with MMPBSA of a simulation between two proteins and got calculations for TDC Complex, TDC Receptor, TDC Ligand, and TDC Delta.

Can someone please explain to me the difference between TDC Complex and TDC Delta, specifically in regards to the heat map differences?

Are there any tutorials that will explain more about data analysis?

Thank you in advance.

[Mario Sergio Valdes]

The TDC (Total Decomposition Contribution) corresponds to SDC (Sidechain Decomposition Contribution) + BDC (Backbone Decomposition Contribution). For ligands other than amino acids, TDC = SDC because BDC is discarded since it does not exist.
Each component of the thermodynamic cycle has TDC, SDC, and BDC in decomp calculation(https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/introduction/). As is the case for the interaction energy calculation, delta = Complex - (Receptor + Ligand). Thus, the deltaTDC = comTDC - recTDC - ligTDC. In terms of charts, the only difference between comTDC and deltaTDC is the values. Commonly, only the delta is used for reports.
You can check more information here (https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA/issues/296)
We are sorry for the lack of information. This is something that we have yet to document better.
Let me know any doubts

[Isaac Silverman]

I am having trouble understanding the Delta compared to the TDC Complex results. For example, in the heat map for the complex, the line across fro each residue is stable and remains an almost solid blue throughout all the frames. However, in the Delta, the same amino acid will now be unstable and vary shades of blue and red. I am confused if the specific residue is stable or not with its energy? How is the specific energy for each residue being calculated using this formula delta = Complex - (Receptor + Ligand)? Is there excess energy that remains over from their interaction? Does the complex or delta show me more valuable data?

Also, when I generate the energy calculation, I only receive value for TDC. In the second link you provided you mentioned that "Unless you are interested in mutating amino acids and want to estimate which ones contribute their side chains, TDC should be sufficient." How would I be able to seee the SDC and BDC on their own? How does having those benefit mutation analysis?

Thank you for you guidance

[Mario Sergio Valdes]

The thermodynamic cycle to compute PB/GB energy is the delta (interaction energy) = Complex (bound state) - [Receptor + Ligand] (unbound state). Then, to compute the interaction energy, either the energy of a residue or the total (sum of all the residues), the energy of the bound and the unbound states must be calculated. In this sense, the complex interaction energy can be used as a stability measure. However, if you want to know whether the receptor and the ligand binding is thermodynamically favorable, you must compute the difference between the bound and unbound state, in this case, it corresponds to the enthalpy variation (delta here is equivalent to ΔH since the entropic component is not calculated).

The chart display depends on the component (com, rec, lig, or delta) and the energy values obtained. Note that in the start dialog, there is a checkbox to check if you want to hide the non-contributing residues. This option is checked by default, so you can hide residues in the delta that are displayed in the complex. In this case, the heatmap can be different from the complex one and asymmetric.
As I have explained, the delta depends on the com, rec, and lig values. In the case of decomposition, the contribution of each residue is calculated as follows:
For residues belonging to the receptor
ΔHres = ΔHres[com] - ΔHres[rec]
For residues belonging to the ligand
ΔHres = ΔHres[com] - ΔHres[lig]

I don't quite understand what you mean regarding the energy values per frame. Can you send me the file COMPACT_MMXSA_RESULTS.mmxsa to examine it?

To compute the SDC and BDC you must set dec_verbose = 1 or dec_verbose = 3 (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/input_file/#decomp-namelist-variables)

An example of SDC and BDC use is for computational peptide design. Say you are computing the interaction energy of the peptide with the receptor and would like to estimate the contribution of each residue to generating affinity-enhancing mutants. In this case, it is useful to know which region of the amino acid has the greatest contribution. If the major contribution is given by the backbone interaction (BDC), then the mutation of this residue would not have such a marked effect on affinity. If, on the other hand, the interaction energy is mainly given by the side chain (SDC), then the mutation would have a greater effect on the affinity. It is important to clarify that these methods are predictive, so the researcher will be responsible for the interpretation of the results and will have to corroborate them experimentally for greater certainty.

[Isaac Silverman]

Thank you for the explanation. So the complex energy is a measure of each residues bound vs unbound energies which can be used to determine stability (a solid colored line across indicating stable). And the delta is a  measurement of the thermodynamic favorability. In that case, can the delta energy be used to rank the strength of the residue interactions? (darker blue meaning more strong and favorable and more red the less favorable and weaker interactions?) Is that all correct?

What I was referring to with the energy values per frame was about difference between the the complex heat map and delta heat map. I have attached the file you requested. I believe that your explanation between complex and delta made me understand why the heat maps are different. The complex heat map residues have a solid color across all the frames of the simulation because it is showing a stable energy at the residue. For example Residue A:SER:189 has a solid blue line across 1000 frames indicating a uniform energy of -100 kcal/mol in the complex energy heat map. Then in the delta heat map Residue A:SER:189 has its energy fluctuating around positive 3 kcal/mol across the 1000 frames with a red line. Would this be indicating that Residue A:SER:189 has a stable energy of interaction as shown by the blue line on the complex heat map, but it is not thermodynamically favorable as indicated by the red line in the delta heat map? I have included images of the two heat maps for reference of this.

Also, regards to the start dialog having the option to hide the non-contributing residues, when I choose to display the non contributing residues, there still seems to be some negative energy values with those residues, although most of the frames show no energy. Is the reason they were left out because even though for a hundred frames or so there was energy, because the rest had very little energy that it was excluded?

[Mario Sergio Valdes]

Thank you for the explanation. So the complex energy is a measure of each residues bound vs unbound energies which can be used to determine stability (a solid colored line across indicating stable). And the delta is a  measurement of the thermodynamic favorability. In that case, can the delta energy be used to rank the strength of the residue interactions? (darker blue meaning more strong and favorable and more red the less favorable and weaker interactions?) Is that all correct?

It is correct. Actually, the delta is the one used to estimate how strong is the interaction.

 

What I was referring to with the energy values per frame was about difference between the the complex heat map and delta heat map. I have attached the file you requested. I believe that your explanation between complex and delta made me understand why the heat maps are different. The complex heat map residues have a solid color across all the frames of the simulation because it is showing a stable energy at the residue. For example Residue A:SER:189 has a solid blue line across 1000 frames indicating a uniform energy of -100 kcal/mol in the complex energy heat map. Then in the delta heat map Residue A:SER:189 has its energy fluctuating around positive 3 kcal/mol across the 1000 frames with a red line. Would this be indicating that Residue A:SER:189 has a stable energy of interaction as shown by the blue line on the complex heat map, but it is not thermodynamically favorable as indicated by the red line in the delta heat map? I have included images of the two heat maps for reference of this.

It is important not to confuse the terms. When we speak of stability here we are referring to the fact that the energy remains in a small fluctuation range (even if it is favorable, i.e. negative). This energy corresponds to the internal energy of the complex as a whole, not to the interaction. To compute the interaction energy it is necessary to take into account the bound state (complex) and the unbound state (receptor or ligand, depending on the selected residue). Then, to better understand what happens you must visualize the energy of the residue of your interest in the complex, in the receptor, and in the delta.
Here I explain the example of the residue R:A:SER:189.
As you can see in the picture attached, the energy value of the delta is generally positive. This is because the energy of this residue in the unbound receptor is more favorable than in the complex. Therefore, this residue will have an unfavorable interaction with the ligand.

 

Also, regards to the start dialog having the option to hide the non-contributing residues, when I choose to display the non contributing residues, there still seems to be some negative energy values with those residues, although most of the frames show no energy. Is the reason they were left out because even though for a hundred frames or so there was energy, because the rest had very little energy that it was excluded?

 No. This algorithm is very simple. Only the residues whose average contribution is less than the set cutoff are hidden.

As far as I can see, the energy of your system is not entirely stable.  Visualizing the structure I see that it is also a very flexible system. Is your trajectory centered? Did you eliminate the PBC? It is important to make sure that the trajectory is consistent

 

[Isaac Silverman]

I think I understand now the equation delta = complex - (receptor+ligand). That the complex energy is of the bound state of the receptor and ligand, and the receptor an ligand energies are of each of their individual unbound states? Then the complex (bound state) minus both of receptor and ligand (unbound states) is the total complex delta. And if the delta in negative the bound state is favored and if its is positive the unbound is favored. Is that right?

Therefore, would the delta values be the most important for analysis?

In regards to the stability of the the system. I believe that I centered my trajectory. After the simulation completed I used these two line to center and fit the system:

gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -center -pbc mol -ur compact

gmx trjconv -s md.tpr -f md_center.xtc -o md_fit.xtc -fit rot+trans

I then used this line to run the gmx_MMPBSA, where 20 was my receptor and 21 was my ligand:

gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 20 21 -ct md_fit.xtc -cp topol.top

Was I correct in using the md_fit.xtc for the energy calculations instead of the md_center.xtc file? Could that be what is making the system unstable? Have I eliminated PBC through these lines?

Thank you so much for your help.

[Mario Sergio Valdes]

I think I understand now the equation delta = complex - (receptor+ligand). That the complex energy is of the bound state of the receptor and ligand, and the receptor an ligand energies are of each of their individual unbound states? Then the complex (bound state) minus both of receptor and ligand (unbound states) is the total complex delta. And if the delta in negative the bound state is favored and if its is positive the unbound is favored. Is that right?

That's right! 

 

Therefore, would the delta values be the most important for analysis?

Generally, yes, at least to estimate the ligand affinity

 

In regards to the stability of the the system. I believe that I centered my trajectory. After the simulation completed I used these two line to center and fit the system:

gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -center -pbc mol -ur compact

gmx trjconv -s md.tpr -f md_center.xtc -o md_fit.xtc -fit rot+trans

I then used this line to run the gmx_MMPBSA, where 20 was my receptor and 21 was my ligand:

gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 20 21 -ct md_fit.xtc -cp topol.top

Was I correct in using the md_fit.xtc for the energy calculations instead of the md_center.xtc file? Could that be what is making the system unstable? Have I eliminated PBC through these lines?

The first command removed the PBC, while de second fit the coordinates to a reference structure. This means that the first one removes the discontinuities of the system (commonly visualized by artificial lines extending from one side to the other on the display). However, if the system rotates or translates, visually, you will observe jumps. The second command aligns the structures to a reference one by removing rotations and translations. Please review this section of our documentation which has visual elements.

As I mentioned before, from what I observed, your structure was very flexible, which can cause this "instability" or jumps in energy. I suggested you remove the PBC because it is common that it is due to this and tends to be forgotten. However, if it was already removed and the energy you get and the jumps, it is probably due to the characteristics of your system. A simple way to find out is to analyze the RMSD. If it is flexible, you should see considerable variations in RMSD and energy.

[Isaac Silverman]

Great! I am happy to understand the delta values now.

Do you think I would receive different energy values if I used the md_center.xtc file instead of the md_fit.xtc file created from the second line? Does changing the structure to a reference one impact the energy calculation. I am having trouble finding that documentation section. Could you please attach a link?

I have calculated the RMSD for the system which ranges between 0.5 and 1.5nm in a square root function manner. Would that signify the system is flexible or that the jumps are just a characteristic of my system as you suggested?

Thank you again for your help.

[Mario Sergio Valdes]

Do you think I would receive different energy values if I used the md_center.xtc file instead of the md_fit.xtc file created from the second line? Does changing the structure to a reference one impact the energy calculation. I am having trouble finding that documentation section. Could you please attach a link?

It should not, as there are no structural changes. The energy is computed on the structure, where each frame is independent. As long as it is consistent, then it does not matter if it is adjusted to the initial structure of the trajectory or not, the energy will be the same. https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/gmx_MMPBSA_running/#before-running-gmx_mmpbsa 

 

I have calculated the RMSD for the system which ranges between 0.5 and 1.5nm in a square root function manner. Would that signify the system is flexible or that the jumps are just a characteristic of my system as you suggested?

Commonly accepted as a good RMSD for proteins simulations should be less than 4A. In this case, 5-15A is a considerable RMSD. This is not necessarily wrong but needs to be reviewed carefully

 

[Isaac Silverman]

Ok, thank you very much for all your help explaining. I tried it and got energy values with slight differences, Ex: total complex energy was -26048 vs -25995 kcal/mol for the fit vs center files. It was similar small differences for almost all energies, including the individual residue decompositions. Is this expected?

Also, when I ran a simulation for a complex and got certain energy values. I then remade the complex with the same two proteins and ran another simulation and got the energy values for that. Similarly to how I described above, I got slightly different values the second time I did it, both using the fit file for the energy calculation. Is this expected that the same proteins will get slightly different values when repeated with another prepared complex? I have found that if an energy calculation is run for the same simulation twice it gives exact values again, but it is this slight difference if I remake the simulation and then run it foor energy calculation.

Some questions about the TDC, BDC, and SDC. I changed my input file to have dec_verbose = 3 but I only got TDC and SDC in my output. idecomp was set to 1 for per residue analysis. Why did I not get the BDC? Also, the graphs for TDC and SDC look identical. What can I conclude about that? The system is two proteins, one with a mutation.

Also, for mutation analysis for effect on binding affinity, is it better to use SDC and BDC of per residue or per wise?

Also, for the line plots of A:SER:189 that you showed me before, how did you get the individual plots for that residue? I only have access to the bar plots for individual residues in the per residue analysis. In the per wise I have access to line, bar, and heat maps for all the individual residues, but not for per residue. I attached image.

Thank you again.

[Mario Sergio Valdes]

Ok, thank you very much for all your help explaining. I tried it and got energy values with slight differences, Ex: total complex energy was -26048 vs -25995 kcal/mol for the fit vs center files. It was similar small differences for almost all energies, including the individual residue decompositions. Is this expected?

Yes, since the no PBC and the fitted trajectory are basically the same. The fitted one is the no PBC trajectory but with no rotation or translations. So, if no structural modifications were made, for energy calculation both work similarly.

 

Also, when I ran a simulation for a complex and got certain energy values. I then remade the complex with the same two proteins and ran another simulation and got the energy values for that. Similarly to how I described above, I got slightly different values the second time I did it, both using the fit file for the energy calculation. Is this expected that the same proteins will get slightly different values when repeated with another prepared complex? I have found that if an energy calculation is run for the same simulation twice it gives exact values again, but it is this slight difference if I remake the simulation and then run it foor energy calculation.

Yes. When you run two different simulations, even starting from the same initial structure, the trajectories may be different. This is because you start with different velocities. Some studies find that MMPB/GBSA works well when you average different short trajectories instead of one long trajectory. If the system is very flexible, one way to "explore" the conformational space, at least in a limited way, is to carry out MD from intermediate states as follows:

Frames                 1                                               100

                               |------------------------------------------------|

 >>> Clustering         |  |               |      |                     |     --> Representative structures

                                   I  II             III    IV                   V

 New MD for each cluster

----

This method is deterministic since energy is computed over a static state with known and invariant variables. That means, if you compute the energy with the same parameters over the same trajectory, you will get the same result.

 

Some questions about the TDC, BDC, and SDC. I changed my input file to have dec_verbose = 3 but I only got TDC and SDC in my output. idecomp was set to 1 for per residue analysis. Why did I not get the BDC? Also, the graphs for TDC and SDC look identical. What can I conclude about that? The system is two proteins, one with a mutation.

I will check this. It should show both the SDC and the BDC for protein

 

Also, for mutation analysis for effect on binding affinity, is it better to use SDC and BDC of per residue or per wise?

They are simply complementary analyses and both are computed with the same option. As I explained in another thread, this allows you to estimate whether the residue contribution comes from the backbone or from the side chain, which would be useful for making point mutations. Note that these are approximations and depend on the structure you are analyzing. It may be that a mutation changes the structure of that region, and this analysis is not sufficient. In principle, if you find that interaction analysis at this level provides some benefit, use it. If you will generate a mutant library because you have the resources, perhaps you can use TDC alone. Keep in mind, that sometimes you have to leave a bit of randomness in these types of studies because the predictions can be considerably biased. For example, in antibody engineering, you estimate which residues are most important, both energetically and structurally, and generate mutant libraries with doped codons (they have a high probability of coding for the native amino acid). This allows the generation of mutants similar to the native Ab but with one or two mutations. The mutated positions, as well as the studies, can be directed to try to find mutants with higher affinity, validate the role of the aa in the interaction, whether structural or energetic, etc. If you use computational predictions to narrow down the search too much, you may leave out solutions that cannot be or have not been explored in your study. So, be careful!

 

Also, for the line plots of A:SER:189 that you showed me before, how did you get the individual plots for that residue? I only have access to the bar plots for individual residues in the per residue analysis. In the per wise I have access to line, bar, and heat maps for all the individual residues, but not for per residue. I attached image.

As you can see in the image, each residue has a bullet on the left. You can click it or double-click on the residue element to display the rest of the charts. This happens for all options. Perhaps this is why you have not quite understood how per-wise works. Please check all the components of each element. Please, keep in mind the description I have given you above of what each graphic represents.