Post-Doctoral • GraduateUndergrad

Postdoctoral Research at the University of Minnesota

University of Minnesota logo

Monte Carlo Simulations Probing the Liquid/Vapor Interface of Water/Hexane Mixtures

Rendering of a simulation of a slab of water with hexane molecules on the top and bottom surfaceKnowledge of the interfacial properties of water/oil mixtures is extremely important for the petrochemical industry and for understanding detergency and hydrophobic effects. In 2010, Javadi et al. observed a significant decrease in the surface tension of water when alkanes, including n-hexane, adsorbed from the vapor phase onto the surface. Mucic et al. found that the surface tension of water decreases with increasing (partial) pressure of n-hexane. Here, we probe the liquid/vapor interface of water/n-hexane mixtures using configurational-bias Monte Carlo simulations in the NpNAT Gibbs ensemble where one box contains a liquid film with constant area and applied normal pressure and the other represents the isotropic vapor phase. We study the effect of n-hexane, at several different partial pressures, on the liquid/vapor surface tension, and analyze the simulation trajectories to determine the n-hexane loading and to provide molecular-level insights on the structure of the interface.

Gibbs ensemble Monte Carlo Simulations for Additive Loading in Surfactant Bilayers

Side by side images of renderings of a simulation of a surfactant bilayer solvated by water, one of which also contains nonane molecules inserted into the center of the bilayerUnderstanding solute uptake into soft micro-structured materials, such as bilayers and micelles, is of specific interest in the pharmaceutical, agricultural, and personal care industries. Currently, there are continuum-solvation methods available that can rapidly estimate the amount of a solute that will load into these soft micro-structured materials. However, these methods do not give any insight into the structural orientation and positioning of the solutes at a molecular level. To obtain this insight, we utilize the Shinoda-Devane-Klein (SDK) force field in conjunction wih configurational-bias Monte Carlo simulations in the osmotic Gibbs ensemble. One simulation box contains a lamellar system comprised of a bilayer formed by triethylene glycol mono-n-decyl ether (C10E3) surfactants surrounded by water with a 50/50 surfactant/water weight ratio. This simulation cell is allowed to undergo volume and shape fluctuation, by changing either the z cell length independently or the x and y cell lengths in a coupled fashion and is found to have a stable lamellar phase in agreement with experiment.

The other simulation box is a vapor-phase reservoir containing solute molecules at a constant applied pressure, and only solute molecules are allowed to transfer between the two simulation boxes. We study both the adsorption isotherm and structural effects of n-nonane, ethyl butyrate, and 1-hexanol on the bilayer at several different reduced pressures with respect to the saturated vapor pressure of the solutes. We observe that both ethyl butyrate and 1-hexanol load closer to the water and hydrophilic end of the surfactant, while n-nonane loads closer to the center of the bilayer, towards the hydrophobic end of the surfactant. These differences in the spatial distribution of the additive molecules lead to distinct changes in the structure of the bilayer. At high reduced pressures, the bilayer ruptures.

Vapor-liquid and Liquid-liquid Equilibria for a Langmuir Monolayer of Pentadecanoic Acid: A Monte Carlo Study

Rendering of a simulation of a monolayer of pentadecanoic acid moleculesA Langmuir monolayer consists of amphiphilic surfactant molecules that adsorb at a water/air interface to form a unimolecular layer. As well as having multiple applications in templating, chemical sensing, and coatings, these monolayers can serve as a model for lipid bilayers. The precise fluid phase behavior of Langmuir monolayers and other two-dimensional systems is difficult to accurately obtain through experiments. Molecular simulations are needed to probe the vapor-liquid equilibria of these complex systems. In 1994, Siepmann, Karaborni, and Klein [J. Phys. Chem. 1994, 98, 6675] used configurational-bias Monte Carlo simulations in the Gibbs ensemble (CBMC-GE) to study the vapor-liquid coexistence curve (VLCC) of a pentadecanoic acid monolayer, but these simulations were limited to a small system comprised of only 60 surfactant chains and relatively short trajectories due to the computational resources available. In this work, we carry out CBMC-GE simulations for much larger systems and using longer trajectories utilizing the TraPPE united-atom force field for the alkyl tails to probe the VLCC of a pentadecanoic monolayer. Furthermore, we investigate modification of the headgroup-headgroup potential to improve the accuracy of the predictions and to allow for extension of the simulations to coexistence of liquid-expanded and liquid-condensed phases.

Monte Carlo Simulations of Surfactant Bilayers

A rendering of a single pentadecanoic acid molecule rotating. Balloon shaped structures around each atom indicate the highest occupied molecule orbitals.

In my postdoctoral research project, I aim to gain a better understanding of the behavior of surfactant molecules through the use of Monte Carlo algorithms.  Due to their amphiphilic properties, surfactants preferentially orient and distribute themselves at water/oil and water/air interfaces, as well as self-assembling into micellar structures. Obtaining physical properties, such as the loading of surfactants, spatial distributions of different compounds, and surface tensions of surfactant systems would be beneficial for innovating new functions for this class of molecules.  Our research group’s highly specialized Monte Carlo algorithms and simulation protocols will allow for the precise study and prediction of the physical properties of surfactant systems.


Graduate Research at the University of Florida

University of Florida logo

Helicobacter pylori Urease Molecular Dynamics Simulations

Helicobactor pylori (HP) are gastrointestinal bacteria that infect the majority of the world’s population with no effective cure.  These bacteria survive stomach acidity by producing HP urease.  My Ph.D. projects entailed the provision of atomic-level understanding of the structure and function of the enzyme in order to pave the way for developing inhibitors for medical treatment.

Molecular Dynamics Simulations of HP Urease in Aqueous Solution:

A slideshow of three rendered images in which the flap covering the active site of helicobacter pylori urease is opening
A demonstration of the opening and closing of the active site-covering flaps.

The flaps that cover the active site cavity of the protein were found to open wider than previously observed, a finding which has significant implications for the design of HP urease inhibitors.

Molecular Dynamics Simulations of HP Urease in 10.5M Aqueous Urea:

A possible shuttling process of urea into the active site of the protein was observed when the behavior of HP urease in 10.5M aqueous urea was studied. These findings allow for the targeting of other regions of the protein that may negatively affect the enzymatic hydrolysis of urea, leading to a greater number of potential inhibitor binding sites.

A full image of the rendering of helicobacter pylori urease.
Horseshoe of high-rmsf regions surrounding the loop (residues 564-568) of the active site-covering flap (D4 and D5 refer to the dimer on which the residues are located): 564-568, red (D4); 771-791, blue (D4); 55-68, green (D5); 100-111, orange (D5); 113-130, black (D5); 177-184, cyan (D5); 222-238, purple (D5); 239-255, ice blue (D5). 

Molecular Dynamics Simulations of HP Urease in 150mM Aqueous NH4Cl:

Interactions between aqueous ammonium ions and HP Urease were monitored. A possible path in which ammonium ions are expelled from the active sites was observed. It was found that blockage of this pathway results in another possible binding site for inhibition.

Urease Docking

We obtained a structure of Klebsiella aerogenes urease (PDB ID: 1EJX) with the flap region added via homology modeling based on 2UBP (Bacillus pasteurii urease). Our group previously identified a hitherto unobserved wide-open flap state, and this was selected as the target for ligand docking studies.  Crystallographic water molecules were stripped from this structure and it was subsequently loaded into Maestro. Next, we prepared a grid centered on the active site, and docked (a) urea, (b) a 36,000 ligand library, (c) several libraries from the ZINC database including biologically relevant compounds, lead-like structures and drug molecules and finally, (d) selected individual compounds from the literature that had been previously identified as promising inhibitors of urease. From the output we generated a list of compounds in their successfully docked poses, sorted by docking score. We are currently analyzing the results with particular focus on the interactions between the binding pocket and the docked ligands.  We hope to identify promising candidates to be submitted for experimental assay in order to determine their efficacy as urease inhibitors.

Six rendering arranged in a two by three grid with the letters A - F labeling each image. The first row is A,B, the second is C,D, and the third is E,F

Depictions of the best-scoring docking poses for each of the 4 potential urease ligands (panels A-D) and the 2 control (panels E-F). Depictions of the best-scoring docking poses for each of the 4 potential urease ligands (panels A-D) and the 2 control (panels E-F). In each figure, the active site flap is shown in a blue cartoon representation, the Ni pair is shown as magenta-colored van der Waals spheres, the coordination sphere is shown in licorice representation colored by atom type, and the docked ligand–as well as the Cys319 and His320 residues–is shown as CPK representations colored by atom type. Shown in order are docked images of caftaric acid (A), desmosine (B), lebetol (C), kanamycin (D), quercetin (E), and eipgallocatechin (F).

Heats of Formation for Transition Metal Complexes

While ab initio calculations and density functional theory have been well validated for a wide variety of organic compounds, there is comparably little supporting work that has been done in the area of transition metal-containing complexes.  We have used fifteen density functionals and the MP2, CCSD, and CCSD(T) methods with Dunning double and triple-z basis sets, both with and without diffuse functions, in order to determine their utility in predicting the heats of formation and bond lengths of third row transition metal complexes.  Preliminary results reveal that the TPSSTPSS and TPSSKCIS density functionals are quite effective in accurately predicting the heat of formation and equilibrium bond lengths for these metal complexes, and that the B3LYP, BLYP, and correlated methods perform much more poorly for the same compounds.  These results suggest levels of computational theory that are suitable for applications such as geometry optimizations of metal-containing enzyme active sites, and for serving as the quantum mechanical segment of QM/MM calculations on larger protein structures.

Amide Hydrolysis

The exact pathway of amide hydrolysis is currently unknown. We use computational methods in order to probe this reaction mechanism. The current experimental and computational data are in some instances inconsistent with one another. We are furthering our theoretical investigations in order to elucidate the hydrolytic pathway.  We are using computational methods in order to further elucidate the various mechanistic pathways of urea hydrolysis by predicting the electronic structures of the ground states, intermediates, and transition states associated with a variety of possible mechanisms. We are using an array of density functionals to characterize these reaction pathways, with the aim of identifying a suitable QM level of theory for studying urease.


Undergraduate Research at Wellesley College

Wellesley College Logo

HIV1 Reverse Transcriptase with  Non-Nucleocide Reverse Transcriptase Inhibitor

Two renderings labeled A and B. A is on the left.
Optimizing RPVs charge distribution to bind promiscuously to multiple variants. (A) The computed optimally promiscuous charge distribution, colored by charge value according to the legend (blue = positive, red = negative). (B) Groups of atoms are colored according to their ΔΔGmoiety,additive, assuming that the drug had the computed optimally promiscuous charge distribution. Red-colored groups contribute unfavorably toward binding and blue-colored groups contribute favorably.

We conducted a systematic analysis of the electrostatic component of binding of three HIV-1 Reverse Transcriptase (RT) inhibitors—nevirapine (NVP), efavirenz (EFZ), and the recently-approved rilpivirine (RPV)—to wild-type (WT) and mutant variants of RT.

Electrostatic charge optimization was applied to determine how suited each molecule’s charge distribution is for binding WT and individual mutants of HIV-1 RT. We also determined the contributions of chemical moieties on each molecule toward the electrostatic component of binding, and show that different regions of a drug molecule may be used for recognition by different RT variants. The electrostatic contribution of certain RT residues toward drug binding was also computed to highlight critical residues for each interaction. Finally, the charge distribution of RPV is optimized to promiscuously bind to three RT variants rather than to each one in turn, with the resulting charge distribution being a compromise between the optimal charge distributions to each individual variant. This work demonstrates that even in a binding site considered quite hydrophobic, electrostatics play a subtle yet varying role that must be considered in designing next-generation molecules that recognize rapidly-mutating targets.

Trypsin/Bovine Pancreatic Trypsin Inhibitor

Rendering of the trypsin-BPTI complex
Rendering of the trypsin-BPTI complex, with trypsin in light blue and BPTI in dark red, with key residues for binding highlighted.                              

This study examines the promiscuity and specificity of the trypsin/bovine pancreatic trypsin inhibitor (BPTI) system through the analysis of the electrostatic component of ΔG. The charge contribution for each residue on the ligand is computationally eliminated and effects on ΔG are observed. Studies are also initiated for residues on the receptor (trypsin) near the binding site. A very positive change in ΔG indicates that that particular residue is important to binding. A parallel study is being conducted to evaluate the efficiency and accuracy of various numerical methods to determine ΔG electrostatic. The most accurate method, the finite difference method, is being compared with boundary element method, generalized Born, boundary-integral-based electrostatic estimation-Coulomb field approximation, and boundary-integral-based electrostatic estimation-preconditioned methods. The goal is to determine the method that optimizes accuracy and efficiency. Another undergraduate at Wellesley is currently completing this project.

Carbon Nano-tubes

Rendering of a carbon nano tube
Carbon Nano-tube

My first research project utilized semi-empirical PM3 and density functional theory (DFT) B3LYP computations in order to determine the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), and band gap energy patterns of increasing length single-walled carbon nanotubes (SWCNTs). Semi-empirical PM3 was first used to determine the band gap energy pattern. Changing the length in a carbon nanotube makes an oscillating logarithmic band gap energy pattern with a period of three points. Each point differed by one ring of carbon hexagons around the circumference. Both the HOMO and LUMO energies also create an oscillating pattern. In these patterns, at a length where HOMO energy is at a local maximum, LUMO and band gap energy is usually at a local minimum. To verify this result, the change in energy of the HOMO, LUMO, and band gap energies were graphed. It was found that in most cases these energies intercept the x-axis at the same point. Although the DFT and semi-empirical computations produced different values for the energies, a similar oscillating pattern was observed.