Advancing Molecular Simulations of Crystal Nucleation: Applications to Clathrate Hydrates

Crystallization is a fundamental physical phenomenon with broad impacts in science and engineering. Nonetheless, mechanisms of crystallization in many systems remain incompletely understood. Molecular dynamics (MD) simulations are a powerful computational technique that, in principle, are well-suited to offer insights into the mechanisms of crystallization. Unfortunately, the waiting time required to observe crystal nucleation in simulated systems often falls far beyond the limits of modern MD simulations. This rare-event problem is the primary barrier to simulation studies of crystallization in complex systems. This dissertation takes a combined approach to advance simulation studies of nucleation in complex systems. First, we apply existing tools to a challenging problem — clathrate hydrate nucleation. We then use methods development, software development, and machine learning to address the specific challenges to simulation studies of crystallization posed by the rare-event problem. Clathrate hydrate formation is an exemplar of crystallization in complex systems. Nucleation of clathrate hydrates generally occurs in systems with interfaces, and even homogeneous hydrate nucleation is inherently a multicomponent process. We address two aspects of clathrate hydrate nucleation which are not well-studied. The first aspect is the effects of interfaces on clathrate hydrate nucleation. Interfaces are common in hydrate systems, yet there are few studies probing the effects of interfaces on clathrate hydrate nucleation. We find that nucleation occurs through a homogeneous mechanism near model hydrophobic and hydrophilic surfaces. The only effect of the surfaces is through a partitioning of guest molecules which results in aggregation of guest molecules at the hydrophobic surface. The second aspect is the effect of guest solubility in water on the homogeneous nucleation mechanism. Experiments show that soluble guests act as strong promoter molecules for hydrate formation, but the molecular mechanisms of this effect are unclear. We apply forward flux sampling (FFS) and a committor analysis to identify good approximations of the reaction coordinate for homogeneous nucleation of hydrates formed from a water-soluble guest molecule. Our results suggest the possibility that the nucleation mechanism for hydrates formed from water-soluble guest molecules is different than the nucleation mechanism for hydrates formed from sparingly soluble guest molecules. FFS studies of crystal nucleation can require hundreds of thousands of individual MD simulations. For complex systems, these simulations easily generate terabytes of intermediate data. Furthermore, each simulation must be completed, analyzed, and individually processed based upon the behavior of the system. The scale of these calculations thus quickly exceeds the practical limits of traditional scripting tools (e.g., bash). In order to apply FFS to study clathrate hydrate nucleation we developed a software package, SAFFIRE. SAFFIRE automates and manages FFS with a user-friendly interface. It is compatible with any simulation software and/or analysis codes. Since SAFFIRE is built on the Hadoop framework, it easily scales to tens or hundreds of nodes. SAFFIRE can be deployed on commodity computing clusters such as the Palmetto cluster at Clemson University or XSEDE resources. Studying crystal nucleation in simulations generally requires selecting an order parameter for advanced sampling a priori. This is particularly challenging since one of the very goals of the study itself may be to elucidate the nucleation mechanism, and thus order parameters that provide a good description of the nucleation process. Furthermore, despite many strengths of FFS, it is somewhat more sensitive to the choice of order parameter than some other advanced sampling methods. To address these challenges, we develop a new method, contour forward flux sampling (cFFS), to perform FFS with multiple order parameters simultaneously. cFFS places nonlinear interfaces on-the-fly from the collective progress of the simulations, without any prior knowledge of the energy landscape or appropriate combination of order parameters. cFFS thus allows testing multiple prospective order parameters on-the-fly. Order parameters clearly play a key role in simulation studies of crystal nucleation. However, developing new order parameters is difficult and time consuming. Using ideas from computer vision, we adapt a specific type of neural network called a PointNet to identify local structural environments (e.g., crystalline environments) in molecular simulations. Our approach requires no system-specific feature engineering and operates on the raw output of the simulations, i.e., atomic positions. We demonstrate the method on crystal structure identification in Lennard-Jones, water, and mesophase systems. The method can even predict the crystal phases of atoms near external interfaces. We demonstrate the versatility of our approach by using our method to identify surface hydrophobicity based solely upon positions and orientations of nearby water molecules. Our results suggest the approach will be broadly applicable to many types of local structure in simulations. We address several interdependent challenges to studying crystallization in molecular simulations by combining software development, method development, and machine learning. While motivated by specific challenges identified during studies of clathrate hydrate nucleation, these contributions help extend the applicability of molecular simulations to crystal nucleation in a broad variety of systems. The next step of the development cycle is to apply these methods on complex systems to motivate further improvements. We believe that continued integration of software, methods, and machine learning will prove a fruitful framework for improving molecular simulations of crystal nucleation

Methods for increasing model accuracy and simulation time scales of biological processes with molecular dynamics

This dissertation presents three research projects on novel methods in computational bio- physics. Each of these projects introduces methodologies to extend the capabilities of molecular dynamics simulations in one way or another. In the first chapter, molecular dynamics simulations and the central role they play in the field of structural biology is introduced to give the reader some background on the common basis of the projects. The second chapter describes the first of these projects, where the molecular dynamics flexible fitting method for refining molecular structures of macromolecules using experimental electron density data is extended to be able to handle high-resolution density data, which are becoming increasingly commonplace. The third chapter focuses on adaptive multilevel splitting, a replica-based sampling technique that was employed in molecular dynamics simulations to measure the rate of drug molecule dissociation, a process that occurs on the order of milliseconds and above, which is out of the reach of typical molecular dynamics simulations. In the final chapter, a kinetic model of diffusion is introduced. This model allows simulation of the diffusion of small molecules in arbitrary potentials, for example, those that characterize the space around and within a membrane protein channel. The adaptive discretization scheme allows simulations between the micro- to millisecond time scales, which are typical of diffusive processes. This collection of projects is a snapshot of the diversity and versatility of current problems in structural biology that can be addressed by molecular dynamics simulations. I hope to instill in the reader a sense of how method development in molecular dynamics will expand the contributions of the field to both scientific and practical pursuits in biology

Computational studies of barrier-crossing in polymer field theory

This dissertation is primarily a survey of the zero-temperature string method, a minimum energy path search algorithm, applied to novel barrier-crossing problems in polymer field theory. I apply the method to both self-consistent field theory (SCFT) and a phase-field model (the Landau-Brazovskii model). In the case of SCFT, the focus is on defect annealing problems in homo+copolymer melts; in the case of the Landau-Brazovskii model, the focus is on finding critical nuclei for the disorder-to-lamellar transition, which is known to be a fluctuation-induced first-order phase transition.In SCFT, applying the string method is computationally demanding in both processing time and memory, especially for fully 3-dimensional simulations at industrially relevant system sizes. I successfully address these challenges on state-of-the-art massively parallel computing architectures (NVIDIA graphics processing units). As a result our group is able to identify free energy barriers and transition mechanisms for a wide variety of defect annealing problems relevant to industrial directed self-assembly (DSA).Nucleation in the Landau-Brazovskii model presents its own challenges. The string method as originally formulated is inefficient for nucleation problems, since many images are wasted tracing out unphysical configurations once the nucleus grows to the edges of the simulation cell. I devise a new truncation-based energy weighting (TBEW) scheme that resolves this issue, and will prove valuable to future researchers using the string method to find critical nuclei.Since the bare Landau-Brazovskii model predicts a second-order transition between disorder and lamellae at a mean-field level, naive application of the zero-temperature string method to this model fails to find a barrier. To circumvent this, I instead apply the string method to a renormalized model that incorporates fluctuations at a mean-field level. Using TBEW and the renormalized model, I investigate nucleation pathways for the disorder-to-lamellar transition, finding anisotropic nuclei in agreement with previous predictions and experimental observations. I also conduct a comprehensive search for experimentally observed nuclei containing various exotic defect structures.Finally, I evaluate the validity of the nucleation pathways obtained from the renormalized model by numerically simulating the bare model with explicit fluctuations. I find that the renormalized model makes good predictions for certain quantities, including the location of the order-disorder transition. However, due to sharp dependence of critical nucleus size on proximity to the order-disorder transition, even slight errors in the predicted ODT lead to large errors in predicted nucleus size. I conclude that the renormalized Landau-Brazovskii model is a poor tool for predicting critical nuclei in the fully fluctuating bare theory at experimentally accessible parameters, and recommend that future studies work with the fluctuating bare theory directly. I recommend several strategies to extract barriers and rates

Physics of Ionic Conduction in Narrow Biological and Artificial Channels

The book reprints a set of important scientific papers applying physics and mathematics to address the problem of selective ionic conduction in narrow water-filled channels and pores. It is a long-standing problem, and an extremely important one. Life in all its forms depends on ion channels and, furthermore, the technological applications of artificial ion channels are already widespread and growing rapidly. They include desalination, DNA sequencing, energy harvesting, molecular sensors, fuel cells, batteries, personalised medicine, and drug design. Further applications are to be anticipated.The book will be helpful to researchers and technologists already working in the area, or planning to enter it. It gives detailed descriptions of a diversity of modern approaches, and shows how they can be particularly effective and mutually reinforcing when used together. It not only provides a snapshot of current cutting-edge scientific activity in the area, but also offers indications of how the subject is likely to evolve in the future

Studying protein-ligand interactions using a Monte Carlo procedure

[eng] Biomolecular simulations have been widely used in the study of protein-ligand interactions; comprehending the mechanisms involved in the prediction of binding affinities would have a significant repercussion in the pharmaceutical industry. Notwithstanding the intrinsic difficulty of sampling the phase space, hardware and methodological developments make computer simulations a promising candidate in the resolution of biophysically relevant problems. In this context, the objective of the thesis is the development of a protocol that permits studying protein-ligand interactions, in view to be applied in drug discovery pipelines. The author contributed to the rewriting PELE, our Monte Carlo sampling procedure, using good practices of software development. These involved testing, improving the readability, modularity, encapsulation, maintenance and version control, just to name a few. Importantly, the recoding resulted in a competitive cutting-edge software that is able to integrate new algorithms and platforms, such as new force fields or a graphical user interface, while being reliable and efficient. The rest of the thesis is built upon this development. At this point, we established a protocol of unbiased all-atom simulations using PELE, often combined with Markov (state) Models (MSM) to characterize the energy landscape exploration. In the thesis, we have shown that PELE is a suitable tool to map complex mechanisms in an accurate and efficient manner. For example, we successfully conducted studies of ligand migration in prolyl oligopeptidases and nuclear hormone receptors (NHRs). Using PELE, we could map the ligand migration and binding pathway in such complex systems in less than 48 hours. On the other hand, with this technique we often run batches of 100s of simulations to reduce the wall-clock time. MSM is a useful technique to join these independent simulations in a unique statistical model, as individual trajectories only need to characterize the energy landscape locally, and the global characterization can be extracted from the model. We successfully applied the combination of these two methodologies to quantify binding mechanisms and estimate the binding free energy in systems involving NHRs and tyorsinases. However, this technique represents a significant computational effort. To reduce the computational load, we developed a new methodology to overcome the sampling limitations caused by the ruggedness of the energy landscape. In particular, we used a procedure of iterative simulations with adaptive spawning points based on reinforcement learning ideas. This permits sampling binding mechanisms at a fraction of the cost, and represents a speedup of an order of magnitude in complex systems. Importantly, we show in a proof-of-concept that it can be used to estimate absolute binding free energies. Overall, we hope that the methodologies presented herein help streamline the drug design process.[spa] Las simulaciones biomoleculares se han usado ampliamente en el estudio de interacciones proteína-ligando. Comprender los mecanismos involucrados en la predicción de afinidades de unión tiene una gran repercusión en la industria farmacéutica. A pesar de las dificultades intrínsecas en el muestreo del espacio de fases, mejoras de hardware y metodológicas hacen de las simulaciones por ordenador un candidato prometedor en la resolución de problemas biofísicos con alta relevancia. En este contexto, el objetivo de la tesis es el desarrollo de un protocolo que introduce un estudio más eficiente de las interacciones proteína-ligando, con vistas a diseminar PELE, un procedimiento de muestreo de Monte Carlo, en el diseño de fármacos. Nuestro principal foco ha sido sobrepasar las limitaciones de muestreo causadas por la rugosidad del paisaje de energías, aplicando nuestro protocolo para hacer analsis detallados a nivel atomístico en receptores nucleares de hormonas, receptores acoplados a proteínas G, tirosinasas y prolil oligopeptidasas, en colaboración con una compañía farmacéutica y de varios laboratorios experimentales. Con todo ello, esperamos que las metodologías presentadas en esta tesis ayuden a mejorar el diseño de fármacos

EXAMINING PROTEIN CONFORMATIONAL DYNAMICS USING COMPUTATIONAL TECHNIQUES: STUDIES ON PHOSPHATIDYLINOSITOL-3-KINASE AND THE SODIUM-IODIDE SYMPORTER

Experimental biophysics techniques used to study proteins, polymers of amino acids that comprise most therapeutic targets of human disease, face limitations in their ability to interrogate the continual structural fluctuations exhibited by these macromolecules in the context of their myriad cellular functions. This dissertation aims to illustrate case studies that demonstrate how protein conformational dynamics can be characterized using computational methods, yielding novel insights into their functional regulation and activity. Towards this end, the work presented here describes two specific membrane proteins of therapeutic relevance: Phosphoinositide 3-kinase (PI3Kα), and the Na+/I- symporter (NIS). The PI3KCA gene, encoding the catalytic subunit of the PI3Kα protein that phosphorylates phosphatidylinositol-4,5-bisphosphate (PIP2) to generate phosphatidylinositol-3,4,5-triphosphate (PIP3), is highly mutated in human cancer. As such, a deeper mechanistic understanding of PI3Kα could facilitate the development of novel chemotherapeutic approaches. The second chapter of this dissertation describes molecular dynamics (MD) simulations that were conducted to determine how PI3Kα conformations are influenced by physiological effectors and the nSH2 domain of a regulatory subunit, p85. The results reported here suggest that dynamic allostery plays a role in populating the catalytically competent conformation of PI3Kα. NIS, a thirteen-helix transmembrane protein found in the thyroid and other tissues, transports iodide, a required constituent of thyroid hormones T3 and T4. Despite extensive experimental information and clinical data, many mechanistic details about NIS remain unresolved. The third chapter of this dissertation describes the results of unbiased and enhanced-sampling MD simulations of inwardly and outwardly open models of bound NIS under an enforced ion gradient. Simulations of NIS in the absence or presence of perchlorate are also described. The work presented in this dissertation aims to add to our mechanistic understanding of NIS ion transport and elucidate conformational states that occur between the inward and outward transitions of NIS in the absence and presence of bound Na+ and I- ions, which can provide valuable insight into its physiological activity and inform therapeutic interventions. Taken together, these case studies demonstrate the ability of computational techniques to provide novel insights into the impact of structural dynamics on the functional regulation of therapeutically important biological macromolecules

Physics of Ionic Conduction in Narrow Biological and Artificial Channels

This is a book about ion channels. It has been written mostly by physical scientists and mathematicians, even though the most widespread and important manifestation of ion channels is in biology, where they are essential to life in all its forms. How do non-biologists get involved in such investigations? Everyone will have their own particular story but, for ourselves, it was the heady combination of scientific curiosity, a wish to contribute to the fundamental understanding of natural phenomena that clearly have crucially important applications, and the realisation that some of our physics knowledge and expertise might be relevant

Computational Approaches to Simulation and Analysis of Large Conformational Transitions in Proteins

abstract: In a typical living cell, millions to billions of proteins—nanomachines that fluctuate and cycle among many conformational states—convert available free energy into mechanochemical work. A fundamental goal of biophysics is to ascertain how 3D protein structures encode specific functions, such as catalyzing chemical reactions or transporting nutrients into a cell. Protein dynamics span femtosecond timescales (i.e., covalent bond oscillations) to large conformational transition timescales in, and beyond, the millisecond regime (e.g., glucose transport across a phospholipid bilayer). Actual transition events are fast but rare, occurring orders of magnitude faster than typical metastable equilibrium waiting times. Equilibrium molecular dynamics (EqMD) can capture atomistic detail and solute-solvent interactions, but even microseconds of sampling attainable nowadays still falls orders of magnitude short of transition timescales, especially for large systems, rendering observations of such "rare events" difficult or effectively impossible. Advanced path-sampling methods exploit reduced physical models or biasing to produce plausible transitions while balancing accuracy and efficiency, but quantifying their accuracy relative to other numerical and experimental data has been challenging. Indeed, new horizons in elucidating protein function necessitate that present methodologies be revised to more seamlessly and quantitatively integrate a spectrum of methods, both numerical and experimental. In this dissertation, experimental and computational methods are put into perspective using the enzyme adenylate kinase (AdK) as an illustrative example. We introduce Path Similarity Analysis (PSA)—an integrative computational framework developed to quantify transition path similarity. PSA not only reliably distinguished AdK transitions by the originating method, but also traced pathway differences between two methods back to charge-charge interactions (neglected by the stereochemical model, but not the all-atom force field) in several conserved salt bridges. Cryo-electron microscopy maps of the transporter Bor1p are directly incorporated into EqMD simulations using MD flexible fitting to produce viable structural models and infer a plausible transport mechanism. Conforming to the theme of integration, a short compendium of an exploratory project—developing a hybrid atomistic-continuum method—is presented, including initial results and a novel fluctuating hydrodynamics model and corresponding numerical code.Dissertation/ThesisDoctoral Dissertation Physics 201

Path Reweighting Methods for underdamped Langevin Dynamics for Molecular Systems

Knowledge about the dynamical properties of biomolecules is essential to understand their function in biological processes. This thesis approaches the task to compute dynamical properties with two different strategies. Part A focuses on Molecular Dynamics (MD) simulations combined with path reweighting. Three of the most widely used underdamped Langevin integrators for MD simulations are the splitting methods BAOAB and BAOA which are available in the MD packages OpenMM and AMBER and the Gromacs Stochastic Dynamics (GSD) integrator implemented in GROMACS. We found that all three integrators are equivalent configurational sampling algorithms and thus yield configurational properties at equivalent accuracy. MD simulations with stochastic integrators such as Langevin integrators offer the possibility to reweight estimated dynamical properties using path reweighting. With path reweighting we can for example recover the original dynamics from MD simulation that have been conducted with enhanced sampling methods. The key component of path reweighting is the path reweighting factor M which strongly depends on the chosen integrator. We derive M_L for underdamped Langevin dynamics propagated by a variant of the Langevin Leapfrog integrator. Additionally, we present two strategies which can be used as blueprints to straightforwardly derive M_L for other Langevin integrators. The previously reported path reweighting factor matches the Euler-Maruyama integrator for overdamped Langevin dynamics and was used as standard reweighting factor even though the MD simulation was conducted with an underdamped Langevin integrator. We prove that this path reweighting factors differs from the exact M_L only by O(ξ^4 ∆t^4) and thus yields highly accurate dynamical reweighting results (∆t is the integration time step, and ξ is the collision rate.). Part B of this thesis combines experimental and theoretical approaches to investigate Multiple Inositol Polyphosphate Phosphatase 1 (MINPP1)-mediated inositol polyphosphate (InsP) networks. We use 13C-labeling experiments combined with nuclear magnetic resonance spectroscopy (NMR) to uncover a novel branch of InsP dephosphorylation in human cells. Additionally, we extract the corresponding reaction rates using a Markovian kinetic scheme as theoretical model to describe the network

Applications of Molecular Dynamics simulations for biomolecular systems and improvements to density-based clustering in the analysis

Molecular Dynamics simulations provide a powerful tool to study biomolecular systems with atomistic detail. The key to better understand the function and behaviour of these molecules can often be found in their structural variability. Simulations can help to expose this information that is otherwise experimentally hard or impossible to attain. This work covers two application examples for which a sampling and a characterisation of the conformational ensemble could reveal the structural basis to answer a topical research question. For the fungal toxin phalloidin—a small bicyclic peptide—observed product ratios in different cyclisation reactions could be rationalised by assessing the conformational pre-organisation of precursor fragments. For the C-type lectin receptor langerin, conformational changes induced by different side-chain protonations could deliver an explanation of the pH-dependency in the protein’s calcium-binding. The investigations were accompanied by the continued development of a density-based clustering protocol into a respective software package, which is generally well applicable for the use case of extracting conformational states from Molecular Dynamics data