What is a Molecular Dynamics simulation?
Molecular Dynamics is the technique of using classical mechanics or quantum mechanics to represent a system and that is using to study as time-dependent, kinetic models of the system for different kinetic and thermodynamic analyses. With the increasing complexity and recent advances, MD has become a part of supercomputer area, and different computer experiments related with biology also done by using that immense machines and complex algorithms.
In our case, MD uses Newtonian mechanics to describe the molecular system through empirical calculations of the potential energy (Leach, 2001). In this methodology, atoms of the system modelled as points with given mass and charge with pre-assigned coordinates of different atoms. The positions and charges are then used to calculate a potential energy for the system. This potential is used to compute the force felt by each of the atoms in the simulation, which is totally known as the force field. By using Newton’s laws of motion over a short time step with the obtained individual forces, a new set of positions and velocities is produced for each of the atoms. The updated values can now be fed back into the first step of the calculation, to update the position of each atom according to classical mechanics, and the process repeated for creating a trajectory of atomic locations and velocities through time. As a result, MM approach on protein systems allows us not just only discover conformations of the proteins, but also their dynamics (Frenkel and Smit, 2002).
Overall, energy function and the calculation of individual forces can be formulated as given below (Note that calculation of energy function will be stated in next sub topic) where U is total energy of the system, Fi is the force applied on each atom (i is the atom number), and mi.(d2ri / dt2) part of the equation is the mass multiplied by acceleration (m.a) part of the classical Newton’s law of motion formula shown below.
A MD simulation contains three stages: setup, equilibration, and production run. The protein structure (a crystal structure, NMR structure or a homology model) is usually not present in its natural environment, and some preparation is therefore necessary. The setup stage involves insertion of water molecules, for example determined by grid calculations (Farantos et al., 2005) into the protein structure and soaking the protein membrane complex into a box of water with included ions. In equillibration, the molecular system is first energy minimized, then gradually heated up to the target temperature of the simulation, and then equilibrated. The equilibration phase is usually a multi-step process, in which the heavy atoms in proteins are first heavily restrained and gradually released until the system is completely unrestrained and stable. After that, system run with pre-defined calculations to produce trajectory data. In this study NAMD molecular dynamics package used throughout these stages of the simulations.
The potential energy of a molecular system is evaluated through few relatively simple mathematical terms accounting for covalent bonded interactions (including bond, bond angle, and torsional terms) and nonbonded interactions (including long-range Coulombic and Van der Waals (VdW) terms). The mathematical terms together with atom/bond parameters, e.g. the radius of individual atoms and the optimal bondlengths, are referred as the force field.
During years, different force fields developed to study all-atom molecular systems. In this study, we use CHARMM22 force field to operate molecular dynamics simulations (Brooks et al., 1983). CHARMM force field interpretes energy as the sum of bonded and non-bonded terms as illustrated in the equation below.
Where bonds term in the energy function accounts for the bond stretches where kb is the bond force constant and b-b0 is the distance from equilibrium that the atom has moved, angles term in the equation accounts for the bond angles where k_θ is the angle force constant and θ-θ_0 is the angle from equilibrium between 3 bonded atoms, dihedrals term is for the dihedrals (a.k.a. torsion angles) where k_ϕ is the dihedral force constant, n is the multiplicity of the function, ϕ is the dihedral angle and ᵹ is the phase shift, the term improper accounts for the impropers, that is out of plane bending, where k_ω is the force constant and ω-ω_0 is the out of plane angle. The Urey-Bradley component (cross-term accounting for angle bending using 1,3 nonbonded interactions) comprises the fifth term, where k_u is the respective force constant and u is the distance between the 1,3 atoms in the harmonic potential. Nonbonded interactions between pairs of atoms (i, j) are represented by the last two terms. By definition, the nonbonded forces are only applied to atom pairs separated by at least three bonds. The van Der Waals (VdW) energy is calculated with a standard 12-6 Lennard-Jones potential and the electrostatic energy with a Coulombic potential. In the Lennard-Jones potential above, Rmin term is not the minimum of the potential, but rather where the Lennard-Jones potential crosses the x-axis (i.e. where the Lennard-Jones potential is zero) (Vanommeslaeghe et al, 2010).
NAMD (Nanoscale Molecular Dynamics) is a parallel molecular dynamics code designed for high-performance simulation of large biomolecular systems. Based on Charm++ parallel objects, NAMD scales the necessary calculations to hundreds of cores for typical simulations and beyond 500,000 cores for the largest simulations. NAMD uses the popular molecular graphics program VMD (Visual Molecular Dynamics) for simulation setup and trajectory analysis, but is also file-compatible with other popular force fields including AMBER, CHARMM, and X-PLOR. NAMD is distributed free of charge with source code and it can be builded to a wide variety of platforms.
VMD is also designed for modeling, visualization, and analysis of biological systems such as proteins, nucleic acids, lipid bilayer assemblies, etc. It may also be used to view more general molecules, as VMD can read standard Protein Data Bank (PDB) files and display the contained structure. VMD provides a wide variety of methods for rendering and coloring a molecule: simple points and lines, CPK spheres and cylinders, licorice bonds, backbone tubes and ribbons, cartoon drawings, and others. VMD can also be used to animate and analyze the trajectory of a molecular dynamics (MD) simulation. In particular, VMD can act as a graphical front end for an external MD program by displaying and animating a molecule that is undergoing simulation on a remote computer.
How to MD with UHEM or TRUBA Supercomputers of Turkey
Getting started;
PDB: Protein Data Base file type, your protein structure presented here.
R value: 0.2 or lower considered as well defined (
Choose complete crystals, otherwise CHARMM GUI will not recognize it
Exclude waters from structure (Might be done by CHARMM-GUI)
Ligands: Another issue, parameters should have presented for every ligand. If not, use cgenff.com. Well-known ligands have pre-defined parameters at CHARMM-GUI files.
SETUP (VMD-VPN CLIENT-WINSCP-PUTTY): You will use VPN client to access the supercomputer, WINSCP to access distant folders at supercomputer and PUTTY to reach control panel to control supercomputer.
CHARMM GUI: An online web-site interfaced database that helps to generate simulation files
Regulatory files:
Run.job: to initiate simulation
Loop.tmp: parameters will be applied for each cycle
Loopnumber.txt: file that counts loop cycle
Colvars.in: collective variables module
*There are also some corrections that needs to be done in the CHARMM-GUI files
1. Choose your crystal structure from the literature (PDB is a very useful database for that job). Make sure it is intact (all amino acids on sequence presented at PDB file) and have a preferable resolution (below 2.5 Angstrom).
*If your structure is not intact, search for online tools.
https://swissmodel.expasy.org/interactive is an online database that can make your protein structure complete based on presented sequence data.
2. Go to the http://www.charmm-gui.org/ , choose INPUT GENERATOR from the left menu, and then, choose the QUICK MD SIMULATOR from again the same menu.
3. Upload your .pdb file and select PDB format to select model/chain, then click NEXT STEP at right-bottom.
4. In that step, choose your chain of interests and ligands, and exclude non-interested molecules, then click next.
5. In that step, you’ll add PDB manipulation options and many different settings can be applied from there. As for standardized MD applications, we use terminal group patching by C-TER and N-TER selections. Besides, you should add disulfide bonds from there.
*If you are not sure about disulfide bonds of your proteins, you might check from
http://clavius.bc.edu/~clotelab/DiANNA/
*If you are not sure about protonated-deprotonated residues, you might also check them from
http://nbcr-222.ucsd.edu/pdb2pqr_2.0.0/
Adjust your desired parameters from this section and then click Next Step for solvating our system of molecules.
6. In that step, water molecules around our molecule will be placed to represent cellular environment. We use standard 10 A edge distance and Monte-Carlo method for ion placing.
*If you want to change the edge distance of water box, use Calculate number of ions button to re-adjust ion concentrations
After the adjustments, click next. Passing this step may take a little longer time, be patient.
7. After step 6 end, pass step 7 by clicking Next Step (nothing to change here for standard simulations).
8. Choose only NAMD from Input Generation Options and then leave the rest.
* We use NVT Ensemble for equilibration input generation and NPT at 30315 K Temperature for Dynamic input generation. These may have considered as the standards of MD simulation with CHARMM force field and NAMD software.
9. As a last step at CHARMM-GUI server, download simulation files from download .tgz link, exlude from the .tgz file and upload charm-gui folder to the necessary folder at your supercomputer inside a folder that is named for your job (ex: your_job).
10. As mentioned above, in addition to these files, you will also need 4 additional file for starting the simulations and some adjustments on CHARMM-GUI files.
11. Upload run.job, loop.tmp, loopnumber.txt, and colvars.in to the ../your_job/charm-gui/namd folder
12. These files have pre-written informations to guide you. Enter your own parameters like project number, simulation type, queue that will be used etc.
13. As a final adjustment, go to the step4_equillibration.inp and step5_production.inp files at ../charm-gui/namd file, open these two and deleted the line
exec python ../checkfft.py $a $b $c >checkfft.str
Then, go to the ../charm-gui folder, copy checkfft.str file from there and paste to the ../your_job/charm-gui/namd
*We use that step to by-pass python execution since its type may change from different computers
14. Finally, all you need to do is typing commands sbatch
15. You might control the process by squeue command for TRUBA or /RS/progs/isler commands for UHEM.
Experimental Design
First, we prepared a flowchart of our experimental procedures as shown in the figure below.
Investigation of Sequences:
After we decided our proteins that will be used in our biosensor system, we used GenBank to obtain DNA sequences of those protein of interests.
IDT Orders:
Once we are convinced with the designs (proper plasmid and restriction site selections) we performed on Benchling, we ordered our DNA fragments in the extent of IDT’s sponsorship.
Deciding Host Strains:
We used Escherichia coli (E. coli) DH5alpha strain for cloning purposes and E. coli Shuffle T7 express cells for protein production host.
Cloning the Parts:
We used pET28a(+) vector for both cloning and expression of the parts we used in this project.
Isolation of Plasmids:
For plasmid isolation, Thermo Fisher GeneJET plasmid mini-prep kit was used.
Transformation:
Obtained plasmids transformed into E. coli Shuffle T7 express strain for protein production.
Measurements:
Fluorescence measurements done via Jasco FP-750 spectrofluorometer.
Excitation and emission parameters for Clover: 505 nm – 515 nm
Excitation and emission parameters for mRuby2: 559 nm – 600 nm
Excitation and emission parameters for FRET effect: 505 nm – 600 nm
