Metadynamics of GK Enzyme and Domain Analysis



In [3]:

    
%pylab inline
from scipy.stats import kde









    



Populating the interactive namespace from numpy and matplotlib



In [10]:

    
HILLS_1EX6 = loadtxt('data/HILLS_1EX6')
HILLS_4F4J = loadtxt('data/HILLS_4F4J')
HILLS_1KJW = loadtxt('data/HILLS_1KJW')
CV_1EX6 = loadtxt('data/COLVAR_1EX6')
CV_4F4J = loadtxt('data/COLVAR_4F4J')
CV_1KJW = loadtxt('data/COLVAR_1KJW')

Ramachandran plots of "hinge" residues

These are the $\psi$ and $\phi$ angles located in the amino acids in the hinge region of the GK family. These are defined in another notebook. In the cell below, GK proteins which demonstrate enzyme and domain behavior are compared.



In [12]:

    
figure(figsize=[20,65])

amino_acids_1EX6 = ["Glycine-30","Phenylalanine-31","Serine-32",
                    "Valine-33","Serine-34","Serine-35","Threonine-36","Threonine-37"]
amino_acids_4F4J = ["Glycine-30","Phenylalanine-31","Serine-32",
                    "Valine-33","Proline-34","Serine-35","Threonine-36","Threonine-37"]

fig = plt.figure(figsize=[20,65])

for i in range(0,8):
    ax = fig.add_subplot(8,2,(2*i+1))
    #subplot(8,2,(2*i+1))
    im = plt.imread('figures/plot.png')
    implot = plt.imshow(im, origin='upper', extent=(-180,180,-180,180), alpha=.4)
    ax.plot(CV_1EX6[:,11+i]*180/pi,CV_1EX6[:,i+2]*180/pi,'+')
    xlabel('$\phi$', fontsize=16)
    ylabel('$\psi$',fontsize=16)
    title(amino_acids_1EX6[i] + ' Ramachandran Plot', fontsize=20)
    ax.axhline(linewidth=.5, color='black')
    ax.axvline(linewidth=.5, color='black')
    axis([-180,180,-180,180])
    
    ax = fig.add_subplot(8,2,(2*i+2))
    #subplot(8,2,(2*i+2))
    implot = plt.imshow(im, origin='upper', extent=(-180,180,-180,180), alpha=.4)
    ax.plot(CV_4F4J[:,11+i]*180/pi, CV_4F4J[:,(i+2)]*180/pi, 'r+')
    xlabel('$\phi$', fontsize=16)
    ylabel('$\psi$',fontsize=16)
    title(amino_acids_4F4J[i] + ' Ramachandran Plot', fontsize=20)
    ax.axhline(linewidth=.5, color='black')
    ax.axvline(linewidth=.5, color='black')
    axis([-180,180,-180,180])









    





<matplotlib.figure.Figure at 0x10834ca50>

Calculating the free energy landscape

The free energy surface of an enzyme and domain are compared using results from a Metadynamics simulation (done by Gromacs/Plumed). Each surface is calculated by summing the Gaussian's added metadynamics.



In [11]:

    
figure(figsize=[20,10])

center = HILLS_1EX6[:,1]
width = HILLS_1EX6[:,2]
height = HILLS_1EX6[:,-2]
x = linspace(min(CV_1EX6[:,1]),max(CV_1EX6[:,1]),len(center))
y1 = zeros(len(center))
Gaussian = 0
for i in range(0,len(center)):
    gaussian = height[i]*exp(-(x-center[i])**2/(2*width[i]**2))
    y1 += gaussian
    Gaussian += gaussian
plot(x*10,-y1*.24, 'b')
title('Free Energy Landscape')
xlabel('Distance between lobes ($\AA$)', fontsize=18)
ylabel('Free Energy (kcal/mol)' , fontsize=18)


center = HILLS_4F4J[:,1]
width = HILLS_4F4J[:,2]
height = HILLS_4F4J[:,-2]
x = linspace(min(CV_4F4J[:,1]),max(CV_4F4J[:,1]),len(center))
y1 = zeros(len(center))
Gaussian = 0
for i in range(0,len(center)):
    gaussian = height[i]*exp(-(x-center[i])**2/(2*width[i]**2))
    y1 += gaussian
    Gaussian += gaussian
