Tutorial 2. Biochemical basics

About      Issues      Tutorials      Documentation </span>

Tutorial 2: Playing with proteins

Here, you'll see how to build, visualize, and simulate a protein structure from the PDB.

In [ ]:
# First, import MDT
import moldesign as mdt
from moldesign import units as u

# This sets up your notebook to draw inline plots:
%matplotlib inline
import numpy as np
from matplotlib.pylab import *

try: import seaborn
except ImportError: pass

1. Download from PDB

In this example, we'll look at 1YU8, a crystal structure of the Villin Headpiece.

In [ ]:
one_yu8 = mdt.read('data/1yu8.pdb')

By evaluating the one_yu8 variable, you can get some basic biochemical information, including metadata about missing residues in this crystal structure (hover over the amino acid sequence to get more information).

In [ ]:

2. Strip water and assign forcefield

Next, we isolate the protein and prepare it using the default Amber forcefield parameters.

In [ ]:
headpiece = mdt.Molecule([res for res in one_yu8.residues if res.type == 'protein'])

In [ ]:
ff = mdt.forcefields.DefaultAmber()
protein = ff.create_prepped_molecule(headpiece)

3. Set up energy model and minimize

Next, we'll set up a full molecular mechanics model using OpenMM, then run a minimization and visualize it.

In [ ]:

In [ ]:

In [ ]:
mintraj = protein.minimize()

In [ ]:

4. Add integrator and run dynamics

In [ ]:

In [ ]:
traj = protein.run(20*u.ps)

In [ ]:

5. Some simple analysis

As in tutorial 1, tutorial objects permit a range of timeseries-based analyses.

In [ ]:
# Plot kinetic energy vs. time
plot(traj.time, traj.kinetic_energy)
xlabel('time / %s' % u.default.time); ylabel('energy / %s' % u.default.energy)

In [ ]:
# Plot time evolution of PHE47's sidechain rotation
residue = protein.chains[0].residues['PHE47']
plot(traj.time, traj.dihedral(residue['CA'], residue['CB']).to(u.degrees))

title('sidechain rotation vs time')
xlabel('time / %s' % u.default.time); ylabel(u'angle / º')

In [ ]:
# Plot distance between C-terminus and N-terminus
chain = protein.chains[0]
plot(traj.time, traj.distance(chain.n_terminal.atoms['N'],

plt.title('bond length vs time')
xlabel('time / %s' % u.default.time); ylabel('distance / %s' % u.default.length)