In [57]:
%pylab inline
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show


Populating the interactive namespace from numpy and matplotlib

Reaction phase plot


In [58]:
StaticAngle=80/360*2*pi # The angle at the top-right of the substrate is 80 degrees.
LowStateDist=0.5        # The distance of the 'reaction center' charge from the right plane,
                        # determining the energy of the low energy state.
HighStateDist=1         # The distance of the 'reaction center' charge from the top plane,
                        # determining the energy of the high energy state.

StaticAngComplement=pi-StaticAngle
MaxAngle=arctan((HighStateDist/LowStateDist-cos(StaticAngComplement))/
                sin(StaticAngComplement))   # The angle of the
                                            # transition state.

ReactionCenterLoc = [0,0]   # We place the 'reaction center' at the origin 
                            # and calculate the location of the other points.

LowStateLoc = [LowStateDist,ReactionCenterLoc[1]]
HighStateLoc = [HighStateDist*cos(StaticAngComplement),
                HighStateDist*sin(StaticAngComplement)]
TransitionStateLoc = [LowStateLoc[0],tan(MaxAngle)*LowStateLoc[0]]
    
def transitionCoords(angle):    # Given an angle relative to 3 O'clock
                                # from the reaction center, the function
                                # returns the xy position of the state point
    global MaxAngle
    global HighStateDist
    global LowStateDist
    global StaticAngComplement
    if angle > MaxAngle:
        dist=HighStateDist/cos(StaticAngComplement-angle)
    if angle <= MaxAngle:
        dist=LowStateDist/cos(angle)
    xcoord=dist*cos(angle)
    ycoord=dist*sin(angle)
    return array([xcoord,ycoord])

AngRange=linspace(StaticAngComplement,0)
EnzymeHeight=TransitionStateLoc[1]

def substrateEnergy(angle):
    coords = transitionCoords(angle)
    return -1/sqrt(sum(coords**2))

def absElectricPotential(loc1,loc2):    # Returns the absolute value of the electric potential energy of two charges
    return 1/sqrt(sum((loc1-loc2)**2))

def enzymeSubstrateEnergy(angle,dist):  # Calulates, given the state angle of the substrate and the distance between
                                        # the substrate and the catalyst the potential energy of the entire system
    coords = transitionCoords(angle)
    Center = array([0,0])
    EnzymeUpper = array([dist,EnzymeHeight])
    EnzymeLower = array([dist,0])
    EnzymeUpperTransition = -1*absElectricPotential(coords,EnzymeUpper)
    EnzymeLowerTransition = absElectricPotential(coords,EnzymeLower)
    EnzymeUpperCenter = absElectricPotential(Center,EnzymeUpper)
    EnzymeLowerCenter = -1 * absElectricPotential(Center,EnzymeLower)
    
    return  substrateEnergy(angle)+EnzymeUpperTransition+EnzymeUpperCenter+EnzymeLowerTransition+EnzymeLowerCenter
figsize(7,5)
plot(-AngRange,[substrateEnergy(x) for x in AngRange],label='No catalyst present')
plot(-AngRange,[enzymeSubstrateEnergy(x,1.45) for x in AngRange],label='Catalyst at 1.45 distance from substrate',color='r')
legend(loc=3)
xlabel('Reaction state [radians]')
ylabel('Potential energy')
savefig('energylandscapemodificationmodel.svg')



In [10]:
LowerLeft = [1.5*HighStateLoc[0],ReactionCenterLoc[1]]
UpperLeft= [1.5*HighStateLoc[0],1.1*TransitionStateLoc[1]]
UpperRight= [1.1*LowStateLoc[0],1.1*TransitionStateLoc[1]]
LowerRight= [1.1*LowStateLoc[0],ReactionCenterLoc[1]]

StateLine = list(zip( * [LowStateLoc,TransitionStateLoc,HighStateLoc] ))
SubstrateBound = list(zip( * [LowerLeft,UpperLeft,UpperRight,LowerRight,LowerLeft]))
Distances = list(zip( * [LowStateLoc,ReactionCenterLoc,HighStateLoc]))

