Using umbrella integration to find minimum free energy paths

Abstract

The analysis of the detailed mechanism of chemical reactions is a key task of computational chemistry. The detailed knowledge may help to improve known processes or even contribute to the development of new ones. In this way ecological and economic demands can be reduced. Furthermore, the course of a reaction path also plays an important role in the action of drugs. Understanding the binding of the active ingredient to receptors, such as proteins, can make it possible to optimize drugs by reducing side effects, or even to find new effects. Unfortunately, chemical reactions are usually too fast to observe intermediate states by experimental methods. This is where theoretical chemistry joins the game.

In theoretical chemistry we combine the coordinates of all atoms with the configuration space. The potential energy forms a hypersurface in this space. Minima represent stable or metastable states. Saddle points represent transition states which are the most unfavourable configurations occurring on the most favourable path between two minima. If the path of a reaction is known, all intermediate states can be observed by theoretical methods. However, these calculations are usually computationally costly. This is the reason why, in contrast to experimental methods, it is in general impossible to sample the whole configuration space. This would in most cases exceed the available computational resources. It is therefore necessary to use techniques enabling a reaction path to be found without sampling the whole configuration space.

In the case of a thermodynamic ensemble, e.g. the contents of a test tube, statistical information has to be included. The corresponding potential is the so-called free energy landscape, for which some degrees of freedom of the configuration space are thermostatistically integrated out. Consequently this function includes statistical and energetic properties. The free energy can in general only be calculated by statistical simulations (Monte--Carlo or molecular dynamics). Unfortunately, the transition states in which we have a special interest are rarely sampled. Special methods have to be applied to get sufficient sampling as well in areas of these rare events.

In this work a non-physical quadratic potential is used to bias the equation of motion of the particles while doing molecular dynamics simulations. In this way unfavourable areas in the configuration space can also be sampled sufficiently. This technique is called umbrella sampling. The bias is applied to one or more coordinates which describe the reaction and therefore are called reaction coordinates. An umbrella sampling run will result in a distribution of the reaction coordinates. The expectation value of this distribution will be located close to the minimum of the bias function without in general corresponding to it. The difference between the minimum of the bias function and the expectation value can be used to calculate the gradient of the underlying free energy surface. Similarly, the covariance of the distribution of the reaction coordinate can be used to calculate the Hessian of the free energy surface. This method of interpreting the data gained by umbrella sampling is called umbrella integration.

At first these values are used for an iterative search of the saddle points, which represent the transition states. From these configurations, free energy paths can be constructed by following the gradient down to the minima. This algorithm was successfully tested for the alanine dipeptide system.

This simple method has the disadvantages that it works serially and that one needs good initial guesses for the saddle points in order to find them. Therefore, in a second part of the work, the established method of nudged elastic band optimization (NEB) is extended for use in the free energy surface. NEB optimization searches for a reaction path. This path is discretized into a number configurations, so-called images. For the sake of equal distribution of the images along the path a non-physical spring force between the images is used. The force, which is actually minimized during the NEB optimisation, consists of the projection of the real force of the underlying potential perpendicular to the path, and the projection of the spring force parallel to the path.

An optimizer is developed which archives quadratic convergence of NEB optimizations in the noise-free potential energy surface of some test systems. This optimiser uses gradients and Hessians at each step. For the free energy surface both values can be calculated by umbrella integration as mentioned above. NEB optimizations within the free energy are performed in this work in the following way: at first a guess path is assumed, e.g. a straight line between two points in the configuration space, usually minima. This path is discretized into a number of images. Molecular dynamics umbrella sampling simulations are performed on each image. The gradient and Hessian from the umbrella integration are fed into the newly developed NEB optimizer. This way one does not need good starting guesses for the saddle points but an interpolation between the much more easily accessible minima is sufficient. Furthermore, the need for independent molecular dynamics runs at each image makes the method intrinsically parallel. The whole method is applied to the well-studied alanine dipeptide system and compared with the results from the serial method. Subsequently the algorithm is applied to a much more costly system of binding a ligand to its receptor in water.

Die genaue Analyse von chemischen Reaktionen ist eine der Schlüsselaufgaben in der theoretischen Chemie. Das Wissen, das dabei gewonnen wird, hilft bekannte Prozesse zu verbessern oder neue Prozesse zu entwickeln und somit die ökologischen und ökonomischen Kosten zu reduzieren. Darüber hinaus spielt der Verlauf von Reaktionspfaden auch eine entscheidende Rolle für die Wirkung von Medikamenten. Durch das genaue Verständnis der Bindungsvorgänge der Wirkstoffe mit Proteinen können Wirkungen optimiert, Nebenwirkungen reduziert oder sogar neue Wirkmechanismen gefunden werden. Leider laufen die meisten Reaktion mit hohen Geschwindigkeiten ab, so dass es oft nicht möglich ist, genaue Informationen über Zwischenzustände aus experimentellen Daten zu erhalten. Diese Lücke wird von chemischen Simulationen geschlossen.

In der theoretischen Chemie fasst man die Gesamtheit aller Kernkoordinaten im so genannten Konfigurationsraum zusammen. Die potentielle Energie bildet eine Hyperfläche in diesem Raum. Minima stellen stabile oder metastabilen Zustände dar, die über längere Zeit existieren können. Sattelpunkte repräsentieren Übergangszustände, die die energetisch ungünstigste Konfiguration auf einem Weg von einem Minimum zu einem anderen darstellen. Kennt man nun den Weg, den eine chemische Reaktion im Konfigurationsraum zurücklegt, können alle Zwischenzustände dargestellt werden. Allerdings sind die nötigen Berechnungen sehr aufwändig. Anders als im Experiment kann in einer Simulationen in der Regel nicht der ganze Konfigurationsraum erkundet werden. Dies würde in den meisten Fällen die vorhandene Computerkapazität überfordern. Es müssen also Techniken zum Einsatz kommen, die den Ablauf einer Reaktion ermitteln können ohne dafür den gesamten Raum aller Freiheitsgrade abtasten zu müssen.

In einem thermodynamischen Ensemble, zum Beispiel dem Inhalt eines Reagenzglases, müssen zusätzlich noch statistische Eigenschaften miteinbezogen werden. Das hierfür relevante Potential ist die Freie-Energie-Landschaft, bei der ein Teil der Freiheitsgrade im Konfigurationsraum statistisch gemittelt wurden. In diese Funktion gehen also statistische und energetische Eigenschaften ein. Die Freie Energie kann in Allgemeinen nur durch statistische Simulationen (molekulardynamische oder Monte-Carlo) berechnet werden. Leider kommen die Übergangszustände, an denen wir besonderes Interesse haben, bei statistischen Simulationen besonders selten vor. Daher müssen Techniken verwendet werden, die eine ausreichende Aufl"osung auch im Bereich dieser seltenen Ereignisse erlauben.

In dieser Arbeit wird ein nicht-physikalisches quadratisches Potential, ein sogenannter "bias", zur Bewegungsgleichung der Teilchen hinzugefügt. Dadurch wird erreicht, dass auch ungünstige Bereiche im Konfigurationsraum ausreichend aufgelöst werden. Dieses Technik wird Umbrella Sampling genannt. Das Bias-Potential wirkt nur entlang einer oder mehrerer Koordinaten, die die Reaktion beschreiben. Sie werden deshalb Reaktionskoordinaten genannt. Eine Umbrella-Sampling-Simulation wird eine Fluktuation der Reaktionskoordinate ergeben, deren Erwartungswert nahe am, aber im Allgemeinen nicht auf dem Minimum des Bias-Potentials liegen wird. Aus der Differenz zwischen Erwartungswert und Minimum des Bias-Potentials kann der Gradient, der zugrundeliegenden Freien Energieoberfläche berechnet werden. Ähnlich kann die Hessematrix der freien Energie aus der Kovarianz-Matrix der Verteilung der Reaktionskoordinate hergleitet werden. Diese Methode zur Interpretation der Daten, die man durch Umbrella Sampling erhält, nennen wir Umbrella-integration.

Zunächst verwenden wir diese Werte um damit die Sattelpunkte, die die Übergangszustände repräsentieren, iterativ zu bestimmen. Ausgehend von diesen Konfigurationen kann der Reaktionspfad durch Verfolgen des Gradienten ermittelt werden. Dieser Algorithmus wurde für das System Alanindipeptid erfolgreich getestet und veröffentlicht.

Diese einfache Methode hat die Nachteile, dass sie zum Einen intrinsisch seriel arbeitet und zum Anderen recht gute Startwert für die Sattelpunkte voraussetzt. In einem zweiten Teil der Arbeit wird deshalb die etablierte Methode der Nudged Elastic Band (NEB)-Optimierungen für die Anwendung in der Freien-Energie-Landschaft erweitert. Bei NEB-Optimierungen beginnt man mit einer ersten Abschätzung für den Reaktionspfad, dieser wird durch eine Anzahl von Konfiguration, sogenannten Bildern, mit gleichem Abstand diskretisiert. Damit die Bilder während der Optimierung gleichmäßig verteilt bleiben, führt man eine nicht physikalische Federkraft ein, die zwischen den Bildern wirkt. Zur eigentlichen Optimierung wird nun die Summe aus der Komponente der Kraft des zugrundeliegenden Potentials, die senkrecht zum Pfad wirkt, und die Komponente der Federkraft die parallel zum Pfad wirkt, minimiert.

In dieser Arbeit wird zuerst ein Optimierer entwickelt, der für bestimmte Testsysteme quadratische Konvergenz von NEB-Optimierungen innerhalb der potentiellen Energie erreicht. Dieser neue Optimierer benutzt für jeden Schritt nicht nur den Gradient sondern auch die Hessematrix des zugrundeliegenden Potentials. Für die Freie-Energie-Oberfläche können beide Werte durch Umbrella-integration ermittelt werden. Um also NEB-Optimierungen in der Freien-Energie Oberfläche durchführen zu können, wird in dieser Arbeit wie folgt vorgegangen. Zuerst wird ein Startpfad z.B. eine gerade Linie zwischen zwei Punkten, meist Minima, als Kleinster-Freier-Energie-Pfad angenommen. Dieser Pfad wird in mehrere Bilder diskretisiert, um dann an jeder dieser Konfigurationen molekulardynamische Simulationen mit Umbrella-integration durchzuführen. Die so erhaltenen Gradienten und Hessematrizen finden nun Eingang in den neu entwickelten NEB-Optimierer. Auf diese Weise benötigt man weniger gute Startwerte für die Sattelpunkt, als Startwerte reicht eine Interpolation zwischen viel leichter zug"anglichen Minima. Darüber hinaus ist die Methode intrinsisch parallel, da zu jedem Bild unabhängige molekulardynamische Simulationen durchgeführt werden müssen. Diese Methode wurde für das System Alanindipeptid erfolgreich getestet und mit den Ergebnissen der seriellen Methode verglichen. Anschließend wird der Algorithmus auch für das deutlich aufwendigere System der Bindung von Ligand und Rezeptor in Wasser eingesetzt.