Search

Molecular dynamics simulations have become an indispensable tool in the study of a wide range of problems from chemistry to physics and biology. However, the presence of kinetic bottlenecks in many problems of interest in conjunction with the timescale limitations of molecular dynamics reduces the scope of such applications. This has catalyzed the development of several enhanced sampling techniques that help to increase the rate of transitions between states of interest. Despite recent advances, many issues remain. In this thesis, we aim to resolve two of these issues. The first concern that we address here is the insufficient sampling of high energy configurations. To rectify this, we utilized the Variationally Enhanced Sampling technique and devised a way to tune the amount of sampling in the different regions of the free energy surface. Besides increasing the configurations visited in the high energy transition regions, this method helps to accelerate the convergence of free energy surfaces. Another problem that is currently present in enhanced sampling simulations is a protocol to accurately predict unbiased dynamics from biased trajectories. In this regard, we propose a new enhanced sampling technique called Gaussian Mixture Based Enhanced Sampling (GAMBES), where a bias potential is constructed using local probability distributions of known metastable states. This bias potential does not act on transition regions between metastable states and only facilitates the transition of a system out of these known states. The trajectories thus generated, can be used not only to calculate static properties, but can also be harvested to calculate transition rates or to understand transition mechanisms. Finally, using GAMBES we examine the role of water in unbinding dynamics in a set of host-guest systems.

The free energy of a process is the fundamental quantity that determines its spontaneity or propensity at a given temperature. In particular, the binding free energy of a drug candidate to its biomolecular target is used as an objective quantity in drug design. Recently, binding kinetics—rates of association (kon) and dissociation (koff)—have also demonstrated utility for their ability to predict efficacy and in some cases have been shown to be more predictive than the binding free energy alone. Some methods exist to calculate binding kinetics from molecular simulations, although these are typically more difficult to calculate than the binding affinity as they depend on details of the transition path ensemble. Assessing these rate constants can be difficult, due to uncertainty in the definition of the bound and unbound states, large error bars and the lack of experimental data. As an additional consistency check, rate constants from simulation can be used to calculate free energies (using the log of their ratio) which can then be compared to free energies obtained experimentally or using alchemical free energy perturbation. However, in this calculation it is not straightforward to account for common, practical details such as the finite simulation volume or the particular definition of the “bound” and “unbound” states. Here we derive a set of correction terms that can be applied to calculations of binding free energies using full reactive trajectories. We apply these correction terms to revisit the calculation of binding free energies from rate constants for a host-guest system that was part of a blind prediction challenge, where significant deviations were observed between free energies calculated with rate ratios and those calculated from alchemical perturbation. The correction terms combine to significantly decrease the error with respect to computational benchmarks, from 3.4 to 0.76 kcal/mol. Although these terms were derived with weighted ensemble simulations in mind, some of the correction terms are generally applicable to free energies calculated using physical pathways via methods such as Markov state modeling, metadynamics, milestoning, or umbrella sampling.

Molecular dynamics simulations represent a powerful tool to gain insights into structural and dynamical features of biomolecular systems. Nevertheless, their recognized limitation in terms of achievable timescales becomes particularly severe when dealing with slow processes. In such cases, the employment of enhanced sampling methods, which allow accelerating the characterization of rare events in a timeframe consistent with conventional computational resources, results as crucial. In particular, such advanced techniques have proven highly valuable in the context of protein folding and, specifically, to explore the conformational ensemble spanned by intrinsically disordered proteins (IDPs). Here, we describe how to set up molecular dynamics simulations with one of these enhanced sampling approaches (namely, Parallel Tempering Metadynamics in the Well-Tempered Ensemble) using the NTAIL peptide as a test case.

Allergies are a widespread problem in western countries, affecting a large part of the population, with levels of prevalence increasingly rising due to reasons still not understood. Evidence accumulated in recent years points to an essential role played by ligands of allergen proteins in the sensitization phase of allergies. In this regard, we recently identified the natural ligand of Pru p 3, a lipid transfer protein, a major allergen from peach fruit and a model of food allergy. The ligand of Pru p 3 has been shown to play a key role in the sensitization to peach and to other plant food sources that provoke cross-reactivity in a large proportion of patients allergic to peach. However, the question of which is the binding pose of this ligand in its carrier protein, and how it can be transferred to receptors of the immune system where it develops its function as a coadjuvant was not elucidated. In this work, different molecular dynamics simulations have been considered as starting points to study the properties of the ligand&ndash, protein system in solution. Besides, an energy landscape based on collective variables that describe the process of ligand motion within the cavity of Pru p 3 was obtained by using well-tempered metadynamics. The simulations revealed the differences between distinct binding modes, and also revealed important aspects of the motion of the ligand throughout its carrier protein, relevant to its binding&ndash, unbinding process. Our findings are potentially interesting for studying protein&ndash, ligand systems beyond the specific case of the allergen protein dealt with here.

This thesis focuses on studies of protein-ligand systems using enhanced sampling methods. In chapter I, I give a brief introduction to the time-scale problem and some enhanced sampling methods. In chapter II, the basics of MD simulation are reviewed. In chapter III, the theoretical backgrounds of umbrella sampling, bias-exchange metadynamics and infrequent metadynamics are presented. In chapter IV, the 5 papers included in this thesis are summarized. In paper 1, we studied the relationship between the antibacterial activities of antimicrobial peptides and their aggregation propensities. We found that an increasing aggregation propensity increases the free energy cost of peptide embedding into the bacterial membrane and decreases antibacterial activity. In paper 2, we employed the umbrella sampling approach to obtain the free energy landscape of Pittsburgh compound-B penetrating into the core binding sites of amyloid βfibrils. Our study suggested that, for the design of probes binding to fibril like proteins, other than the binding affinity, the dynamics of probes in the fibrils should also be considered. In paper 3, we studied the coupled folding and binding process of the intrinsically disordered protein p53 to MDM2 with bias-exchange metadynamics and infrequent metadynamics. We reconstructed the free energy landscape and built a kinetic network for this process. In paper 4, we studied the binding modes of ASEM with a chimera structure of α7 nicotinic acetylcholine receptor with well-tempered metadynamics. We found that an important residue, Trp53, can significantly affect the stabilities of the binding modes. In paper 5, we proposed an efficient method to estimate the transition times of rare events in biomolecular systems. In chapter V, I present a conclusion of this thesis and propose an outlook related to the selection of collective variables for enhanced sampling methods. QC 2019-05-10

Ligand-protein interactions play a key role in protein functions. Ligand binding to biological targets can be effectively studied using molecular dynamics (MD) simulations with enhanced sampling technique to explore the long time-scales of such processes [1]⁠. The study-case here investigated is the Hypoxia Inducible Factor-2α (HIF-2α), a transcription factor that responds to decreases in available oxygen in the cellular environment (hypoxia). A domain of this protein (the HIF-2α PAS-B) contains a preformed cavity inaccessible to the solvent and able to bind artificial ligands which role as potential anticancer drugs is attracting great interest. Several crystallographic structures of HIF-2α in complex with different ligands have become available, and both thermodynamic and kinetic data associated to the binding/unbinding process of some ligands have been experimentally determined [2]. The aim of this study is to get an accurate mechanistic description of ligand binding to HIF-2α using computational methods able to simulate the dynamics of the process. Among the MD-based enhanced sampling techniques, we tested and compared: Steered MD [3] to identify the exit paths of the ligands from the pocket of HIF-2α PAS-B; Umbrella Sampling [4]⁠ to predict the free-energy profiles on the previously identified pathways; Path Collective Variables [5] to optimize the paths and obtain a more accurate description of the process. After an initial test of the different methods on a single ligand, the study will be extended to all the HIF-2α ligands for which experimental information is available to set up a reliable computational protocol for predictive studies. [1] Rydzewski, J. & Nowak, W., Phys. Life Rev. 22–23, 58–74 (2017) [2] Key, J., Scheuermann, T. H., Anderson, P. C., Daggett, V. & Gardner, K. H., J. Am. Chem. Soc. 131, 17647– 54 (2009) [3] Patel, J. S., Berteotti, A., Ronsisvalle, S., Rocchia, W. & Cavalli, A., J. Chem. Inf. Model. 54, 470–480 (2014) [4] Kästner, J., Wiley Interdiscip. Rev. Comput. Mol. Sci. 1, 932–942 (2011) [5] Tiwary, P., Limongelli, V., Salvalaglio, M. & Parrinello, M., Proc. Natl. Acad. Sci. 112, E386–E391 (2015)

In this thesis I focused my efforts in developing and applying computational approaches aimed to characterize the complex phenomena of molecular recognition, including both kinetic and thermodynamic aspects. In details, I worked on i) prediction of the correct binding mode of small molecules within biological targets of pharmaceutical interest, ii) prediction of the kinetics of protein-ligand unbinding and iii) estimation of covalent drug affinity and reactivity in mutated forms of the targeted protein. Procedures based on the application of enhanced sampling methods are described and employed to solve these challenging tasks.

Die archaeale Protonenpumpe Bacteriorhodopsin nutzt Sonnenlicht um einen Protonengradienten über der Zellmembran aufzubauen. Der Mechanismus der Energieumwandlung basiert auf einem gut erforschten Photozyklus. Obwohl schon über 45 Jahre an diesem bioenergetischen Modellsystem geforscht wurde, sind noch immer wenig Informationen über den letzten Schritt dieses Photozykluses vorhanden: den Übergang vom O zum Grundzustand (bR). Es wird jedoch angenommen, dass ein langreichweitiger Protonentransfer den bR Zustand wiederherstellt. Diese Arbeit zielt darauf ab, weiteren Aufschluss über die mechanistischen Details und thermody- namischen/kinetischen Eigenschaften des O→bR Übergangs zu geben. Um dies zu bewerkstelligen, wurden computergestütze Rechenmethoden als brauchbare Ergänzung zum Experiment verwendet, da die Handhabung von Membranproteinen und die gesetzten Ziele für experimentelle Methoden eine schwierige Aufgabe darstellen. Zunächst wurden Strukturmodelle mit Hilfe von enhanced sampling Methoden für den O, bR und den O* Intermediatszustand erstellt. Für die Simulation der (langreichweitigen) Protonentransfer- reaktionen wurde ein hybrides QM/MM Set-up basierend auf der semiempirischen Quantenchemie- Methode DFTB3 und dem CHARMM36 Kraftfeld verwendet. Die mit dem Protonentransfer asso- ziierte Änderung der Gibbs’schen Freien Energie wurde mit Hilfe von Freie-Energie-Techniken in Kombination mit einer fortgeschrittenen Reaktionskoordinate, die auf der center of excess charge- Darstellung basiert, aufgeklärt. Die Ergebnisse weisen darauf hin, dass es sich bei der O→bR Konversion um einen leicht exer- gonen Prozess handelt, in dem der O* Zustand ein metastabiler Zwischenzustand bildet. Bezüglich des Protonentransfer-Pfads wurde eine interessante Feststellung gemacht: Der langreichweitige Pro- tonentransfer wird über einen Protonen-Loch-/inversen Grotthuss-Mechanismus vollzogen. Zudem wird dieser Ladungstransfer begleitet von der Umorientierung einer funktionell wichtigen, positiv geladenen Arginin Seitenkette sowie der Ausbildung der proton release group. Im Zuge dieser Studien, wurde des Weiteren die Genauigkeit von DFTB3 für Protonentransfer- reaktionen evaluiert. Die Methode bewies sich als leistungsstark für die Beschreibung dieser Reak- tionen mit einer geringfügigen Tendenz für die Unterschätzung von Reaktionsbarrieren. Um diesen Fehler zu korrigieren, der aus der Unterschätzung der kurzreichweitigen Pauli Repulsion in DFTB3 resultiert, wurde die empirische Delta-Pauli Korrektur entwickelt. Allerdings gab es bisher keine Pa- rameter für dieses Model. Daher wurde in dieser Arbeit der erste Parametersatz für CHNO-basierte Molekülsysteme entworfen und an relevanten organischen/biochemischen Systemen evaluiert. Die Ergebnisse zeigen, dass Delta-Pauli erfolgreich kurzreichweitige Pauli Repulsion in DFTB3 einführt. Dies hat zur Folge, dass Delta-Pauli die Beschreibung von Molekülsystemen bezüglich nichtkovalenter Wechselwirkungen, inter- und intramolekularen Reaktionsbarrieren und Gleichge- wichtsgeometrien verbessert.

We introduce a computational method to discover polymorphs in molecular crystals at finite temperature. The method is based on reproducing the crystallization process starting from the liquid and letting the system discover the relevant polymorphs. This idea, however, conflicts with the fact that crystallization has a time scale much longer than that of molecular simulations. In order to bring the process within affordable simulation time, we enhance the fluctuations of a collective variable by constructing a bias potential with well tempered metadynamics. We use as collective variable an entropy surrogate based on an extended pair correlation function that includes the correlation between the orientation of pairs of molecules. We also propose a similarity metric between configurations based on the extended pair correlation function and a generalized Kullback-Leibler divergence. In this way, we automatically classify the configurations as belonging to a given polymorph using our metric and a hierarchical clustering algorithm. We find all relevant polymorphs for both substances and we predict new polymorphs. One of them is stabilized at finite temperature by entropic effects., Comment: 7 pages, 4 figures

Molecular dynamics simulations have played a fundamental role in numerous fields of science by providing insights into the structure and dynamics of complex systems at the atomistic level. However, exhaustive sampling by standard molecular dynamics is in most cases computationally prohibitive, and the time scales accessible remain significantly shorter than many biological processes of interest. In particular, in the study of ion channels, realistic models to describe permeation and gating require accounting for large numbers of particles and accurate interaction potentials, which severely limits the length of the simulations. To overcome such limitations, several advanced methods have been proposed among which is metadynamics. In this algorithm, an external bias potential to accelerate sampling along selected collective variables is introduced. This bias potential discourages visiting regions of the configurational space already explored. In addition, the bias potential provides an estimate of the free energy as a function of the collective variables chosen once the simulation has converged. In this review, recent contributions of metadynamics to the field of ion channels are discussed, including how metadynamics has been used to search for transition states, predict permeation pathways, treat conformational flexibility that underlies the coupling between gating and permeation, or compute free energy of permeation profiles. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

In the realm of protein-protein interactions, the assembly process of homooligomers plays a fundamental role because the majority of proteins fall into this category. A comprehensive understanding of this multistep process requires the characterization of the driving molecular interactions and the transient intermediate species. The latter are often short-lived and thus remain elusive to most experimental investigations. Molecular simulations provide a unique tool to shed light onto these complex processes complementing experimental data. Here we combine advanced sampling techniques, such as metadynamics and parallel tempering, to characterize the oligomerization landscape of fibritin foldon domain. This system is an evolutionarily optimized trimerization motif that represents an ideal model for experimental and computational mechanistic studies. Our results are fully consistent with previous experimental nuclear magnetic resonance and kinetic data, but they provide a unique insight into fibritin foldon assembly. In particular, our simulations unveil the role of nonspecific interactions and suggest that an interplay between thermodynamic bias toward native structure and residual conformational disorder may provide a kinetic advantage.

The use of free-energy landscapes rationalizes a wide range of aspects of protein behavior by providing a clear illustration of the different states accessible to these molecules, as well as of their populations and pathways of interconversion. The determination of the free-energy landscapes of proteins by computational methods is, however, very challenging as it requires an extensive sampling of their conformational spaces. We describe here a technique to achieve this goal with relatively limited computational resources by incorporating nuclear magnetic resonance (NMR) chemical shifts as collective variables in metadynamics simulations. As in this approach the chemical shifts are not used as structural restraints, the resulting free-energy landscapes correspond to the force fields used in the simulations. We illustrate this approach in the case of the third Ig-binding domain of protein G from streptococcal bacteria (GB3). Our calculations reveal the existence of a folding intermediate of GB3 with nonnative structural elements. Furthermore, the availability of the free-energy landscape enables the folding mechanism of GB3 to be elucidated by analyzing the conformational ensembles corresponding to the native, intermediate, and unfolded states, as well as the transition states between them. Taken together, these results show that, by incorporating experimental data as collective variables in metadynamics simulations, it is possible to enhance the sampling efficiency by two or more orders of magnitude with respect to standard molecular dynamics simulations, and thus to estimate free-energy differences among the different states of a protein with a k B T accuracy by generating trajectories of just a few microseconds.

This thesis is devoted to the study of activated processes at the nanoscale by means of state-of-the-art computer simulations. On the one hand, we are interested in addressing specific aspects of chemico-physical processes of great practical importance, namely the early steps of wet carbon dioxide capture, and the diffusion of an impurity in cadmium telluride. On the other hand, we assess the limitations of current methodologies and propose new tools and protocols to extend their range of applicability. The simulation of the dynamics of activated processes at the atomistic level is a delicate issue facing two major difficulties: (i) an accurate description of the changes happening in the electronic structure has to be obtained ; and (ii) a statistically meaningful portion of the configuration space has to be explored, including scarcely populated regions like transition states. By employing finite temperature ab initio molecular dynamics empowered by metadynamics, it is possible to successfully overcome these challenges. We apply these methods and develop general robust protocols for these simulations and for the analysis of the results. Applications refer to CO2 chemical reactions, which are relevant to the science of carbon capture, and sulfur diffusion in CdTe, a photovoltaic material of current great interest. An unprecedented characterization of these processes is obtained, leading to new understanding of the experimentally observed behavior, well beyond the reach of state-of-the-art calculations. For example, our simulations disclose the active role played by water in the capture of carbon dioxide, and reveal the easiness by which sulfur is trapped in a CdTe lattice.