Hairpin formation

From OxDNA
Revision as of 16:36, 13 April 2012 by Romano (talk | contribs)
Jump to navigation Jump to search

Introduction

In this example, you will simulate a single strand of length 18 and sequence GCGTTGCTTCTCCAACGC at 334 K (~61 C) in three different ways:

  • with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is inputMD.
  • with an MD simulation of the sequence-dependent (SD) model. The input file is inputMD_seq_dep.
  • with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is inputTRAP.

The traps act between the pairs depicted in blue and red in the sequence GCGTTGCTTCTCCAACGC. The details of the interaction associated to the traps can be changed in the file hairpin\_forces.dat.

This strand, if T is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop. The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the \textit{hairpin} example (Section~\ref{sec:setup}) and how to run it (Section~\ref{sec:run}). Section~\ref{sec:results} contains results and plots extracted from the simulation output. In the following, $EXEC refers to the oxDNA executable.

Preparation

The script run.sh generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 10^8} steps by default, takes approximately one hour on a modern CPU. The default run.sh expects $EXEC to be in the ../.. directory. If this is not the case, open run.sh and change the variable CODEDIR accordingly.

If you only want to generate the initial configuration, you can issue ./run.sh --generate-only. Then you can run the simulations by yourself. The generated initial configuration files are initial.top (which contains the topology) and initial.conf (which contains positions and orientations of the nucleotides).

Running

In order to run the whole example, it is sufficient to issue the command ./run.sh (or bash run.sh). As described in Section~\ref{sec:setup}, you can generate the initial configuration and then run the simulations by hand. The three simulations, described in Section~\ref{sec:introduction}, can be performed issuing $EXEC input, where input is a text file that specifies the simulation configuration. For this example, three files have been prepared: inputMD, inputMD_seq_dep and inputTRAP. Table~\ref{tbl:sim} report all the files associated to each simulations.

Type Input energy file trajectory file last configuration file log file
SA model inputMD energy.dat trajectory.dat last_conf.dat log.dat
SD model inputMD_seq_dep energy_seq_dep.dat trajectory_seq_dep.dat last_conf_seq_dep.dat log_seq_dep.dat
SA model with traps inputTRAP energy_trap.dat trajectory_trap.dat last_conf_trap.dat log_trap.dat

Results

As mentioned in Section~\ref{sec:introduction}, temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure~\ref{fig:histo_MD} shows the probability distribution histogram $P(U_{HB})$ for the hydrogen-bonding (HB) energy $U_{HB}$ of the system (in simulation units~\cite{ouldridge_jcp}). Both the minimum in $P(U_{HB})$ and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since $U_{HB}$ oscillates between $0$ (no HBs) and $\approx -3.5$, which corresponds to $4-5$ hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed~\cite{ouldridge_jcp}. This is clearly visible in the snapshot at the beginning of this document, in which $4$ base-pairs are formed.