plot(x*10,-y1*.24, 'r')
title('Free Energy Landscape')
xlabel('Distance between lobes ($\AA$)', fontsize=18)
ylabel('Free Energy (kcal/mol)' , fontsize=18)


center = HILLS_1KJW[:,1]
width = HILLS_1KJW[:,2]
height = HILLS_1KJW[:,-2]
x = linspace(min(CV_1KJW[:,1]),max(CV_1KJW[:,1]),len(center))
y1 = zeros(len(center))
Gaussian = 0
for i in range(0,len(center)):
    gaussian = height[i]*exp(-(x-center[i])**2/(2*width[i]**2))
    y1 += gaussian
    Gaussian += gaussian
plot(x*10,-y1*.24, 'g')
title('Free Energy Landscape')
xlabel('Distance between lobes ($\AA$)', fontsize=18)
ylabel('Free Energy (kcal/mol)' , fontsize=18)









    Out[11]:





<matplotlib.text.Text at 0x10807f750>



In [5]:

    
figure(figsize=[10,7])

time = linspace(0,5000,2505)
time2 = linspace(0,6000,3006)
print len(time)
distance = CV_1EX6[:,1]
bias = CV_1EX6[:,-1]
#time = CV_1EX6[:,0]
print len(distance)
plot(time,distance, 'b')

distance = CV_4F4J[:,1]
bias = CV_4F4J[:,-1]
#time = CV_4F4J[:,0]
plot(time2,distance, 'r')

distance = CV_1KJW[:,1]
bias = CV_1KJW[:,-1]
#time = CV_1KJW[:,0]
plot(time,distance, 'g')

title('Distance of between two lobes versus time')
xlabel('time (ns)')
ylabel('Distance (angtroms)')









    



2505
2505






    Out[5]:





<matplotlib.text.Text at 0x108279950>

Energy Landscape from Ramachandran Plots



In [13]:

    
aa = ["Glycine-30","Phenylalanine-31","Serine-32",
                    "Valine-33","Serine-34","Serine-35","Threonine-36","Threonine-37"]

fig, axes = plt.subplots(ncols=2, nrows=1)
fig.set_figwidth(16)
fig.set_figheight(8)

bb = 1
Phi = CV_4F4J[:,bb+1]*180/pi
Psi = CV_4F4J[:,bb+10]*180/pi
data = array([Psi,Phi])

y = Phi
x = Psi
nbins = 200

# Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
k = kde.gaussian_kde(data)
xi, yi = np.mgrid[-180:180:nbins*1j, -180:180:nbins*1j]
zi = k(np.vstack([xi.flatten(), yi.flatten()]))

ax1 = fig.add_subplot(1,2,1)
ax1.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap='hot')
xlabel('$\phi$', fontsize=20)
ylabel('$\psi$',fontsize=20)
title(aa[bb-1] + ' Ramachandran Plot', fontsize=20)
ax1.axhline(linewidth=.5, color='black')
ax1.axvline(linewidth=.5, color='black')
axis([-180,180,-180,180])

subplot(1,2,2)
im = imread('figures/plot.png')
implot = imshow(im, origin='upper', extent=(-180,180,-180,180), alpha=.4)
plot(Psi,Phi,'+')
xlabel('$\phi$', fontsize=20)
ylabel('$\psi$',fontsize=20)
title(aa[bb-1] +' Ramachandran Plot', fontsize=20)
axhline(linewidth=.5, color='black')
axvline(linewidth=.5, color='black')
axis([-180,180,-180,180])









    Out[13]:





[-180, 180, -180, 180]



In [7]:









    Out[7]:





array([ 0.        ,  2.0183727 , -0.75083694, -2.89619331, -0.77088044,
        2.54351772, -1.20559529, -1.24322182, -3.11587292, -1.65795392,
        0.86154933,  2.89907795,  0.15724387,  0.82028822, -2.71803007,
       -2.14082749,  1.28061089, -1.53569532, -1.77332954, -2.1920222 ,  0.        ])



In [ ]: