In this example, we'll compute the solvent accessible surface area of one of the residues in our protien accross each frame in a MD trajectory. We're going to use our trustly alanine dipeptide trajectory for this calculation, but in a real-world situtation you'll probably want to use a more interesting peptide or protein, especially one with a Trp residue.
In [ ]:
%matplotlib inline
from __future__ import print_function
import numpy as np
import mdtraj as md
We'll use the algorithm from Shrake and Rupley for computing the SASA. Here's the function in MDTraj:
In [ ]:
help(md.shrake_rupley)
In [ ]:
trajectory = md.load('ala2.h5')
sasa = md.shrake_rupley(trajectory)
print(trajectory)
print('sasa data shape', sasa.shape)
The computed sasa array contains the solvent accessible surface area for each atom in each frame of the trajectory. Let's sum over all of the atoms to get the total SASA from all of the atoms in each frame.
In [ ]:
total_sasa = sasa.sum(axis=1)
print(total_sasa.shape)
In [ ]:
from matplotlib.pylab import *
plot(trajectory.time, total_sasa)
xlabel('Time [ps]', size=16)
ylabel('Total SASA (nm)^2', size=16)
show()
We probably don't really have enough data do compute a meaningful autocorrelation, but for more realistic datasets, this might be something that you want to do.
In [ ]:
def autocorr(x):
"Compute an autocorrelation with numpy"
x = x - np.mean(x)
result = np.correlate(x, x, mode='full')
result = result[result.size//2:]
return result / result[0]
semilogx(trajectory.time, autocorr(total_sasa))
xlabel('Time [ps]', size=16)
ylabel('SASA autocorrelation', size=16)
show()
In [ ]: