Monte Carlo Methods for Simulation of Protein Folding and Titration

Abstract

The present work is divided into two parts: The first part presents the folding of peptides using a Monte Carlo method for the long-term dynamics of proteins. The second part contains studies of the titration and redox behavior of the bacterial photosynthetic reaction center and of myoglobin. By combining both parts, I suggest a new method to calculate titration behavior taking into account full conformational flexibility.

The Knapp group has been working on a Monte Carlo method for the long-term dynamics of proteins for a long time. I improved this method for the first part of this work. After that, it was possible to fold not only a model peptide to a helix-turn-helix motif, as already done before, but also to simulate the folding of a ß-hairpin. However, so far these calculations all lacked a solvent model. Therefore, the Analytical Continuum Solvent model of Michael Schaefer was integrated into our method. Using this model, the formation and melting of a polyalanine helix could be simulated in agreement with results from a molecular dynamics simulation with an explicit solvent model. Also, the folding of a fragment of ribonuclease A was successful, but the folding of the ß-hairpin forming peptide BH8 failed. This unveiled so far unknown problems of our protein model with rigid peptide planes and special effective torsion potentials for the backbone torsion angles. The folding of BH8 was successful using a fully flexible protein model. However, using the flexible model the folding is much less efficient as expected for the rigid model.

The calculation of protonation and redox states of titratable and redox-active groups of proteins is done on the basis of continuum electrostatics, similar to the solvent model used in the folding simulation described above. In continuum electrostatics, the solvent is represented by a medium of a high dielectric constant. By solving the Poisson-Boltzmann equation numerically on a grid, the resulting electrostatic potentials (and thus also the interaction energies) can be calculated. In this way, the total electrostatic energy of a specific protonation and redox state of a molecular system can be determined. A protein contains usually a lot of titratable or redox-active groups. For n of such groups, the total number of possible states is extremely large (2n) so that the calculation of thermodynamical averages, which are necessary to determine titration curves and redox potentials, by exact summation is computationally infeasible. Therefore, I apply in the present work a Monte Carlo method where an importance sampling of the huge number of possible states is performed according to the Metropolis criterion. The method yields good results after relatively short sampling times. The statistical error of these results can be probed by evaluating a correlation function.

By applying these methods, I could gain a lot of valuable insights. For the first time, the energy of the electron transfer from QA.- to QB in the bacterial reaction center was calculated correctly. The calculations also yield important hints for the so far unsettled sequence of the following electron transfer and protonation events involving the quinones of the reaction center. In addition, the conformational gating hypothesis of the electron transfer from QA.- to QB could be supported from the viewpoint of theory, and the insight into the related processes was deepened. I present in this work several approaches to include conformational ensembles and conformational flexibility in the calculation of titration behavior. By explicit consideration of conformational relaxation, the dielectric constant of the protein interior could be decreased remarkably in accordance with fundamental principles. By including an ensemble of x-ray structures, the pH induced conformational changes of myoglobin could be reproduced. The pKa values calculated in this study are in better agreement with experimental values than any other result obtained by non-trivial methods before.

In the concluding outlook, I discuss the future improvements of the Monte Carlo method for the long-term dynamics of proteins, and the assets and drawbacks of existing methods to calculate titration behavior considering conformational flexibility. I show how to avoid most of the drawbacks by a new method on the basis of a combination of the Monte Carlo dynamics with the titration calculation. This method realizes the unrestricted and unbiased sampling of the conformational space and the space of titration states at the same time.

Die vorliegende Arbeit gliedert sich in zwei Teile. Der erste Teil behandelt die Faltung von Peptiden mit einer Monte-Carlo-Methode zur Langzeitdynamik von Proteinen. Der zweite Teil enthält Studien zum Titrations- und Redoxverhalten des bakteriellen photosynthetischen Reaktionszentrums und von Myoglobin. Durch die Kombination beider Teile skizziere ich am Schluß der Arbeit eine neue Methode zur Titrationsberechnung von Proteinen mit vollständiger Konformationsflexibilität.

Die Arbeitsgruppe Knapp arbeitet bereits seit langer Zeit an der Entwicklung einer Monte-Carlo-Methode zur Langzeitdynamik von Proteinen. Diese Methode wurde von mir für den ersten Teil meiner Arbeit weiterentwickelt. Es gelang daraufhin nicht nur, wie bereits zuvor, ein Modellpeptid zu einem Helix-Turn- Helix-Motif zu falten, sondern auch die Faltung eines ß-Hairpins zu simulieren. Diese Berechnungen beinhalteten jedoch noch kein Modell für das Lösungsmittel. In einem weiteren Schritt wurde daher das Analytical Continuum Solvent-Modell von Michael Schaefer in unsere Methode integriert. Mit diesem Modell konnte der Aufbau und das Schmelzen einer alpha-Helix aus Polyalanin in Übereinstimmung mit Ergebnissen einer Molekulardynamik mit explizitem Lösungsmittelmodell simuliert werden. Auch die Faltung eines Fragmentes der Ribonuclease A war erfolgreich, wohingegen die Faltung des ß-Hairpin-bildenden Peptids BH8 mißlang. Hier zeigten sich bislang nicht erkannte Schwächen unseres Proteinmodells mit starrer Geometrie der Peptidebene und speziellen effektiven Potentialen für die Torsionswinkel. Die Faltung von BH8 konnte erfolgreich mit einem voll flexiblen Proteinmodell simuliert werden, allerdings weitaus weniger effizient als man es bei Verwendung des starren Proteinmodells erwarten würde.

Die Berechnung der Protonierungs- und Redoxzustände von titrierbaren bzw. redoxaktiven Gruppen in Proteinen erfolgte, wie auch schon beim für die Simulation der Proteinfaltung verwendeten Lösungsmittelmodell, auf Basis der Kontinuumselektrostatik, bei der das Lösungsmittel durch einen Bereich erhöhter Dielektrizitätskonstante simuliert wird. Durch numerisches Lösen der Poisson-Boltzmann-Gleichung auf einem Gitter lassen sich die resultierenden elektrostatischen Potentiale und damit auch die elektrostatischen Interaktionsenergien berechnen. Damit kann man die elektrostatische Energie eines bestimmten Protonierungs- und Redoxzustandes des Gesamtproteins bestimmen. Ein Protein verfügt in der Regel über eine große Zahl n von titrierbaren bzw. redoxaktiven Gruppen. Die Gesamtzahl an möglichen Zuständen (2n) ist extrem groß, so daß sich thermodynamische Mittelwerte, die zur Bestimmung von Titrationskurven und Redoxpotentialen nötig sind, praktisch nicht exakt berechnen lassen. Daher verwende ich in der vorliegenden Arbeit eine Monte-Carlo-Methode, mit der eine Teilmenge der Vielzahl möglicher Zustände gemäß des Metropoliskriteriums durchmustert wird. Durch diese Methode erlangt man nach kurzer Zeit Ergebnisse hoher statistischer Genauigkeit, was durch Auswerten einer Korrelationsfunktion gezeigt werden kann.

Mit der genanten Vorgehensweise war es mir möglich, eine Vielzahl wertvoller Erkenntnisse zu sammeln. So konnte erstmalig die Energie des Elektronentransfers von QA.- zu QB im bakteriellen Reaktionszentrum korrekt berechnet werden. Die Berechnungen ergaben auch wichtige Hinweise für die zuvor ungeklärte weitere Reihenfolge der Elektronentransfer- und Protonierungsschritte der Chinone im Reaktionszentrum. Desweiteren konnte auch von theoretischer Seite die conformational gating-Hypothese betreffend den Elektronentransfers von QA.- zu QB untermauert und das Verständnis der zugehörigen Vorgänge im Reaktionszentrum vertieft werden. Ich entwickele in der vorliegenden Arbeit mehrere Ansätze zur Einführung multipler Konformationen und von Konformationsflexibilität in die Titrationsberechnungen. In einem Fall konnte durch die explizit berücksichtigte strukturelle Relaxiation die Dielektrizitätskonstante für das Innere des Proteins entsprechend allgemeiner Prinzipien stark verringert werden. In einem anderen Fall wurde die Energetik pH-induzierter struktureller Änderungen von Myoglobin durch Einbeziehen eines Ensembles von Kristallstrukturen verstanden. Die dabei berechneten pKa-Werte zeigen eine Übereinstimmung mit dem Experiment, die in dieser Genauigkeit niemals zuvor mit nicht-trivialen Methoden zur pKa-Wertberechnung erzielt worden sind. In einem Ausblick erörtere ich die in Zukunft möglichen Verbesserungen der Monte- Carlo-Methode zur Langzeitdynamik von Proteinen und die Vor- und Nachteile der existierenden Methoden zur Titrationsberechnung mit Konformationsflexibilität. Ich zeige auf, wie sich durch eine zukünftige Kombination der Monte-Carlo- Dynamik mit der Titrationsberechnung die meisten Nachteile der existierenden Methoden umgehen lassen. Diese Methode verwirklicht das uneingeschränkte und verzerrungsfreie gleichzeitige Durchmustern des Konformationsraumes und des Raumes der Titrationszustände.