Research‎ > ‎

Research Blog

The Standard State in Binding Free Energy Calculations

posted Jul 30, 2019, 4:10 PM by Emilio Gallicchio   [ updated Jul 30, 2019, 7:01 PM ]

Phys. Chem. Intro

Q: What is the standard free energy of binding?

It is the free energy of binding between a ligand and a receptor when they are in an ideal solution at the standard concentration (C=1 M). The standard free energy of binding is usually denoted by ΔGb

Q: how is ΔGmeasured in practice?

There are many ways to measure a standard binding free energy, the most straightforward is to measure the binding constant Kb, or equivalently the dissociation constant Kd= 1/Kb

Q: But wait, the relation above does not make sense because it takes the log of a quantity in concentration units.

In Physical Chemistry equilibrium reaction constants are dimensionless. The dimensionless nature of the equilibrium constant for binding is clear from its definition:

where R is the receptor, L is the ligand and RL is the complex. Changing the units of concentration does not change the binding constant.

Q: yeah ... I am not sure, my friend works in this lab and they measured the Kd of a drug and he said it's 1.8 nM

In fields such as Biochemistry, Medicinal Chemistry, and Biology, equilibrium constants are very confusingly reported with units, sometimes very strange units. Unless they are doing something really odd, such as redefining the standard state, read the above as Kd = 1.8 x 10-9, no units, and you'll be fine.

Q: really, I can change the standard state?

If you want, sure, but maybe don't call it "standard". A standard it's just a reference that people have agreed on, a bit like the origin of a coordinate system. You can certainly define your own reference solution state, I don't know ... 34,512 molecules per cubic millimeter. (As an example, for acid/base reactions, biologists like to use 10-7 M as the standard concentration of the hydronium ion, causing endless confusion.) Anyway, if you do use your own standard do not expect to get the same "standard" free energy of binding and equilibrium constant as the rest of the world, which uses molar units, for the same reason that you can't expect to get the same distance from the origin if you change the origin.

Let's stick to the 1 M standard state, please ...

The Standard State in Binding Free Energy Calculations

Q: Why do I need to report the standard free energy of binding?

The binding free energy depends on the concentrations of receptor and ligand. For example, the binding free energy is zero at equilibrium concentrations. It is meaningless to report a free energy of binding without specifying the concentrations it corresponds to. By stating that a free energy of binding is "standard", we tell people that it corresponds to the standard state at 1 M concentrations.

Q: Do I need to do a simulation at the standard state? That is, do I need to place enough ligand and receptor molecules in the simulation box so that they are at 1 M concentration?

No. You can do a calculation at any concentration you want. As the simulation progresses, monitor the number of complexes, RL, and the number of free R and L to measure their equilibrium concentrations. That will give you the equilibrium constant and standard free energy of binding using the relation above.

However, the binding constant is rarely computed by "counting" in this way. We like to do computational alchemy. It is faster and more reliable.

Q: Okay, I did an alchemical calculation. I decoupled the ligand from the solution then I coupled it to the receptor. I fed the data into MBAR etc. which spit out a free energy. Can I call it a standard free energy of binding?

Generally, no. What you have calculated is the excess component of the free energy of binding. To turn it into a standard free energy of binding, add the ideal term:

where Vsite is the volume of the binding "region". 

Q: Okay, that sounds easy. Where does the quantity ΔGideal come from?

 ΔGideal  is the reversible work for transferring a ligand from an ideal solution at concentration C to the binding site of volume Vsite . Read these papers:
    1. MK Gilson, JA Given, BL Bush, JA McCammon. The statistical thermodynamic basis for computation of binding affinities: a critical review, Biophysical Journal, 72, 1047 (1997).
    2. M Mihailescu, MK Gilson. On the theory of noncovalent binding. Biophysical Journal 87 (1), 23-36 (2004).
    3. E Gallicchio and RM Levy, Recent Theoretical and Computational Advances for Modeling Protein-Ligand Binding Affinities. Advances in Protein Chemistry and Structural Biology, 85, 27-80, (2011). 

Q: But how do I get  Vsite?

Getting Vsite is easy. You have set it yourself. It is whatever volume the ligand is allowed to explore during the process of alchemically coupling the ligand to the receptor.

Q: What do you mean? I did not restrain the ligand.

Then Vsite is the volume of the simulation box. Did you run the simulation long enough so that the ligand visited the whole box? Probably not, eh!? Then the simulation is not converged. Even more troubling is that you also have implicitly defined the complex RL as any conformation in which R and L are in a region of solution of the same volume as the simulation box. This is probably not what you wanted. You probably wanted to measure the standard free energy of binding of the ligand to a specific region of the receptor, not for the whole receptor and even including the solvent. You should have restrained the ligand to that region. Vsite  is the volume of whatever that region is that you meant to consider.

Q: Right, right ... I forgot about it. I did apply a restraint potential to avoid the "wandering ligand" problem when the ligand and receptor are decoupled.

Ok, great. Than you have a Vsite after all. Is the restraint potential based on the CM-CM distance between ligand and receptor atoms? If so, simply integrate the Boltzmann factor of the restraint potential over the three coordinates of space to get Vsite. Or is the restraining potential something more complicated?

Q: Yeah ... more complicated. The orientation of the ligand was also restrained.

No problem. Integrate over the orientational angles as well and get the angular binding site volume Ωsite and add -kB T ln Ωsite⁄8 π 2 to the ideal standard state factor. Read the papers I told you about.

Q: No, no. I restrained the orientation of the ligand by tethering two atoms of the ligand to two receptor atoms.

That is not allowed. It would perturb the intramolecular conformational distributions of receptor and ligand when they are not coupled. Read this paper:

Boresch S, Tettinger F, Leitgeb M, Karplus M. Absolute binding free energies: A quantitative approach for their calculation. J. Phys. Chem. B. 107, 9535–9551 (2003)

Q: I see. But after coupling the ligand to the receptor I turned off the restraints. Does that change  Vsite?

The binding restraints are not meant to be turned off. They are on during the whole coupling alchemical leg and they stay on.

Q: But then I am not simulating the real thing! The real complex does not have restraints!

I never said that you are simulating the real thing. We are doing computational alchemy! The binding restraints define the complex. With the binding site restraints, you are essentially specifying the set of conformations of the receptor-ligand system that form what you call the "complexed" state. You can't change the definition of the complex depending on the alchemical state of the system.

Q: Okay, now I know you finally went off the deep end. You cannot be right because then the standard free energy of binding would depend on how I set the restraints!

Precisely. The free energy of the complex depends on how you defined it. This is not different from, say, having to define the ranges of φ and ψ backbone angles of the alpha helix state in order to compute the alpha-helix population of a peptide. Change the range of the angles, and the population, as well as the corresponding free energy, changes. You need to define the complexed state before you can measure its free energy!

Q: But then I do compare the calculated free energy to the experiment? The experimental system does not have restraints. What does the experiment measure?

Good question and a very interesting topic. Read paper no. 2 above.

Q: Okay, but now I am worried that I have used harmonic restraints centered on a specific pose of the ligand which I got from docking and I am not sure it is the correct bound pose. Is it a problem?

It could be. It is safer to use flat-bottom harmonic restraints for the binding site volume. That way you have more tolerance and a greater chance of including the "correct" pose.

Q: I am starting to get this. But I heard that imposing restraints and then releasing them can improve convergence. Is that wrong?

It is absolutely fine to impose and then release additional restraining potentials along the alchemical path to improve convergence. These do not need to obey the same strict rules for the binding site restraints discussed above. Just make sure that the binding site restraints stay on throughout the alchemical simulation.

Q: Now I understand why many people prefer to compute relative binding free energies (FEP) rather than the standard binding free energies! They do not need to worry about standard states, the binding site volume, and all this stuff!

Actually, all of this applies to relative binding free energies as well. The relative binding free energy is the difference between two standard binding free energies. Because these depend on the definition of the binding site, so does the relative binding free energy. In principle, without restraints, a ligand undergoing an FEP transformation will eventually dissociate and wander around the simulation box adversely affecting the relative free energy estimate. The longer the simulation the higher the chance that this will happen. You want your estimate to improve, not worsen, as you make the simulation longer. Binding site restraints should be used in FEP as well ... but people rarely do it. 

Don't ask. I really don't know why.

Binding Free Energy of Posaconazole in Two Protonation States to a Model of Captisol

posted May 17, 2019, 6:15 PM by Sheenam   [ updated May 19, 2019, 10:50 AM ]

Posaconazole is a triazole antifungal drug that is used to treat invasive infections by Candida species and Aspergillus species in severely immunocompromised patients. However, its property as a weakly basic and poorly aqueous soluble drug results in poor bioavailability and variable absorption. Hence, for a better intravenous administration of the drug under aqueous conditions, the drug comes in a composition which consists of a solubilizing agent, a modified β-cyclodextrin such as Captisol which supposedly solubilizes with the drug in acidic medium. To test this hypothesis, we calculate the free binding energy of the drug Posaconazole to Captisol in different protonated states to understand the pH dependence of the observed solubility of posaconazole in the presence of Captisol via computational methods. The results of the computational experiments conducted here confirm the hypothesis that this effect is due to the higher affinity of the protonated form of Posaconazole relative to the unprotonated one.


Fig 1: Captisol-Posaconazole Complex (Posaconazole in Green)



Conformational Propensities of FXR Drug Inhibitors in Water Solution

posted May 7, 2019, 7:18 AM by Brenda Neuman   [ updated May 8, 2019, 5:59 AM ]

Hypercholesterolemia is a major cause of heart disease. A potential drug therapy for hypercholesterolemia is FXR inhibition because of the role of FXR receptors in bile acid and cholesterol metabolism. Drug molecules need to often reorganize conformations to bind protein receptors. The free energy penalty of reorganization is related to how often in solution the drug molecule adopts the same conformation as in the protein. By comparing the conformations of the FXR inhibitors in an aqueous environment to the conformation of these inhibitors in the enzyme substrate complex the probable effectiveness of the inhibitor as a drug can be determined. Different force fields were used to evaluate the distribution of conformations of two FXR inhibitors in an aqueous environment in order to analyze their reorganization penalty. The results from the different force fields were compared and it was determined that there was similarities and also significant differences between the conformations the different force fields produced and the reorganization penalties associated with them.

                                                                                            Figure 3. FXR inhibitor in Water Box


Research Poster that was presented at Brooklyn College Science Day on 05/03/19 can be viewed through the following link 

Heating and Cooling Strategy to Improve Estimates of CYFIP1p-Derived Peptides’ Helicity Probabilities in Biased Molecular Dynamics Simulations

posted Feb 24, 2019, 2:31 PM by Irene R   [ updated Feb 26, 2019, 6:32 PM ]

A CYFIP1p-derived peptide in VMD: bond lengths of 2.5 Å (units used by VMD) or less

for the middle three bonds would indicate a helical structure.



By Irene Rostovsky


Since the binding region of CYFIP1p can be used as a model for a precursor molecule for the development of a drug to treat Fragile X syndrome and cancer, we wanted to calculate the probability of finding the peptide corresponding to this binding region, and another peptide differing by one residue, in a helical conformation. Due to the difficulty of attaining convergence of random and helical starting conformations for one of our peptides, we performed a series of simulations at 6.25 kJ/mol bias force (previously determined as the optimal force by Megan Wang), where we manipulated temperatures in a manner that is similar conceptually to replica exchange, but is simpler to perform. We alternated between 400 K and 300 K for succeeding simulations to disturb the initial states and speed up exploration of different conformations. Our results so far suggest that alternating heating and cooling is an effective way to help peptides break out of their initial conformations, and form/ fully break alpha helices at a reasonable rate over 1 ms. The next step is to perform statistical analysis on our results and to test the possible convergence we have seen through our new heating/ cooling method.


Water energy contributions to structural inhomogeneity of curcubit[8]uril

posted Jun 15, 2018, 4:54 PM by Fatlum Hajredini   [ updated Jun 17, 2018, 7:09 AM ]

Expulsion of water from cavities on host binding sites can contribute a great deal of energy to binding of ligands. This can be attributed to the gain in entropy as the waters become unstructured. Accounting for such effects in binding free energy calculations greatly improves the prediction accuracy. This effect is likely to also contribute to the structural reorganization of macromolecules. Using molecular dynamics simulation and the recently developed SSTMap suite we investigate the contribution of water thermodynamics to the tendency of curcubit[8]uril to assuming a compacted structure. Hydration site analysis shows that structuring of waters at the cavity of cb8 carries an entropic penalty which was expected to be relieved upon compaction of cb8. A cb8ligand complex had overall more favorable hydration upon assuming a more compacted structure compared to its more extended counterpart, suggesting that the compaction was due to water expulsion. To investigate whether the same phenomenon would be observed in the cb8 complex alone, simulations were carried in an explicit solvation model, and an implicit solvation model which lacks water structure information. Contrary to expectation, when simulated in the implicit solvation model cb8 is predominantly compacted, and also assumes an overall more compact structure compared to any of the conformations in the explicit water simulation. Simulations in explicit solvent show a wide distribution of states, with the extended conformation being much more populated then the counterpart in the implicit solvent simulation. Taken together, these findings suggest that structural networks of water, when treated explicitly can have a significant impact on the structural reorganization of macromolecules. 


                    State Distribution of cb8 when simulated in implicit solvent (blue), or explicit solvent (red) as determined by differences in their radii of gyration. 

Tripeptide Stabilized Nanoemulsions for Cancer Therapy

posted May 24, 2018, 9:18 PM by Marium

A carboplatin derivative consisting of two oleic acids attached to a platinum head, used to kill cancer cells, is investigated in this project through chemical simulation software. One of the challenges of this class of anticancerous drugs is its high level of toxicity to the human body, which results in low clinical dosings of the drug and ultimately ineffective therapy. Thus, for drug delivery to the human body, the carboplatin derivative requires an emulsifier to stabilize it and shield its toxic agents. The emulsifiers used in this project are tripeptides because of their biodegradability and it has been shown that they are promising self assembly candidates in a previous study. In the first phase of this project, general features of oleic acid aggregates in water solution are investigated. In the second, tripeptides, specifically KYF (Lysine- Tyrosine - Phenylalanine) and DFF (Aspartic Acid- Diphenylalanine), are added to these aggregates in order to study the tripeptide stabilized nanoemulsions. The simulations gave insight about the interactions present in the nanoemulsions. 
Structure of KYF and oleic acid aggregate
Structure of KYF(green) and oleic acid aggregate
(carbons are gray, carboxyl groups are in red, amine groups are in blue)


Use of Umbrella Sampling Methods to Estimate Probability of CYFIP1p Helical Conformation

posted May 20, 2018, 11:58 AM by Megan Wang

Previously, we made modifications to the wild-type CYFIP1p to stabilize its α-helical structure and increase its binding efficiency to EIF4E. In this study, we utilized umbrella sampling methods to accelerate the sampling of possible peptide conformations, such that the probability of these conformations could be readily calculated from an MD simulation that would otherwise have taken months to collect a sufficient number of relevant samples. Using two Python scripts, we were able to apply a bias force onto CYFIP1p and measure the proportion of instances that CYFIP1p maintained a helical conformation. We effectively estimated the probability of such an outcome to be 0.18 percent if the bias force were nonexistent.

Potential bias force applied to interacting atoms in molecular system.  

Research Report

Possible Computational Evidence for Enhanced α-Helix in Modified CYFIP1-derived Peptides

posted Dec 7, 2017, 11:03 AM by Megan Wang   [ updated Mar 17, 2018, 9:29 PM by Emilio Gallicchio ]


Dysregulation of the eukaryotic translation initiation factor (eIf4E) has been shown to exist in Fragile X Syndrome, a leading cause of intellectual disabilities such as autism. The cytoplasmic FMRP-interacting protein 1 (CYFIP1) plays a key regulatory role in repressing associated mRNA translation by binding to eIF4E. A crucial secondary structure element in the interactions between a CYFIP1-derived peptide (CYFIP1p) and eIF4E is the α-helix, and as a consequence, improving the persistence of α-helicity of the peptide could lead to improved binding efficiency. In our study, we made computer-aided chemical modifications in an effort to further stabilize the α-helix structures of peptides derived from wild-type CYFIP1p. Our findings suggest the addition of a staple comprised of alkyls or an aromatic ring has no significant impact on the secondary structure elements of the CYFIP1-derived peptides. However, modifications comprised of a long chemical staple combined with a mutation from serine to alanine resulted in improved α-helix stability and thus, exhibited a potential for enhanced binding efficiency. 

Snapshot of the trajectory for a CYFIP1-derived peptide exhibiting α-helix fold. 



Gaussian-based Volume and Surface Area algorithm for GPUs

posted Jan 8, 2017, 9:26 AM by Emilio Gallicchio   [ updated Jun 5, 2018, 11:07 AM ]


The model proposed by Grant & Pickup [J. Phys. Chem. 99, 3503-3510 (1995)] is a well-established method to estimate accurately volume and surface areas of molecules. The method is based on the inclusion-exclusion formula of statistical physics (also known as the Poincaré formula). It states that the volume of an object composed of multiple overlapping bodies is given by the sum of the volumes of the bodies, minus the sum of the overlap volumes between pairs of bodies, plus the sum of the overlap volumes between triplets of bodies, and so on:
Atomic volumes are represented by Gaussian densities and overlap volumes are computed using standard Gaussian overlap integrals. The model is fully analytic. For example, the solvent-accessible surface area of an atom is obtained as the derivative of the molecular volume with respect to the atomic radius of an atom.

The model leads to a tree-based algorithm in which overlap volumes are recursively evaluated. 
overlap tree
We have used this model extensively to implement the non-polar component and the self-adjusting pair-wise descreening method for the Generalized Born solvation electrostatic component of the of the AGBNP solvation model. Our CPU implementation is based on a dept-first traversal of the tree.

Lately we have been working on a GPU implementation of the Gaussian volumetric model. As one can imagine, tree-based recursive algorithms do not easily lend themselves to GPUs. However, we finally found an efficient solution based on a flat arrays representation of the tree and breadth-first traversal (see paper in JCC).
GPU tree traversal
The resulting GPU implementation is 50 to 100 times faster than our best CPU implementation. With the help with the development team at Stanford, the algorithm is now available as a plugin of the OpenMM molecular mechanics package.





Dopamine D3 receptor antagonists

posted Oct 31, 2016, 1:16 PM by Pierpaolo Cordone   [ updated Feb 1, 2018, 9:10 PM by Emilio Gallicchio ]

According to Newman et al., in all D3 receptor antagonists a salt bridge between the protonated amine in the primary pharmcophore and Asp 110 is observed. Previous antagonists candidates having a stepholidine ring as a primary pharmacophore have been synthetized in our lab. The protonated nitrogen of the pharmacophore makes the same interaction with Asp 110 observed in Newman et al.10 Some of the molecules used in this study like 216F (figure 4 left) have the same feature. However, molecules like 217F don't do the same interaction (figure 4 right). In order to find out the reason why those molecules interact differently 216F and 217F were superimposed (figure 5). From the superimposition it looks that there is electronic repulsion between a cyano group in para of the aromatic ring and the carbonyl of Val 189. This repulsion would cause the molecule to rearrange in another way in order to fit in the receptor. This rearrangement would prevent the protonated amine of 217F to make the salt bridge with Asp 110 in the primary binding site (OBS). This phenomenon occurs in all molecules with a substitution in para.

Figure 4. left, SG-216F makes the salt bridge with Asp 110. Righ, SG-217F does not make the salt bridge interaction with Asp 110. 

Figure 5. SG-216F in blue and SG-217F in pink. SG-217F cannot make the salt bridge interaction because of the clash between the cyano group and valine 189.

1-10 of 27