Team:ZJU-China/model s1 3.htm

From 2012.igem.org

Latest revision as of 01:34, 27 October 2012

Molecular Dynamics

Introduction

Molecular dynamics (MD) is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms by numerically solving the Newton's equations. The results of molecular dynamics simulations includes the trajectory of molecules and atoms in the system and other parameters such as energy potential, temperature and pressure fluctuations, number of hydrogen bonds of the system, which may be used to find the properties of such complex systems analytically.

Methods

Here in our project, we use GROMACS 4.5[5] to run MD and analyze the state of theophylline-riboscaffold complex. GROMACS is a versatile package to perform molecular dynamics on macromolecules like proteins, lipids and nucleic acids. It is fast, has multiple choices for force field, and comes with a large selection of flexible tools for trajectory analysis.

Although GROMACS provides step by step tutorials for beginners, they are all for proteins or protein-ligand complex. Our riboscaffold-theophylline system is similar to protein-ligand complex to some degree, but they have their own characteristics. We, as beginners, are trying to develop a protocol for the MD simulation of ligand sensible RNA devices based on canonical tutorials for protein-ligand complex [6, 7]. It may be coarse now, and we will modify it and add more details after more experiments, we hope it is helpful for people who are also beginners and interested in RNA molecular dynamics.

Step1: Separate the RNA and ligand from a complex file.

We used the complex file generated by AutoDock and separate them into two PDB files with the coordinate of atoms conserved.

Step2: Prepare the protein topology with pdb2gmx.

Choose AMBER force fields because they are more suitable for RNA MD. In our cases, we chose AMBER99SB.

Step3: Prepare the ligand topology using external tools.

Proper treatment of ligands is one of the most challenging tasks in molecular simulation. The ligand should be treated with the force field that is consistent with the receptor. Unlike the tutorials for protein-ligand complex, in which the PRODRG web server with GROMOS force field is used, we prepared our ligand with Antechamber [8] and tleap [9], free tools in the AmberTools package. The output files are a top and a gro file for the ligand. The top file should be converted into itp file for reusability and included in the receptor top file.

Step4: Solvation

At this point, it is necessary to define the unit cell and fill it with water, just as described in other tutorials for any other MD simulation.

Step5: Adding ions

A .tpr file must be generated for the input. Therefore, we will setup energy minimization, which means running grompp before genion. Energy minimization in this step is different form that in the following step, here we just want to know the charges of our system. Our system has a net charge of -157e, so we add 157 Na+ to simulate the status in cytoplasm.

Step6: Energy Minimization

Having added ions, we need to rerun energy minimization. In the .mdp file, remember to set energygrps to the right names of the elements in the system. It is recommended to checked the [molecule] part of the .top file to make sure all kinds of molecules in the system are correctly included.

Step7: Equilibration

Similar to the protein-ligand complex, in the stage of equilibration, we should first restrain our ligand (but it is not fixed!) around the binding site during equilibration, control the temperature coupling, and run two phases of equilibration, NVT (canonical ensemble, moles N, volume V and temperature T are conserved.) and NPT (isothermal–isobaric ensemble, in which moles N, pressure P and temperature T are conserved.)

Step8: Production MD

After two phases of equilibration, the system has reached the desired temperature and pressure. So we will run MD again and get the products.

View the source code and more details here.

Results

Figure 7. The potential energy of the system.

The potential energy equilibrates about an average of ~ –3274000 kJ/mol.

Figure 8. Temperature fluctuation.

Figure 9. Pressure fluctuation.

The temperature and pressure plot shows the normal oscillation behavior. The average temperature is about 310 K, which is what we desired.

Figure 10. Hydrogen bonds.

The number of hydrogen bonds ranges mainly between 3 and 4. This figure shows the interaction is rather stable.

Figure 11. RMSD of RNA.

Figure 12. RMSD of theophylline.

The two figures about RMSD indicates that the RNA, as a receptor, equilibrates rather quickly. But the small ligand, theophylline, equilibrates slower, with a larger range of fluctuation (about 0.02, but in RNA, the range is within 0.05) in conformation. This is normal in ligand-receptor systems.

In conclusion, the preliminary analysis results show that the simulation system generated by our workflow acts similar to other ligand-receptor systems under the environment of molecular dynamics. Actually, with this method, we can conduct more analysis on the level of molecules and even atoms. We will continue our way to explore this magic microworld!