plot(StateLine[0],StateLine[1],'b')
plot(SubstrateBound[0],SubstrateBound[1],'g',linewidth=2)
plot(Distances[0],Distances[1],'k--')
annotate("Low energy state",color='b', xy=(LowStateLoc[0],LowStateLoc[1]),
         xytext=(LowStateLoc[0]+0.3,LowStateLoc[1]+0.3),
        arrowprops=dict(facecolor='black', shrink=0.1,width=1),size='x-large')
annotate("Transition state",color='b', 
         xy=(TransitionStateLoc[0],TransitionStateLoc[1]),
         xytext=(TransitionStateLoc[0]+0.3,TransitionStateLoc[1]),
        arrowprops=dict(facecolor='black', shrink=0.1,width=1),size='x-large')
annotate("High energy state",color='b', xy=(HighStateLoc[0],HighStateLoc[1]),
         xytext=(HighStateLoc[0]+0.1,HighStateLoc[1]+0.3),
        arrowprops=dict(facecolor='black', shrink=0.1,width=1),size='x-large')

StateAngle=pi/2.5
StateLoc = transitionCoords(StateAngle)

figsize(8,6)
plot(StateLoc[0],StateLoc[1],'ko')

annotate("State point", xy=(StateLoc[0],StateLoc[1]),
         xytext=(StateLoc[0]-0.5,StateLoc[1]-0.2),
        arrowprops=dict(facecolor='black', shrink=0.1,width=1),size='x-large')

annotate("Reaction center", xy=(0,0),
         xytext=(-0.7,0.3),color='red',
        arrowprops=dict(facecolor='black', shrink=0.1,width=1),size='x-large')

plot(0,0,'ro')
StateAngleLine=list(zip(*[ReactionCenterLoc,StateLoc]))
plot(StateAngleLine[0],StateAngleLine[1],'m--')

## Draw the state angle
ax=[]
ay=[]
r=0.2
for a in linspace(0,StateAngle,10):
    ax.append(r*cos(a))
    ay.append(r*sin(a))
plot(ax,ay,'m')    

annotate("State angle",color='m',
         xy=(0.12,0.12),
         xytext=(0.7,0.7),
        arrowprops=dict(facecolor='black', shrink=0.1,width=1),size='x-large')

plt.axes().set_aspect('equal', 'datalim')
ylim(-0.1,1.4)
axis('off')
savefig('subsbody5.pdf')
print(LowerRight)


[0.55, 0]

The catalyst

Our catalyst consists of two opposing charges fixed in a rectangle (see figure below). The upper charge has the same sign as the reaction center and the bottom charge has the same sign as the state point.


In [5]:
CatLowerLeft = [0,0]
CatLowerRight = [1,0]
CatUpperLeft = [0,EnzymeHeight*1.1]
CatUpperRight = [1,EnzymeHeight*1.1]
enzymeLine=list(zip(*[CatLowerLeft,CatLowerRight,
                      CatUpperRight,CatUpperLeft,CatLowerLeft]))
figsize(8,6)
plot(enzymeLine[0],enzymeLine[1],'g',linewidth=2)
plot(0.1,EnzymeHeight,'ro')
plot(0.1,0,'ko')
plt.axes().set_aspect('equal', 'datalim')
ylim(-0.1,1.4)
axis('off')
title('Catalyst sketch')


Out[5]:
<matplotlib.text.Text at 0x7f1ea332f550>

Catalyst effect on the potential energy landscape

In our model, the catalyst location with respect to the substrate is confined to a horizontal line , such that a single number represents the distance between the reaction center and the lower charge of the catalyst. Because the upper charge of the catalyst attracts the state point charge and the lower charge of the catalyst attracts the reaction center, when the substrate is far from the catalyst and the state point charge is in the high energy state, the substrate will be attracted to the catalyst. This can be seen in the figure below as the potential energy of the free substrate (blue) at its high energy state (angle=-1.8) is higher than the potential energy of the substrate at this state but in proximity to the catalyst (green). Furthermore, as the figure below shows, at a distance of 1.45 (arbitrary units), the catalyst charges interact with the substrate charges to produce a modified energy landscape that eliminates the activation energy barrier (green line vs. blue line). Finally, when the substrate is in the low energy state, because the state point charge is closer to the lower catalyst charge than the reaction center charge, the substrate is repelled from the catalyst. This is reflected in the figure as at angle=0 the potential energy of the substrate in proximity to the catalyst (green) is higher than the potential energy of the free substrate (blue).


In [ ]:


In [6]:
figsize(8,6)
plot(-AngRange,[substrateEnergy(x) for x in AngRange],label='No catalyst present')
plot(-AngRange,[enzymeSubstrateEnergy(x,1.45) for x in AngRange],label='Catalyst at 1.45 distance from substrate')
legend(loc=3)
xlabel('State angle')
ylabel('Potential energy')
title('Bound vs free potential energy landscape of the substrate')


Out[6]:
<matplotlib.text.Text at 0x7f1ea322cdd8>

The complete energy landscape and expected dynamics

The complete energy landscape of the interaction between the substrate and the catalyst is shown below. It depicts the potential energy of the system as a function of both the state angle and the distance between the catalyst and the substrate. When starting at a state of high energy free substrate, the system is expected to proceed down the energy landscape to minimize its potential energy. While not necessarily exact, the magenta line represents the expected trace of gradient decsent down the energy landscape tracing the transition of the system until it reaches a minimal potential energy point which is the unbound, low energy state, substrate. The two energy landscapes depicted above (free substrate and bound substrate) are plotted in their respective colors in this 3D plot and represent the bounds of the energy landscape (theoretically the free enzyme is at an infinite distance from the substrate but the figure only shows the energy landscape when the distance between the two is 2, which is large enough to make the interaction negligible).


In [7]:
EnzymeXs=linspace(5*LowStateDist,1.45)  # The bounds of the distance between the
                                        # substrate and the catalyst

X,Y = meshgrid(AngRange,EnzymeXs)
Z=X.copy()
for x in range(X.shape[0]):
    for y in range(X.shape[1]):
        Z[x,y]=enzymeSubstrateEnergy(X[x,y],Y[x,y])
X=-X

In [9]:
fig = plt.figure()
figsize(10,8)
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, 
                     cmap=cm.RdBu,linewidth=0)
ax.set_xlabel('angle')
ax.set_ylabel('dist')
ax.set_zlabel('energy')
x=0
y=0
gradx = [X[x,y]]
grady = [Y[x,y]]
gradz = [Z[x,y]]

def nextIndex(x,y):
    global X
    global Y
    global Z
    current = Z[x,y]
    for i in range(max(x-1,0),min(x+2,X.shape[0]-1)):
        for j in range(max(y-1,0),min(y+2,X.shape[1]-1)):
            if Z[i,j]<=current:
                current = Z[i,j]
                retx = i
                rety = j
    return (retx,rety)


while nextIndex(x,y) != (x,y):
    (x,y) = nextIndex(x,y)
    gradx.append(X[x,y])
    grady.append(Y[x,y])
    gradz.append(Z[x,y])
    
ax.plot(gradx,grady,gradz,'y',linewidth=2)

ax.set_xlabel('State angle')
ax.set_ylabel('Distance')
ax.set_zlabel('Potentialenergy')

ax.plot(-AngRange,[EnzymeXs[0]] * len(AngRange),
        [substrateEnergy(x) for x in AngRange],'C0',linewidth=2)
ax.plot(-AngRange,[EnzymeXs[-1]] * len(AngRange),
        [enzymeSubstrateEnergy(x,1.45) for x in AngRange],'C3',linewidth=2)

savefig('catlandscape.pdf')



In [52]:
Xs=linspace(3*LowStateDist,LowerLeft[0])  # 
Ys=linspace(0,UpperLeft[1]*2)
BoundaryX=array([UpperLeft[0],UpperRight[0],LowerRight[0]])
BoundaryY=array([UpperLeft[1],UpperRight[1],LowerRight[1]])
BoundaryZ=zeros(3)

