Biomodeling: Folding of TrpCage

TrpCage protein

https://www.rcsb.org/structure/1L2Y

Processing PDB for Gromacs using gmx pdb2gmx

Gromacs flow chart taken from the official documentation.

topol.top: topology file

Determine the periodic boundary condition

TITLE     TC5B

REMARK    THIS IS A SIMULATION BOX

CRYST1    1.000    1.000    1.000  90.00  90.00  90.00 P 1           1

MODEL        1

ATOM      1  N   ASN A   1      -8.901   4.127  -0.555  1.00  0.00           N

ATOM      2  H1  ASN A   1      -8.345   3.925   0.252  1.00  0.00            

ATOM      3  H2  ASN A   1      -8.685   5.045  -0.887  1.00  0.00            

https://en.wikipedia.org/wiki/Lattice_constant

Periodic images of TrpCage under truncated dodecaheron-type PBC

Solvate

[ molecules ]

; Compound        #mols

Protein_chain_A     1

SOL              2157

mdp: MD Option file

Add ions to the system and Minimize

MD simulation

MD for 10 ns


check box size and temperature

Analysis


Analysis Methods

Q, RMSD, and Rg


import numpy as np

import matplotlib

matplotlib.use('Agg')

import matplotlib.pyplot as plt


#######################################

# Load trajectory

traj = md.load('all.1ns.xtc', top='native.pdb')

native = md.load_pdb('native.pdb')


#######################################

# Calculate Q, the fraction of native contacts

# https://mdtraj.org/1.9.4/examples/native-contact.html

q = best_hummer_q(traj, native)

plt.plot(q)

plt.xlabel('Frame', fontsize=14)

plt.ylabel('Q(X)', fontsize=14)

plt.show()


#######################################

# Radius of gyration

rg = md.compute_rg(traj)

plt.plot(rg)

plt.xlabel('Frame', fontsize=14)

plt.ylabel('Rg', fontsize=14)

plt.show()

2D Free Energy Landscape 


import numpy as np

import matplotlib

matplotlib.use('Agg')

import matplotlib.pyplot as plt


D = <2D histogram matrix of Q, Rg>

D = -0.5962*np.log(D)      # conversion to free energy in kcal/mol


x = np.linspace(0, 1, nx)

y = np.linspace(0.5, 1.5, ny)

X, Y = np.meshgrid(x, y)


plt.figure(figsize=(8,6))

plt.axes(aspect='equal')

plt.axis([0,1,0.5,1.5])

plt.xlabel('Q')

plt.ylabel('Radius of gyration (nm)')


cmap = plt.get_cmap('terrain')

CS = plt.contourf(X, Y, D, cmap=cmap, alpha=0.8, levels=[-7,-6,-5,-4,-3,-2,-1,0]) 

plt.colorbar(CS, boundaries=[-7,-6,-5,-4,-3,-2,-1,0], label='Free energy (kcal/mol)')


#plt.show()

plt.savefig('output.png', transparent=True, dpi=600)


Ramachandran Plot and DSSP


import numpy as np

import matplotlib

matplotlib.use('Agg')

import matplotlib.pyplot as plt


#######################################

# Secondary structure by DSSP

dssp = md.compute_dssp(traj, simplified=True)

print(dssp)


#######################################

# Phi & Psi angles

phi = md.compute_phi(traj)

print(phi)


Phi and Psi angles in a protein backbone

Ramachandran plot

Homework