Molecular dynamics simulations of a small synthetic protein

1. Objective

The objectives of this practical are

•    to  perform. a series of  molecular dynamics (MD) simulation of a protein in vacuum and in water

•    to  obtain  trajectory  files  to  be  analysed  using  readily  available  tools  and  by  developing Python-based tools using available libraries.

•    an introduction to running calculations on shared high performance computing systems.

In this session, your aim is to perform. molecular dynamics on a small synthetic protein (chignolin) that is known to fold rapidly. You will use a widely-used simulation package, GROMACS, for  molecular dynamics. Due to the size of the simulations, these will be run on the University’s high performance computer, Lyceum (this is a part of a larger computing system; Lyceum is used for undergraduate teaching and research).  You should have received an email about your account setup on Lyceum; you can log in with your university credentials.

Some of your  molecular  dynamics simulations will  be  run  during this  session. To  have  molecular dynamics that will run quickly enough, you will first simulate the protein in vacuum. Your aims will be to run molecular dynamics on a starting geometry taken from the Protein Data Bank (PDB) - these are experimentally determined structures - and on a geometry of the protein that is generated on the computer in an extended (unfolded) state.

We will use a python library MDanalysis to analyse the results of the molecular dynamics, giving you scripts that you can use to analyse later calculations that you will set up to run over the time between this session and the 3rd  molecular dynamics practical session.

So the second aim of the practical is to set up calculations of the protein in water that you can set running by the end of the session. For this, we need to get a few preliminary steps performed during today’s session: energy minimisation and ‘equilibration’.

In total, we are aiming for the following simulations to compare their behaviour:

PDB starting geometry: in vacuum at 280 K, in vacuum at 310 K,  in water at 280 K, in water at 310 K

Linear (unfolded) geometry: in vacuum at 280 K, in vacuum at 310 K, in water at 280 K, in water at 310 K

Vacuum runs are to be run for 30 ns in length. Simulations in water should be 160 ns in length.

2. Practical report

Your report for this week should include a few components.

You should have a number of moIecuIar dynamics simuIations compIete by the end of today’s practicaI session. These are: simulations of the protein in vacuum at two temperatures, as well as energy minimisation and (equiIibration’moIecuIar dynamics of the protein in water. These two stages of the simuIations in water precede the (production’simuIations, which you wiII anaIyse in the third practicaI session (week 25). Please submit in zip or tar format the contents of your working directory on Lyceum. You should have a well-organised directory structure for your calculations. We can help in making the zip or tar file.

Submit a Jupyter notebook of your analysis scripts. You will use a python package called MDanalysis  to perform. analysis of your molecular dynamics. You should have code that performs a few types of analysis: i) monitoring the root mean squared deviation in structure between the simulated structure and the reference structure of the protein over the course of the molecular dynamics. You should submit this code, along with results from analysis of the protein in vacuum (at two temperatures) and equilibration of the protein in water; ii) code to monitor the distance between the two ends of the protein through the course of the simulations, along with plots from the vacuum simulations and equilibration of the protein in water and iii) code to calculate a histogram of the distance between the two ends of the protein, followed by calculation of the free energy. You will use these scripts to analyse the (production’runs that are performed over the time between this session and the 3rd  molecular dynamics practical.

You will be asked to comment on some of the results. You may include these comments in the jupyter notebook or as a separate pdf.

The zip (or tar) and pdf might not be able to be submitted on eassignments. These can be emailed to

[email protected]

3. Requirements

This exercise requires access to a machine running any Linux flavour with GROMACS installed. Other useful software packages are:

-     Visual Molecular Dynamics5, 6  (VMD) - www.ks.uiuc.edu/Research/vmd

-     PyMOL7 一 pymol.org

-     MDAnalysis8 - mdanalysis.org using your university login details

-     Access to the  Lyceum  HPC. From a Windows machine, MobaXTerm is recommended. See

https://hpc.soton.ac.uk/redmine/projects/iridis-5-support/wiki/Connecting_to_Iridis_5

The theory behind MD simulations is not covered in these notes. See section 3.3 of Essentials of computational chemistry theories and models 2nd ed. (Christopher Cramer), available electronically on the university library website.

You will need to gain some basic skills in using the HPC system Lyceum, which uses the Linux operating system. The following resources might be helpful (and please ask the demonstrator!):

https://cheatography.com/davechild/cheat-sheets/linux-command-line/pdf/

http://www.ee.surrey.ac.uk/Teaching/Unix/

For visualising the protein, you can use PyMol or VMD and the following should be helpful. https://chryswoods.com/dynamics/visualisation/README.html(VMD)

https://bioquest.org/nimbios2010/wp-

content/blogs.dir/files/2010/07/pymol_tutorial3.pdf(PyMol)

4. PDB structure files

The first step to perform. a simulation is to obtain a structure PDB file for the protein of interest. In this practical, we will study the small synthetic protein chignolin, you can access a structure at the RCSB10  Protein Data Bank (rcsb.org) using PDB ID 1UAO.

This file describes the 3D structure of the protein using a PDB formatting scheme. In a PDB file, the first 8 columns describe the type of information records in each line, as laid out in Table 1, following the specifications  available at wwpdb.org12. The  description  of the  3D  structure of a  molecule  is essentially condensed in the COORDINATE section, in ATOM lines, where all atom info is included. In the case of biomolecules, atoms are organized in residues, either protein or nucleic acid residues.

Table 1. Sections allowed in a PDB file and their content.

ATOM lines constitution is laid out inTable 212. In the case of proteins, other molecules present in the structure file (water molecules, ions, etc.) that are not part of the protein are described in

HETATM lines. HETATM and ATOM lines are organized in the same way and are exemplified in Figure 1.

Table 2. ATOM lines in a PDB file.

Figure 1. ATOM/HETATM line organisation in a PDB file.

5. Topology preparation

The pdb file contains an ensemble of 18 conformers of chignolin determined by NMR. First we need to remove the unwanted NMR structures from the PDB file. We can do this using any text editor. You may keep any of the 18 structures.

The system topology can be created using the pdb2gmx module. Before running the command for GROMACS, you will need to load the module, which sets up the software for you. You will need to do this each time you login with the command:

module load gromacs/2022.3/intel

Now that you have the module loaded, a simple interactive version can be run with the command line

gmx_mpi pdb2gmx -f 1UAO_clean.pdb -o 1UAO_processed.gro

and a more complete and user-independent command would use more options, such as

gmx_mpi pdb2gmx -f 1UAO_clean.pdb -o 1UAO.gro -p topology.top -ff amber99sb -water tip3p -ignh

This command invokes GROMACS (gmx_mpi) module pdb2gmx with the option -f to read a file and -o to write a file, invoking the system's topology with -p and specifying the forcefield (-ff) and water model (-water). Another option frequently used is the -ignh one, in which all H atoms in the input file are disregarded and H atoms are de novo assigned to the molecule based on the force field files.

A full list of options for these and all other GROMACS module can be found in the online documentation (manual.gromacs.org), under the correct version – just search for command line reference.

The pdb2gmx command generates various ouput files, in particular the topology.top one, that describes the system and will be updated in several steps, and a number of gro extension files, a simpler version of PDB files for internal GROMACS use together with a topology.

A sample gro file is shown in

Figure 2. The first line of the files is a title for the system, and the second one lists the total number of atoms in the system. The last one includes the vectors that define the simulation box, in the sequence v1(x) v2(y) v3(z) v1(y) v1(z) v2(x) v2(z) v3(x) v3(y), in which only the three first values are required to describe a parallelepipedal simulation box, while the others are considered zero.

Topology files contain all the necessary information about the system structure, organized by atoms, bonds, angles, dihedrals, impropers, etc, allow the inclusion of other topologies (other molecules) in the system, and allow applying restrictions – applying a force to a given set of atoms to prevent those atoms from moving too much.

In a topology file, comments are marked with a semicolon ; - lines starting with # are not comments.

Figure 2. A sample gro file – initial and final sections.

LYSOZYME

1323

1LYS N 1 3.536 2.234 -1.198

1LYS H1 2 3.612 2.288 -1.236

1LYS H2 3 3.470 2.214 -1.270

[…]

129LEU C 1321 4.212 1.238 0.737

129LEU O1 1322 4.173 1.228 0.856

129LEU O2 1323 4.323 1.268 0.691

5.90620 6.84510 3.05170

The steps described in this section are the necessary ones, and are enough if you are in fact using the 1UAO model. If you use other proteins, you will need to be sure the file you retrieve can be used. Search for anisotropy lines, and other alternative locations. Many experimental models also miss some atoms and even residues. While one or two atoms are easily added using Jmol, adding whole residues is more demanding. Browse the internet for your problem – someone is likely to have had it before and found a solution!