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')
one_yu8.draw()

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 [ ]:
one_yu8

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 [ ]:
protein.set_energy_model(mdt.models.OpenMMPotential)

In [ ]:
protein.configure_methods()

In [ ]:
mintraj = protein.minimize()

In [ ]:
mintraj.draw()

4. Add integrator and run dynamics


In [ ]:
protein.set_integrator(mdt.integrators.OpenMMLangevin,
                       temperature=300*u.kelvin,
                       timestep=2.0*u.fs,
                       frame_interval=2.0*u.ps)

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

In [ ]:
traj.draw()

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'],
                              chain.c_terminal.atoms['C']))

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