ReactionLineX=array([HighStateLoc[0],TransitionStateLoc[0],LowStateLoc[0]])
ReactionLineY=array([HighStateLoc[1],TransitionStateLoc[1],LowStateLoc[1]])

def subsHighPot(x,y):
    return 1/sqrt(x**2+y**2)-1/sqrt((x-HighStateLoc[0])**2+(y-HighStateLoc[1])**2)

X,Y = meshgrid(Xs,Ys)
Z=X.copy()
for x in range(X.shape[0]):
    for y in range(X.shape[1]):
        if X[x,y]>LowerRight[0] or Y[x,y] > UpperRight[1]:
            Z[x,y]=subsHighPot(X[x,y],Y[x,y])
        else:
            Z[x,y]=0

fig = plt.figure()
figsize(10,8)
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, 
                     cmap=cm.RdBu,linewidth=0)

plot(BoundaryX,BoundaryY,BoundaryZ,'g')
plot(ReactionLineX,ReactionLineY,BoundaryZ,'b')
plot([HighStateLoc[0]],[HighStateLoc[1]],[0],'ko')
plot([0],[0],[0],'ro')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('electric potential')
ax.set_xticks([0,0.5,1,1.5])
ax.set_yticks([0,0.5,1,1.5,2,2.5])
ax.set_zticks([-3,-2,-1,0])
figsize(5.5,4)
savefig('subsstartpot.pdf')



In [55]:
def subsTransPot(x,y):
    return 1/sqrt(x**2+y**2)-1/sqrt((x-TransitionStateLoc[0])**2+(y-TransitionStateLoc[1])**2)

X,Y = meshgrid(Xs,Ys)
Z=X.copy()
for x in range(X.shape[0]):
    for y in range(X.shape[1]):
        if X[x,y]>LowerRight[0] or Y[x,y] > UpperRight[1]:
            Z[x,y]=subsTransPot(X[x,y],Y[x,y])
        else:
            Z[x,y]=0

fig = plt.figure()
figsize(10,8)
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, 
                     cmap=cm.RdBu,linewidth=0)
plot(BoundaryX,BoundaryY,BoundaryZ,'g')
plot(ReactionLineX,ReactionLineY,BoundaryZ,'b')
plot([TransitionStateLoc[0]],[TransitionStateLoc[1]],[0],'ko')
plot([0],[0],[0],'ro')

ax.set_zlim(-14,2)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('electric potential')
ax.set_xticks([0,0.5,1,1.5])
ax.set_yticks([0,0.5,1,1.5,2,2.5])
ax.set_zticks([-12,-8,-4,0])
figsize(5.5,4)
savefig('substranspot.pdf')



In [56]:
X,Y = meshgrid(Xs,Ys)

Z=X.copy()
for x in range(X.shape[0]):
    for y in range(X.shape[1]):
        if X[x,y]>LowerRight[0] or Y[x,y] > UpperRight[1]:
            Z[x,y]=subsTransPot(X[x,y],Y[x,y])-subsHighPot(X[x,y],Y[x,y])
        else:
            Z[x,y]=0

fig = plt.figure()
figsize(10,8)
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, 
                     cmap=cm.RdBu,linewidth=0)
plot(BoundaryX,BoundaryY,BoundaryZ,'g')
plot(ReactionLineX,ReactionLineY,BoundaryZ,'b')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('electric potential')
ax.set_xticks([0,0.5,1,1.5])
ax.set_yticks([0,0.5,1,1.5,2,2.5])
ax.set_zticks([-12,-8,-4,0])
figsize(5.5,4)

savefig('subsdiffpot.pdf')



In [46]:
AngRange=linspace(StaticAngComplement,0)
EnzymeHeight=TransitionStateLoc[1]

figsize(7,5)
plot(-AngRange,[substrateEnergy(x) for x in AngRange],label='No catalyst present')
plot(-AngRange,[enzymeSubstrateEnergy(x,LowerRight[0]) for x in AngRange],label='Catalyst at minimal distance from substrate',color='r')
legend(loc=0)
xlabel('Reaction state [radians]')
ylabel('Potential energy')
savefig('energylandscapemodificationmodelMax.pdf')



In [ ]: