Research Blog
Maximum Likelihood Inference of the Symmetric DoubleWell Potential
The doublewell potential serves as an optimal, onedimensional model for exploring physical phenomena. This project aimed to estimate the probability distribution of the doublewell potential fitted to a Gaussian mixture model by maximum likelihood inference. Hypothetical datasets of the quadratic function were generated using smartdarting MonteCarlo simulations under the Metropolis acceptance criterion. A major challenge in studying free energy landscapes is sampling efficiency, therefore a larger displacement was implemented in the MonteCarlo simulations that are generally done at smaller displacements. Although a conventional MonteCarlo approach can accomplish displacement from one free energy minimum to the next, this is done at the expense of accuracy, such that the acceptance rate of the simulation deviates from an optimal 50%. The generated datasets demonstrated close resemblance to normalized Boltzmann Distribution functions at two temperatures. TensorFlow was utilized to obtain the most optimal parameters for observing the generated datasets based on the Gaussian Mixture Model. Our results demonstrated that TensorFlow approximates a set of optimal parameters for the Gaussian Mixture Model that accurately resemble the Boltzmann Distribution at 300 K, however the method is unable to do so with similar accuracy at 2000 K. 
The Standard State in Binding Free Energy Calculations
Phys. Chem. IntroQ: What is the standard free energy of binding?
Q: how is ΔG^{◦}_{b }measured in practice?
Q: But wait, the relation above does not make sense because it takes the log of a quantity in concentration units.
Q: yeah ... I am not sure, my friend works in this lab and they measured the K_{d} of a drug and he said it's 1.8 nM
Q: really, I can change the standard state?
The Standard State in Binding Free Energy CalculationsQ: Why do I need to report the standard free energy of binding?
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?
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?
Q: Okay, that sounds easy. Where does the quantity ΔG^{◦}_{ideal }come from? ΔG^{◦}_{ideal } is the reversible work for transferring a ligand from an ideal solution at concentration C^{◦} to the binding site of volume V_{site} . Read these papers:
Q: But how do I get V_{site}?
Q: What do you mean? I did not restrain the ligand.
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.
Q: Yeah ... more complicated. The orientation of the ligand was also restrained.
Q: No, no. I restrained the orientation of the ligand by tethering two atoms of the ligand to two receptor atoms.
Q: I see. But after coupling the ligand to the receptor I turned off the restraints. Does that change V_{site}?
Q: But then I am not simulating the real thing! The real complex does not have restraints!
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!
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?
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?
Q: I am starting to get this. But I heard that imposing restraints and then releasing them can improve convergence. Is that wrong?
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!

Binding Free Energy of Posaconazole in Two Protonation States to a Model of Captisol
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: CaptisolPosaconazole Complex (Posaconazole in Green) 
Conformational Propensities of FXR Drug Inhibitors in Water Solution
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. 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 CYFIP1pDerived Peptides’ Helicity Probabilities in Biased Molecular Dynamics Simulations
A CYFIP1pderived 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
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 cb8•ligand 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
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. 
Use of Umbrella Sampling Methods to Estimate Probability of CYFIP1p Helical Conformation
Previously, we made modifications to the wildtype 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 CYFIP1derived Peptides
Gaussianbased Volume and Surface Area algorithm for GPUs
The model proposed by Grant & Pickup [J. Phys. Chem. 99, 35033510 (1995)] is a wellestablished method to estimate accurately volume and surface areas of molecules. The method is based on the inclusionexclusion 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 solventaccessible 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 treebased algorithm in which overlap volumes are recursively evaluated. We have used this model extensively to implement the nonpolar component and the selfadjusting pairwise descreening method for the Generalized Born solvation electrostatic component of the of the AGBNP solvation model. Our CPU implementation is based on a deptfirst traversal of the tree. Lately we have been working on a GPU implementation of the Gaussian volumetric model. As one can imagine, treebased 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 breadthfirst traversal (see paper in JCC). 