Coarse-grained models for molecular dynamics of RNA. Evolutionary algorithm optimization of the potential energy function parameters.

Zrozumienie dynamiki cząstek kwasów nukleinowych jest zadaniem fascynującym, choć trudnym w wykonaniu. Obecnie dostępne są liczne struktury kwasów nukleinowych otrzymane przy użyciu dyfrakcji promieniowania X, mikroskopii elektronowej lub jądrowego rezonansu magnetycznego. Jednakże są to struktury statyczne, uśrednione w czasie. Niezwykle trudno jest z nich odczytać, jak wygląda dynamika tych cząstek. Symulacje komputerowe przy użyciu dynamiki molekularnej pozwalają na uzupełnienie statycznych struktur molekuł o ruch. Od lat 70-tych metoda ta jest stale ulepszana i stosowana do licznych zagadnień z pogranicza biologii, fizyki i chemii. Ograniczenia symulacji przy użyciu dynamiki molekularnej zmniejszają się każdego roku dzięki stale rosnącej mocy obliczeniowej dostępnych komputerów. Jednakże nadal znaleźć można cząsteczki zbyt duże lub skale czasowe zbyt długie dla dynamiki molekularnej. Na dzień dzisiejszy najdłuższe symulacje przy użyciu pełnoatomowej (każdy atom w cząsteczce symulowany jest niezależnie) dynamiki molekularnej obejmują czas milisekund dla średniej wielkości białek. Jednakże istotne procesy biologiczne, jak np. translacja białek na rybosomie, trwają sekundy lub minuty. Symulacja tych problemów przy użyciu metod pełnoatomowych możliwa będzie najprawdopodobniej dopiero za kilka lat. Rozwiązaniem na chwilę obecną może być użycie modeli zredukowanych (lub gruboziarnistych). Cząstka w takim modelu symulowana jest nie jako zbiór atomów, ale zbiór centrów interakcji (kulek) odpowiadających grupom chemicznym lub nawet całym nukleotydom i aminokwasom. Liczba symulowanych obiektów i oddziaływań pomiędzy nimi znacznie maleje, a zatem ograniczenia symulacji ulegają redukcji o kilka rzędów wielkości. Jednakże uproszczenie symulacji pociąga za sobą pogorszenie zgodności symulacji z danymi eksperymentalnymi. Jak silny jest ten spadek, zależy od tego, jak dobrze przybliżona jest energia potencjalna układu przez model zredukowany. Jeżeli przybliżenie jest dobrze dopasowane do potrzeb konkretnego problemu biologicznego, to wyniki otrzymane przy użyciu modelu zredukowanego będą satysfakcjonujące. W przeciwnym przypadku mogą okazać się zupełnie bezużyteczne. Przygotowanie dobrego przybliżenia energii potencjalnej (pola siłowego) nie jest jednak zadaniem prostym. Parametry pola siłowego dla modelu zredukowanego mogą być otrzymane z analizy zbioru znanych struktur eksperymentalnych lub z analizy symulacji pełnoatomowej dynamiki molekularnej. Ze względu na niedoskonałości danych wejściowych parametry otrzymane wymienionymi sposobami często wymagają dalszej optymalizacji za pomocą licznych prób i~błędów. W tej rozprawie zaprezentowano metodę pozwalająca na optymalizację parametrów gruboziarnistego pola siłowego przy użyciu algorytmów metaheurystycznych. Dzięki temu drugi etap tworzenia pola, optymalizacja za pomocą prób i błędów, jest zautomatyzowany i mniej arbitralny. Program komputerowy sprawdza tysiące kombinacji parametrów gruboziarnistego pola siłowego, by zaproponować optymalne rozwiązanie. Algorytm ocenia zestawy parametrów, wykonując przy ich użyciu testowe symulacje dynamiki molekularnej, a następnie weryfikując zgodność trajektorii z tych dynamik z danymi odniesienia. Razem z opracowaniem metody optymalizacji parametrów przygotowane zostało oprogramowanie, RedMDStream, implementujące metodę. Jest ono dostępne pod adresem: http://bionano.cent.uw.edu.pl/Software/RedMD Program pozwala nie tylko na optymalizację parametrów, ale także można przy jego pomocy przygotować struktury gruboziarniste na podstawie współrzędnych atomów i opisu topologii modelu gruboziarnistego. RedMDStream pozwala również na prowadzenie gruboziarnistych dynamik molekularnych i ich ocenę pod kątem stabilności, mobilności atomów oraz dystrybucji odległości w trakcie symulacji. W rozprawie, poza opisem metody i oprogramowania, przedstawiono też przykład przygotowania gruboziarnistego pola siłowego dla termometru RNA. Opisano całą procedurę: (1) wykonanie pełnoatomowej dynamiki molekularnej, jako punktu odniesienia dla dalszej parametryzacji, (2) projektowanie modelu zredukowanego z jedną kulką opisującą nukleotyd, (3) optymalizacja parametrów gruboziarnistego pola siłowego. Pokazano również, w jaki sposób powyższa metoda pozwala na wybór najlepszego zestawu równań dla funkcji energii potencjalnej.

Understanding dynamics of nucleic acid structures is a fascinating although not a simple task. Currently many X--ray, electron microscopy or nuclear magnetic resonance structures of nucleic acids are available, however these are still static images averaged over time and captured in at most few conformational states. Internal motions of these molecules can be hardly deduced from these data. The movement into the structures can be introduced using a computer simulation with the molecular dynamics (MD) method. From its first applications in 1970s the MD method was constantly improved and is currently applied to a wide range of problems. MD applications extend each year thanks to a constant rise in available computational power. However, there are still molecules too large or timescales too long for efficient MD simulations. Currently the longest MD simulations for moderate size proteins achieve millisecond timescales with the state of the art computer clusters and atomistic scale models (with motions of each atom simulated independently). However, there are biological processes that last over seconds or minutes, like protein translation on the ribosome. Such tasks will not be accessible for full--atomistic models for at least a few more years. Using reduced, coarse--grained (CG) models can work as a remedy. A molecule in a reduced model is represented as a set of interaction centers (beads) that represent functional groups or even whole residues, nucleotides or amino acids. The number of simulated entities is largely reduced once beads are used instead of atoms. Such reduction makes simulations more efficient by orders of magnitude, increasing limits on the system size or timescale. However, compliance with experimental data is compromised by the reduction. How severe is the loss in the quality is dependent on a potential energy function approximation (force field) used by a reduced model. If a force field is well adapted for the aim of the simulation, the results are satisfactory, otherwise useless data might be produced. The task is unfortunately not trivial. Force field parameters can be deduced by analyzing structural databases or trajectories from atomistic force fields. For some cases, like structure prediction, such approach is sufficient. However, models for dynamics usually require additional trial and error optimization due to imperfections of the first guess. This thesis presents a method for optimization of the force field parameters for coarse--grained molecular dynamics models using a metaheuristic algorithm. The method provides an objective, automated way to explore thousands of combinations of these parameters to find the optimal set. The algorithm assesses the quality of the parameters, by running and analyzing decoy simulations and checking how well the force field with these parameters models the dynamics of biomolecules. A RedMDStream application implementing the algorithm was prepared together with this thesis and is made available to download from a following webpage: http://bionano.cent.uw.edu.pl/Software/RedMD The application provides not only the force field optimization routines, but also allows preparing coarse--grained structures from atomistic data using a provided coarse--grain topology. The application enables running molecular dynamics simulations and analyzing the outcome of the dynamics in terms of the stability, mobility and distance distributions. Along with the description of the method and software, an example of coarse--grained model parameterization for an RNA thermometer hairpin is presented. This example provides a proof of concept of the method. All steps of the procedure are presented: (1) obtaining a reference full--atomistic molecular dynamics trajectory, (2) designing a reduced model with one bead representing a whole nucleotide and (3) optimizing the model parameters. It is also shown how the procedure allows the user to make choices about the potential functional form.