\documentclass[a4paper,12pt,twoside]{article} \usepackage[utf8]{inputenc} \usepackage{ngerman} \usepackage{graphicx} \usepackage{subfigure} \usepackage{threeparttable} \usepackage{a4wide} \usepackage{esvect} \usepackage{pstricks} \usepackage{pst-plot} \usepackage{pstricks-add} \usepackage{multicol,epsf} \usepackage{ae,aecompl} \usepackage{helvet} \usepackage{amsmath} \usepackage{amssymb} \usepackage{bm} \usepackage{bbm} \usepackage[mathcal]{euscript} \usepackage[version=3]{mhchem} \usepackage[T1]{fontenc} \usepackage{hyperref} \definecolor{linkcolor}{rgb}{0,0,0} \hypersetup{colorlinks=true, breaklinks=true, linkcolor=linkcolor, menucolor=linkcolor, urlcolor=linkcolor, bookmarksopen=true, bookmarksnumbered=false} \numberwithin{equation}{subsection} \numberwithin{figure}{subsection} \pagestyle{headings} \begin{document} \newcommand{\diff}[1]{ \, \mathrm{d} #1 \, } \newcommand{\diffm}[2]{ \, \mathrm{d}^{#1} \mathbf{#2} \, } \newcommand{\db}{ \, \mathrm{dB} } \newcommand{\comb}{ \mathrm{comb} } \newcommand{\F}{ \mathcal{F} } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Titelseite % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \normalsize \rule{0mm}{0,4cm} \thispagestyle{empty} \begin{figure}[htbp] \centering \includegraphics[scale=0.1]{UniGoe.png} \label{Georg-August-Universität Göttingen} \end{figure} \normalsize \begin{center} \Huge Versuch B-2 \\[5mm] {\bf Molekulardynamiksimulationen \\ von Proteinen} \end{center} \normalsize \vspace{5mm} \centerline{\epsfysize=2cm \hfill {\bf Praktikum für Fortgeschrittene} \hfill \epsfxsize=2cm} \centerline{\epsfysize=2cm \hfill am Dritten Physikalischen Institut \hfill \epsfxsize=2cm} \centerline{\epsfysize=2cm \hfill der Universität Göttingen\hfill \epsfxsize=2cm} \vspace{10mm} \begin{center} 18. Januar 2009 \end{center} \vspace{7mm} \begin{center} \begin{tabular}{rl} \rule{5cm}{0mm} & \rule{5,5cm}{0mm} \\ {\it Praktikant} \ \ & {\bf Johannes Dörr} \\ & \href{mailto:mail@johannesdoerr.de}{mail@johannesdoerr.de} \\ & \href{http://physik.johannesdoerr.de}{physik.johannesdoerr.de} \\ \rule{5cm}{0mm} & \rule{5,5cm}{0mm} \\ {\it Durchführung } \ \ & am 16.12.2008 \\ & zusammen mit \\ & Jan Schumann-Bischoff \\ \rule{5cm}{0mm} & \rule{5,5cm}{0mm} \\ {\it Betreuer} \ \ & Christian Kappel \\ \rule{5cm}{0mm} & \rule{5,5cm}{0mm} \\ \end{tabular} \end{center} \newpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Unterschrift % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \thispagestyle{empty} \begin{center} \rule{0cm}{15cm} \begin{tabular}{l} {\it Unterschrift des Praktikanten:} \\ \rule{0cm}{1.5cm} \\ \rule{8cm}{0.3pt} \\ \small{Johannes Dörr - Göttingen, den 18.01.2009} \end{tabular} \end{center} \newpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Inhaltsverzeichnis % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \tableofcontents \newpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (Dokument) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Einleitung} Proteine sind Makromoleküle, welche in jeglichen Lebewesen vorkommen. Sie erfüllen die verschiedensten Aufgaben (z.B. Transport von Molekülen usw.). Die Kenntnis ihres Aufbaus und ihrer Funktion sind daher wichtig für das Verstehen von biologischen Systemen. Neben Experimenten trägt mittlerweile auch die Simulation der Proteinen zum Verständnis bei. Bei einer solchen Molekulardynamiksimulationen wird der Ort und die Geschwindigkeit der einzelnen Atome numerisch unter Berücksichtigung verschiedenster Vereinfachungen berechnet. Dieser Versuch soll einen Einblick in diese Technik geben und durch das Simulieren einfacher Systeme die Funktionsweise einer solchen Simulation erklären. \section{Theorie} \subsection{Molekulare Wechselwirkungen} Die Coulomb-Kraft zwischen zwei Teilchen der Ladung $q_1$ und $q_2$ im Abstand $r$ ist gegeben durch \begin{align*} F \propto \frac{q_1\cdot q_2}{r^2} \ \ . \end{align*} Sie ist neben dem Aufbau von Atomen entscheidend für molekulare Wechselwirkungen. In diesem Abschnitt werden die wichtigsten Wechselwirkungen beschrieben, die bei der Betrachtung von Proteinen eine entscheidende Rolle spielen. \subsubsection{Kovalente Atombindungen} Edelgase sind sehr reaktionsträge und kommen daher atomar vor. Die {\it Oktettregel} besagt, dass auch andere Stoffe eine {\it Edelgaskonfiguration} anstreben, was durch kovalente Autombindungen geschehen kann. Dabei teilt sich ein Atom mit einem anderen ein oder mehrere Elektronen, um 8 Elektronen in der Außenschale zu haben. Beispielsweise Sauerstoff tritt praktisch immer als Molekül $O_2$ auf, da das einzelne Atom nur 6 Elektronen in der Außenschale besitzt, durch Austauschen von zwei Elektronen mit einem anderen Sauerstoff-Atom beide aber eine Edelgaskonfiguration annehmen. \begin{figure}[h] \centering \includegraphics[width=0.5\textwidth]{bilder/sauerstoff_kov_bind.png} % sauerstoff_kov_bind.png: 323x51 pixel, 72dpi, 11.39x1.80 cm, bb=0 0 323 51 \caption{Kovaltene Bindung bei Sauerstoff.} \label{fig:sauerstoff_kov_bind} \end{figure} Da es sich bei den Bindungspartnern hierbei um dieselbe Atomsorte handelt, liegt der Ladungsschwerpunkt genau zwischen beiden Atomen. \subsubsection{Wasserstoff-Brücken} Handelt es sich jedoch um Atomone verschiedener Elektronegativität, entstehen Partialladungen. Wichtiges Beispiel hierfür ist Wasser. Hierbei teilen sich zwei Wasserstoffmoleküle jeweils ein Elektron mit dem Sauerstoff. Da Wasserstoff sehr {\it elektropositiv} ist, halten sich die Elektronen sehr nah beim Sauerstoff auf, sodass die positive Ladung der Wasserstoff-Protonen kaum kompensiert wird. Diese positive Ladung wirkt anziehend auf die negative Seite eines anderen Wasssermoleküls. Diese {\it Wasserstoff-Brücke} ist verantwortlich für die besonderen Eigenschaften von Wasser wie zum Beispiel die sehr hohe Oberflächenspannung oder der späte Siedepunkt. \begin{figure}[h] \centering \includegraphics[width=0.4\textwidth]{bilder/Watermolecule.png} % Watermolecule.png: 729x662 pixel, 80dpi, 23.15x21.02 cm, bb=0 0 656 596 \caption{Geometrie und Partialladung eines Wassermoleküls.} \label{fig:wassermolekuel} \end{figure} \subsubsection{Hydrophober Effekt} Mischt man nun in die polaren Wassermoleküle unpolare Atome/Moleküle, so können diese keine Bindung mit dem Wasser eingehen. Als Folge der Entropiemaximierung (2. Hauptsatz der Thermodynamik) trennen sich polare und unpolare Moleküle unter Minimierung der Oberfläche. Dies nennt man den {\it hydrophoben Effekt}. Durch die verringerte Zahl von polaren Molekülen direkt neben unpolaren steigt die Entropie. In der Biologie ist dies entscheidend für die Faltung von Proteinen im Wasser. Im Alltag beobachtet man, dass sich in der Suppe an der Oberfläche kleine Fettäugchen bilden. \subsubsection{Van-der-Waals-Wechselwirkung} Auch in einer Ansammlung von unpolaren Molekülen, die keine kovalenten Bindungen bilden, treten Wechselwirkungen auf. Diese liegen darin begründet, dass die Ladungsverteilung im Atom nicht statisch ist. Es treten fluktuierende Partialladungen auf, die in den Nachbaratomen wiederum Dipolmomente induzieren können. Die dadurch entstehende anziehende Wechselwirkung kann angenähert werden durch das {\it Lennard-Jones-Potential}: \begin{align*} \phi(r) = -\frac{A}{r^6} + \frac{B}{r^{12}} \ \ . \end{align*} Sein Verlauf ist in Abbildung \ref{fig:lennardjones} dargestellt. Die Van-der-Waals-Wechselwirkung tritt immer auf, oft ist sie jedoch im Vergleich mit den anderen Bindungen sehr schwach. \begin{figure}[h] \centering \includegraphics[width=0.7\textwidth]{bilder/lennardjones.png} \caption{Verlauf des Lennard-Jones-Potentials} \label{fig:lennardjones} \end{figure} \subsection{Aufbau von Proteinen} Proteine zählen zu den wichtigsten Bausteinen in der Biologie. Sie funktionieren beispielsweise als Katalysator für chemische Reaktionen oder als Transportmittel für andere Moleküle. Sie sind Makromoleküle, also vergleichsweise groß, und bestehen aus einer Abfolge von {\it Aminosäuren}. Im Menschen findet man 22 verschiedene Aminosäuren, obwohl es noch wesentlich mehr gibt. Die Aminosäuren sind über die s.g. {\it Rückgratatome}, nämlich zwei Kohlenstoff- und einem Stickstoffatom, über eine kovalente Atombindung miteinander verkettet. An dem einen Kohlenstoffatom sitzt die {\it Seitengruppe}, die die chemische Eigenschaft der Aminosäure bestimmt und sie damit von einer anderen unterscheidet. Abbildung \ref{fig:aminokette} zeigt ein Beispiel für eine Kette aus fünf Aminosäuren mit Rückgrat und Seitenkette. Echte Proteine haben hingegen eine Länge von über 100 Aminosäuren. \begin{figure}[h] \centering \includegraphics[width=0.7\textwidth]{bilder/protein.png} \caption{Die Rückgratatome der dargestellten Aminosäurenkette sind als Kugeln, die Seitenkette als Stäbe dargestellt.} \label{fig:aminokette} \end{figure} Die für ein Protein charakteristische Abfolge der Aminosäuren nennt man {\it Primärstruktur}. Sie ist codiert in der DNA eines Organismus, mit der dieser sie vererben kann. Bedingt durch den hydrophoben Effekt faltet sich ein Protein augenblicklich im Wasser. Dabei bilden sich zwischen Atomen der Seitenketten und denen der Rückgratatome Wasserstoffbrücken aus, wodurch räumliche Strukturen wie $\alpha$-Helices und $\beta$-Faltblätter (siehe Abbildung \ref{fig:aminokette}) entstehen. Sie bilden die {\it Sekundärstruktur} des Proteins. \begin{figure}[h] \centering \includegraphics[width=0.5\textwidth]{bilder/sekundaer.png} \caption{{\it links:} $\beta$-Faltblatt; {\it rechts:} $\alpha$-Helix.} \label{fig:aminokette} \end{figure} Diese Struktur faltet sich nun noch einmal und bildet die {\it Tertiärstruktur}, die sowohl durch Wasserstoffbrücken aber auch durch einfache Coulomb-Anziehung stabilisiert wird. Ein fertig gefaltetes Protein ist in Abbildung \ref{fig:protein2} dargestellt. Anordnungen von Proteinen nennt man schließlich {\it Quartärstruktur}. \begin{figure}[h] \centering \includegraphics[width=0.45\textwidth]{bilder/protein2.png} \caption{Gefaltetes Protein.} \label{fig:protein2} \end{figure} Die Faltung eines Proteins ist nicht nur ein Nebeneffekt der Minimierung der potentiellen Energie. Sie ist vielmehr Ausschlaggebend für seine Funktionsweise. Zum Beispiel das innerhalb dieses Versuchs betrachtete Protein Lysozym faltet sich zu einer Art schwingenden Schere, die Zellwände von Bakterien zerstören kann. Auf Grund der Komplexität der Faltung ist es bis dato nicht möglich, diese zu simulieren, also vorherzusagen. Im Laufe des Versuchs sind wir also auf Datenbestände angewiesen, in denen Proteinstrukturen, die mit Röntgenbeugungsuntersuchungen gewonnen wurden, gespeichert sind. Deren Dynamik kann dann mit Methoden simuliert werden, die im Folgenden dargestellt werden. \subsection{Simulation von Molekulardynamik} Die Molekulardynamiksimulation (MD-Simulation) erlaubt es, die Wechselwirkungen von Atomen bzw. Molekülen für einen gewissen Zeitraum unter vorgegebenen Bedingungen zu berechnen und grafisch darzustellen. Für eine exakte Berechnung der Vorgänge müsste das System quantenmechanisch betrachtet und seine Schrödingergleichung gelöst werden. Dies ist mit der heute verfügbaren Rechenleistung jedoch nur bis zu einer Molekülanzahl von 10 möglich. Da Proteine aus einigen 100 Aminosäuren bestehen, ist man bei ihrer Simulation auf Näherungen angewiesen. Die MD-Simulation macht deshalb drei Vereinfachungen: (1) Das ungleiche Massenverhältnis von Elektronen und Atomen wird ausgenutzt, um ihre Freiheitsgrade zu separieren. Dabei wird angenommen, dass sich die sehr schnellen Elektronen in einem effektiven Potential der ruhenden Atomkerne bewegen. Dies ist die {\it Born-Oppenheimer-Näherung}. (2) Dieses effektive Potential wird durch eine Summe analytischer Funktionen genähert, die leicht (schnell) zu berechnen ist. Hierbei gehen sowohl die Bindungen, die die chemische Struktur bestimmen, als auch langreichweitigere Wechselwirkungen ein. Die Parameter stammen hierbei aus Experimenten sowie aus quantenmechanischen Berechnungen. Während des Versuchs wird das OPLS-Kraftfeld (Optimized Potential for Liquid Simulations) verwendet, das sich besonders gut für die Simulation von Lösungsmitteleffekten eignet. Es gibt noch weitere, deren Verwendung je nach Komplexität der betrachteten Moleküle zu mehr oder weniger abweichenden Ergebnissen führt. (3) Die Bewegung der Atomkerne wird mit den {\it Newtonschen Bewegungsgleichungen}, also klassisch, gelöst, wobei es sich selbstverständlich auch hierbei nur um eine Näherung handelt. Dabei kommt der {\it leapfrog-Algorithmus} zum Einsatz. Bei der dabei durchgeführten numerischen Integration ist darauf zu achten, dass die Schrittweite $\Delta t$, in der Geschwindigkeit und Position berechnet werden, klein genug ist, um die stattfindenden Bewegungen zu erfassen. Bei den von uns untersuchten Proteinen ist hierfür eine Femtosekunde angemessen. Um die Simulation trotz Näherungen möglichst genau zu halten, ist es ratsam, vor dem Start eine {\it Energieminimierung} durchzuführen. Dabei wird die vorgegebene Ausgangslage der Moleküle, beispielsweise ein Protein, so lange leicht variiert, bis ein Minimum der Energie vorliegt. Dies verhindet Konstellationen (z.B. Überlappungen von Atomen), in denen sehr starke Kräfte wirken und die deshalb zu unrealistischen Ergebnissen führen würden. Bei kleinen Systemen mit abzählbar vielen Molekülen wie in unserem Fall sind {\it periodische Randbedingungen} sehr praktikabel, da hierdurch keine Teilchen verloren gehen. Man umgibt das System im Prinzip mit belibig vielen Kopien von sich selbst. Verlässt ein Teilchen das System auf der einen Seite, trifft sofort ein neues von der gegenüberliegenden Seite ein. \subsection{Analysen} Neben der bildlichen Anschauung von molekularen Prozessen erlaubt die MD-Simulation noch rechnerische Analysen. \subsubsection{Root Mean Square Deviation (RMSD)} Der RMSD-Wert sagt aus, wie weit eine Stuktur im zeitlichen Verlauf von einer Referenzstruktur abweicht: \begin{align*} \mathrm{RMSD} &= \sqrt{\frac{1}{N}\sum_{i=1}^N(\mathbf{x}_i-\mathbf{x}_{i,\mathrm{ref}})^2} \ \ . \end{align*} Dabei ist $N$ die Anzahl der Elemene der Struktur, $\mathbf{x}_i$ die Position des $i$-ten Elements und $\mathbf{x}_{i,\mathrm{ref}}$ sein Referenzwert. Dabei hängen die $\mathbf{x}_i$ und damit der RMSD-Wert von der Zeit ab. Es handelt sich um eine Mittelung über die Elemente. \subsubsection{Root Mean Square Fluctuation (RMSF)} Der RMSF-Wert hingegen gibt die über die Simulationszeit gemittelte Abweichung von der Referenzstruktur für jedes Element aus. Hiermit lässt sich beispielsweise erkennen, welche Teile des Systems stark schwingen. \subsubsection{Principal Component Analysis (PCA)} Manchmal scheint die Bewehung eines simulierten Moleküls unstrukturiert. Vermutet man dennoch ein Bewegungsmuster, so kann eine PCA nützlich sein. Dazu bestimmt man für alle Zeitschritte die Kovarianzmatrix aller Atome um ihren Mittelwert, wobei $x_i$ für jede Raumkoordinate steht, also $i=1...3N$: \begin{align*} C_{ij} &= \left\langle M_{ii}^{1/2}(x_i-\langle x_i\rangle)M_{jj}^{1/2}(x_j-\langle x_j\rangle)\right\rangle \ \ . \end{align*} Dabei enthält die Matrix $M$ die Massen der Atome. Die Matrix enthält, wie stark die Atombewegungen miteinander korrelieren. Wird sie diagonalisiert, so erhält man sie zur Basis ihrer Eigenvektoren. Die größten Eigenwerte zeigen dann die größte kollektive Bewegung des Systems. \subsubsection{Mean Square Deviation (MSD)} Zur Bestimmung der Diffusionskonstanten $D$ wird die Formel \begin{align*} \lim_{t\rightarrow \infty}\left\langle \left| \mathbf{x}(t)-\mathbf{x}(0)\right|^2\right\rangle &=6Dt \end{align*} verwendet. Beim Auftragen der mittleren quadratischen Entfernung des Teilchens von seinem Ausgangspunkt in Abhängigkeit von der Zeit ergibt sich aus der Steigung die Diffusionskonstante. \section{Durchführung} \subsection{Teil 1 - Argon} In diesem Versuchsteil soll man sich mit der Software \verb|Gromacs| vertraut machen und Grenzen herausgefunden werden. Dazu wird das Van-der-Waals Gas Argon bei einer Temperatur von 100k simuliert und die Trajektorie eines einzelnen Atoms beobachtet. Danach wird das Gas in folgenden Schritten abgekühlt und aufgeheizt: \begin{itemize} \item Abkühlung von 100K auf 25K in 1ns \item Abkühlung von 100K auf 25K in 5ns \item Aufheizen von 50K auf 125K in 5ns \end{itemize} Die Phasenübergänge sollen näher untersucht werden. Im nächsten Schritt wird flüssiges Argon untersucht. Folgende Systeme sollen analysiert werden: \begin{itemize} \item Eine Simulation bei 94.4K für 1ns \item Eine Simulation für 1ns, bei der das Argon von 100K auf 0K abgekühlt wird \end{itemize} Aus den Ergebnissen werden dann die MSD und die RDF bestimmt. \subsection{Teil 2 - Kleines Protein} Im zweiten Versuchsteil wird das Protein G bei 300K und konstantem Druck für 50ps simuliert. Die Struktur des Proteins entnehmen wir der RCSB Protein Datenbank. Es werden die folgenden Schritte durchgeführt. \begin{itemize} \item Topologie erstellen \item Größe des Simulationssystems festlegen \item Hinzufügen von Wassermolekülen in das Simulationssystem, da die Röntgenbeudung, aus der die Struktur ermittelt wurde, diese nicht erkennt. \item Hinzufügen von Ionen \item Energieminimierung des Systems \end{itemize} Aus der Simulation werden die folgenden Observablen ermittelt: potentielle Energie, RMSD-Werte, RMSF-Werte und B-Faktoren. Außerdem werden die Häufigkeit von bestimmten Sekundärstrukturen (Ramachandran-Plot) sowie die Größe des {\it Solvent Accessible Surface} bestimmt. \subsection{Teil 3 - Principal Component Analysis (PCA) von Lysozym} Weil eine Simulation des Protein Lysozym zu zeitaufwändig wäre, wird das Ergebnis einer bereits durchgeführten Simulation verwendet. Es werden die Eigenwerte (der Kovarianzmatrix) der Proteinbewegung bestimmt. Ebenfalls wird eine Projektion der ersten 8 Eigenvektoren analysiert. Die größte Hauptbewegungsrichtung wird näher betrachtet. \section{Auswertung} \subsection{Teil 1 - Argon} Es wurde eine Simulation von Argon bei einer Temperatur von 100K über die Dauer von 1ns laufen gelassen. Beobachtet man nur ein Atom, so ist zu erkennen, dass dieses Atom durch Stöße mit anderen Atomen die Geschwindigkeit und die Richtung ändert. Die für die Simulation eingestellten periodischen Randbedingungen zeigen sich darin, dass das Atom, wenn es an einen Rand stößt, verschwindet und am gegenüberliegenden Rand mit gleicher Richtung und Geschwindigkeit wieder auftaucht. \subsubsection{Phasenübergänge} In der Simulation wurde das Argon jeweils in 1ns und 5ns von 100K auf 25K angekühlt (Abb. \ref{fig:abkuehlung_argon}). \begin{figure}[htbp] \subfigure[1ns]{ \label{fig:abkuehlen_100k_25k_1ns} \includegraphics[width=0.50\textwidth]{bilder/teil1/abkuehlen_100k_25k_1ns.png}} \subfigure[5ns]{ \label{fig:abkuehlen_100k_25k_5ns} \includegraphics[width=0.50\textwidth]{bilder/teil1/abkuehlen_100k_25k_5ns.png}} \caption{Abkühlung von Argon von 100k auf 25k in 1ns und 5ns.} \label{fig:abkuehlung_argon} \end{figure} Es ist zu beobachten, dass sich bei 1ns mehrere kleine Grüppchen von Argon\-atomen gebildet haben und bei 5ns sich alle (bis auf ein) Atome zu einer großen Gruppe zusammenschließen. 1ns reicht als für eine Simulation der Kondensation nicht aus. Bei 5ns kann man hingegen davon ausgehen, dass das Gas sehr gut kondensiert ist. Beim Aufheitzen von 25k auf 100k in 5ns sieht man, dass sich aus der einen großen Gruppe von Atomen mit zunehmender Temperatur immer mehr herausfliegen, bis sich das Argon wieder vollständig im gasförmigen Zustand befindet. Plottet man die potentielle Energie für alle drei Simulationen über die Temperatur (Abb. \ref{fig:pot_energy_argon}), lässt sich jeweils der Siedepunkt ablesen. \begin{figure}[htbp] \subfigure[Abkühlung 1ns]{ \label{fig:abkuehlen_100k_25k_1ns_pot_energy} \includegraphics[width=0.50\textwidth]{bilder/teil1/abuehlen_100k_25k_1ns.pdf}} \subfigure[Abkühlung 5ns]{ \label{fig:abkuehlen_100k_25k_5ns_pot_energy} \includegraphics[width=0.50\textwidth]{bilder/teil1/abuehlen_100k_25k_5ns.pdf}} \subfigure[Aufheizen 5ns]{ \label{fig:aufheitzen_25k_100k_5ns_pot_energy} \includegraphics[width=0.50\textwidth]{bilder/teil1/aufheitzen_25k_100k_5ns.pdf}} \caption{Potentielle Energie in Abh. der Temperatur zur Bestimmung des Siedepunktes.} \label{fig:pot_energy_argon} \end{figure} Man erkennt den Siedepunkt daran, dass beim Abkühlen aufgrund der Grüppchenbildung die potentielle Energie anfängt abzunehmen. Beim Aufheizen nimmt ab dem Siedepunkt die potentielle Energie zu. Die Ergebnisse sind in Tab. \ref{tab:siedetemp} dargestellt. \begin{table}[htpb] \centering \begin{tabular}{|c|c|} \hline \textbf{Plot} & \textbf{Siedetemperatur} \\ \hline Abkühlung von 100k auf 25k in 1ns & 50K \\ Abkühlung von 100k auf 25k in 5ns & 55K \\ Aufheizen von 25k auf 100k in 5ns & 88K \\ \hline \end{tabular} \caption{Siedetemperatur von Argon aus den Simulationen.} \label{tab:siedetemp} \end{table} Der Literaturwert für die Siedetemperatur liegt bei 87.3K. Die Werte aus der Abkühlung weichen deutlich von diesem Wert ab. Lediglich der Wert aus der Simulation des Aufheizens liegt ziemlich genau beim Literaturwert.\\ \\ \textbf{Diffusionskonstante von flüssigem Argon}\\ Wir haben 94.4k kaltes flüßiges Argon 1ns lang simuliert. Mit dem MSD-Verfahren haben eine Diffusionskonstante von \begin{align*} D &= 2.3254(2084)\cdot 10^{-5}\frac{\textrm{cm}}{\textrm{s}} \end{align*} berechnet. Der Literaturwert liegt bei $2.43\cdot 10^{-5}\frac{\textrm{cm}}{\textrm{s}}$ und weicht somit 4\% vom Literaturwert ab.\\ \\ \textbf{RDF}\\ Die RDF (\textit{radial distribution function}) beschreibt die Verteilung von Materie um einen bestimmten Punkt herum. In unserem Fall wird die RDF um alle Atome gemittelt. Wir haben die RDF aus zwei Simulationen gewonnen: zum einen Argon bei 94.4K und zum anderen Argon, welches innerhalb von 1ns von 100K auf 0K abgekühlt wird. Für letztere Simulation haben wir die RDF nur für die zweite Hälfte der Erstarrungssimulation bestimmt. Die Ergebnisse sind in Abb. \ref{fig:abkuehlung_argon_rdf} dargestellt. \begin{figure}[htbp] \subfigure[94.4K]{ \label{fig:rdf_94k} \includegraphics[width=0.50\textwidth]{bilder/teil1/RDF_fluessiges_Argon_94k.pdf}} \subfigure[Abkühlung 100K auf 0K]{ \label{fig:rdf_abkuehlen_100k_0k_1ns} \includegraphics[width=0.50\textwidth]{bilder/teil1/RDF_fluessiges_Argon_100k_0k_abkuehlen.pdf}} \caption{RDF für Argon bei 94.4K und für die zweite Hälfte der Erstarrungssimulation von Argon, welches von 100K auf 0K in 1ns abgekühlt wird.} \label{fig:abkuehlung_argon_rdf} \end{figure} Bei beiden Plots ist zu erkennen, dass sich innerhalb eines Radius von ca. 0.4nm keine weiteren Atome befinden. Dies ist gerade der dopellte Atomradius. Desweiteren sieht man in Abb. \ref{fig:rdf_94k} bei 0.4nm einen relativ schmalen Peak, welcher von einem deutlich breiteren gefolgt wird. Die RDF fällt oberhalb von 0.4nm nicht mehr auf 0 ab, was auf thermische Fluktuationen des Argons zurück zu führen ist. In Abb. 4.1.\ref{fig:rdf_abkuehlen_100k_0k_1ns} ist zu erkennen, dass für einen Radius größer als 0.4nm die RDF neben scharfen Peaks auch wieder auf Null abfällt, was zeigt, dass hier die Atome in der kristallinen Struktur dichter und regelmäßiger gepackt sind. Für das Abkühlen von 100K auf 0K in 1ns kann man, wenn man die pot. Energie über die Temperatur aufträgt (Abb. \ref{fig:rdf_abkuehlen_siedepunkt}), den Schmelzpunkt bestimmen. \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil1/abkuehlen_100k_0k_1ns_rdf.pdf} % abkuehlen_100k_0k_1ns_rdf.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{Potentielle Energie in Abhängigkeit von der Temperatur zur Bestimmung des Schmelzpunktes bei der Abkühlung von 100K auf 0K in 1ns.} \label{fig:rdf_abkuehlen_siedepunkt} \end{figure} Der Schmelzpunkt liegt etwa bei 48K. Dies weicht deutlich vom Literaturwert von 83.8K ab. \subsection{Teil 2 - Simulation eines kleinen Proteins} In diesem Versuchsteil wird ein kleines Protein bei 300K bei konstantem Druck simuliert. Die Simulationsdauer beträgt 50ps. \subsubsection{Potentielle Energie} Abbildung \ref{fig:pot_en_teil2} zeigt die potentielle Energie in Abhängigkeit von der Zeit. Der Wert geht schnell gegen einen bestimmten Wert, die thermischen Fluktuationen, die das Protein ausführt, sind dabei aber deutlich zu erkennen. \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil2/pot_energy_time.pdf} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{Potentielle Energie in Abhängigkeit von der Zeit.} \label{fig:pot_en_teil2} \end{figure} \subsubsection{RMSD-Wert} Der RMSD-Wert in Abhängigkeit von der Zeit ist in Abbildung \ref{fig:rmsd_teil2} aufgetragen. Er scheint nicht (zumindest innerhalb der gemessenen Zeit) gegen einen Grenzwert zu gehen. Das Protein erreicht also nicht die Referenzstruktur. Auch hier sind thermischen Fluktuationen zu sehen. \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil2/rmsd.pdf} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{RMSD-Wert in Abhängigkeit von der Zeit.} \label{fig:rmsd_teil2} \end{figure} \subsubsection{RMSF-Wert und B-Faktor} Abbildung \ref{fig:rmsf_teil2} zeigt den RMSF-Wert für alle Atome des Proteins. Es ist zu erkennen, dass einige Atome besonders stark fluktuieren. Die B-Faktoren sind in Abbildung \ref{fig:brac_teil2} farblich dargestellt. Er zeigt, wie stark ein Röntgenreflex auf Grund der starken Fluktuationen des jeweiligen Atoms abgeschwächt wird. \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil2/rmsf.pdf} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{RMSF-Wert in Abhängigkeit vom Atom-Index.} \label{fig:rmsf_teil2} \end{figure} \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil2/bfac.png} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{B-Faktor in farblicher Darstellung.} \label{fig:brac_teil2} \end{figure} \subsubsection{Solvent Accessible Surface} Die Größe der Oberfläche, die von den Wassermolekülen (generell auch von einem anderen Lösungsmittel) zugänglich ist, nennt man {\it Solvent Accessible Surface}. Diese kann bestimmt werden, indem eine Kugel, deren Radius der Größe eines Moleküls entspricht (1.4\AA{} bei Wasser) auf dem Protein abgerollt wird. Wichtig hierbei ist auch, ob die Stellen, die das Lösungsmittel berühren, hydrophil oder hydrophob sind. Letzterer Fall ist energetisch ungünstiger, weshalb die hydrophobe Fläche des Proteins kleiner ausfällt als die hydrophile (siehe Abb. \ref{fig:solv_acc_surface}). \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil2/solvent_accessible_surface.pdf} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{Solvend Accessible Surface im Verlauf der Simulation.} \label{fig:solv_acc_surface} \end{figure} \subsubsection{Ramachandran-Plot} Am {\it Ramachandran-Plot} (siehe Abb. \ref{fig:ramachandran_teil2}) ist zu erkennen, wie oft bestimmte Formen im Protein auftreten. Er zeigt das statistische Vorkommen von Kombinationen zweier Diederwinkel. In unserem Fall sind die Alpha-Helix und Beta-Sheets am häufigsten. \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil2/ramachandran.pdf} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{Ramachandran-Plot des untersuchten Proteins.} \label{fig:ramachandran_teil2} \end{figure} \subsection{Teil 3 - PCA von Lysozym} Mit einer PCA (Principal Component Analysis) lassen sich die Hauptbewegungsrichtungen des Lysozym bestimmen. Die größten Eigenwerte der Kovarianzmatrix sind in Abb. \ref{fig:kovmatrix_teil3} dargestellt. Bis auf die ersten 80 sind alle anderen ca. 1500 Eigenwerte nahezu Null. \begin{figure} \centering \includegraphics[width=0.7\textwidth]{bilder/teil3/eigenwerte.pdf} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{Eigenwerte der Kovarianzmatrix.} \label{fig:kovmatrix_teil3} \end{figure} Signifikant von Null verschieden sind nur die ersten 8. Dies bedeutet, dass es ca. 8 verschiedene Richtungen gibt, in die das Protein merklich fluktuiert. In Abb. \ref{fig:projektion} ist die Projektion der ersten 8 Eigenwerte dargestellt. \begin{figure} \centering \includegraphics[width=0.5\textwidth]{bilder/teil3/eigenvec_1-8.pdf} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{Projektion der ersten 8 Eigenvektoren.} \label{fig:projektion} \end{figure} In der Zeitentwicklung der Eigenvektoren ist eine ungenaue Periodizität zu erkennen. Das Protein führt also periodische Bewegungen aus. Dabei führt der größte Eigenvektor die langsamste Bewegung aus. Diese wäre nicht erkennbar, wenn die Simulationsdauer kürzer gewählt werden würde (z.B. 1ns bis 10ns). Mit Hilfe des Programms \verb|g_anaeig| lassen sich die Bewegungen der einzelnen Eigenvektoren isoliert berechnen. Tun wir dies für den größten Eigenvektor, so ist eine deutliche Schneidbewegung zu erkennen (Abb. \ref{fig:protein_schneid}). \begin{figure} \centering \includegraphics[width=1\textwidth]{bilder/teil3/schneidbew.png} % pot_energy_time.pdf: 792x612 pixel, 72dpi, 27.94x21.59 cm, bb=0 0 792 612 \caption{Schneidebewegung des Lysozym. \textit{links}: offen, \textit{mitte}: halboffen, \textit{rechts}: geschlossen} \label{fig:protein_schneid} \end{figure} Das Protein ist in der Lage, Zellwände abzubauen. Dies steht in Einklang mit der von uns beobachteten Schneidefunktion. \section{Diskussion} Der Versuch liefert einen guten Einblick in die Simulation von Molekülen sowie in molekulare Vorgänge. Es stellt sich jedoch heraus, dass eine solche Simulation nicht immer korrekte Ergebnisse liefert, wie die Untersuchung des Argongases gezeigt hat. Ebenso ist der nötige (Rechen-)Aufwand sehr deutlich geworden. Unsere Eingangsgfrage, ob die gesamte Struktur von Proteinen, also auch ihre räumliche Faltung simuliert werden kann, hat sich somit sehr deutlich beantwortet. Für die qualitative Betrachtung von Schmelz- und Siedevorgängen eignet sich der Versuch sehr gut und ist deshalb durchaus empfehlenswert. \end{document}